1 |
|
|
/* |
2 |
|
|
Teem: Tools to process and visualize scientific data and images . |
3 |
|
|
Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago |
4 |
|
|
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann |
5 |
|
|
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah |
6 |
|
|
|
7 |
|
|
This library is free software; you can redistribute it and/or |
8 |
|
|
modify it under the terms of the GNU Lesser General Public License |
9 |
|
|
(LGPL) as published by the Free Software Foundation; either |
10 |
|
|
version 2.1 of the License, or (at your option) any later version. |
11 |
|
|
The terms of redistributing and/or modifying this software also |
12 |
|
|
include exceptions to the LGPL that facilitate static linking. |
13 |
|
|
|
14 |
|
|
This library 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 GNU |
17 |
|
|
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 this library; if not, write to Free Software Foundation, Inc., |
21 |
|
|
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
22 |
|
|
*/ |
23 |
|
|
|
24 |
|
|
|
25 |
|
|
#include "pull.h" |
26 |
|
|
#include "privatePull.h" |
27 |
|
|
|
28 |
|
|
int |
29 |
|
|
_pullPointProcessNeighLearn(pullTask *task, pullBin *bin, pullPoint *point) { |
30 |
|
|
|
31 |
|
|
/* sneaky: we just learn (and discard) the energy, and this function |
32 |
|
|
will do the work of learning the neighbors */ |
33 |
|
|
_pullEnergyFromPoints(task, bin, point, NULL); |
34 |
|
|
return 0; |
35 |
|
|
} |
36 |
|
|
|
37 |
|
|
static double |
38 |
|
|
_pointEnergyOfNeighbors(pullTask *task, pullBin *bin, pullPoint *point, |
39 |
|
|
double *fracNixed) { |
40 |
|
|
double enr; |
41 |
|
|
unsigned int ii, xx; |
42 |
|
|
pullPoint *her; |
43 |
|
|
|
44 |
|
|
enr = 0; |
45 |
|
|
xx = 0; |
46 |
|
|
for (ii=0; ii<point->neighPointNum; ii++) { |
47 |
|
|
her = point->neighPoint[ii]; |
48 |
|
|
if (her->status & PULL_STATUS_NIXME_BIT) { |
49 |
|
|
xx += 1; |
50 |
|
|
} else { |
51 |
|
|
enr += _pullEnergyFromPoints(task, bin, her, NULL); |
52 |
|
|
} |
53 |
|
|
} |
54 |
|
|
*fracNixed = (point->neighPointNum |
55 |
|
|
? AIR_CAST(double, xx)/point->neighPointNum |
56 |
|
|
: 0); |
57 |
|
|
return enr; |
58 |
|
|
} |
59 |
|
|
|
60 |
|
|
int |
61 |
|
|
_pullPointProcessAdding(pullTask *task, pullBin *bin, pullPoint *point) { |
62 |
|
|
static const char me[]="_pullPointProcessAdding"; |
63 |
|
|
unsigned int npi, iter, api; |
64 |
|
|
double noffavg[4], npos[4], enrWith, enrWithout, |
65 |
|
|
fracNixed, newSpcDist, tmp; |
66 |
|
|
pullPoint *newpnt; |
67 |
|
|
int E; |
68 |
|
|
|
69 |
|
|
task->pctx->count[pullCountAdding] += 1; |
70 |
|
|
if (point->neighPointNum && task->pctx->targetDim |
71 |
|
|
&& task->pctx->flag.popCntlEnoughTest) { |
72 |
|
|
unsigned int plenty, tardim; |
73 |
|
|
tardim = task->pctx->targetDim; |
74 |
|
|
if (task->pctx->flag.zeroZ && tardim > 1) { |
75 |
|
|
/* GLK unsure of tardim == 1 logic here */ |
76 |
|
|
tardim -= 1; |
77 |
|
|
} |
78 |
|
|
plenty = (1 == tardim |
79 |
|
|
? 3 |
80 |
|
|
: (2 == tardim |
81 |
|
|
? 7 |
82 |
|
|
: (3 == tardim |
83 |
|
|
? 13 /* = 1 + 12 |
84 |
|
|
= 1 + coordination number of 3D sphere packing */ |
85 |
|
|
: 0 /* shouldn't get here */))); |
86 |
|
|
/* |
87 |
|
|
if (0 == (point->idtag % 100)) { |
88 |
|
|
printf("!%s: #num %d >?= plenty %d\n", me, point->neighPointNum, plenty); |
89 |
|
|
} |
90 |
|
|
*/ |
91 |
|
|
if (point->neighPointNum >= plenty) { |
92 |
|
|
/* there's little chance that adding points will reduce energy */ |
93 |
|
|
return 0; |
94 |
|
|
} |
95 |
|
|
} |
96 |
|
|
/* |
97 |
|
|
printf("!%s: point->pos = (%g,%g,%g,%g)\n", me, |
98 |
|
|
point->pos[0], point->pos[1], point->pos[2], point->pos[3]); |
99 |
|
|
printf("!%s: point->neighPointNum = %u\n", me, |
100 |
|
|
point->neighPointNum); |
101 |
|
|
printf("!%s: task->pctx->targetDim = %u\n", me, |
102 |
|
|
task->pctx->targetDim); |
103 |
|
|
printf("!%s: task->pctx->flag.popCntlEnoughTest = %d\n", me, |
104 |
|
|
task->pctx->flag.popCntlEnoughTest); |
105 |
|
|
*/ |
106 |
|
|
ELL_4V_SET(noffavg, 0, 0, 0, 0); |
107 |
|
|
for (npi=0; npi<point->neighPointNum; npi++) { |
108 |
|
|
double off[4]; |
109 |
|
|
ELL_4V_SUB(off, point->pos, point->neighPoint[npi]->pos); |
110 |
|
|
/* normalize the offset */ |
111 |
|
|
ELL_3V_SCALE(off, 1/task->pctx->sysParm.radiusSpace, off); |
112 |
|
|
if (task->pctx->haveScale) { |
113 |
|
|
off[3] /= task->pctx->sysParm.radiusScale; |
114 |
|
|
} |
115 |
|
|
ELL_4V_INCR(noffavg, off); |
116 |
|
|
} |
117 |
|
|
if (point->neighPointNum) { |
118 |
|
|
/* |
119 |
|
|
if (0 == (point->idtag % 100)) { |
120 |
|
|
printf("%s: len(offavg) %g >?= thresh %g\n", me, |
121 |
|
|
ELL_4V_LEN(noffavg)/point->neighPointNum, |
122 |
|
|
_PULL_NEIGH_OFFSET_SUM_THRESH); |
123 |
|
|
} |
124 |
|
|
*/ |
125 |
|
|
if (ELL_4V_LEN(noffavg)/point->neighPointNum |
126 |
|
|
< _PULL_NEIGH_OFFSET_SUM_THRESH) { |
127 |
|
|
/* we have neighbors, and they seem to be balanced well enough; |
128 |
|
|
don't try to add */ |
129 |
|
|
return 0; |
130 |
|
|
} |
131 |
|
|
} |
132 |
|
|
if (task->pctx->energySpecR->energy->well(&newSpcDist, |
133 |
|
|
task->pctx->energySpecR->parm)) { |
134 |
|
|
/* HEY: if we don't actually have a well, what is the point of |
135 |
|
|
guessing a new distance? */ |
136 |
|
|
newSpcDist = _PULL_NEWPNT_DIST; |
137 |
|
|
} |
138 |
|
|
/* compute offset (normalized) direction from current point location */ |
139 |
|
|
if (!point->neighPointNum) { |
140 |
|
|
/* we had no neighbors, have to pretend like we did */ |
141 |
|
|
airNormalRand_r(noffavg + 0, noffavg + 1, task->rng); |
142 |
|
|
airNormalRand_r(noffavg + 2, noffavg + 3, task->rng); |
143 |
|
|
if (!task->pctx->flag.zeroZ) { |
144 |
|
|
noffavg[2] = 0; |
145 |
|
|
} |
146 |
|
|
if (!task->pctx->haveScale) { |
147 |
|
|
noffavg[3] = 0; |
148 |
|
|
} |
149 |
|
|
} |
150 |
|
|
if (task->pctx->constraint) { |
151 |
|
|
double proj[9], tmpvec[3]; |
152 |
|
|
_pullConstraintTangent(task, point, proj); |
153 |
|
|
ELL_3MV_MUL(tmpvec, proj, noffavg); |
154 |
|
|
ELL_3V_COPY(noffavg, tmpvec); |
155 |
|
|
} |
156 |
|
|
ELL_4V_NORM(noffavg, noffavg, tmp); |
157 |
|
|
ELL_3V_SCALE(noffavg, task->pctx->sysParm.radiusSpace, noffavg); |
158 |
|
|
noffavg[3] *= task->pctx->sysParm.radiusScale; |
159 |
|
|
ELL_4V_SCALE(noffavg, newSpcDist, noffavg); |
160 |
|
|
/* set new point location */ |
161 |
|
|
ELL_4V_ADD2(npos, noffavg, point->pos); |
162 |
|
|
/* |
163 |
|
|
printf("!%s: new test pos @ (%g,%g,%g,%g)\n", me, |
164 |
|
|
npos[0], npos[1], npos[2], npos[3]); |
165 |
|
|
*/ |
166 |
|
|
if (!_pullInsideBBox(task->pctx, npos)) { |
167 |
|
|
if (task->pctx->verbose > 2) { |
168 |
|
|
printf("%s: new pnt would start (%g,%g,%g,%g) outside bbox, nope\n", |
169 |
|
|
me, npos[0], npos[1], npos[2], npos[3]); |
170 |
|
|
} |
171 |
|
|
return 0; |
172 |
|
|
} |
173 |
|
|
/* initial pos is good, now we start getting serious */ |
174 |
|
|
newpnt = pullPointNew(task->pctx); |
175 |
|
|
if (!newpnt) { |
176 |
|
|
biffAddf(PULL, "%s: couldn't spawn new point from %u", me, point->idtag); |
177 |
|
|
return 1; |
178 |
|
|
} |
179 |
|
|
ELL_4V_COPY(newpnt->pos, npos); |
180 |
|
|
/* set status to indicate this is an unbinned point, with no |
181 |
|
|
knowledge of its neighbors */ |
182 |
|
|
newpnt->status |= PULL_STATUS_NEWBIE_BIT; |
183 |
|
|
/* satisfy constraint if needed */ |
184 |
|
|
if (task->pctx->constraint) { |
185 |
|
|
int constrFail; |
186 |
|
|
if (_pullConstraintSatisfy(task, newpnt, |
187 |
|
|
_PULL_CONSTRAINT_TRAVEL_MAX, |
188 |
|
|
&constrFail)) { |
189 |
|
|
biffAddf(PULL, "%s: on newbie point %u (spawned from %u)", me, |
190 |
|
|
newpnt->idtag, point->idtag); |
191 |
|
|
pullPointNix(newpnt); return 1; |
192 |
|
|
} |
193 |
|
|
if (constrFail) { |
194 |
|
|
/* constraint satisfaction failed, which isn't an error for us, |
195 |
|
|
we just don't try to add this point. Can do immediate nix |
196 |
|
|
because no neighbors know about this point. */ |
197 |
|
|
pullPointNix(newpnt); |
198 |
|
|
return 0; |
199 |
|
|
} |
200 |
|
|
if (!_pullInsideBBox(task->pctx, newpnt->pos)) { |
201 |
|
|
if (task->pctx->verbose > 2) { |
202 |
|
|
printf("%s: post constr newpnt %u (%g,%g,%g,%g) outside bbox; nope\n", |
203 |
|
|
me, newpnt->idtag, newpnt->pos[0], newpnt->pos[1], |
204 |
|
|
newpnt->pos[2], newpnt->pos[3]); |
205 |
|
|
} |
206 |
|
|
newpnt = pullPointNix(newpnt); |
207 |
|
|
return 0; |
208 |
|
|
} |
209 |
|
|
} |
210 |
|
|
/* |
211 |
|
|
printf("!%s: new test pos @ (%g,%g,%g,%g)\n", me, |
212 |
|
|
newpnt->pos[0], newpnt->pos[1], newpnt->pos[2], newpnt->pos[3]); |
213 |
|
|
*/ |
214 |
|
|
/* do some descent, on this point only, which (HACK!) we do by |
215 |
|
|
changing the per-task process mode . . . */ |
216 |
|
|
task->processMode = pullProcessModeDescent; |
217 |
|
|
E = 0; |
218 |
|
|
for (iter=0; iter<task->pctx->iterParm.addDescent; iter++) { |
219 |
|
|
double diff[4]; |
220 |
|
|
if (!E) E |= _pullPointProcessDescent(task, bin, newpnt, |
221 |
|
|
/* ignoreImage; actually it is |
222 |
|
|
this use which motivated the |
223 |
|
|
creation of ignoreImage */ |
224 |
|
|
AIR_TRUE); |
225 |
|
|
if (newpnt->status & PULL_STATUS_STUCK_BIT) { |
226 |
|
|
if (task->pctx->verbose > 2) { |
227 |
|
|
printf("%s: possible newpnt %u stuck @ iter %u; nope\n", me, |
228 |
|
|
newpnt->idtag, iter); |
229 |
|
|
} |
230 |
|
|
newpnt = pullPointNix(newpnt); |
231 |
|
|
/* if we don't change the mode back, then pullBinProcess() won't |
232 |
|
|
know to try adding for the rest of the bins it sees, bad HACK */ |
233 |
|
|
task->processMode = pullProcessModeAdding; |
234 |
|
|
return 0; |
235 |
|
|
} |
236 |
|
|
if (!_pullInsideBBox(task->pctx, newpnt->pos)) { |
237 |
|
|
if (task->pctx->verbose > 2) { |
238 |
|
|
printf("%s: newpnt %u went (%g,%g,%g,%g) outside bbox; nope\n", |
239 |
|
|
me, newpnt->idtag, newpnt->pos[0], newpnt->pos[1], |
240 |
|
|
newpnt->pos[2], newpnt->pos[3]); |
241 |
|
|
} |
242 |
|
|
newpnt = pullPointNix(newpnt); |
243 |
|
|
task->processMode = pullProcessModeAdding; |
244 |
|
|
return 0; |
245 |
|
|
} |
246 |
|
|
ELL_4V_SUB(diff, newpnt->pos, point->pos); |
247 |
|
|
ELL_3V_SCALE(diff, 1/task->pctx->sysParm.radiusSpace, diff); |
248 |
|
|
diff[3] /= task->pctx->sysParm.radiusScale; |
249 |
|
|
if (ELL_4V_LEN(diff) > _PULL_NEWPNT_STRAY_DIST) { |
250 |
|
|
if (task->pctx->verbose > 2) { |
251 |
|
|
printf("%s: newpnt %u went too far %g from old point %u; nope\n", |
252 |
|
|
me, newpnt->idtag, ELL_4V_LEN(diff), point->idtag); |
253 |
|
|
} |
254 |
|
|
newpnt = pullPointNix(newpnt); |
255 |
|
|
task->processMode = pullProcessModeAdding; |
256 |
|
|
return 0; |
257 |
|
|
} |
258 |
|
|
/* still okay to continue descending */ |
259 |
|
|
newpnt->stepEnergy *= task->pctx->sysParm.opporStepScale; |
260 |
|
|
} |
261 |
|
|
/* |
262 |
|
|
printf("!%s: new test pos @ (%g,%g,%g,%g)\n", me, |
263 |
|
|
newpnt->pos[0], newpnt->pos[1], newpnt->pos[2], newpnt->pos[3]); |
264 |
|
|
{ |
265 |
|
|
double posdiff[4]; |
266 |
|
|
ELL_4V_SUB(posdiff, newpnt->pos, point->pos); |
267 |
|
|
printf("!%s: at dist %g\n", me, ELL_4V_LEN(posdiff)); |
268 |
|
|
} |
269 |
|
|
*/ |
270 |
|
|
/* now that newbie point is final test location, see if it meets |
271 |
|
|
the live thresh, if there is one */ |
272 |
|
|
if (task->pctx->ispec[pullInfoLiveThresh] |
273 |
|
|
&& 0 > pullPointScalar(task->pctx, newpnt, pullInfoLiveThresh, |
274 |
|
|
NULL, NULL)) { |
275 |
|
|
/* didn't meet threshold */ |
276 |
|
|
newpnt = pullPointNix(newpnt); |
277 |
|
|
task->processMode = pullProcessModeAdding; |
278 |
|
|
return 0; |
279 |
|
|
} |
280 |
|
|
if (task->pctx->ispec[pullInfoLiveThresh2] /* HEY: copy & paste */ |
281 |
|
|
&& 0 > pullPointScalar(task->pctx, newpnt, pullInfoLiveThresh2, |
282 |
|
|
NULL, NULL)) { |
283 |
|
|
/* didn't meet threshold */ |
284 |
|
|
newpnt = pullPointNix(newpnt); |
285 |
|
|
task->processMode = pullProcessModeAdding; |
286 |
|
|
return 0; |
287 |
|
|
} |
288 |
|
|
if (task->pctx->ispec[pullInfoLiveThresh3] /* HEY: copy & paste */ |
289 |
|
|
&& 0 > pullPointScalar(task->pctx, newpnt, pullInfoLiveThresh3, |
290 |
|
|
NULL, NULL)) { |
291 |
|
|
/* didn't meet threshold */ |
292 |
|
|
newpnt = pullPointNix(newpnt); |
293 |
|
|
task->processMode = pullProcessModeAdding; |
294 |
|
|
return 0; |
295 |
|
|
} |
296 |
|
|
/* see if the new point should be nixed because its at a volume edge */ |
297 |
|
|
if (task->pctx->flag.nixAtVolumeEdgeSpace |
298 |
|
|
&& (newpnt->status & PULL_STATUS_EDGE_BIT)) { |
299 |
|
|
newpnt = pullPointNix(newpnt); |
300 |
|
|
task->processMode = pullProcessModeAdding; |
301 |
|
|
return 0; |
302 |
|
|
} |
303 |
|
|
/* no problem with live thresh, test energy by first learn neighbors */ |
304 |
|
|
/* we have to add newpnt to task's add queue, just so that its neighbors |
305 |
|
|
can see it as a possible neighbor */ |
306 |
|
|
api = airArrayLenIncr(task->addPointArr, 1); |
307 |
|
|
task->addPoint[api] = newpnt; |
308 |
|
|
task->processMode = pullProcessModeNeighLearn; |
309 |
|
|
_pullEnergyFromPoints(task, bin, newpnt, NULL); |
310 |
|
|
/* and teach the neighbors their neighbors (possibly including newpnt) */ |
311 |
|
|
for (npi=0; npi<newpnt->neighPointNum; npi++) { |
312 |
|
|
_pullEnergyFromPoints(task, bin, newpnt->neighPoint[npi], NULL); |
313 |
|
|
} |
314 |
|
|
/* now back to the actual mode we came here with */ |
315 |
|
|
task->processMode = pullProcessModeAdding; |
316 |
|
|
enrWith = (_pullEnergyFromPoints(task, bin, newpnt, NULL) |
317 |
|
|
+ _pointEnergyOfNeighbors(task, bin, newpnt, &fracNixed)); |
318 |
|
|
newpnt->status |= PULL_STATUS_NIXME_BIT; /* turn nixme on */ |
319 |
|
|
enrWithout = _pointEnergyOfNeighbors(task, bin, newpnt, &fracNixed); |
320 |
|
|
newpnt->status &= ~PULL_STATUS_NIXME_BIT; /* turn nixme off */ |
321 |
|
|
if (enrWith < enrWithout) { |
322 |
|
|
/* energy is clearly lower *with* newpnt, so we want to add it, which |
323 |
|
|
means keeping it in the add queue where it already is */ |
324 |
|
|
if (task->pctx->logAdd) { |
325 |
|
|
double posdiff[4]; |
326 |
|
|
ELL_4V_SUB(posdiff, newpnt->pos, point->pos); |
327 |
|
|
fprintf(task->pctx->logAdd, "%u %g %g %g %g %g %g\n", newpnt->idtag, |
328 |
|
|
ELL_3V_LEN(posdiff)/task->pctx->sysParm.radiusSpace, |
329 |
|
|
AIR_ABS(posdiff[3])/task->pctx->sysParm.radiusScale, |
330 |
|
|
newpnt->pos[0], newpnt->pos[1], newpnt->pos[2], newpnt->pos[3]); |
331 |
|
|
} |
332 |
|
|
} else { |
333 |
|
|
/* adding point is not an improvement, remove it from the add queue */ |
334 |
|
|
airArrayLenIncr(task->addPointArr, -1); |
335 |
|
|
/* ugh, have to signal to neighbors that its no longer their neighbor. |
336 |
|
|
HEY this is the part that is totally screwed with multiple threads */ |
337 |
|
|
task->processMode = pullProcessModeNeighLearn; |
338 |
|
|
newpnt->status |= PULL_STATUS_NIXME_BIT; |
339 |
|
|
for (npi=0; npi<newpnt->neighPointNum; npi++) { |
340 |
|
|
_pullEnergyFromPoints(task, bin, newpnt->neighPoint[npi], NULL); |
341 |
|
|
} |
342 |
|
|
task->processMode = pullProcessModeAdding; |
343 |
|
|
newpnt = pullPointNix(newpnt); |
344 |
|
|
} |
345 |
|
|
return 0; |
346 |
|
|
} |
347 |
|
|
|
348 |
|
|
int |
349 |
|
|
_pullPointProcessNixing(pullTask *task, pullBin *bin, pullPoint *point) { |
350 |
|
|
static const char me[]="_pullPointProcessNixing"; |
351 |
|
|
double enrWith, enrNeigh, enrWithout, fracNixed; |
352 |
|
|
|
353 |
|
|
task->pctx->count[pullCountNixing] += 1; |
354 |
|
|
|
355 |
|
|
if (0 && 289 == task->pctx->iter) { |
356 |
|
|
fprintf(stderr, "!%s(%04u): hello lthr %p -> %g %g %g\n", me, point->idtag, |
357 |
|
|
task->pctx->ispec[pullInfoLiveThresh], |
358 |
|
|
pullPointScalar(task->pctx, point, pullInfoLiveThresh, NULL, NULL), |
359 |
|
|
point->pos[0], point->pos[1] |
360 |
|
|
); |
361 |
|
|
} |
362 |
|
|
/* if there's a live thresh, do we meet it? */ |
363 |
|
|
if (task->pctx->ispec[pullInfoLiveThresh] |
364 |
|
|
&& 0 > pullPointScalar(task->pctx, point, pullInfoLiveThresh, |
365 |
|
|
NULL, NULL)) { |
366 |
|
|
point->status |= PULL_STATUS_NIXME_BIT; |
367 |
|
|
return 0; |
368 |
|
|
} |
369 |
|
|
/* HEY copy & paste */ |
370 |
|
|
if (task->pctx->ispec[pullInfoLiveThresh2] |
371 |
|
|
&& 0 > pullPointScalar(task->pctx, point, pullInfoLiveThresh2, |
372 |
|
|
NULL, NULL)) { |
373 |
|
|
point->status |= PULL_STATUS_NIXME_BIT; |
374 |
|
|
return 0; |
375 |
|
|
} |
376 |
|
|
/* HEY copy & paste */ |
377 |
|
|
if (task->pctx->ispec[pullInfoLiveThresh3] |
378 |
|
|
&& 0 > pullPointScalar(task->pctx, point, pullInfoLiveThresh3, |
379 |
|
|
NULL, NULL)) { |
380 |
|
|
point->status |= PULL_STATUS_NIXME_BIT; |
381 |
|
|
return 0; |
382 |
|
|
} |
383 |
|
|
|
384 |
|
|
/* if many neighbors have been nixed, then system is far from convergence, |
385 |
|
|
so energy is not a very meaningful guide to whether to nix this point |
386 |
|
|
NOTE that we use this function to *learn* fracNixed */ |
387 |
|
|
enrNeigh = _pointEnergyOfNeighbors(task, bin, point, &fracNixed); |
388 |
|
|
if (fracNixed < task->pctx->sysParm.fracNeighNixedMax) { |
389 |
|
|
/* is energy lower without us around? */ |
390 |
|
|
enrWith = enrNeigh + _pullEnergyFromPoints(task, bin, point, NULL); |
391 |
|
|
point->status |= PULL_STATUS_NIXME_BIT; /* turn nixme on */ |
392 |
|
|
enrWithout = _pointEnergyOfNeighbors(task, bin, point, &fracNixed); |
393 |
|
|
if (enrWith <= enrWithout) { |
394 |
|
|
/* Energy isn't distinctly lowered without the point, so keep it; |
395 |
|
|
turn off nixing. If enrWith == enrWithout == 0, as happens to |
396 |
|
|
isolated points, then the difference between "<=" and "<" |
397 |
|
|
keeps the isolated points from getting nixed */ |
398 |
|
|
point->status &= ~PULL_STATUS_NIXME_BIT; /* turn nixme off */ |
399 |
|
|
} |
400 |
|
|
/* else energy is certainly higher with the point, do nix it */ |
401 |
|
|
} |
402 |
|
|
|
403 |
|
|
return 0; |
404 |
|
|
} |
405 |
|
|
|
406 |
|
|
int |
407 |
|
|
_pullIterFinishNeighLearn(pullContext *pctx) { |
408 |
|
|
static const char me[]="_pullIterFinishNeighLearn"; |
409 |
|
|
|
410 |
|
|
/* a no-op for now */ |
411 |
|
|
AIR_UNUSED(pctx); |
412 |
|
|
AIR_UNUSED(me); |
413 |
|
|
|
414 |
|
|
return 0; |
415 |
|
|
} |
416 |
|
|
|
417 |
|
|
int |
418 |
|
|
_pullIterFinishAdding(pullContext *pctx) { |
419 |
|
|
static const char me[]="_pullIterFinishAdding"; |
420 |
|
|
unsigned int taskIdx; |
421 |
|
|
|
422 |
|
|
pctx->addNum = 0; |
423 |
|
|
for (taskIdx=0; taskIdx<pctx->threadNum; taskIdx++) { |
424 |
|
|
pullTask *task; |
425 |
|
|
task = pctx->task[taskIdx]; |
426 |
|
|
if (task->addPointNum) { |
427 |
|
|
unsigned int pointIdx; |
428 |
|
|
int added; |
429 |
|
|
for (pointIdx=0; pointIdx<task->addPointNum; pointIdx++) { |
430 |
|
|
pullPoint *point; |
431 |
|
|
pullBin *bin; |
432 |
|
|
point = task->addPoint[pointIdx]; |
433 |
|
|
point->status &= ~PULL_STATUS_NEWBIE_BIT; |
434 |
|
|
if (pullBinsPointMaybeAdd(pctx, point, &bin, &added)) { |
435 |
|
|
biffAddf(PULL, "%s: trouble binning new point %u", me, point->idtag); |
436 |
|
|
return 1; |
437 |
|
|
} |
438 |
|
|
if (added) { |
439 |
|
|
pctx->addNum++; |
440 |
|
|
} else { |
441 |
|
|
unsigned int npi, xpi; |
442 |
|
|
if (pctx->verbose) { |
443 |
|
|
printf("%s: decided NOT to add new point %u\n", me, point->idtag); |
444 |
|
|
} |
445 |
|
|
/* HEY: copied from above */ |
446 |
|
|
/* ugh, have to signal to neigs that its no longer their neighbor */ |
447 |
|
|
task->processMode = pullProcessModeNeighLearn; |
448 |
|
|
point->status |= PULL_STATUS_NIXME_BIT; |
449 |
|
|
for (npi=0; npi<point->neighPointNum; npi++) { |
450 |
|
|
_pullEnergyFromPoints(task, bin, point->neighPoint[npi], NULL); |
451 |
|
|
} |
452 |
|
|
task->processMode = pullProcessModeAdding; |
453 |
|
|
/* can't do immediate nix for reasons GLK doesn't quite understand */ |
454 |
|
|
xpi = airArrayLenIncr(task->nixPointArr, 1); |
455 |
|
|
task->nixPoint[xpi] = point; |
456 |
|
|
} |
457 |
|
|
} |
458 |
|
|
airArrayLenSet(task->addPointArr, 0); |
459 |
|
|
} |
460 |
|
|
} |
461 |
|
|
if (pctx->verbose && pctx->addNum) { |
462 |
|
|
printf("%s: ADDED %u\n", me, pctx->addNum); |
463 |
|
|
} |
464 |
|
|
return 0; |
465 |
|
|
} |
466 |
|
|
|
467 |
|
|
void |
468 |
|
|
_pullNixTheNixed(pullContext *pctx) { |
469 |
|
|
unsigned int binIdx; |
470 |
|
|
|
471 |
|
|
pctx->nixNum = 0; |
472 |
|
|
for (binIdx=0; binIdx<pctx->binNum; binIdx++) { |
473 |
|
|
pullBin *bin; |
474 |
|
|
unsigned int pointIdx; |
475 |
|
|
bin = pctx->bin + binIdx; |
476 |
|
|
pointIdx = 0; |
477 |
|
|
while (pointIdx < bin->pointNum) { |
478 |
|
|
pullPoint *point; |
479 |
|
|
point = bin->point[pointIdx]; |
480 |
|
|
if (pctx->flag.nixAtVolumeEdgeSpace |
481 |
|
|
&& (point->status & PULL_STATUS_EDGE_BIT)) { |
482 |
|
|
point->status |= PULL_STATUS_NIXME_BIT; |
483 |
|
|
} |
484 |
|
|
if (point->status & PULL_STATUS_NIXME_BIT) { |
485 |
|
|
pullPointNix(point); |
486 |
|
|
/* copy last point pointer to this slot */ |
487 |
|
|
bin->point[pointIdx] = bin->point[bin->pointNum-1]; |
488 |
|
|
airArrayLenIncr(bin->pointArr, -1); /* will decrement bin->pointNum */ |
489 |
|
|
pctx->nixNum++; |
490 |
|
|
} else { |
491 |
|
|
pointIdx++; |
492 |
|
|
} |
493 |
|
|
} |
494 |
|
|
} |
495 |
|
|
return; |
496 |
|
|
} |
497 |
|
|
|
498 |
|
|
int |
499 |
|
|
_pullIterFinishNixing(pullContext *pctx) { |
500 |
|
|
static const char me[]="_pullIterFinishNixing"; |
501 |
|
|
unsigned int taskIdx; |
502 |
|
|
|
503 |
|
|
_pullNixTheNixed(pctx); |
504 |
|
|
/* finish nixing the things that we decided not to add */ |
505 |
|
|
for (taskIdx=0; taskIdx<pctx->threadNum; taskIdx++) { |
506 |
|
|
pullTask *task; |
507 |
|
|
task = pctx->task[taskIdx]; |
508 |
|
|
if (task->nixPointNum) { |
509 |
|
|
unsigned int xpi; |
510 |
|
|
for (xpi=0; xpi<task->nixPointNum; xpi++) { |
511 |
|
|
pullPointNix(task->nixPoint[xpi]); |
512 |
|
|
} |
513 |
|
|
airArrayLenSet(task->nixPointArr, 0); |
514 |
|
|
} |
515 |
|
|
} |
516 |
|
|
if (pctx->verbose && pctx->nixNum) { |
517 |
|
|
printf("%s: NIXED %u\n", me, pctx->nixNum); |
518 |
|
|
} |
519 |
|
|
return 0; |
520 |
|
|
} |
521 |
|
|
|