ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/BioCocoa/BioSwarm/trunk/Models/BioSwarmModel.m
Revision: 350
Committed: Wed Aug 19 16:45:07 2015 UTC (4 years, 9 months ago) by schristley
File size: 34726 byte(s)
Log Message:
MPI timing
Line File contents
1 /*
2 BioSwarmModel.m
3
4 Copyright (C) 2012 Scott Christley
5
6 Author: Scott Christley <schristley@mac.com>
7 Date: April 2012
8
9 This file is part of the BioSwarm Framework.
10
11 This library is free software; you can redistribute it and/or
12 modify it under the terms of the GNU Lesser General Public
13 License as published by the Free Software Foundation; either
14 version 2 of the License, or (at your option) any later version.
15
16 This library is distributed in the hope that it will be useful,
17 but WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 Library General Public License for more details.
20
21 You should have received a copy of the GNU Lesser General Public
22 License along with this library; if not, write to the Free
23 Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
24 Boston, MA 02111 USA.
25 */
26
27 #import "BioSwarmModel.h"
28 #import "BioSwarmController.h"
29 #import "Grid2DArray.h"
30
31 #import "internal.h"
32
33 // encode strings
34 static NSString *idEncode = nil;
35 static NSString *intEncode = nil;
36 static NSString *floatEncode = nil;
37 static NSString *doubleEncode = nil;
38 static NSString *longEncode = nil;
39 static NSString *boolEncode = nil;
40
41 @implementation BioSwarmModel
42
43 - initWithModel: (NSDictionary *)aModel andKey: (NSString *)aKey andController: anObj
44 {
45 return [self initWithModel: aModel andKey: aKey andController: anObj andParent: nil];
46 }
47
48 - initWithModel: (NSDictionary *)aModel andKey: (NSString *)aKey andController: anObj andParent: aParent
49 {
50 id s;
51
52 [super init];
53
54 delegate = nil;
55 modelController = [anObj retain];
56 model = [aModel retain];
57 keyName = [aKey retain];
58 parentModel = aParent;
59 childModels = [NSMutableArray new];
60 spatialContext = [aModel objectForKey: @"spatialContext"];
61 //printf("parent child spatial: %p %p %p\n", parentModel, childModels, spatialContext);
62
63 // setup species
64 species = [aModel objectForKey: @"species"];
65 numOfSpecies = [species count];
66 //printf("%d species\n", numOfSpecies);
67 speciesData = NULL;
68
69 isDebug = NO;
70 s = [aModel objectForKey: @"debug"];
71 if (s) isDebug = [s boolValue];
72
73 isBinaryFile = [modelController isBinaryFile];
74 s = [aModel objectForKey: @"binaryFile"];
75 if (s) isBinaryFile = [s boolValue];
76
77 name = [model objectForKey: @"name"];
78 if (!name) name = [modelController modelName];
79 CHECK_PARAM(name, "name");
80
81 // numerical scheme
82 s = [aModel objectForKey: @"numericalScheme"];
83 //CHECK_PARAM(s, "numericalScheme")
84 numericalScheme = s;
85
86 // Multiplicity is number of systems to manage.
87 // Separate from multiple models which is managed externally.
88 // Useful for models which require multiple copies of submodels.
89 multiplicity = 1;
90
91 dt = 1.0;
92 s = [aModel objectForKey: @"dt"];
93 //CHECK_PARAM(s, "dt")
94 if (s) dt = [s doubleValue];
95
96 timeScale = 1;
97 s = [aModel objectForKey: @"timeScale"];
98 if (s) timeScale = [s doubleValue];
99
100 runNumber = [modelController outputRun];
101 outputSuffix = [modelController outputSuffix];
102
103 paramMin = NULL;
104 paramMax = NULL;
105
106 return self;
107 }
108
109 - newFromModelIndividual
110 {
111 return [[[self class] alloc] initWithModel: model];
112 }
113
114 - (void)dealloc
115 {
116 //printf("BioSwarmModel dealloc\n");
117 if (paramMin) free(paramMin);
118 if (paramMax) free(paramMax);
119 [model release];
120 [modelController release];
121 //[parentModel release];
122 [childModels release];
123 [keyName release];
124 if (speciesData) free(speciesData);
125
126 [super dealloc];
127 }
128
129 - (NSDictionary *)model { return model; }
130 - (NSString *)modelName { return name; }
131 - (id)modelController { return modelController; }
132 - (id)parentModel { return parentModel; }
133 - (NSMutableArray *)childModels { return childModels; }
134 - (void)addChildModel: aModel { [childModels addObject: aModel]; }
135 - (int)modelNumber { return modelNumber; }
136 - (void)setModelNumber: (int)aNum { modelNumber = aNum; }
137
138 - modelWithNumber: (int)aNum
139 {
140 if (modelNumber == aNum) return self;
141
142 int i;
143 for (i = 0; i < [childModels count]; ++i) {
144 id aModel = [[childModels objectAtIndex: i] modelWithNumber: aNum];
145 if (aModel) return aModel;
146 }
147
148 return nil;
149 }
150
151 - (NSString *)keyName { return keyName; }
152
153 - modelWithKey: (NSString *)aKey
154 {
155 if ([keyName isEqualToString: aKey]) return self;
156
157 int i;
158 for (i = 0; i < [childModels count]; ++i) {
159 id aModel = [[childModels objectAtIndex: i] modelWithKey: aKey];
160 if (aModel) return aModel;
161 }
162
163 return nil;
164 }
165
166 - (int)numOfModelObjects { return 1; }
167
168 - (SpatialModel *)spatialModel
169 {
170 // by default, ask our parent, if no parent then ask controller
171 if (parentModel) return [parentModel spatialModel];
172 else return [modelController spatialModel];
173 }
174 - (SpatialModel *)spatialModelForContext: aModel
175 {
176 if (!spatialContext) return [self spatialModel];
177 else {
178 NSUInteger i = [childModels indexOfObject: aModel];
179 if (i == NSNotFound) {
180 printf("ERROR: attempt to get spatial model for invalid context.\n");
181 return nil;
182 }
183 NSString *s = [spatialContext objectAtIndex: i];
184 printf("spatial context: %ld %s\n", (unsigned long)i, [s UTF8String]);
185 if ([s isEqualToString: @"parent"]) return [parentModel spatialModel];
186 else if ([s isEqualToString: @"nil"]) return nil;
187 }
188 return nil;
189 }
190
191 - (id)modelOfTypeClass: (Class)aClass
192 {
193 int i;
194
195 if ([self isKindOfClass: aClass]) return self;
196
197 for (i = 0; i < [childModels count]; ++i) {
198 id m = [[childModels objectAtIndex: i] modelOfTypeClass: aClass];
199 if (m) return m;
200 }
201
202 return nil;
203 }
204
205 - (void)setMultiplicity: (int)aNum { multiplicity = aNum; }
206 - (int)multiplicity { return multiplicity; }
207
208 - delegate { return delegate; }
209 - (void)setDelegate: anObj { delegate = anObj; }
210
211 - (int)numOfSpecies { return numOfSpecies; }
212 - (NSString *)nameOfSpecies: (int)theSpecies { return [species objectAtIndex: theSpecies]; }
213 - (id)valueOfSpecies: (int)theSpecies withMultiplicity: (int)aNum
214 {
215 if (!numOfSpecies) return nil;
216
217 switch (encodeTag) {
218 case 1: {
219 int *sData = speciesData;
220 return [NSNumber numberWithInt: sData[theSpecies]];
221 break;
222 }
223 case 2: {
224 long *sData = speciesData;
225 return [NSNumber numberWithLong: sData[theSpecies]];
226 break;
227 }
228 case 3: {
229 float *sData = speciesData;
230 return [NSNumber numberWithFloat: sData[theSpecies]];
231 break;
232 }
233 case 4: {
234 double *sData = speciesData;
235 return [NSNumber numberWithDouble: sData[theSpecies]];
236 break;
237 }
238 }
239 return nil;
240 }
241
242 - (BOOL)isDebug
243 {
244 // the controller may globally turn on debugging
245 //if (!isDebug) return [modelController isDebug];
246 return isDebug;
247 }
248 - (void)setDebug: (BOOL)aFlag { isDebug = aFlag; }
249
250 - (BOOL)isBinaryFile { return isBinaryFile; }
251 - (void)setBinaryFile: (BOOL)aFlag { isBinaryFile = aFlag; }
252
253 - (void)setRunNumber: (int)aNum { runNumber = aNum; }
254 - (int)runNumber { return runNumber; }
255
256 - (double)dt { return dt; }
257
258 #if 0
259 #pragma mark == MANAGE MODELS AND SIMULATIONS ==
260 #endif
261
262 //
263 // pass messages down the model hierarchy
264 //
265
266 - (int)numModelsInHierarchy
267 {
268 int i, cnt = 1;
269 for (i = 0; i < [childModels count]; ++i)
270 cnt += [[childModels objectAtIndex: i] numModelsInHierarchy];
271 return cnt;
272 }
273
274 static int gcd(int a, int b)
275 {
276 if (b == 0) return a;
277 if (a == b) return a;
278 return gcd(b, a % b);
279 }
280
281 - (double)smallestDTInHierarchy
282 {
283 int i;
284
285 double sdt = dt;
286 for (i = 0; i < [childModels count]; ++i) {
287 double cdt = [[childModels objectAtIndex: i] smallestDTInHierarchy];
288 if (cdt < sdt) sdt = cdt;
289 }
290
291 return sdt;
292 }
293
294 - (int)largestTimeScaleInHierarchy
295 {
296 int i;
297
298 int lts = 0;
299 for (i = 0; i < [childModels count]; ++i) {
300 int ts = [[childModels objectAtIndex: i] largestTimeScaleInHierarchy];
301 if (ts > lts) lts = ts;
302 }
303 if (timeScale > lts) lts = timeScale;
304
305 return lts;
306 }
307
308 - (int)actionTimeScale
309 {
310 int lts = [modelController largestTimeScale];
311 int ats = gcd(lts, timeScale);
312 printf("gcd(%d,%d) = %d, %d\n", lts, timeScale, ats, lts / ats);
313 return lts / ats;
314 }
315
316 - (void)performHierarchyInvocation: (NSInvocation *)anInvocation
317 {
318 int i;
319 if ([self respondsToSelector: [anInvocation selector]]) [anInvocation invokeWithTarget: self];
320 for (i = 0; i < [childModels count]; ++i)
321 [[childModels objectAtIndex: i] performHierarchyInvocation: anInvocation];
322 }
323
324 - (void)setupHierarchyDataFilesWithState: (int)aState
325 {
326 int i;
327 [self setupDataFilesWithState: aState];
328 for (i = 0; i < [childModels count]; ++i)
329 [[childModels objectAtIndex: i] setupHierarchyDataFilesWithState: aState];
330 }
331
332 - (void)allocateHierarchyWithEncode:(NSString *)anEncode
333 {
334 int i;
335
336 if ([anEncode isEqual: [BioSwarmModel intEncode]]) {
337 encodeTag = BCENCODE_INT;
338 } else if ([anEncode isEqual: [BioSwarmModel longEncode]]) {
339 encodeTag = BCENCODE_LONG;
340 } else if ([anEncode isEqual: [BioSwarmModel floatEncode]]) {
341 encodeTag = BCENCODE_FLOAT;
342 } else if ([anEncode isEqual: [BioSwarmModel doubleEncode]]) {
343 encodeTag = BCENCODE_DOUBLE;
344 } else {
345 // throw exception or something
346 NSLog(@"ERROR: BioSwarmModel unsupported encoding %@\n", anEncode);
347 return;
348 }
349 dataEncode = [anEncode retain];
350
351 if (numOfSpecies) {
352 switch (encodeTag) {
353 case BCENCODE_INT:
354 speciesData = malloc(sizeof(int) * numOfSpecies);
355 for (i = 0; i < numOfSpecies; ++i) ((int *)speciesData)[i] = 0;
356 break;
357 case BCENCODE_LONG:
358 speciesData = malloc(sizeof(long) * numOfSpecies);
359 for (i = 0; i < numOfSpecies; ++i) ((long *)speciesData)[i] = 0;
360 break;
361 case BCENCODE_FLOAT:
362 speciesData = malloc(sizeof(float) * numOfSpecies);
363 for (i = 0; i < numOfSpecies; ++i) ((float *)speciesData)[i] = 0;
364 break;
365 case BCENCODE_DOUBLE:
366 speciesData = malloc(sizeof(double) * numOfSpecies);
367 for (i = 0; i < numOfSpecies; ++i) ((double *)speciesData)[i] = 0;
368 break;
369 }
370 }
371
372 [self allocateSimulationWithEncode: anEncode];
373 for (i = 0; i < [childModels count]; ++i)
374 [[childModels objectAtIndex: i] allocateHierarchyWithEncode: anEncode];
375 }
376
377 - (void)initializeHierarchy
378 {
379 int i;
380 [self initializeSimulation];
381 for (i = 0; i < [childModels count]; ++i)
382 [[childModels objectAtIndex: i] initializeHierarchy];
383 }
384
385 - (BOOL)writeHierarchyDataWithCheck: (BOOL)aFlag
386 {
387 int i;
388 BOOL res = YES;
389 res &= [self writeDataWithCheck: aFlag];
390 for (i = 0; i < [childModels count]; ++i)
391 res &= [[childModels objectAtIndex: i] writeHierarchyDataWithCheck: aFlag];
392 return res;
393 }
394
395 - (BOOL)readHierarchyData
396 {
397 int i;
398 BOOL res = YES;
399 res &= [self readData];
400 for (i = 0; i < [childModels count]; ++i)
401 res &= [[childModels objectAtIndex: i] readHierarchyData];
402 return res;
403 }
404
405 - (void)rewindHierarchyFiles
406 {
407 int i;
408 [self rewindFiles];
409 for (i = 0; i < [childModels count]; ++i)
410 [[childModels objectAtIndex: i] rewindHierarchyFiles];
411 }
412
413 - (void)cleanupHierarchy
414 {
415 int i;
416 [self cleanupSimulation];
417 for (i = 0; i < [childModels count]; ++i)
418 [[childModels objectAtIndex: i] cleanupHierarchy];
419 }
420
421 //
422 // These should be implemented by subclasses
423 //
424
425 - (void)setupDataFilesWithState: (int)aState
426 {
427 }
428
429 - (void)allocateSimulationWithEncode: (NSString *)anEncode
430 {
431 }
432
433 - (void)initializeSimulation
434 {
435 }
436
437 - (void)initializeFileCompleted
438 {
439 }
440
441 - (BOOL)writeDataWithCheck: (BOOL)aFlag
442 {
443 return YES;
444 }
445
446 - (BOOL)readData
447 {
448 return YES;
449 }
450
451 - (void)rewindFiles
452 {
453 }
454
455 - (void)stepSimulation
456 {
457 }
458
459 - (void)cleanupSimulation
460 {
461 }
462
463 //
464 // Parameter operations
465 //
466 - (int)numOfParameters { return numOfParameters; }
467 - (NSString *)nameOfParameter: (int)theParam { return [parameterNames objectAtIndex: theParam]; }
468
469 - (BOOL)hierarchyParameterList:(NSArray *)includeList exclude:(NSArray *)excludeList
470 names:(NSMutableArray *)pNames indexes:(NSMutableArray *)pIndexes
471 modelNumbers:(NSMutableArray *)pModelNumbers
472 modelKeys:(NSMutableArray *)pModelKeys
473 {
474 int i, k;
475 int numParams = [self numOfParameters];
476 BOOL flag = YES;
477
478 if (numParams) {
479 for (k = 0; k < numParams; ++k) {
480 NSString *pn = [self nameOfParameter: k];
481 if (includeList && ([includeList indexOfObject: pn] == NSNotFound)) continue;
482 if (excludeList && ([excludeList indexOfObject: pn] != NSNotFound)) continue;
483
484 // save parameter info
485 if (pNames) [pNames addObject: pn];
486 if (pIndexes) [pIndexes addObject: [NSNumber numberWithInt: k]];
487 if (pModelNumbers) [pModelNumbers addObject: [NSNumber numberWithInt: modelNumber]];
488 if (pModelKeys) [pModelKeys addObject: keyName];
489 }
490
491 // make sure parameter ranges are setup
492 //if (![self setupParameterRanges: nil checkNames: parameterNames]) flag = NO;
493 }
494
495 for (i = 0; i < [childModels count]; ++i)
496 flag &= [[childModels objectAtIndex: i] hierarchyParameterList: includeList
497 exclude: excludeList
498 names: pNames
499 indexes: pIndexes
500 modelNumbers: pModelNumbers
501 modelKeys: pModelKeys];
502
503 return flag;
504 }
505
506 - (BOOL)setupParameterRanges: (NSDictionary *)ranges doCheck: (BOOL)aFlag
507 {
508 id paramRanges;
509 BOOL good = YES;
510 int i, numParams = [self numOfParameters];
511
512 if (numParams) {
513 if (!ranges) {
514 paramRanges = [model objectForKey: @"parameterRanges"];
515 CHECK_PARAMF(paramRanges, "parameterRanges", good);
516 } else paramRanges = ranges;
517
518 if (paramMin) free(paramMin);
519 if (paramMax) free(paramMax);
520 paramMin = malloc(sizeof(double) * numParams);
521 paramMax = malloc(sizeof(double) * numParams);
522 BOOL flag = YES;
523 for (i = 0; i < numParams; ++i) {
524 NSString *s = [self nameOfParameter: i];
525
526 NSString *paramName = [NSString stringWithFormat: @"%@_min", s];
527 id param = [paramRanges objectForKey: paramName];
528 if (aFlag) CHECK_PARAMF(param, [paramName UTF8String], flag);
529 paramMin[i] = [param doubleValue];
530
531 paramName = [NSString stringWithFormat: @"%@_max", s];
532 param = [paramRanges objectForKey: paramName];
533 if (aFlag) CHECK_PARAMF(param, [paramName UTF8String], flag);
534 paramMax[i] = [param doubleValue];
535 }
536 good &= flag;
537 }
538
539 return good;
540 }
541
542 - (BOOL)setupParameterRanges: (NSDictionary *)ranges checkNames: (NSArray *)names
543 {
544 id paramRanges;
545 BOOL good = YES;
546 int i, numParams = [self numOfParameters];
547
548 if (numParams) {
549 if (!ranges) {
550 paramRanges = [model objectForKey: @"parameterRanges"];
551 CHECK_PARAMF(paramRanges, "parameterRanges", good);
552 } else paramRanges = ranges;
553
554 if (paramMin) free(paramMin);
555 if (paramMax) free(paramMax);
556 paramMin = malloc(sizeof(double) * numParams);
557 paramMax = malloc(sizeof(double) * numParams);
558 BOOL flag = YES;
559 for (i = 0; i < numParams; ++i) {
560 NSString *s = [self nameOfParameter: i];
561
562 NSString *paramName = [NSString stringWithFormat: @"%@_min", s];
563 id param = [paramRanges objectForKey: paramName];
564 if (names) {
565 NSUInteger idx = [names indexOfObject: s];
566 if (idx != NSNotFound) CHECK_PARAMF(param, [paramName UTF8String], flag);
567 }
568 paramMin[i] = [param doubleValue];
569
570 paramName = [NSString stringWithFormat: @"%@_max", s];
571 param = [paramRanges objectForKey: paramName];
572 if (names) {
573 NSUInteger idx = [names indexOfObject: s];
574 if (idx != NSNotFound) CHECK_PARAMF(param, [paramName UTF8String], flag);
575 }
576 paramMax[i] = [param doubleValue];
577 }
578 good &= flag;
579 }
580
581 for (i = 0; i < [childModels count]; ++i)
582 good &= [[childModels objectAtIndex: i] setupParameterRanges: ranges checkNames: names];
583
584 return good;
585 }
586
587 - (double *)parameterMinRange { return paramMin; }
588 - (double *)parameterMaxRange { return paramMax; }
589
590 - (double)randomValueforParameter: (int)aParam
591 {
592 if ((aParam < 0) || (aParam >= [self numOfParameters])) {
593 printf("ERROR: Invalid parameter index: %d\n", aParam);
594 return 0.0;
595 }
596
597 // check if parameter ranges are setup
598 if (!paramMin)
599 if (![self setupParameterRanges: nil doCheck: NO]) {
600 printf("ERROR: Could not setup parameter ranges.\n");
601 return 0.0;
602 }
603 if (paramMin[aParam] == paramMax[aParam]) {
604 printf("ERROR: Invalid parameter range defined for %s.\n", [[self nameOfParameter: aParam] UTF8String]);
605 return 0.0;
606 }
607
608 // determine scales of range
609 int minScale = determine_scale_double(paramMin[aParam]);
610 int maxScale = determine_scale_double(paramMax[aParam]);
611
612 int signCheck = 0;
613 if ((paramMin[aParam] < 0) && (paramMax[aParam] < 0)) signCheck = 1;
614 if (((paramMin[aParam] < 0) && (paramMax[aParam] >= 0))
615 || ((paramMin[aParam] >= 0) && (paramMax[aParam] < 0))) signCheck = 2;
616
617 double newValue = 0.0;
618 if (minScale == maxScale) {
619 double delta = paramMax[aParam] - paramMin[aParam];
620 newValue = RAND_NUM * delta + paramMin[aParam];
621 } else {
622 do {
623 double deltaScale = maxScale - minScale;
624 double newScale = RAND_NUM * deltaScale + minScale;
625 newScale = round(newScale);
626 //printf("newScale: %lf %d %d %lf\n", deltaScale, minScale, maxScale, newScale);
627 //double minVal = paramMin[aParam] / pow(10.0, minScale);
628 //double maxVal = paramMax[aParam] / pow(10.0, maxScale);
629 //double delta = maxVal - minVal;
630 //newValue = RAND_NUM * delta + minVal;
631 newValue = RAND_NUM;
632 if (signCheck == 1) newValue = -newValue;
633 if ((signCheck == 2) && (RAND_NUM < 0.5)) newValue = -newValue;
634 //printf("newVal: %lf %lf %lf %lf\n", delta, minVal, maxVal, newValue);
635 newValue = newValue * pow(10.0, newScale);
636 } while ((newValue < paramMin[aParam]) || (newValue > paramMax[aParam]));
637 }
638
639 return newValue;
640 }
641
642 - (double)perturbValueForParameter: (int)aParam withinRange: (BOOL)aFlag standardDeviation: (double)SD
643 {
644 if ((aParam < 0) || (aParam >= [self numOfParameters])) {
645 printf("ERROR: Invalid parameter index: %d\n", aParam);
646 return 0.0;
647 }
648
649 double origValue = [[model objectForKey: [self nameOfParameter: aParam]] doubleValue];
650 //NSLog(@"origValue = %lf", origValue);
651 double newValue = origValue;
652
653 if (aFlag) {
654 // check if parameter ranges are setup
655 if (!paramMin)
656 if (![self setupParameterRanges: nil doCheck: NO]) {
657 printf("ERROR: Could not setup parameter ranges.\n");
658 return 0.0;
659 }
660
661 do {
662 newValue = origValue + normal_dist_rand(0, SD);
663 } while ((newValue < paramMin[aParam]) || (newValue > paramMax[aParam]));
664 } else {
665 // don't worry about range, just keep same sign
666 BOOL isPositive = YES;
667 if (origValue < 0) isPositive = NO;
668
669 do {
670 newValue = origValue + normal_dist_rand(0, SD);
671 } while (((newValue < 0) && isPositive) || ((newValue >= 0) && !isPositive));
672 }
673
674 return newValue;
675 }
676
677 - (int)modifyRandomParameter: (NSMutableDictionary *)params
678 {
679 // check if parameter ranges are setup
680 if (!paramMin)
681 if (![self setupParameterRanges: nil doCheck: NO]) {
682 printf("ERROR: Could not setup parameter ranges.\n");
683 return -1;
684 }
685
686 // how many parameters
687 int numParams = [self numOfParameters];
688
689 int rp = (int)(RAND_NUM * numParams);
690 NSString *ps = [self nameOfParameter: rp];
691
692 if (!ps) printf("ERROR: Could not randomly pick parameter: %d %d\n", numParams, rp);
693 else {
694 double delta = paramMax[rp] - paramMin[rp];
695 double pv = [[params objectForKey: ps] doubleValue];
696 double newpv = RAND_NUM * delta + paramMin[rp];
697
698 #if 0
699 double scale = determine_scale(pv);
700 BOOL done = NO;
701 while (!done) {
702 newpv = normal_dist_rand(pv, scale);
703 // maintain the same sign
704 if ((pv < 0) && (newpv < 0)) done = YES;
705 if ((pv >= 0) && (newpv >= 0)) done = YES;
706 }
707 #endif
708 #if 0
709 if (fabs(pv) < 0.001) {
710 if (pv >= 0) newpv = RAND_NUM;
711 else newpv = - RAND_NUM;
712 } else {
713 if (RAND_NUM < 0.5) newpv = pv + (RAND_NUM * pv);
714 else newpv = pv - (RAND_NUM * pv);
715 }
716 #endif
717
718 [params setObject: [NSNumber numberWithDouble: newpv] forKey: ps];
719 printf("change parameter %s(%d): %lf -> %lf\n", [ps UTF8String], rp, pv, newpv);
720 }
721
722 return rp;
723 }
724
725 // species operations
726 - (BOOL)hierarchySpeciesList:(NSMutableArray *)sNames indexes:(NSMutableArray *)sIndexes
727 modelNumbers:(NSMutableArray *)sModelNumbers
728 modelKeys:(NSMutableArray *)sModelKeys;
729 {
730 int i, k;
731 int numSpecies = [self numOfSpecies];
732 BOOL flag = YES;
733
734 if (numSpecies) {
735 for (k = 0; k < numSpecies; ++k) {
736 NSString *pn = [self nameOfSpecies: k];
737
738 // save species info
739 if (sNames) [sNames addObject: pn];
740 if (sIndexes) [sIndexes addObject: [NSNumber numberWithInt: k]];
741 if (sModelNumbers) [sModelNumbers addObject: [NSNumber numberWithInt: modelNumber]];
742 if (sModelKeys) [sModelKeys addObject: keyName];
743 }
744 }
745
746 for (i = 0; i < [childModels count]; ++i)
747 flag &= [[childModels objectAtIndex: i] hierarchySpeciesList: sNames
748 indexes: sIndexes
749 modelNumbers: sModelNumbers
750 modelKeys: sModelKeys];
751
752 return flag;
753 }
754
755 // perturbation operations
756 - (BOOL)setValueForAllSpecies: (double)aValue { return NO; }
757 - (BOOL)setValue: (double)aValue forSpecies: (NSString *)speciesName
758 {
759 if (!numOfSpecies) return NO;
760
761 NSUInteger idx = [species indexOfObject: speciesName];
762 if (idx == NSNotFound) return NO;
763
764 switch (encodeTag) {
765 case BCENCODE_INT: {
766 int *sData = speciesData;
767 sData[idx] = (int)aValue;
768 break;
769 }
770 case BCENCODE_LONG: {
771 long *sData = speciesData;
772 sData[idx] = (long)aValue;
773 break;
774 }
775 case BCENCODE_FLOAT: {
776 float *sData = speciesData;
777 sData[idx] = (float)aValue;
778 break;
779 }
780 case BCENCODE_DOUBLE: {
781 double *sData = speciesData;
782 sData[idx] = aValue;
783 break;
784 }
785 }
786
787 return YES;
788 }
789
790 - (BOOL)adjustValue: (double)aValue forSpecies: (NSString *)speciesName { return NO; }
791 - (BOOL)knockoutSpecies: (NSString *)speciesName { return NO; }
792 - (BOOL)setValue: (double)aValue forParameter: (NSString *)pName { return NO; }
793
794 // Singleton encode strings
795 + (NSString *)idEncode
796 {
797 if (!idEncode) idEncode = [[NSString alloc] initWithUTF8String: @encode(id)];
798 return idEncode;
799 }
800
801 + (NSString *)intEncode
802 {
803 if (!intEncode) intEncode = [[NSString alloc] initWithUTF8String: @encode(int)];
804 return intEncode;
805 }
806
807 + (NSString *)longEncode
808 {
809 if (!longEncode) longEncode = [[NSString alloc] initWithUTF8String: @encode(long)];
810 return longEncode;
811 }
812
813 + (NSString *)floatEncode
814 {
815 if (!floatEncode) floatEncode = [[NSString alloc] initWithUTF8String: @encode(float)];
816 return floatEncode;
817 }
818
819 + (NSString *)doubleEncode
820 {
821 if (!doubleEncode) doubleEncode = [[NSString alloc] initWithUTF8String: @encode(double)];
822 return doubleEncode;
823 }
824
825 + (NSString *)boolEncode
826 {
827 if (!boolEncode) boolEncode = [[NSString alloc] initWithUTF8String: @encode(BOOL)];
828 return boolEncode;
829 }
830
831 @end
832
833 //
834 // GPU
835 //
836
837 @implementation BioSwarmModel (GPU)
838 - (NSMutableDictionary *)controllerGPUCode: (NSString *)prefixName { return nil; }
839 - (NSMutableString *)definitionsGPUCode: (NSString *)prefixName { return nil; }
840 - (NSMutableString *)kernelGPUCode: (NSString *)prefixName { return nil; }
841 - (NSMutableString *)cleanupGPUCode: (NSString *)prefixName { return nil; }
842
843 @end
844
845 //
846 // Run GPU
847 //
848
849 @implementation BioSwarmModel (GPURun)
850
851 - (void)allocateHierarchyGPUData: (ModelGPUData **)gpuData modelInstances: (int)numInstances
852 {
853 // allocate GPU data for ourself
854 getGPUFunctions *func = [modelController gpuFunctions];
855 void *gpuFunctions = (func)(modelNumber, 0);
856 gpuData[modelNumber] = [self allocGPUData: numInstances withGPUFunctions: gpuFunctions];
857 [self allocGPUKernel: gpuData[modelNumber]];
858
859 // propogate message down model hierarchy
860 int i;
861 for (i = 0; i < [childModels count]; ++i)
862 [[childModels objectAtIndex: i] allocateHierarchyGPUData: gpuData modelInstances: numInstances];
863 }
864
865 - (void)freeHierarchyGPUData: (ModelGPUData **)gpuData
866 {
867 // free GPU data for ourself
868 [self freeGPUData: gpuData[modelNumber]];
869
870 // propogate message down model hierarchy
871 int i;
872 for (i = 0; i < [childModels count]; ++i)
873 [[childModels objectAtIndex: i] freeHierarchyGPUData: gpuData];
874 }
875
876 - (void)assignHierarchyDataGPUData: (ModelGPUData **)gpuData ofInstance: (int)aNum toGPU: (BOOL)aFlag
877 {
878 // assign GPU data for ourself
879 [self assignData: gpuData[modelNumber] ofNumber: aNum toGPU: aFlag];
880
881 // propogate message down model hierarchy
882 int i;
883 for (i = 0; i < [childModels count]; ++i)
884 [[childModels objectAtIndex: i] assignHierarchyDataGPUData: gpuData ofInstance: aNum toGPU: aFlag];
885 }
886
887 - (void)assignHierarchyParametersGPUData: (ModelGPUData **)gpuData ofInstance: (int)aNum toGPU: (BOOL)aFlag
888 {
889 // assign GPU data for ourself
890 [self assignParameters: gpuData[modelNumber] ofNumber: aNum toGPU: aFlag];
891
892 // propogate message down model hierarchy
893 int i;
894 for (i = 0; i < [childModels count]; ++i)
895 [[childModels objectAtIndex: i] assignHierarchyParametersGPUData: gpuData ofInstance: aNum toGPU: aFlag];
896 }
897
898 - (void)transferHierarchyGPUData: (ModelGPUData **)gpuData toGPU: (BOOL)aFlag
899 {
900 // transfer GPU data for ourself
901 [self transferGPUKernel: gpuData[modelNumber] toGPU: aFlag];
902
903 // propogate message down model hierarchy
904 int i;
905 for (i = 0; i < [childModels count]; ++i)
906 [[childModels objectAtIndex: i] transferHierarchyGPUData: gpuData toGPU: aFlag];
907 }
908
909 - (void)invokeHierarchyGPUData: (ModelGPUData **)gpuData currentTime: (double)currentTime endTime: (double)endTime timeScale: (int)currentTimeScale
910 {
911 if (currentTimeScale == timeScale) {
912 // invoke GPU kernel for ourself
913 [self invokeGPUKernel: gpuData[modelNumber] currentTime: currentTime endTime: endTime];
914 }
915
916 // propogate message down model hierarchy
917 int i;
918 for (i = 0; i < [childModels count]; ++i)
919 [[childModels objectAtIndex: i] invokeHierarchyGPUData: gpuData currentTime: currentTime endTime: endTime timeScale: currentTimeScale];
920 }
921
922 - (void)invokeHierarchyGPUData: (ModelGPUData **)gpuData currentTime: (double)currentTime endTime: (double)endTime
923 {
924 // Need to break up the total simulation time into smaller time chunks
925 // appropriate for the submodels
926 double deltaTime = [modelController smallestDT];
927 double newTime = currentTime;
928 double newEndTime = newTime + deltaTime;
929
930 BOOL done = NO;
931 while (!done) {
932 // TODO: need to handle the different time scales for each model
933 // right now assume each is on same time scale
934 [self invokeHierarchyGPUData: gpuData currentTime: newTime endTime: newEndTime timeScale: timeScale];
935
936 newTime += deltaTime;
937 newEndTime = newTime + deltaTime;
938 if (TIME_EQUAL(newTime, endTime) || (newTime >= endTime)) done = YES;
939 }
940 }
941
942 - (void)releaseHierarchyGPUData: (ModelGPUData **)gpuData
943 {
944 // call release kernel for ourself
945 [self releaseGPUKernel: gpuData[modelNumber]];
946
947 // propogate message down model hierarchy
948 int i;
949 for (i = 0; i < [childModels count]; ++i)
950 [[childModels objectAtIndex: i] releaseHierarchyGPUData: gpuData];
951 }
952
953 - (void *)allocGPUData: (int)numModels withGPUFunctions: (void *)gpuFunctions { return NULL; }
954 - (void)freeGPUData: (void *)data { }
955 - (void)assignData: (void *)data ofNumber: (int)aNum toGPU: (BOOL)aFlag { }
956 - (int)assignParameters: (void *)data ofNumber: (int)aNum toGPU: (BOOL)aFlag { return 0; }
957 - (void)allocGPUKernel: (void *)data { }
958 - (void)transferGPUKernel: (void *)data toGPU: (BOOL)aFlag { }
959 - (void)invokeGPUKernel: (void *)data currentTime: (double)currentTime endTime: (double)endTime { }
960 - (void)releaseGPUKernel: (void *)data { }
961
962 @end
963
964 //
965 // Visualization
966 //
967
968 @implementation BioSwarmModel (Povray)
969
970 - (void)povrayHierarchyCode: (NSMutableString *)povrayCode
971 {
972 int i;
973 NSString *povray = [self povrayCode];
974 if (povray) [povrayCode appendString: povray];
975
976 for (i = 0; i < [childModels count]; ++i)
977 [[childModels objectAtIndex: i] povrayHierarchyCode: povrayCode];
978 }
979
980 - (NSMutableString *)povrayCode
981 {
982 return nil;
983 }
984
985 @end
986
987 //
988 // GPU and MPI
989 //
990
991 @implementation BioSwarmModel (GPUMPI)
992 - (NSMutableDictionary *)controllerGPUMPICode: (NSString *)prefixName { return [self controllerGPUCode: prefixName]; }
993 - (NSMutableString *)definitionsGPUMPICode: (NSString *)prefixName { return [self definitionsGPUCode: prefixName]; }
994 - (NSMutableString *)kernelGPUMPICode: (NSString *)prefixName { return [self kernelGPUCode: prefixName]; }
995 - (NSMutableString *)cleanupGPUMPICode: (NSString *)prefixName { return [self cleanupGPUCode: prefixName]; }
996 @end
997
998 @implementation BioSwarmModel (GPUMPIRun)
999
1000 - (void)allocateHierarchyGPUMPIData: (ModelGPUData **)gpuData modelInstances: (int)numInstances
1001 {
1002 // we assume GPU data already allocated
1003 [self allocGPUMPIKernel: gpuData[modelNumber]];
1004
1005 // propogate message down model hierarchy
1006 int i;
1007 for (i = 0; i < [childModels count]; ++i)
1008 [[childModels objectAtIndex: i] allocateHierarchyGPUMPIData: gpuData modelInstances: numInstances];
1009 }
1010
1011 - (void)freeHierarchyGPUMPIData: (ModelGPUData **)gpuData
1012 {
1013 // free GPU data for ourself
1014 [self freeGPUMPIData: gpuData[modelNumber]];
1015
1016 // propogate message down model hierarchy
1017 int i;
1018 for (i = 0; i < [childModels count]; ++i)
1019 [[childModels objectAtIndex: i] freeHierarchyGPUMPIData: gpuData];
1020 }
1021
1022 - (void)transferHierarchyGPUMPIData: (ModelGPUData **)gpuData toGPU: (BOOL)aFlag
1023 {
1024 // transfer GPU data for ourself
1025 [self transferGPUMPIKernel: gpuData[modelNumber] toGPU: aFlag];
1026
1027 // propogate message down model hierarchy
1028 int i;
1029 for (i = 0; i < [childModels count]; ++i)
1030 [[childModels objectAtIndex: i] transferHierarchyGPUMPIData: gpuData toGPU: aFlag];
1031 }
1032
1033 - (void)releaseHierarchyGPUMPIData: (ModelGPUData **)gpuData
1034 {
1035 // call release kernel for ourself
1036 [self releaseGPUMPIKernel: gpuData[modelNumber]];
1037
1038 // propogate message down model hierarchy
1039 int i;
1040 for (i = 0; i < [childModels count]; ++i)
1041 [[childModels objectAtIndex: i] releaseHierarchyGPUMPIData: gpuData];
1042 }
1043
1044 - (void)allocGPUMPIKernel: (void *)data { return; }
1045 - (void)transferGPUMPIKernel: (void *)data toGPU: (BOOL)aFlag { return; }
1046 - (void)invokeGPUMPIKernel: (void *)data currentTime: (double)currentTime endTime: (double)endTime { return; }
1047 - (void)releaseGPUMPIKernel: (void *)data { return; }
1048
1049 - (void)initiateGPUMPIData: (void *)data mode: (int)aMode { return; }
1050 - (void)finalizeGPUMPIData: (void *)data mode: (int)aMode { return; }
1051 @end
1052
1053 //
1054 // Common functions
1055 //
1056
1057 // Generate normally distributed random number
1058 double normal_dist_rand(double mean, double sd)
1059 {
1060 double fac, radius, v1, v2;
1061 double rd1Value, rd2Value;
1062
1063 do {
1064 rd1Value = RAND_NUM;
1065 rd2Value = RAND_NUM;
1066 v1 = (2.0 * rd1Value) - 1.0;
1067 v2 = (2.0 * rd2Value) - 1.0;
1068 radius = v1*v1 + v2*v2;
1069 } while (radius >= 1.0);
1070 fac = sqrt (-2.0 * log (radius) / radius);
1071 return v2 * fac * sd + mean; // use fixed params
1072 }
1073
1074 int determine_scale_double(double value)
1075 {
1076 double v;
1077 int scale = 0;
1078
1079 if (value < 0) v = -value;
1080 else v = value;
1081
1082 if (v == 0.0) return scale;
1083
1084 if (v >= 10) {
1085 while (v >= 10) {
1086 v = v / 10;
1087 ++scale;
1088 }
1089 } else {
1090 while (v < 1) {
1091 v = v * 10;
1092 --scale;
1093 }
1094 }
1095
1096 //NSLog(@"%lf %d\n", value, scale);
1097 return scale;
1098 }
1099
1100 int determine_scale_float(float value)
1101 {
1102 float v;
1103 int scale = 0;
1104
1105 if (value < 0) v = -value;
1106 else v = value;
1107
1108 if (v == 0.0) return scale;
1109
1110 if (v >= 10) {
1111 while (v >= 10) {
1112 v = v / 10;
1113 ++scale;
1114 }
1115 } else {
1116 while (v < 1) {
1117 v = v * 10;
1118 --scale;
1119 }
1120 }
1121
1122 return scale;
1123 }
1124
1125 // which way to round error correction based upon direction
1126 // and the coordinate dimension
1127 BOOL direction_round(int direction, int dim)
1128 {
1129 BOOL roundUp = NO; // round down for default
1130
1131 switch (direction) {
1132 case S_LEFT:
1133 roundUp = NO;
1134 break;
1135 case S_RIGHT:
1136 if (dim == 0) roundUp = YES;
1137 break;
1138 case S_UP:
1139 if (dim == 1) roundUp = YES;
1140 break;
1141 case S_DOWN:
1142 roundUp = NO;
1143 break;
1144 case S_LEFT_UP:
1145 if (dim == 1) roundUp = YES;
1146 break;
1147 case S_RIGHT_UP:
1148 if (dim == 0) roundUp = YES;
1149 if (dim == 1) roundUp = YES;
1150 break;
1151 case S_LEFT_DOWN:
1152 roundUp = NO;
1153 break;
1154 case S_RIGHT_DOWN:
1155 if (dim == 0) roundUp = YES;
1156 break;
1157 case S_FRONT:
1158 roundUp = NO;
1159 break;
1160 case S_BACK:
1161 if (dim == 2) roundUp = YES;
1162 break;
1163 case S_LEFT_FRONT:
1164 roundUp = NO;
1165 break;
1166 case S_RIGHT_FRONT:
1167 if (dim == 0) roundUp = YES;
1168 break;
1169 case S_LEFT_BACK:
1170 if (dim == 2) roundUp = YES;
1171 break;
1172 case S_RIGHT_BACK:
1173 if (dim == 0) roundUp = YES;
1174 if (dim == 2) roundUp = YES;
1175 break;
1176 case S_UP_FRONT:
1177 if (dim == 1) roundUp = YES;
1178 break;
1179 case S_DOWN_FRONT:
1180 roundUp = NO;
1181 break;
1182 case S_UP_BACK:
1183 if (dim == 1) roundUp = YES;
1184 if (dim == 2) roundUp = YES;
1185 break;
1186 case S_DOWN_BACK:
1187 if (dim == 2) roundUp = YES;
1188 break;
1189 case S_LEFT_UP_FRONT:
1190 if (dim == 1) roundUp = YES;
1191 break;
1192 case S_RIGHT_UP_FRONT:
1193 if (dim == 0) roundUp = YES;
1194 if (dim == 1) roundUp = YES;
1195 break;
1196 case S_LEFT_DOWN_FRONT:
1197 roundUp = NO;
1198 break;
1199 case S_RIGHT_DOWN_FRONT:
1200 if (dim == 0) roundUp = YES;
1201 break;
1202 case S_LEFT_UP_BACK:
1203 if (dim == 1) roundUp = YES;
1204 if (dim == 2) roundUp = YES;
1205 break;
1206 case S_RIGHT_UP_BACK:
1207 if (dim == 0) roundUp = YES;
1208 if (dim == 1) roundUp = YES;
1209 if (dim == 2) roundUp = YES;
1210 break;
1211 case S_LEFT_DOWN_BACK:
1212 if (dim == 2) roundUp = YES;
1213 break;
1214 case S_RIGHT_DOWN_BACK:
1215 if (dim == 0) roundUp = YES;
1216 if (dim == 2) roundUp = YES;
1217 break;
1218 }
1219
1220 return roundUp;
1221 }
1222
1223 //
1224 // Simple profiling function
1225 //
1226
1227 static time_t bc_profile_time = 0;
1228
1229 void bc_record_time(const char *msg)
1230 {
1231 time_t delta = 0;
1232 time_t current = time(NULL);
1233
1234 if (!bc_profile_time) {
1235 bc_profile_time = time(NULL);
1236 current = bc_profile_time;
1237 } else {
1238 delta = current - bc_profile_time;
1239 bc_profile_time = current;
1240 }
1241
1242 printf("TIME: (%d secs elapsed) %s, %s", delta, msg, asctime(localtime(&bc_profile_time)));
1243 }