SPEED
Data Types | Functions/Subroutines | Variables
kdtree2_module Module Reference

Data Types

type  interval
 
type  kdtree2
 
type  tree_node
 

Functions/Subroutines

type(kdtree2) function, pointer, public kdtree2_create (input_data, dim, sort, rearrange)
 
subroutine, public kdtree2_destroy (tp)
 
subroutine, public kdtree2_n_nearest (tp, qv, nn, results)
 
subroutine, public kdtree2_n_nearest_around_point (tp, idxin, correltime, nn, results)
 
subroutine, public kdtree2_r_nearest (tp, qv, r2, nfound, nalloc, results)
 
subroutine, public kdtree2_r_nearest_around_point (tp, idxin, correltime, r2, nfound, nalloc, results)
 
integer function, public kdtree2_r_count (tp, qv, r2)
 
integer function, public kdtree2_r_count_around_point (tp, idxin, correltime, r2)
 
subroutine, public kdtree2_n_nearest_brute_force (tp, qv, nn, results)
 
subroutine, public kdtree2_r_nearest_brute_force (tp, qv, r2, nfound, results)
 
subroutine, public kdtree2_sort_results (nfound, results)
 

Variables

integer, parameter bucket_size = 12
 

Function/Subroutine Documentation

◆ kdtree2_create()

type (kdtree2) function, pointer, public kdtree2_module::kdtree2_create ( real(kdkind), dimension(:,:), target  input_data,
integer, intent(in), optional  dim,
logical, intent(in), optional  sort,
logical, intent(in), optional  rearrange 
)

Definition at line 611 of file kdtree2.f90.

612 !
613 ! create the actual tree structure, given an input array of data.
614 !
615 ! Note, input data is input_data(1:d,1:N), NOT the other way around.
616 ! THIS IS THE REVERSE OF THE PREVIOUS VERSION OF THIS MODULE.
617 ! The reason for it is cache friendliness, improving performance.
618 !
619 ! Optional arguments: If 'dim' is specified, then the tree
620 ! will only search the first 'dim' components
621 ! of input_data, otherwise, dim is inferred
622 ! from SIZE(input_data,1).
623 !
624 ! if sort .eqv. .true. then output results
625 ! will be sorted by increasing distance.
626 ! default=.false., as it is faster to not sort.
627 !
628 ! if rearrange .eqv. .true. then an internal
629 ! copy of the data, rearranged by terminal node,
630 ! will be made for cache friendliness.
631 ! default=.true., as it speeds searches, but
632 ! building takes longer, and extra memory is used.
633 !
634 ! .. Function Return Cut_value ..
635 type (kdtree2), pointer :: mr
636 integer, intent(in), optional :: dim
637 logical, intent(in), optional :: sort
638 logical, intent(in), optional :: rearrange
639 ! ..
640 ! .. Array Arguments ..
641 real(kdkind), target :: input_data(:,:)
642 !
643 integer :: i
644 ! ..
645 allocate (mr)
646 mr%the_data => input_data
647 ! pointer assignment
648
649 if (present(dim)) then
650 mr%dimen = dim
651 else
652 mr%dimen = size(input_data,1)
653 end if
654 mr%n = size(input_data,2)
655
656 if (mr%dimen > mr%n) then
657 ! unlikely to be correct
658 write (*,*) 'KD_TREE_TRANS: likely user error.'
659 write (*,*) 'KD_TREE_TRANS: You passed in matrix with D=',mr%dimen
660 write (*,*) 'KD_TREE_TRANS: and N=',mr%n
661 write (*,*) 'KD_TREE_TRANS: note, that new format is data(1:D,1:N)'
662 write (*,*) 'KD_TREE_TRANS: with usually N >> D. If N =approx= D, then a k-d tree'
663 write (*,*) 'KD_TREE_TRANS: is not an appropriate data structure.'
664 stop
665 end if
666
667 call build_tree(mr)
668
669 if (present(sort)) then
670 mr%sort = sort
671 else
672 mr%sort = .false.
673 endif
674
675 if (present(rearrange)) then
676 mr%rearrange = rearrange
677 else
678 mr%rearrange = .true.
679 endif
680
681 if (mr%rearrange) then
682 allocate(mr%rearranged_data(mr%dimen,mr%n))
683 do i=1,mr%n
684 mr%rearranged_data(:,i) = mr%the_data(:, &
685 mr%ind(i))
686 enddo
687 else
688 nullify(mr%rearranged_data)
689 endif
690

