file:/local_home/local_home/hugo/neurospaces_project/heccer/source/c/snapshots/0/heccer.c (Mon Jul 14 10:50:33 2008
) HOME
1: static char *pcVersionTime="(08/07/13) Sunday, July 13, 2008 18:47:45 hugo";
2:
3: //
4: // Heccer : a compartmental solver that implements efficient Crank-Nicolson
5: // integration for neuronal models.
6: //
7:
8: //////////////////////////////////////////////////////////////////////////////
9: //'
10: //' Heccer : testbed C implementation
11: //'
12: //' Copyright (C) 2006-2008 Hugo Cornelis
13: //'
14: //' functional ideas .. Hugo Cornelis, hugo.cornelis@gmail.com
15: //'
16: //' coding ............ Hugo Cornelis, hugo.cornelis@gmail.com
17: //'
18: //////////////////////////////////////////////////////////////////////////////
19:
20:
21: #include <stdarg.h>
22: #include <stdlib.h>
23: #include <stdio.h>
24: #include <string.h>
25:
26: #include "heccer/heccer.h"
27:
28:
29: static int HeccerAggregatorsInitialize(struct Heccer *pheccer);
30:
31:
32: /// **************************************************************************
33: ///
34: /// SHORT: HeccerAggregatorsInitialize()
35: ///
36: /// ARGS.:
37: ///
38: /// pheccer...: a heccer.
39: ///
40: /// RTN..: int
41: ///
42: /// Result of operation.
43: ///
44: /// DESCR: Zero out result arrays for aggregation operators.
45: ///
46: /// **************************************************************************
47:
48: static int HeccerAggregatorsInitialize(struct Heccer *pheccer)
49: {
50: //- set default result : ok
51:
52: int iResult = TRUE;
53:
54: //- zero out ...
55:
56: int i;
57:
58: for (i = 0 ; i < pheccer->vm.iAggregators ; i++)
59: {
60: pheccer->vm.pdAggregators[i] = 0.0;
61: }
62:
63: //- return result
64:
65: return(iResult);
66: }
67:
68:
69: /// **************************************************************************
70: ///
71: /// SHORT: HeccerAggregatorsCompile()
72: ///
73: /// ARGS.:
74: ///
75: /// pheccer...: a heccer.
76: ///
77: /// RTN..: int
78: ///
79: /// Result of operation.
80: ///
81: /// DESCR: Allocate result arrays for aggregation operators.
82: ///
83: /// **************************************************************************
84:
85: int HeccerAggregatorsCompile(struct Heccer *pheccer)
86: {
87: //- check for errors
88:
89: if (pheccer->iErrorCount)
90: {
91: return(FALSE);
92: }
93:
94: //- set default result : ok
95:
96: int iResult = TRUE;
97:
98: //- allocate memory
99:
100: pheccer->vm.pdAggregators = calloc(pheccer->vm.iAggregators, sizeof(*pheccer->vm.pdAggregators));
101:
102: if (!pheccer->vm.pdAggregators)
103: {
104: return(FALSE);
105: }
106:
107: //- return result
108:
109: return(iResult);
110: }
111:
112:
113: /// **************************************************************************
114: ///
115: /// SHORT: HeccerCanCompile()
116: ///
117: /// ARGS.:
118: ///
119: /// pheccer...: a heccer.
120: ///
121: /// RTN..: int
122: ///
123: /// Model can be compiled.
124: ///
125: /// DESCR: Can the model be compiled, given the current options ?
126: ///
127: /// **************************************************************************
128:
129: int HeccerCanCompile(struct Heccer *pheccer)
130: {
131: //- check for errors
132:
133: if (pheccer->iErrorCount)
134: {
135: return(FALSE);
136: }
137:
138: //- set default result : ok
139:
140: int iResult = TRUE;
141:
142: #define MINIMAL_TIME_STEP 1e-30
143:
144: if (pheccer->dStep < MINIMAL_TIME_STEP)
145: {
146: HeccerError
147: (pheccer,
148: NULL,
149: "illegal time step (smallest is %g), cannot compile\n", MINIMAL_TIME_STEP);
150:
151: return(FALSE);
152: }
153:
154: if (pheccer->inter.iCompartments == 0)
155: {
156: HeccerError
157: (pheccer,
158: NULL,
159: "no compartments found in the intermediary, cannot compile this model\n");
160:
161: return(FALSE);
162: }
163:
164: //- return result
165:
166: return(iResult);
167: }
168:
169:
170: /// **************************************************************************
171: ///
172: /// SHORT: HeccerCompileP1()
173: ///
174: /// ARGS.:
175: ///
176: /// pheccer...: a heccer.
177: ///
178: /// RTN..: int
179: ///
180: /// success of operation.
181: ///
182: /// DESCR: Compile a model into an intermediary format.
183: ///
184: /// **************************************************************************
185:
186: int HeccerCompileP1(struct Heccer *pheccer)
187: {
188: //- check for errors
189:
190: if (pheccer->iErrorCount)
191: {
192: return(FALSE);
193: }
194:
195: //- if the intermediary is already available
196:
197: if (pheccer->iStatus >= HECCER_STATUS_PHASE_2)
198: {
199: //- return success
200:
201: return(TRUE);
202: }
203:
204: //- set default result : ok
205:
206: int iResult = TRUE;
207:
208: //- get access to the translation service
209:
210: struct TranslationService *pts = pheccer->pts;
211:
212: /* struct TranslationServiceData *ptsd = pts->ptsd; */
213:
214: //- build up intermediary using the available service
215:
216: ComponentInspector si = pts->segments_inspector;
217:
218: iResult = iResult && si(pheccer, pts);
219:
220: ComponentInspector mi = pts->mechanisms_inspector;
221:
222: iResult = iResult && mi(pheccer, pts);
223:
224: //- set new status
225:
226: if (iResult)
227: {
228: pheccer->iStatus = HECCER_STATUS_PHASE_2;
229: }
230: else
231: {
232: //t do something sensible here
233:
234: //t HeccerError()
235: }
236:
237: //- return result
238:
239: return(iResult);
240: }
241:
242:
243: /// **************************************************************************
244: ///
245: /// SHORT: HeccerCompileP2()
246: ///
247: /// ARGS.:
248: ///
249: /// pheccer...: a heccer.
250: ///
251: /// RTN..: int
252: ///
253: /// success of operation.
254: ///
255: /// DESCR: Analyze the model, build indices for optimization.
256: ///
257: /// Internally, Heccer addresses mechanisms using their
258: /// compartment's schedule number. So the minimum degree
259: /// algorithm must run first before the mechanisms can be
260: /// compiled.
261: ///
262: /// NOTE.:
263: ///
264: /// This function can be used for testing internals of heccer,
265: /// just be sure to provide a consistent intermediary image.
266: ///
267: /// **************************************************************************
268:
269: int HeccerCompileP2(struct Heccer *pheccer)
270: {
271: //- check for errors
272:
273: if (pheccer->iErrorCount)
274: {
275: return(FALSE);
276: }
277:
278: //- set default result : ok
279:
280: int iResult = TRUE;
281:
282: //- first sanity
283:
284: if (!HeccerCanCompile(pheccer))
285: {
286: iResult = FALSE;
287: }
288:
289: //- do minimum degree
290:
291: iResult = iResult && HeccerMinimumDegree(pheccer);
292:
293: //- index & sort mechanisms
294:
295: iResult = iResult && HeccerMechanismSort(pheccer);
296:
297: //- return result
298:
299: return(iResult);
300: }
301:
302:
303: /// **************************************************************************
304: ///
305: /// SHORT: HeccerCompileP3()
306: ///
307: /// ARGS.:
308: ///
309: /// pheccer...: a heccer.
310: ///
311: /// RTN..: int
312: ///
313: /// success of operation.
314: ///
315: /// DESCR: Compile the intermediary format into byte code.
316: ///
317: /// Uses indices, initialized with HeccerCompileP2().
318: ///
319: /// NOTE.:
320: ///
321: /// This function can be used for testing internals of heccer,
322: /// just be sure to provide a consistent intermediary image.
323: ///
324: /// **************************************************************************
325:
326: int HeccerCompileP3(struct Heccer *pheccer)
327: {
328: //- check for errors
329:
330: if (pheccer->iErrorCount)
331: {
332: return(FALSE);
333: }
334:
335: //- set default result : ok
336:
337: int iResult = TRUE;
338:
339: //- compile compartments to byte code
340:
341: iResult = iResult && HeccerCompartmentCompile(pheccer);
342:
343: //- compile mechanisms to byte code
344:
345: iResult = iResult && HeccerMechanismCompile(pheccer);
346:
347: //- link mechanisms
348:
349: iResult = iResult && HeccerMechanismLink(pheccer);
350:
351: //t perhaps should discretize channels overhere ?
352:
353: //- allocate memory for aggragate results
354:
355: iResult = iResult && HeccerAggregatorsCompile(pheccer);
356:
357: //- return result
358:
359: return(iResult);
360: }
361:
362:
363: /// **************************************************************************
364: ///
365: /// SHORT: HeccerDump()
366: ///
367: /// ARGS.:
368: ///
369: /// pheccer...: a heccer.
370: /// pfile.....: stdio file.
371: /// iSelection: selection to dump.
372: ///
373: /// RTN..: int
374: ///
375: /// success of operation.
376: ///
377: /// DESCR: Call the dump functions, with the given selection.
378: ///
379: /// The selection is the boolean or of zero or more of the following :
380: ///
381: ///
382: ///
383: /// The shorthand HECCER_DUMP_ALL selects everything.
384: ///
385: /// **************************************************************************
386:
387: int HeccerDumpV(struct Heccer *pheccer)
388: {
389: return(HeccerDump(pheccer, stdout, HECCER_DUMP_ALL));
390: }
391:
392: int HeccerDump(struct Heccer *pheccer, FILE *pfile, int iSelection)
393: {
394: //- check for errors
395:
396: if (pheccer->iErrorCount)
397: {
398: return(FALSE);
399: }
400:
401: //- set default result : ok
402:
403: int iResult = TRUE;
404:
405: if (!pfile)
406: {
407: pfile = stdout;
408: }
409:
410: //- name: always present
411:
412: fprintf(pfile, "Heccer (pcName) : (%s)\n", pheccer->pcName ? pheccer->pcName : "(null)");
413:
414: //- status: reflects phases of compilation.
415:
416: fprintf(pfile, "Heccer (iStatus) : (%i)\n", pheccer->iStatus);
417:
418: fprintf(pfile, "Heccer (iErrorCount) : (%i)\n", pheccer->iErrorCount);
419:
420: //- options and operation mode.
421:
422: fprintf(pfile, "Heccer Options (iOptions) : (%i)\n", pheccer->ho.iOptions);
423:
424: fprintf(pfile, "Heccer Options (dIntervalStart) : (%g)\n", pheccer->ho.dIntervalStart);
425:
426: fprintf(pfile, "Heccer Options (dIntervalEnd) : (%g)\n", pheccer->ho.dIntervalEnd);
427:
428: fprintf(pfile, "Heccer Options (dConcentrationGateStart) : (%g)\n", pheccer->ho.dConcentrationGateStart);
429:
430: fprintf(pfile, "Heccer Options (dConcentrationGateEnd) : (%g)\n", pheccer->ho.dConcentrationGateEnd);
431:
432: fprintf(pfile, "Heccer Options (iIntervalEntries) : (%i)\n", pheccer->ho.iIntervalEntries);
433:
434: fprintf(pfile, "Heccer Options (iSmallTableSize) : (%i)\n", pheccer->ho.iSmallTableSize);
435:
436: //- simulation time
437:
438: fprintf(pfile, "Heccer (dTime) : (%g)\n", pheccer->dTime);
439:
440: //- time step
441:
442: fprintf(pfile, "Heccer (dStep) : (%g)\n", pheccer->dStep);
443:
444: //- identification service : translates serials to math components
445:
446: if (iSelection & HECCER_DUMP_SERVICE)
447: {
448: fprintf(pfile, "Heccer (pts) : (%p)\n", pheccer->pts);
449: fprintf(pfile, "Heccer (ped) : (%p)\n", pheccer->ped);
450: }
451:
452: //- dump intermediary
453:
454: iResult = iResult && HeccerIntermediaryDump(&pheccer->inter, pfile, iSelection);
455:
456: //- dump the indices
457:
458: iResult = iResult && HeccerIndexersDump(&pheccer->indexers, pfile, iSelection);
459:
460: //- dump the tables
461:
462: iResult = iResult && HeccerTablesDump(&pheccer->tgt, pfile, iSelection);
463:
464: //- dump compartment arrays
465:
466: iResult = iResult && HeccerVMDump(&pheccer->vm, pfile, iSelection);
467:
468: //- flush output
469:
470: fflush(pfile);
471:
472: //- return result
473:
474: return(iResult);
475: }
476:
477:
478: /// **************************************************************************
479: ///
480: /// SHORT: HeccerError()
481: ///
482: /// ARGS.:
483: ///
484: /// pheccer..: a heccer.
485: /// pcContext: context of error.
486: /// pcError..: error string.
487: ///
488: /// RTN..: int
489: ///
490: /// success of operation.
491: ///
492: /// DESCR: Register an error, print to stderr.
493: ///
494: /// **************************************************************************
495:
496: int HeccerError(struct Heccer *pheccer, char *pcContext, char *pcError, ...)
497: {
498: //- set default result: ok
499:
500: int iResult = 1;
501:
502: //- give diagnostics
503:
504: //- print to stderr
505:
506: fprintf(stderr, "Heccer the hecc :");
507:
508: /* if (pcContext) */
509: /* { */
510: /* fprintf */
511: /* (stderr, */
512: /* "%s: *** Error: ", */
513: /* pcContext); */
514: /* } */
515:
516: //v stdargs list
517:
518: va_list vaList;
519:
520: //- get start of stdargs
521:
522: va_start(vaList, pcError);
523:
524: //- give diagnostics
525:
526: vfprintf(stderr, pcError, vaList);
527:
528: fprintf(stderr, "\n");
529:
530: //- end stdargs
531:
532: va_end(vaList);
533:
534: //- negate status: this heccer is in error
535:
536: pheccer->iStatus = - pheccer->iStatus;
537:
538: //- increment total error count
539:
540: pheccer->iErrorCount++;
541:
542: //- return result
543:
544: return(iResult);
545: }
546:
547:
548: /// **************************************************************************
549: ///
550: /// SHORT: HeccerGetVersion()
551: ///
552: /// ARGS.:
553: ///
554: /// RTN..: char *
555: ///
556: /// Version identifier.
557: ///
558: /// DESCR: Obtain version identifier.
559: ///
560: /// **************************************************************************
561:
562: char * HeccerGetVersion(void)
563: {
564: // $Format: " static char *pcVersion=\"${label}\";"$
565: static char *pcVersion="network-12";
566:
567: return(pcVersion);
568: }
569:
570:
571: /// **************************************************************************
572: ///
573: /// SHORT: HeccerHecc()
574: ///
575: /// ARGS.:
576: ///
577: /// pheccer...: a heccer.
578: ///
579: /// RTN..: int
580: ///
581: /// success of operation.
582: ///
583: /// DESCR: Compute one step of simulation time.
584: ///
585: /// **************************************************************************
586:
587: static int HeccerHecc(struct Heccer *pheccer)
588: {
589: //- check for errors
590:
591: if (pheccer->iErrorCount)
592: {
593: return(FALSE);
594: }
595:
596: //- set default result : ok
597:
598: int iResult = TRUE;
599:
600: //- update the simulation time
601:
602: pheccer->dTime += pheccer->dStep;
603:
604: //- initialize the aggregate results
605:
606: if (pheccer->vm.pdAggregators)
607: {
608: iResult = iResult && HeccerAggregatorsInitialize(pheccer);
609: }
610:
611: //! I am undecided where to make the difference between CN and BE.
612: //! From cosmetic viewpoint, here we should only delegate to
613: //! HeccerMechanismHecc() and HeccerCompartmentHecc().
614: //! Then the difference can be made by looking at the options
615: //! (bad~?), or by having byte-code generated for the two specific
616: //! cases (slower~?).
617:
618: //! On the other hand, this is less work for the moment and delays
619: //! the decisions involved in the above (perhaps even unnecessary
620: //! to be made~?).
621:
622: //- perform mechanism operations (including concentration pools)
623:
624: //! doing it like this, assumes that the mechanism values are
625: //! correctly given to HeccerInitiate() at time (now - dt/2). Does
626: //! not make a practical difference.
627:
628: iResult = iResult && HeccerMechanismSolveCN(pheccer);
629:
630: //- perform the compartment operations
631:
632: if (pheccer->vm.iVms != 1)
633: {
634: iResult = iResult && HeccerCompartmentSolveCN(pheccer);
635: }
636:
637: //- return result
638:
639: return(iResult);
640: }
641:
642:
643: /// **************************************************************************
644: ///
645: /// SHORT: HeccerHeccs()
646: ///
647: /// ARGS.:
648: ///
649: /// pheccer...: a heccer.
650: /// dTime.....: current time.
651: ///
652: /// RTN..: int
653: ///
654: /// success of operation.
655: ///
656: /// DESCR: Call HeccerHecc() until dTime.
657: ///
658: /// **************************************************************************
659:
660: int HeccerHeccs(struct Heccer *pheccer, double dTime)
661: {
662: //- check for errors
663:
664: if (pheccer->iErrorCount)
665: {
666: return(FALSE);
667: }
668:
669: //- set default result : ok
670:
671: int iResult = TRUE;
672:
673: //- while not reached simulation time
674:
675: while (pheccer->dTime < dTime)
676: {
677: //- do a single hecc
678:
679: iResult = iResult && HeccerHecc(pheccer);
680:
681: //! perhaps should move the advance of the local time to this
682: //! point ? Would allow to remove this test ...
683: //! I don't care for the moment.
684:
685: if (iResult == FALSE)
686: {
687: return(iResult);
688: }
689: }
690:
691: //- return result
692:
693: return(iResult);
694: }
695:
696:
697: /// **************************************************************************
698: ///
699: /// SHORT: HeccerInitiate()
700: ///
701: /// ARGS.:
702: ///
703: /// pheccer...: a heccer.
704: ///
705: /// RTN..: int
706: ///
707: /// success of operation.
708: ///
709: /// DESCR: Fill the data arrays with initial values.
710: ///
711: /// TODO.:
712: ///
713: /// I guess it should be possible to override the initial value
714: /// array using an additional argument. Providing NULL means not
715: /// overriden. Probably the initial values should be separated
716: /// out from the intermediary and have their own data structure.
717: /// In that case, we can also annotate the initial values with
718: /// flags, e.g. indicating that they are calculated or not (for
719: /// channel equilibrium state fi.).
720: ///
721: /// **************************************************************************
722:
723: int HeccerInitiate(struct Heccer *pheccer)
724: {
725: //- check for errors
726:
727: if (pheccer->iErrorCount)
728: {
729: return(FALSE);
730: }
731:
732: //- set default result : ok
733:
734: int iResult = TRUE;
735:
736: //- initiate compartments
737:
738: iResult = iResult && HeccerCompartmentInitiate(pheccer);
739:
740: /* //- initiate mechanisms */
741:
742: /* iResult = iResult && HeccerMechanismInitiate(pheccer); */
743:
744: //- return result
745:
746: return(iResult);
747: }
748:
749:
750: /// **************************************************************************
751: ///
752: /// SHORT: HeccerNew()
753: ///
754: /// ARGS.:
755: ///
756: /// pc.....: name of this heccer, may be NULL.
757: /// pts....: translation service.
758: /// ped....: event distribution service.
759: /// peq....: event queuing service.
760: ///
761: /// RTN..: struct Heccer *
762: ///
763: /// Instantiated heccer, NULL for failure.
764: ///
765: /// DESCR: Create a new heccer using defaults.
766: ///
767: /// Defaults include option values and time step. Look at the
768: /// code to see what they really are.
769: ///
770: /// **************************************************************************
771:
772: struct Heccer *
773: HeccerNew
774: (char *pc,
775: struct TranslationService *pts,
776: struct EventDistributor *ped,
777: struct EventQueuer *peq)
778: {
779: //- set result : initialized heccer
780:
781: struct Heccer *pheccerResult
782: = HeccerNewP1
783: (
784: pc,
785: pts,
786: ped,
787: peq,
788: 0, // HECCER_OPTION_LOGICAL_BRANCH_SCHEDULING,
789: 2e-5
790: );
791:
792: //- return result
793:
794: return(pheccerResult);
795: }
796:
797:
798: /// **************************************************************************
799: ///
800: /// SHORT: HeccerNewP1()
801: ///
802: /// ARGS.:
803: ///
804: /// pc.......: name of this heccer, may be NULL.
805: /// pts......: identification service.
806: /// ped......: event distribution service.
807: /// peq......: event queuing service.
808: /// iOptions.: see heccer.h.
809: /// dStep....: required time step (from the time constants of the model).
810: ///
811: /// RTN..: struct Heccer *
812: ///
813: /// Instantiated heccer, NULL for failure.
814: ///
815: /// DESCR: Create a new heccer.
816: ///
817: /// **************************************************************************
818:
819: struct Heccer *
820: HeccerNewP1
821: (char *pc,
822: struct TranslationService *pts,
823: struct EventDistributor *ped,
824: struct EventQueuer *peq,
825: int iOptions,
826: double dStep)
827: {
828: //- set result : a new heccer
829:
830: struct Heccer *pheccerResult
831: = (struct Heccer *)calloc(1, sizeof(struct Heccer));
832:
833: if (!pheccerResult)
834: {
835: return(NULL);
836: }
837:
838: //- set name
839:
840: pheccerResult->pcName = pc;
841:
842: //- set naming service
843:
844: pheccerResult->pts = pts;
845:
846: //- set event distribution service
847:
848: pheccerResult->ped = ped;
849:
850: //- set event queuing service
851:
852: pheccerResult->peq = peq;
853:
854: //- set options and time step
855:
856: pheccerResult->ho.iOptions = iOptions;
857:
858: pheccerResult->dStep = dStep;
859:
860: //- initialize table discretization options
861:
862: pheccerResult->ho.dIntervalStart = HECCER_INTERVAL_DEFAULT_START;
863: pheccerResult->ho.dIntervalEnd = HECCER_INTERVAL_DEFAULT_END;
864:
865: pheccerResult->ho.dConcentrationGateStart = HECCER_INTERVAL_CONCENTRATION_GATE_DEFAULT_START;
866: pheccerResult->ho.dConcentrationGateEnd = HECCER_INTERVAL_CONCENTRATION_GATE_DEFAULT_END;
867:
868: pheccerResult->ho.iIntervalEntries = HECCER_INTERVAL_DEFAULT_ENTRIES;
869:
870: pheccerResult->ho.iSmallTableSize = HECCER_INTERPOL_INTERVAL_DEFAULT_ENTRIES;
871:
872: //- set new status
873:
874: pheccerResult->iStatus = HECCER_STATUS_PHASE_1;
875:
876: //- return result
877:
878: return(pheccerResult);
879: }
880:
881:
882: /// **************************************************************************
883: ///
884: /// SHORT: HeccerNewP2()
885: ///
886: /// ARGS.:
887: ///
888: /// pc......: name of this heccer, may be NULL.
889: /// pinter..: intermediary with a complete numerical model definition.
890: ///
891: /// RTN..: struct Heccer *
892: ///
893: /// Instantiated heccer, NULL for failure.
894: ///
895: /// DESCR: Create a new heccer.
896: ///
897: /// **************************************************************************
898:
899: struct Heccer *HeccerNewP2(char *pc, struct Intermediary *pinter)
900: {
901: //- set result : initialized heccer
902:
903: struct Heccer *pheccerResult = HeccerNew(pc, NULL, NULL, NULL);
904:
905: //- link in intermediary
906:
907: memcpy(&pheccerResult->inter, pinter, sizeof(*pinter));
908:
909: //- set new status
910:
911: pheccerResult->iStatus = HECCER_STATUS_PHASE_2;
912:
913: //- return result
914:
915: return(pheccerResult);
916: }
917:
918:
919: // Local variables:
920: // eval: (add-hook 'write-file-hooks 'time-stamp)
921: // time-stamp-start: "pcVersionTime="
922: // time-stamp-format: "\"(%02y/%02m/%02d) %a, %b %d, %y %02H:%02M:%02S %u\";"
923: // time-stamp-end: "$"
924: // End:
925:
Generated by Xrefactory version 2.0.14 on Thu Jul 24 22:41:20 2008