Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "GridCoordinatesObject.h"
23 : #include "tools/Random.h"
24 : #include "tools/Matrix.h"
25 :
26 : namespace PLMD {
27 : namespace gridtools {
28 :
29 175 : void GridCoordinatesObject::setup( const std::string& geom, const std::vector<bool>& ipbc,
30 : const unsigned& np, const double& fib_cutoff ) {
31 175 : if( geom=="flat" ) {
32 143 : gtype=flat;
33 143 : dimension = ipbc.size();
34 32 : } else if( geom=="fibonacci" ) {
35 32 : gtype=fibonacci;
36 32 : dimension = 3;
37 : } else {
38 0 : plumed_merror( geom + " is invalid geometry type");
39 : }
40 :
41 175 : if( gtype==flat ) {
42 143 : bounds_set=false;
43 143 : npoints=0;
44 143 : pbc.resize( ipbc.size() );
45 328 : for(unsigned i=0; i<ipbc.size(); ++i) {
46 185 : pbc[i]=ipbc[i];
47 : }
48 32 : } else if( gtype==fibonacci ) {
49 32 : if( np==0 ) {
50 4 : plumed_merror("the number of points on the sphere should be at least 1");
51 : }
52 30 : bounds_set=true;
53 30 : root5 = sqrt(5);
54 30 : npoints = np;
55 30 : golden = ( 1 + sqrt(5) ) / 2.0;
56 30 : igolden = golden - 1;
57 30 : fib_increment = 2*pi*igolden;
58 30 : log_golden2 = std::log( golden*golden );
59 30 : fib_offset = 2 / static_cast<double>( npoints );
60 30 : fib_shift = fib_offset/2 - 1;
61 :
62 30 : std::vector<double> icoord( dimension ), jcoord( dimension );
63 : // Find minimum distance between each pair of points
64 30 : std::vector<unsigned> tindices( dimension );
65 30 : std::vector<double> mindists( npoints );
66 7792 : for(unsigned i=0; i<npoints; ++i) {
67 7762 : getFibonacciCoordinates( i, icoord );
68 7762 : mindists[i] = 0;
69 2633750 : for(unsigned j=0; j<npoints; ++j) {
70 2625988 : if( i==j ) {
71 7762 : continue ; // Points are not neighbors to themselves
72 : }
73 2618226 : getFibonacciCoordinates( j, jcoord );
74 : // Calculate the dot product
75 : double dot=0;
76 10472904 : for(unsigned k=0; k<dimension; ++k) {
77 7854678 : dot += icoord[k]*jcoord[k];
78 : }
79 2618226 : if( dot>mindists[i] ) {
80 133331 : mindists[i]=dot;
81 : }
82 : }
83 : }
84 : // And now take minimum of dot products
85 30 : const double minDist=*std::min_element(mindists.begin(),mindists.end());
86 : double final_cutoff;
87 30 : if( fib_cutoff<-1 ) {
88 : final_cutoff=-1;
89 : } else {
90 27 : final_cutoff = cos( acos( fib_cutoff ) + acos( minDist ) );
91 : }
92 :
93 : // And now construct the neighbor list
94 30 : fib_nlist.resize( npoints );
95 7792 : for(unsigned i=0; i<npoints; ++i) {
96 7762 : fib_nlist[i].resize(0);
97 7762 : getFibonacciCoordinates( i, icoord );
98 2633750 : for(unsigned j=0; j<npoints; ++j) {
99 2625988 : if( i==j ) {
100 7762 : continue ; // Points are not neighbors to themselves
101 : }
102 2618226 : getFibonacciCoordinates( j, jcoord );
103 : // Calculate the dot product
104 : double dot=0;
105 10472904 : for(unsigned k=0; k<dimension; ++k) {
106 7854678 : dot += icoord[k]*jcoord[k];
107 : }
108 2618226 : if( dot>final_cutoff ) {
109 1347586 : fib_nlist[i].push_back(j);
110 : }
111 : }
112 : }
113 : } else {
114 0 : plumed_error();
115 : }
116 173 : }
117 :
118 139 : void GridCoordinatesObject::setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax,
119 : const std::vector<std::size_t>& binsin, std::vector<double>& spacing ) {
120 : plumed_dbg_assert( smin.size()==dimension && smax.size()==dimension );
121 139 : plumed_assert( gtype==flat && (spacing.size()==dimension || binsin.size()==dimension) );
122 139 : str_min.resize( dimension );
123 139 : str_max.resize( dimension );
124 139 : nbin.resize( dimension );
125 139 : min.resize( dimension );
126 139 : max.resize( dimension );
127 139 : dx.resize( dimension );
128 139 : stride.resize( dimension );
129 :
130 139 : npoints=1;
131 139 : bounds_set=(smin[0]!="auto" && smax[0]!="auto");
132 139 : if( bounds_set ) {
133 177 : for(unsigned i=1; i<dimension; ++i) {
134 38 : if( smin[i]=="auto" || smax[i]=="auto" ) {
135 0 : bounds_set=false;
136 0 : break;
137 : }
138 : }
139 : }
140 312 : for(unsigned i=0; i<dimension; ++i) {
141 177 : str_min[i]=smin[i];
142 : str_max[i]=smax[i];
143 177 : if( bounds_set ) {
144 177 : Tools::convert( str_min[i], min[i] );
145 177 : Tools::convert( str_max[i], max[i] );
146 177 : if( max[i]<=min[i] ) {
147 2 : plumed_merror("GRID_MAX should be larger than GRID_MIN");
148 : }
149 : }
150 176 : if( spacing.size()==dimension && (spacing[i]<0 || (spacing[i]==0 && binsin.size()!=dimension)) ) {
151 4 : plumed_merror("GRID_SPACING should be positive");
152 : }
153 174 : if( spacing.size()==dimension && binsin.size()==dimension ) {
154 0 : if( spacing[i]==0 ) {
155 0 : nbin[i] = binsin[i];
156 0 : } else if( bounds_set ) {
157 0 : double range = max[i] - min[i];
158 0 : nbin[i] = std::round( range / spacing[i]);
159 0 : dx[i]=spacing[i];
160 : // This check ensures that nbins is set correctly if spacing is set the same as the number of bins
161 0 : if( nbin[i]!=binsin[i] ) {
162 0 : plumed_merror("mismatch between input spacing and input number of bins");
163 : }
164 : }
165 174 : } else if( binsin.size()==dimension ) {
166 171 : nbin[i]=binsin[i];
167 171 : dx[i] = ( max[i] - min[i] ) / static_cast<double>( nbin[i] );
168 3 : } else if( spacing.size()==dimension && bounds_set ) {
169 : // as in METAD: the grid ends at GRID_MAX, with a spacing not larger than GRID_SPACING
170 3 : nbin[i] = std::ceil(( max[i] - min[i] ) / spacing[i]);
171 3 : dx[i] = ( max[i] - min[i] ) / static_cast<double>( nbin[i] );
172 0 : } else if( bounds_set ) {
173 0 : plumed_error();
174 : }
175 174 : if( bounds_set && nbin[i]==0 ) {
176 2 : plumed_merror("GRID_BIN should be at least 1");
177 : }
178 173 : if( !pbc[i] ) {
179 93 : max[i] +=dx[i];
180 93 : nbin[i]+=1;
181 : }
182 173 : stride[i]=npoints;
183 173 : npoints*=nbin[i];
184 : }
185 135 : if( spacing.size()!=dimension && bounds_set ) {
186 132 : spacing.resize(dimension);
187 302 : for(unsigned i=0; i<dimension; ++i) {
188 170 : spacing[i]=dx[i];
189 : }
190 : }
191 135 : }
192 :
193 42387699 : unsigned GridCoordinatesObject::getIndex( const std::vector<unsigned>& indices ) const {
194 : plumed_dbg_assert( gtype==flat && indices.size()==dimension );
195 : // indices are flattended using a column-major order
196 42387699 : unsigned index=indices[dimension-1];
197 122020144 : for(unsigned i=dimension-1; i>0; --i) {
198 79632445 : index=index*nbin[i-1]+indices[i-1];
199 : }
200 42387699 : return index;
201 : }
202 :
203 0 : bool GridCoordinatesObject::inbounds( const std::vector<double>& point ) const {
204 0 : return inbounds( View<const double>( point.data(), point.size() ) );
205 : }
206 :
207 0 : bool GridCoordinatesObject::inbounds( const View<const double> point ) const {
208 0 : if( gtype==fibonacci ) {
209 : return true;
210 : }
211 : plumed_dbg_assert( bounds_set && point.size()==dimension );
212 0 : for(unsigned i=0; i<dimension; ++i) {
213 0 : if( pbc[i] ) {
214 0 : continue;
215 : }
216 0 : if( point[i]<min[i] || point[i]>(max[i]-dx[i]) ) {
217 : return false;
218 : }
219 : }
220 : return true;
221 : }
222 :
223 0 : void GridCoordinatesObject::getIndices( const std::vector<double>& point, std::vector<unsigned>& indices ) const {
224 0 : getIndices( View<const double>( point.data(), point.size() ), indices );
225 0 : }
226 :
227 1169232 : void GridCoordinatesObject::getIndices( const View<const double> point, std::vector<unsigned>& indices ) const {
228 : plumed_dbg_assert( gtype==flat && bounds_set && point.size()==dimension && indices.size()==dimension );
229 3453599 : for(unsigned i=0; i<dimension; ++i) {
230 2284368 : indices[i]=std::floor( (point[i] - min[i])/dx[i] );
231 2284368 : if( pbc[i] ) {
232 88672 : indices[i]=indices[i]%nbin[i];
233 2195696 : } else if( indices[i]>=nbin[i] ) {
234 2 : plumed_merror("point is outside grid range");
235 : }
236 : }
237 1169231 : }
238 :
239 2520 : unsigned GridCoordinatesObject::getIndex( const std::vector<double>& point ) const {
240 2520 : return getIndex( View<const double>(point.data(),point.size()) );
241 : }
242 :
243 594217 : unsigned GridCoordinatesObject::getIndex( const View<const double> point ) const {
244 : plumed_dbg_assert( bounds_set && point.size()==dimension );
245 594217 : if( gtype==flat ) {
246 592117 : std::vector<unsigned> indices(dimension);
247 592117 : getIndices( point, indices );
248 592116 : return getIndex( indices );
249 2100 : } else if( gtype==fibonacci ) {
250 2100 : return getFibonacciIndex( point );
251 : } else {
252 0 : plumed_error();
253 : }
254 : }
255 :
256 134280 : unsigned GridCoordinatesObject::getFibonacciIndex( const View<const double> p ) const {
257 : plumed_dbg_assert( gtype==fibonacci );
258 : // Convert input point to coordinates on cylinder
259 : int k=2;
260 134280 : double phi = atan2( p[2], p[0] ), sinthet2 = 1 - p[1]*p[1];
261 : // Calculate power to raise golden ratio
262 134280 : if( sinthet2<epsilon ) {
263 : k = 2;
264 : } else {
265 134280 : k = std::floor( std::log( npoints*pi*root5*sinthet2 ) / log_golden2 );
266 134280 : if( k<2 ) {
267 : k = 2;
268 : }
269 : }
270 134280 : double Fk = pow( golden, k ) / root5, F0 = std::round(Fk), F1 = std::round(Fk*golden);
271 134280 : Matrix<double> B(2,2), invB(2,2);
272 134280 : std::vector<double> thisp(3);
273 134280 : B(0,0) = 2*pi*((F0+1)*igolden - std::floor((F0+1)*igolden)) - fib_increment;
274 134280 : B(0,1) = 2*pi*((F1+1)*igolden - std::floor((F1+1)*igolden)) - fib_increment;
275 134280 : B(1,0) = -2*F0/npoints;
276 134280 : B(1,1) = -2*F1/npoints;
277 134280 : Invert( B, invB );
278 134280 : std::vector<double> vv(2), rc(2);
279 134280 : vv[0]=-phi;
280 134280 : vv[1] = p[1] - fib_shift;
281 134280 : mult( invB, vv, rc );
282 134280 : std::vector<int> c(2);
283 134280 : c[0]=std::floor(rc[0]);
284 134280 : c[1]=std::floor(rc[1]);
285 : unsigned outind=0;
286 : double mindist = 10000000.;
287 671400 : for(int s=0; s<4; ++s) {
288 537120 : double ttt, costheta = B(1,0)*( c[0] + s%2 ) + B(1,1)*( c[1] + s/2 ) + fib_shift;
289 537120 : if( costheta>1 ) {
290 : ttt=1;
291 530862 : } else if( costheta<-1 ) {
292 : ttt=-1;
293 : } else {
294 : ttt=costheta;
295 : }
296 537120 : costheta = 2*ttt - costheta;
297 537120 : unsigned i = std::floor( 0.5*npoints*(1+costheta) );
298 537120 : getFibonacciCoordinates( i, thisp );
299 : double dist=0;
300 2148480 : for(unsigned j=0; j<3; ++j) {
301 1611360 : double tmp=thisp[j]-p[j];
302 1611360 : dist += tmp*tmp;
303 : }
304 537120 : if( dist<mindist ) {
305 294295 : outind = i;
306 : mindist = dist;
307 : }
308 : }
309 134280 : return outind;
310 134280 : }
311 :
312 56848157 : void GridCoordinatesObject::convertIndexToIndices( const unsigned& index, const std::vector<unsigned>& nnbin, std::vector<unsigned>& indices ) const {
313 : plumed_dbg_assert( gtype==flat );
314 56848157 : unsigned kk=index;
315 56848157 : indices[0]=index%nnbin[0];
316 107861389 : for(unsigned i=1; i<dimension-1; ++i) {
317 51013232 : kk=(kk-indices[i-1])/nnbin[i-1];
318 51013232 : indices[i]=kk%nnbin[i];
319 : }
320 56848157 : if(dimension>=2) { // I think this is wrong
321 55883264 : indices[dimension-1]=(kk-indices[dimension-2])/nnbin[dimension-2];
322 : }
323 56848157 : }
324 :
325 5177699 : void GridCoordinatesObject::getIndices( const unsigned& index, std::vector<unsigned>& indices ) const {
326 : plumed_dbg_assert( gtype==flat );
327 5177699 : convertIndexToIndices( index, nbin, indices );
328 5177699 : }
329 :
330 1175350 : void GridCoordinatesObject::getGridPointCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
331 1175350 : std::vector<unsigned> tindices( dimension );
332 1175350 : getGridPointCoordinates( ipoint, tindices, x );
333 1175350 : }
334 :
335 1985967 : void GridCoordinatesObject::getGridPointCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
336 : plumed_dbg_assert( bounds_set && x.size()==dimension && tindices.size()==dimension && ipoint<npoints );
337 1985967 : if( gtype==flat ) {
338 1975071 : getFlatGridCoordinates( ipoint, tindices, x );
339 10896 : } else if( gtype==fibonacci ) {
340 10896 : getFibonacciCoordinates( ipoint, x );
341 : } else {
342 0 : plumed_error();
343 : }
344 1985967 : }
345 :
346 0 : void GridCoordinatesObject::putCoordinateAtValue( const unsigned& ind, const double& val, std::vector<double>& coords ) const {
347 0 : std::vector<double> point( dimension );
348 0 : getGridPointCoordinates( ind, point );
349 0 : if( gtype==flat ) {
350 0 : if( coords.size()!=(dimension+1) ) {
351 0 : coords.resize( (dimension+1) );
352 : }
353 0 : for(unsigned i=0; i<dimension; ++i) {
354 0 : coords[i]=point[i];
355 : }
356 0 : coords[point.size()]=val;
357 0 : } else if( gtype==fibonacci ) {
358 0 : if( coords.size()!=3 ) {
359 0 : coords.resize(3);
360 : }
361 0 : for(unsigned i=0; i<3; ++i) {
362 0 : coords[i] = val*point[i];
363 : }
364 : } else {
365 0 : plumed_error();
366 : }
367 0 : }
368 :
369 1975071 : void GridCoordinatesObject::getFlatGridCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
370 : plumed_dbg_assert( gtype==flat );
371 1975071 : getIndices( ipoint, tindices );
372 6193793 : for(unsigned i=0; i<dimension; ++i) {
373 4218722 : x[i] = min[i] + dx[i]*tindices[i];
374 : }
375 1975071 : }
376 :
377 5799992 : void GridCoordinatesObject::getFibonacciCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
378 : plumed_dbg_assert( gtype==fibonacci );
379 5799992 : x[1] = (ipoint*fib_offset) + fib_shift;
380 5799992 : double r = sqrt( 1 - x[1]*x[1] );
381 5799992 : double phi = ipoint*fib_increment;
382 5799992 : x[0] = r*cos(phi);
383 5799992 : x[2] = r*sin(phi);
384 : double norm=0;
385 23199968 : for(unsigned j=0; j<3; ++j) {
386 17399976 : norm+=x[j]*x[j];
387 : }
388 5799992 : norm = sqrt(norm);
389 23199968 : for(unsigned j=0; j<3; ++j) {
390 17399976 : x[j] = x[j] / norm;
391 : }
392 5799992 : }
393 :
394 552562 : void GridCoordinatesObject::getSplineNeighbors( const unsigned& mybox, unsigned& nneighbors, std::vector<unsigned>& mysneigh ) const {
395 : plumed_dbg_assert( gtype==flat );
396 552562 : mysneigh.resize( static_cast<unsigned>(pow(2.,dimension)) );
397 :
398 : unsigned inind;
399 552562 : nneighbors = 0;
400 552562 : std::vector<unsigned> tmp_indices( dimension );
401 552562 : std::vector<unsigned> my_indices( dimension );
402 552562 : getIndices( mybox, my_indices );
403 2779388 : for(unsigned i=0; i<mysneigh.size(); ++i) {
404 : unsigned tmp=i;
405 : inind=0;
406 6750996 : for(unsigned j=0; j<dimension; ++j) {
407 4524170 : unsigned i0=tmp%2+my_indices[j];
408 4524170 : tmp/=2;
409 4524170 : if(!pbc[j] && i0==nbin[j]) {
410 40800 : continue;
411 : }
412 4483370 : if( pbc[j] && i0==nbin[j]) {
413 : i0=0;
414 : }
415 4483370 : tmp_indices[inind++]=i0;
416 : }
417 2226826 : if( inind==dimension ) {
418 2186226 : mysneigh[nneighbors++]=getIndex( tmp_indices );
419 : }
420 : }
421 552562 : }
422 :
423 561694 : std::vector<std::string> GridCoordinatesObject::getMin() const {
424 : plumed_dbg_assert( gtype==flat );
425 561694 : return str_min;
426 : }
427 :
428 561694 : std::vector<std::string> GridCoordinatesObject::getMax() const {
429 : plumed_dbg_assert( gtype==flat );
430 561694 : return str_max;
431 : }
432 :
433 573146 : std::vector<std::size_t> GridCoordinatesObject::getNbin( const bool& shape ) const {
434 : plumed_dbg_assert( gtype==flat && nbin.size()==dimension );
435 573146 : std::vector<std::size_t> ngrid( dimension );
436 1951806 : for(unsigned i=0; i<dimension; ++i) {
437 1378660 : if( !pbc[i] && !shape ) {
438 665970 : ngrid[i]=nbin[i] - 1;
439 : } else {
440 712690 : ngrid[i]=nbin[i];
441 : }
442 : }
443 573146 : return ngrid;
444 : }
445 :
446 0 : void GridCoordinatesObject::getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
447 : unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
448 0 : getNeighbors( View( pp.data(), pp.size() ), nneigh, num_neighbors, neighbors );
449 0 : }
450 :
451 343692 : void GridCoordinatesObject::getNeighbors( const View<const double> pp, const std::vector<unsigned>& nneigh,
452 : unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
453 : plumed_dbg_assert( bounds_set );
454 :
455 343692 : if( gtype == flat ) {
456 : plumed_dbg_assert( nneigh.size()==dimension );
457 211512 : std::vector<unsigned> indices( dimension );
458 490397 : for(unsigned i=0; i<dimension; ++i) {
459 278885 : indices[i] = std::floor( (pp[i]-min[i])/dx[i] );
460 : }
461 211512 : getNeighbors( indices, nneigh, num_neighbors, neighbors );
462 132180 : } else if( gtype == fibonacci ) {
463 132180 : unsigned find = getFibonacciIndex( View<const double>(pp.data(), pp.size()) );
464 132180 : num_neighbors = 1 + fib_nlist[find].size();
465 132180 : if( neighbors.size()<num_neighbors ) {
466 132180 : neighbors.resize( num_neighbors );
467 : }
468 132180 : neighbors[0]=find;
469 18175217 : for(unsigned i=0; i<fib_nlist[find].size(); ++i) {
470 18043037 : neighbors[1+i] = fib_nlist[find][i];
471 : }
472 : } else {
473 0 : plumed_error();
474 : }
475 343692 : }
476 :
477 227976 : void GridCoordinatesObject::getNeighbors( const std::vector<unsigned>& indices, const std::vector<unsigned>& nneigh,
478 : unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
479 : plumed_dbg_assert( gtype==flat && bounds_set && nneigh.size()==dimension );
480 :
481 : unsigned num_neigh=1;
482 227976 : std::vector<unsigned> small_bin( dimension );
483 556253 : for(unsigned i=0; i<dimension; ++i) {
484 328277 : small_bin[i]=(2*nneigh[i]+1);
485 328277 : if( pbc[i] && small_bin[i]>nbin[i] ) {
486 104400 : small_bin[i]=nbin[i];
487 : }
488 328277 : num_neigh *=small_bin[i];
489 : }
490 227976 : if( neighbors.size()!=num_neigh ) {
491 211515 : neighbors.resize( num_neigh );
492 : }
493 :
494 227976 : num_neighbors=0;
495 227976 : std::vector<unsigned> s_indices(dimension);
496 227976 : std::vector<unsigned> t_indices(dimension);
497 51898434 : for(unsigned index=0; index<num_neigh; ++index) {
498 : bool found=true;
499 51670458 : convertIndexToIndices( index, small_bin, s_indices );
500 204692505 : for(unsigned i=0; i<dimension; ++i) {
501 153022047 : int i0=s_indices[i]-nneigh[i]+indices[i];
502 153022047 : if (i0<0) {
503 14398084 : if (!pbc[i]) {
504 : found = false;
505 : } else {
506 8410766 : i0=nbin[i]-(-i0)%nbin[i];
507 : }
508 138623963 : } else if(static_cast<unsigned>(i0)>=nbin[i]) { // i0 is >=0
509 15504616 : if (!pbc[i]) {
510 : found = false;
511 : } else {
512 8653206 : i0%=nbin[i];
513 : }
514 : }
515 153022047 : if(!pbc[i]) {
516 35502337 : if( i0<0 || static_cast<unsigned>(i0)>=nbin[i] ) { // i0 is >=0
517 : found=false;
518 : }
519 : } else {
520 117519710 : if(i0<0) {
521 0 : i0=nbin[i]-(-i0)%nbin[i];
522 : } else {
523 117519710 : i0%=nbin[i];
524 : }
525 : }
526 153022047 : t_indices[i]=static_cast<unsigned>(i0);
527 : }
528 51670458 : if( found ) {
529 38975870 : neighbors[num_neighbors]=getIndex( t_indices );
530 38975870 : num_neighbors++;
531 : }
532 : }
533 227976 : }
534 :
535 : }
536 : }
537 :
|