References bucket_size, and kdtree2_precision_module::kdkind.

Referenced by make_nh_enhanced_nnsearch().

Here is the caller graph for this function:

◆ kdtree2_destroy()

subroutine, public kdtree2_module::kdtree2_destroy ( type (kdtree2), pointer  tp)

Definition at line 988 of file kdtree2.f90.

989 ! Deallocates all memory for the tree, except input data matrix
990 ! .. Structure Arguments ..
991 type (kdtree2), pointer :: tp
992 ! ..
993 call destroy_node(tp%root)
994
995 deallocate (tp%ind)
996 nullify (tp%ind)
997
998 if (tp%rearrange) then
999 deallocate(tp%rearranged_data)
1000 nullify(tp%rearranged_data)
1001 endif
1002
1003 deallocate(tp)
1004 return
1005
1006 contains
1007 recursive subroutine destroy_node(np)
1008 ! .. Structure Arguments ..
1009 type (tree_node), pointer :: np
1010 ! ..
1011 ! .. Intrinsic Functions ..
1012 intrinsic ASSOCIATED
1013 ! ..
1014 if (associated(np%left)) then
1015 call destroy_node(np%left)
1016 nullify (np%left)
1017 end if
1018 if (associated(np%right)) then
1019 call destroy_node(np%right)
1020 nullify (np%right)
1021 end if
1022 if (associated(np%box)) deallocate(np%box)
1023 deallocate(np)
1024 return
1025
1026 end subroutine destroy_node
1027

Referenced by make_nh_enhanced_nnsearch().

Here is the caller graph for this function:

◆ kdtree2_n_nearest()

subroutine, public kdtree2_module::kdtree2_n_nearest ( type (kdtree2), pointer  tp,
real(kdkind), dimension(:), intent(in), target  qv,
integer, intent(in)  nn,
type(kdtree2_result), dimension(:), target  results 
)

Definition at line 1030 of file kdtree2.f90.

1031 ! Find the 'nn' vectors in the tree nearest to 'qv' in euclidean norm
1032 ! returning their indexes and distances in 'indexes' and 'distances'
1033 ! arrays already allocated passed to this subroutine.
1034 type (kdtree2), pointer :: tp
1035 real(kdkind), target, intent (In) :: qv(:)
1036 integer, intent (In) :: nn
1037 type(kdtree2_result), target :: results(:)
1038
1039
1040 sr%ballsize = huge(1.0)
1041 sr%qv => qv
1042 sr%nn = nn
1043 sr%nfound = 0
1044 sr%centeridx = -1
1045 sr%correltime = 0
1046 sr%overflow = .false.
1047
1048 sr%results => results
1049
1050 sr%nalloc = nn ! will be checked
1051
1052 sr%ind => tp%ind
1053 sr%rearrange = tp%rearrange
1054 if (tp%rearrange) then
1055 sr%Data => tp%rearranged_data
1056 else
1057 sr%Data => tp%the_data
1058 endif
1059 sr%dimen = tp%dimen
1060
1061 call validate_query_storage(nn)
1062 sr%pq = pq_create(results)
1063
1064 call search(tp%root)
1065
1066 if (tp%sort) then
1067 call kdtree2_sort_results(nn, results)
1068 endif
1069! deallocate(sr%pqp)
1070 return

References kdtree2_sort_results(), and kdtree2_priority_queue_module::pq_create().

Referenced by make_nh_enhanced_nnsearch().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ kdtree2_n_nearest_around_point()

