1 /*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 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 /* the Init() functions are up here for easier reference */
29
30 void
_pullIterParmInit(pullIterParm * iterParm)31 _pullIterParmInit(pullIterParm *iterParm) {
32
33 iterParm->min = 0;
34 iterParm->max = 0;
35 iterParm->stuckMax = 4;
36 iterParm->constraintMax = 15;
37 iterParm->popCntlPeriod = 10;
38 iterParm->addDescent = 10;
39 iterParm->callback = 1;
40 iterParm->snap = 0;
41 iterParm->energyIncreasePermitHalfLife = 0;
42 return;
43 }
44
45 void
_pullSysParmInit(pullSysParm * sysParm)46 _pullSysParmInit(pullSysParm *sysParm) {
47
48 sysParm->alpha = 0.5;
49 sysParm->beta = 0.5;
50 sysParm->gamma = 1;
51 sysParm->separableGammaLearnRescale = 8;
52 sysParm->theta = 0.0;
53 sysParm->wall = 1;
54 sysParm->radiusSpace = 1;
55 sysParm->radiusScale = 1;
56 sysParm->binWidthSpace = 1.001; /* supersititious */
57 sysParm->neighborTrueProb = 1.0;
58 sysParm->probeProb = 1.0;
59 sysParm->stepInitial = 1;
60 sysParm->opporStepScale = 1.0;
61 sysParm->backStepScale = 0.5;
62 sysParm->constraintStepMin = 0.0001;
63 sysParm->energyDecreaseMin = 0.001;
64 sysParm->energyDecreasePopCntlMin = 0.02;
65 sysParm->energyIncreasePermit = 0.0;
66 sysParm->fracNeighNixedMax = 0.25;
67 return;
68 }
69
70 void
_pullFlagInit(pullFlag * flag)71 _pullFlagInit(pullFlag *flag) {
72
73 flag->permuteOnRebin = AIR_FALSE;
74 flag->noPopCntlWithZeroAlpha = AIR_FALSE;
75 flag->useBetaForGammaLearn = AIR_FALSE;
76 flag->restrictiveAddToBins = AIR_TRUE;
77 flag->energyFromStrength = AIR_FALSE;
78 flag->nixAtVolumeEdgeSpace = AIR_FALSE;
79 flag->constraintBeforeSeedThresh = AIR_FALSE;
80 flag->noAdd = AIR_FALSE;
81 flag->popCntlEnoughTest = AIR_TRUE; /* really needs to be true by default */
82 flag->binSingle = AIR_FALSE;
83 flag->allowCodimension3Constraints = AIR_FALSE;
84 flag->scaleIsTau = AIR_FALSE;
85 flag->startSkipsPoints = AIR_FALSE; /* must be false by default */
86 return;
87 }
88
89 int
_pullIterParmCheck(pullIterParm * iterParm)90 _pullIterParmCheck(pullIterParm *iterParm) {
91 static const char me[]="_pullIterParmCheck";
92
93 if (!( 1 <= iterParm->constraintMax
94 && iterParm->constraintMax <= 500 )) {
95 biffAddf(PULL, "%s: iterParm->constraintMax %u not in range [%u,%u]",
96 me, iterParm->constraintMax, 1, _PULL_CONSTRAINT_ITER_MAX);
97 return 1;
98 }
99 return 0;
100 }
101
102 int
pullIterParmSet(pullContext * pctx,int which,unsigned int pval)103 pullIterParmSet(pullContext *pctx, int which, unsigned int pval) {
104 static const char me[]="pullIterParmSet";
105
106 if (!pctx) {
107 biffAddf(PULL, "%s: got NULL pointer", me);
108 return 1;
109 }
110 if (!AIR_IN_OP(pullIterParmUnknown, which, pullIterParmLast)) {
111 biffAddf(PULL, "%s: iter parm %d not valid", me, which);
112 return 1;
113 }
114 switch(which) {
115 case pullIterParmMin:
116 pctx->iterParm.min = pval;
117 break;
118 case pullIterParmMax:
119 pctx->iterParm.max = pval;
120 break;
121 case pullIterParmStuckMax:
122 pctx->iterParm.stuckMax = pval;
123 break;
124 case pullIterParmConstraintMax:
125 pctx->iterParm.constraintMax = pval;
126 break;
127 case pullIterParmPopCntlPeriod:
128 pctx->iterParm.popCntlPeriod = pval;
129 break;
130 case pullIterParmAddDescent:
131 pctx->iterParm.addDescent = pval;
132 break;
133 case pullIterParmCallback:
134 pctx->iterParm.callback = pval;
135 break;
136 case pullIterParmSnap:
137 pctx->iterParm.snap = pval;
138 break;
139 case pullIterParmEnergyIncreasePermitHalfLife:
140 pctx->iterParm.energyIncreasePermitHalfLife = pval;
141 if (pval) {
142 pctx->eipScale = pow(0.5, 1.0/pval);
143 } else {
144 pctx->eipScale = 1;
145 }
146 break;
147 default:
148 biffAddf(me, "%s: sorry, iter parm %d valid but not handled?", me, which);
149 return 1;
150 }
151 return 0;
152 }
153
154 #define CHECK(thing, min, max) \
155 if (!( AIR_EXISTS(sysParm->thing) \
156 && min <= sysParm->thing && sysParm->thing <= max )) { \
157 biffAddf(PULL, "%s: sysParm->" #thing " %g not in range [%g,%g]", \
158 me, sysParm->thing, min, max); \
159 return 1; \
160 }
161
162 int
_pullSysParmCheck(pullSysParm * sysParm)163 _pullSysParmCheck(pullSysParm *sysParm) {
164 static const char me[]="_pullSysParmCheck";
165
166 /* these reality-check bounds are somewhat arbitrary */
167 CHECK(alpha, 0.0, 1.0);
168 CHECK(beta, 0.0, 1.0);
169 /* HEY: no check on gamma? */
170 /* no check on theta */
171 CHECK(wall, 0.0, 100.0);
172 CHECK(radiusSpace, 0.000001, 80.0);
173 CHECK(radiusScale, 0.000001, 80.0);
174 CHECK(binWidthSpace, 1.0, 15.0);
175 CHECK(neighborTrueProb, 0.02, 1.0);
176 CHECK(probeProb, 0.02, 1.0);
177 if (!( AIR_EXISTS(sysParm->stepInitial)
178 && sysParm->stepInitial > 0 )) {
179 biffAddf(PULL, "%s: sysParm->stepInitial %g not > 0", me,
180 sysParm->stepInitial);
181 return 1;
182 }
183 CHECK(opporStepScale, 1.0, 5.0);
184 CHECK(backStepScale, 0.01, 0.99);
185 CHECK(constraintStepMin, 0.00000000000000001, 0.1);
186 CHECK(energyDecreaseMin, -0.2, 1.0);
187 CHECK(energyDecreasePopCntlMin, -1.0, 1.0);
188 CHECK(energyIncreasePermit, 0.0, 1.0);
189 CHECK(fracNeighNixedMax, 0.01, 0.99);
190 return 0;
191 }
192 #undef CHECK
193
194 int
pullSysParmSet(pullContext * pctx,int which,double pval)195 pullSysParmSet(pullContext *pctx, int which, double pval) {
196 static const char me[]="pullSysParmSet";
197
198 if (!pctx) {
199 biffAddf(PULL, "%s: got NULL pointer", me);
200 return 1;
201 }
202 if (!AIR_IN_OP(pullSysParmUnknown, which, pullSysParmLast)) {
203 biffAddf(PULL, "%s: sys parm %d not valid", me, which);
204 return 1;
205 }
206 switch(which) {
207 case pullSysParmAlpha:
208 pctx->sysParm.alpha = pval;
209 break;
210 case pullSysParmBeta:
211 pctx->sysParm.beta = pval;
212 break;
213 case pullSysParmGamma:
214 pctx->sysParm.gamma = pval;
215 break;
216 case pullSysParmSeparableGammaLearnRescale:
217 pctx->sysParm.separableGammaLearnRescale = pval;
218 break;
219 case pullSysParmTheta:
220 pctx->sysParm.theta = pval;
221 break;
222 case pullSysParmStepInitial:
223 pctx->sysParm.stepInitial = pval;
224 break;
225 case pullSysParmRadiusSpace:
226 pctx->sysParm.radiusSpace = pval;
227 break;
228 case pullSysParmRadiusScale:
229 pctx->sysParm.radiusScale = pval;
230 break;
231 case pullSysParmBinWidthSpace:
232 pctx->sysParm.binWidthSpace = pval;
233 break;
234 case pullSysParmNeighborTrueProb:
235 pctx->sysParm.neighborTrueProb = pval;
236 break;
237 case pullSysParmProbeProb:
238 pctx->sysParm.probeProb = pval;
239 break;
240 case pullSysParmOpporStepScale:
241 pctx->sysParm.opporStepScale = pval;
242 break;
243 case pullSysParmBackStepScale:
244 pctx->sysParm.backStepScale = pval;
245 break;
246 case pullSysParmEnergyDecreasePopCntlMin:
247 pctx->sysParm.energyDecreasePopCntlMin = pval;
248 break;
249 case pullSysParmEnergyDecreaseMin:
250 pctx->sysParm.energyDecreaseMin = pval;
251 break;
252 case pullSysParmConstraintStepMin:
253 pctx->sysParm.constraintStepMin = pval;
254 break;
255 case pullSysParmEnergyIncreasePermit:
256 pctx->sysParm.energyIncreasePermit = pval;
257 break;
258 case pullSysParmFracNeighNixedMax:
259 pctx->sysParm.fracNeighNixedMax = pval;
260 break;
261 case pullSysParmWall:
262 pctx->sysParm.wall = pval;
263 break;
264 default:
265 biffAddf(me, "%s: sorry, sys parm %d valid but not handled?", me, which);
266 return 1;
267 }
268 return 0;
269 }
270
271 /*
272 ******** pullFlagSet
273 **
274 ** uniform way of setting all the boolean-ish flags
275 */
276 int
pullFlagSet(pullContext * pctx,int which,int flag)277 pullFlagSet(pullContext *pctx, int which, int flag) {
278 static const char me[]="pullFlagSet";
279
280 if (!pctx) {
281 biffAddf(PULL, "%s: got NULL pointer", me);
282 return 1;
283 }
284 if (!AIR_IN_OP(pullFlagUnknown, which, pullFlagLast)) {
285 biffAddf(PULL, "%s: flag %d not valid", me, which);
286 return 1;
287 }
288 switch (which) {
289 case pullFlagPermuteOnRebin:
290 pctx->flag.permuteOnRebin = flag;
291 break;
292 case pullFlagNoPopCntlWithZeroAlpha:
293 pctx->flag.noPopCntlWithZeroAlpha = flag;
294 break;
295 case pullFlagUseBetaForGammaLearn:
296 pctx->flag.useBetaForGammaLearn = flag;
297 break;
298 case pullFlagRestrictiveAddToBins:
299 pctx->flag.restrictiveAddToBins = flag;
300 break;
301 case pullFlagEnergyFromStrength:
302 pctx->flag.energyFromStrength = flag;
303 break;
304 case pullFlagNixAtVolumeEdgeSpace:
305 pctx->flag.nixAtVolumeEdgeSpace = flag;
306 break;
307 case pullFlagConstraintBeforeSeedThresh:
308 pctx->flag.constraintBeforeSeedThresh = flag;
309 break;
310 case pullFlagNoAdd:
311 pctx->flag.noAdd = flag;
312 break;
313 case pullFlagPopCntlEnoughTest:
314 pctx->flag.popCntlEnoughTest = flag;
315 break;
316 case pullFlagBinSingle:
317 pctx->flag.binSingle = flag;
318 break;
319 case pullFlagAllowCodimension3Constraints:
320 pctx->flag.allowCodimension3Constraints = flag;
321 break;
322 case pullFlagScaleIsTau:
323 pctx->flag.scaleIsTau = flag;
324 break;
325 case pullFlagStartSkipsPoints:
326 pctx->flag.startSkipsPoints = flag;
327 break;
328 default:
329 biffAddf(me, "%s: sorry, flag %d valid but not handled?", me, which);
330 return 1;
331 }
332 return 0;
333 }
334
335 /*
336 ** HEY: its really confusing to have the array of per-CONTEXT volumes.
337 ** I know they're there to be copied upon task creation to create the
338 ** per-TASK volumes, but its easy to think that one is supposed to be
339 ** doing something with them, or that changes to them will have some
340 ** effect . . .
341 */
342 int
pullVerboseSet(pullContext * pctx,int verbose)343 pullVerboseSet(pullContext *pctx, int verbose) {
344 static const char me[]="pullVerboseSet";
345 unsigned int volIdx, taskIdx;
346
347 if (!pctx) {
348 biffAddf(PULL, "%s: got NULL pointer", me);
349 return 1;
350 }
351 pctx->verbose = verbose;
352 for (volIdx=0; volIdx<pctx->volNum; volIdx++) {
353 int v;
354 v = verbose > 0 ? verbose - 1 : 0;
355 gageParmSet(pctx->vol[volIdx]->gctx, gageParmVerbose, v);
356 }
357 for (taskIdx=0; taskIdx<pctx->threadNum; taskIdx++) {
358 for (volIdx=0; volIdx<pctx->volNum; volIdx++) {
359 int v;
360 v = verbose > 0 ? verbose - 1 : 0;
361 gageParmSet(pctx->task[taskIdx]->vol[volIdx]->gctx, gageParmVerbose, v);
362 }
363 }
364 return 0;
365 }
366
367 int
pullThreadNumSet(pullContext * pctx,unsigned int threadNum)368 pullThreadNumSet(pullContext *pctx, unsigned int threadNum) {
369 static const char me[]="pullThreadNumSet";
370
371 if (!pctx) {
372 biffAddf(PULL, "%s: got NULL pointer", me);
373 return 1;
374 }
375 pctx->threadNum = threadNum;
376 return 0;
377 }
378
379 int
pullRngSeedSet(pullContext * pctx,unsigned int rngSeed)380 pullRngSeedSet(pullContext *pctx, unsigned int rngSeed) {
381 static const char me[]="pullRngSeedSet";
382
383 if (!pctx) {
384 biffAddf(PULL, "%s: got NULL pointer", me);
385 return 1;
386 }
387 pctx->rngSeed = rngSeed;
388 return 0;
389 }
390
391 int
pullProgressBinModSet(pullContext * pctx,unsigned int bmod)392 pullProgressBinModSet(pullContext *pctx, unsigned int bmod) {
393 static const char me[]="pullProgressBinModSet";
394
395 if (!pctx) {
396 biffAddf(PULL, "%s: got NULL pointer", me);
397 return 1;
398 }
399 pctx->progressBinMod = bmod;
400 return 0;
401 }
402
403 int
pullCallbackSet(pullContext * pctx,void (* iter_cb)(void * data_cb),void * data_cb)404 pullCallbackSet(pullContext *pctx,
405 void (*iter_cb)(void *data_cb),
406 void *data_cb) {
407 static const char me[]="pullCallbackSet";
408
409 if (!pctx) {
410 biffAddf(PULL, "%s: got NULL pointer", me);
411 return 1;
412 }
413 pctx->iter_cb = iter_cb;
414 pctx->data_cb = data_cb;
415 return 0;
416 }
417
418 /*
419 ******** pullInterEnergySet
420 **
421 ** This is the single function for setting the inter-particle energy,
422 ** which is clumsy because which pullEnergySpecs are necessary is
423 ** different depending on interType. Pass NULL for those not needed.
424 **
425 ** Note that all the pullEnergySpecs inside the pctx are created
426 ** by pullContextNew, so they should never be NULL. When a pullEnergySpec
427 ** is not needed for a given interType, we set it to pullEnergyZero
428 ** and a vector of NaNs.
429 */
430 int
pullInterEnergySet(pullContext * pctx,int interType,const pullEnergySpec * enspR,const pullEnergySpec * enspS,const pullEnergySpec * enspWin)431 pullInterEnergySet(pullContext *pctx, int interType,
432 const pullEnergySpec *enspR,
433 const pullEnergySpec *enspS,
434 const pullEnergySpec *enspWin) {
435 static const char me[]="pullInterEnergySet";
436 unsigned int zpi;
437 double zeroParm[PULL_ENERGY_PARM_NUM];
438
439 if (!pctx) {
440 biffAddf(PULL, "%s: got NULL pointer", me);
441 return 1;
442 }
443 if (!AIR_IN_OP(pullInterTypeUnknown, interType, pullInterTypeLast)) {
444 biffAddf(PULL, "%s: interType %d not valid", me, interType);
445 return 1;
446 }
447 for (zpi=0; zpi<PULL_ENERGY_PARM_NUM; zpi++) {
448 zeroParm[zpi] = AIR_NAN;
449 }
450
451 #define CHECK_N_COPY(X) \
452 if (!ensp##X) { \
453 biffAddf(PULL, "%s: need non-NULL ensp" #X " for interType %s", me, \
454 airEnumStr(pullInterType, interType)); \
455 return 1; \
456 } \
457 pullEnergySpecCopy(pctx->energySpec##X, ensp##X)
458
459 switch (interType) {
460 case pullInterTypeJustR:
461 /* 1: phi(r,s) = phi_r(r) */
462 case pullInterTypeUnivariate:
463 /* 2: phi(r,s) = phi_r(u); u = sqrt(r*r+s*s) */
464 CHECK_N_COPY(R);
465 pullEnergySpecSet(pctx->energySpecS, pullEnergyZero, zeroParm);
466 pullEnergySpecSet(pctx->energySpecWin, pullEnergyZero, zeroParm);
467 break;
468 case pullInterTypeSeparable:
469 /* 3: phi(r,s) = phi_r(r)*phi_s(s) */
470 CHECK_N_COPY(R);
471 CHECK_N_COPY(S);
472 pullEnergySpecSet(pctx->energySpecWin, pullEnergyZero, zeroParm);
473 break;
474 case pullInterTypeAdditive:
475 /* 4: phi(r,s) = beta*phi_r(r)*win(s) + (1-beta)*win(r)*phi_s(s) */
476 CHECK_N_COPY(R);
477 CHECK_N_COPY(S);
478 CHECK_N_COPY(Win);
479 break;
480 default:
481 biffAddf(PULL, "%s: interType %d valid but no handled?", me, interType);
482 return 1;
483 }
484 #undef CHECK_N_COPY
485
486 pctx->interType = interType;
487 return 0;
488 }
489
490 /*
491 ** you can pass in a NULL FILE* if you want
492 */
493 int
pullLogAddSet(pullContext * pctx,FILE * flog)494 pullLogAddSet(pullContext *pctx, FILE *flog) {
495 static const char me[]="pullLogAddSet";
496
497 if (!(pctx)) {
498 biffAddf(PULL, "%s: got NULL pointer", me);
499 return 1;
500 }
501
502 pctx->logAdd = flog;
503 return 0;
504 }
505