subroutine, public kdtree2_module::kdtree2_n_nearest_around_point ( type (kdtree2), pointer  tp,
integer, intent(in)  idxin,
integer, intent(in)  correltime,
integer, intent(in)  nn,
type(kdtree2_result), dimension(:), target  results 
)

Definition at line 1073 of file kdtree2.f90.

1074 ! Find the 'nn' vectors in the tree nearest to point 'idxin',
1075 ! with correlation window 'correltime', returing results in
1076 ! results(:), which must be pre-allocated upon entry.
1077 type (kdtree2), pointer :: tp
1078 integer, intent (In) :: idxin, correltime, nn
1079 type(kdtree2_result), target :: results(:)
1080
1081 allocate (sr%qv(tp%dimen))
1082 sr%qv = tp%the_data(:,idxin) ! copy the vector
1083 sr%ballsize = huge(1.0) ! the largest real(kdkind) number
1084 sr%centeridx = idxin
1085 sr%correltime = correltime
1086
1087 sr%nn = nn
1088 sr%nfound = 0
1089
1090 sr%dimen = tp%dimen
1091 sr%nalloc = nn
1092
1093 sr%results => results
1094
1095 sr%ind => tp%ind
1096 sr%rearrange = tp%rearrange
1097
1098 if (sr%rearrange) then
1099 sr%Data => tp%rearranged_data
1100 else
1101 sr%Data => tp%the_data
1102 endif
1103
1104 call validate_query_storage(nn)
1105 sr%pq = pq_create(results)
1106
1107 call search(tp%root)
1108
1109 if (tp%sort) then
1110 call kdtree2_sort_results(nn, results)
1111 endif
1112 deallocate (sr%qv)
1113 return

References kdtree2_sort_results(), and kdtree2_priority_queue_module::pq_create().

Here is the call graph for this function:

◆ kdtree2_n_nearest_brute_force()

subroutine, public kdtree2_module::kdtree2_n_nearest_brute_force ( type (kdtree2), pointer  tp,
real(kdkind), dimension(:), intent(in)  qv,
integer, intent(in)  nn,
type(kdtree2_result), dimension(:)  results 
)

Definition at line 1694 of file kdtree2.f90.

1695 ! find the 'n' nearest neighbors to 'qv' by exhaustive search.
1696 ! only use this subroutine for testing, as it is SLOW! The
1697 ! whole point of a k-d tree is to avoid doing what this subroutine
1698 ! does.
1699 type (kdtree2), pointer :: tp
1700 real(kdkind), intent (In) :: qv(:)
1701 integer, intent (In) :: nn
1702 type(kdtree2_result) :: results(:)
1703
1704 integer :: i, j, k
1705 real(kdkind), allocatable :: all_distances(:)
1706 ! ..
1707 allocate (all_distances(tp%n))
1708 do i = 1, tp%n
1709 all_distances(i) = square_distance(tp%dimen,qv,tp%the_data(:,i))
1710 end do
1711 ! now find 'n' smallest distances
1712 do i = 1, nn
1713 results(i)%dis = huge(1.0)
1714 results(i)%idx = -1
1715 end do
1716 do i = 1, tp%n
1717 if (all_distances(i)<results(nn)%dis) then
1718 ! insert it somewhere on the list
1719 do j = 1, nn
1720 if (all_distances(i)<results(j)%dis) exit
1721 end do
1722 ! now we know 'j'
1723 do k = nn - 1, j, -1
1724 results(k+1) = results(k)
1725 end do
1726 results(j)%dis = all_distances(i)
1727 results(j)%idx = i
1728 end if
1729 end do
1730 deallocate (all_distances)

◆ kdtree2_r_count()

integer function, public kdtree2_module::kdtree2_r_count ( type (kdtree2), pointer  tp,
real(kdkind), dimension(:), intent(in), target  qv,
real(kdkind), intent(in)  r2 
)

Definition at line 1244 of file kdtree2.f90.

1245 ! Count the number of neighbors within square distance 'r2'.
1246 type (kdtree2), pointer :: tp
1247 real(kdkind), target, intent (In) :: qv(:)
1248 real(kdkind), intent(in) :: r2
1249 integer :: nfound
1250 ! ..
1251 ! .. Intrinsic Functions ..
1252 intrinsic huge
1253 ! ..
1254 sr%qv => qv
1255 sr%ballsize = r2
1256
1257 sr%nn = 0 ! flag for fixed r search
1258 sr%nfound = 0
1259 sr%centeridx = -1
1260 sr%correltime = 0
1261
1262 nullify(sr%results) ! for some reason, FTN 95 chokes on '=> null()'
1263
1264 sr%nalloc = 0 ! we do not allocate any storage but that's OK
1265 ! for counting.
1266 sr%ind => tp%ind
1267 sr%rearrange = tp%rearrange
1268 if (tp%rearrange) then
1269 sr%Data => tp%rearranged_data
1270 else
1271 sr%Data => tp%the_data
1272 endif
1273 sr%dimen = tp%dimen
1274
1275 !
1276 !sr%dsl = Huge(sr%dsl) ! set to huge positive values
1277 !sr%il = -1 ! set to invalid indexes
1278 !
1279 sr%overflow = .false.
1280
1281 call search(tp%root)
1282
1283 nfound = sr%nfound
1284
1285 return

◆ kdtree2_r_count_around_point()

integer function, public kdtree2_module::kdtree2_r_count_around_point ( type (kdtree2), pointer  tp,
integer, intent(in)  idxin,
integer, intent(in)  correltime,
real(kdkind), intent(in)  r2 
)

Definition at line 1288 of file kdtree2.f90.

1290 ! Count the number of neighbors within square distance 'r2' around
1291 ! point 'idxin' with decorrelation time 'correltime'.
1292 !
1293 type (kdtree2), pointer :: tp
1294 integer, intent (In) :: correltime, idxin
1295 real(kdkind), intent(in) :: r2
1296 integer :: nfound
1297 ! ..
1298 ! ..
1299 ! .. Intrinsic Functions ..
1300 intrinsic huge
1301 ! ..
1302 allocate (sr%qv(tp%dimen))
1303 sr%qv = tp%the_data(:,idxin)
1304 sr%ballsize = r2
1305
1306 sr%nn = 0 ! flag for fixed r search
1307 sr%nfound = 0
1308 sr%centeridx = idxin
1309 sr%correltime = correltime
1310 nullify(sr%results)
1311
1312 sr%nalloc = 0 ! we do not allocate any storage but that's OK
1313 ! for counting.
1314
1315 sr%ind => tp%ind
1316 sr%rearrange = tp%rearrange
1317
1318 if (sr%rearrange) then
1319 sr%Data => tp%rearranged_data
1320 else
1321 sr%Data => tp%the_data
1322 endif
1323 sr%dimen = tp%dimen
1324
1325 !
1326 !sr%dsl = Huge(sr%dsl) ! set to huge positive values
1327 !sr%il = -1 ! set to invalid indexes
1328 !
1329 sr%overflow = .false.
1330
1331 call search(tp%root)
1332
1333 nfound = sr%nfound
1334
1335 return

◆ kdtree2_r_nearest()

subroutine, public kdtree2_module::kdtree2_r_nearest ( type (kdtree2), pointer  tp,
real(kdkind), dimension(:), intent(in), target  qv,
real(kdkind), intent(in)  r2,
integer, intent(out)  nfound,
integer, intent(in)  nalloc,
type(kdtree2_result), dimension(:), target  results 
)

Definition at line 1116 of file kdtree2.f90.

1117 ! find the nearest neighbors to point 'idxin', within SQUARED
1118 ! Euclidean distance 'r2'. Upon ENTRY, nalloc must be the
1119 ! size of memory allocated for results(1:nalloc). Upon
1120 ! EXIT, nfound is the number actually found within the ball.
1121 !
1122 ! Note that if nfound .gt. nalloc then more neighbors were found
1123 ! than there were storage to store. The resulting list is NOT
1124 ! the smallest ball inside norm r^2
1125 !
1126 ! Results are NOT sorted unless tree was created with sort option.
1127 type (kdtree2), pointer :: tp
1128 real(kdkind), target, intent (In) :: qv(:)
1129 real(kdkind), intent(in) :: r2
1130 integer, intent(out) :: nfound
1131 integer, intent (In) :: nalloc
1132 type(kdtree2_result), target :: results(:)
1133
1134 !
1135 sr%qv => qv
1136 sr%ballsize = r2
1137 sr%nn = 0 ! flag for fixed ball search
1138 sr%nfound = 0
1139 sr%centeridx = -1
1140 sr%correltime = 0
1141
1142 sr%results => results
1143
1144 call validate_query_storage(nalloc)
1145 sr%nalloc = nalloc
1146 sr%overflow = .false.
1147 sr%ind => tp%ind
1148 sr%rearrange= tp%rearrange
1149
1150 if (tp%rearrange) then
1151 sr%Data => tp%rearranged_data
1152 else
1153 sr%Data => tp%the_data
1154 endif
1155 sr%dimen = tp%dimen
1156
1157 !
1158 !sr%dsl = Huge(sr%dsl) ! set to huge positive values
1159 !sr%il = -1 ! set to invalid indexes
1160 !
1161
1162 call search(tp%root)
1163 nfound = sr%nfound
1164 if (tp%sort) then
1165 call kdtree2_sort_results(nfound, results)
1166 endif
1167
1168 if (sr%overflow) then
1169 write (*,*) 'KD_TREE_TRANS: warning! return from kdtree2_r_nearest found more neighbors'
1170 write (*,*) 'KD_TREE_TRANS: than storage was provided for. Answer is NOT smallest ball'
1171 write (*,*) 'KD_TREE_TRANS: with that number of neighbors! I.e. it is wrong.'
1172 endif
1173
1174 return

References kdtree2_sort_results().

Here is the call graph for this function:

◆ kdtree2_r_nearest_around_point()

subroutine, public kdtree2_module::kdtree2_r_nearest_around_point ( type (kdtree2), pointer  tp,
integer, intent(in)  idxin,
integer, intent(in)  correltime,
real(kdkind), intent(in)  r2,
integer, intent(out)  nfound,
integer, intent(in)  nalloc,
type(kdtree2_result), dimension(:), target  results 
)

Definition at line 1177 of file kdtree2.f90.

1179 !
1180 ! Like kdtree2_r_nearest, but around a point 'idxin' already existing
1181 ! in the data set.
1182 !
1183 ! Results are NOT sorted unless tree was created with sort option.
1184 !
1185 type (kdtree2), pointer :: tp
1186 integer, intent (In) :: idxin, correltime, nalloc
1187 real(kdkind), intent(in) :: r2
1188 integer, intent(out) :: nfound
1189 type(kdtree2_result), target :: results(:)
1190 ! ..
1191 ! .. Intrinsic Functions ..
1192 intrinsic huge
1193 ! ..
1194 allocate (sr%qv(tp%dimen))
1195 sr%qv = tp%the_data(:,idxin) ! copy the vector
1196 sr%ballsize = r2
1197 sr%nn = 0 ! flag for fixed r search
1198 sr%nfound = 0
1199 sr%centeridx = idxin
1200 sr%correltime = correltime
1201
1202 sr%results => results
1203
1204 sr%nalloc = nalloc
1205 sr%overflow = .false.
1206
1207 call validate_query_storage(nalloc)
1208
1209 ! sr%dsl = HUGE(sr%dsl) ! set to huge positive values
1210 ! sr%il = -1 ! set to invalid indexes
1211
1212 sr%ind => tp%ind
1213 sr%rearrange = tp%rearrange
1214
1215 if (tp%rearrange) then
1216 sr%Data => tp%rearranged_data
1217 else
1218 sr%Data => tp%the_data
1219 endif
1220 sr%rearrange = tp%rearrange
1221 sr%dimen = tp%dimen
1222
1223 !
1224 !sr%dsl = Huge(sr%dsl) ! set to huge positive values
1225 !sr%il = -1 ! set to invalid indexes
1226 !
1227
1228 call search(tp%root)
1229 nfound = sr%nfound
1230 if (tp%sort) then
1231 call kdtree2_sort_results(nfound,results)
1232 endif
1233
1234 if (sr%overflow) then
1235 write (*,*) 'KD_TREE_TRANS: warning! return from kdtree2_r_nearest found more neighbors'
1236 write (*,*) 'KD_TREE_TRANS: than storage was provided for. Answer is NOT smallest ball'
1237 write (*,*) 'KD_TREE_TRANS: with that number of neighbors! I.e. it is wrong.'
1238 endif
1239
1240 deallocate (sr%qv)
1241 return

References kdtree2_sort_results().

Here is the call graph for this function:

◆ kdtree2_r_nearest_brute_force()

subroutine, public kdtree2_module::kdtree2_r_nearest_brute_force ( type (kdtree2), pointer  tp,
real(kdkind), dimension(:), intent(in)  qv,
real(kdkind), intent(in)  r2,
integer, intent(out)  nfound,
type(kdtree2_result), dimension(:)  results 
)

Definition at line 1734 of file kdtree2.f90.

1735 ! find the nearest neighbors to 'qv' with distance**2 <= r2 by exhaustive search.
1736 ! only use this subroutine for testing, as it is SLOW! The
1737 ! whole point of a k-d tree is to avoid doing what this subroutine
1738 ! does.
1739 type (kdtree2), pointer :: tp
1740 real(kdkind), intent (In) :: qv(:)
1741 real(kdkind), intent (In) :: r2
1742 integer, intent(out) :: nfound
1743 type(kdtree2_result) :: results(:)
1744
1745 integer :: i, nalloc
1746 real(kdkind), allocatable :: all_distances(:)
1747 ! ..
1748 allocate (all_distances(tp%n))
1749 do i = 1, tp%n
1750 all_distances(i) = square_distance(tp%dimen,qv,tp%the_data(:,i))
1751 end do
1752
1753 nfound = 0
1754 nalloc = size(results,1)
1755
1756 do i = 1, tp%n
1757 if (all_distances(i)< r2) then
1758 ! insert it somewhere on the list
1759 if (nfound .lt. nalloc) then
1760 nfound = nfound+1
1761 results(nfound)%dis = all_distances(i)
1762 results(nfound)%idx = i
1763 endif
1764 end if
1765 enddo
1766 deallocate (all_distances)
1767
1768 call kdtree2_sort_results(nfound,results)
1769
1770

References kdtree2_sort_results().

Here is the call graph for this function:

◆ kdtree2_sort_results()

subroutine, public kdtree2_module::kdtree2_sort_results ( integer, intent(in)  nfound,
type(kdtree2_result), dimension(:), target  results 
)

Definition at line 1773 of file kdtree2.f90.

1774 ! Use after search to sort results(1:nfound) in order of increasing
1775 ! distance.
1776 integer, intent(in) :: nfound
1777 type(kdtree2_result), target :: results(:)
1778 !
1779 !
1780
1781 !THIS IS BUGGY WITH INTEL FORTRAN
1782 ! If (nfound .Gt. 1) Call heapsort(results(1:nfound)%dis,results(1:nfound)%ind,nfound)
1783 !
1784 if (nfound .gt. 1) call heapsort_struct(results,nfound)
1785
1786 return

Referenced by kdtree2_n_nearest(), kdtree2_n_nearest_around_point(), kdtree2_r_nearest(), kdtree2_r_nearest_around_point(), and kdtree2_r_nearest_brute_force().

Here is the caller graph for this function:

Variable Documentation

◆ bucket_size

integer, parameter kdtree2_module::bucket_size = 12

Definition at line 519 of file kdtree2.f90.

519 integer, parameter :: bucket_size = 12

Referenced by kdtree2_create().