file:/local_home/local_home/hugo/neurospaces_project/heccer/source/c/snapshots/0/compartment.c (Mon Jun 16 00:02:57 2008
) HOME
1: //
2: // Heccer : a compartmental solver that implements efficient Crank-Nicolson
3: // integration for neuronal models.
4: //
5:
6: //////////////////////////////////////////////////////////////////////////////
7: //'
8: //' Heccer : testbed C implementation
9: //'
10: //' Copyright (C) 2006-2008 Hugo Cornelis
11: //'
12: //' functional ideas .. Hugo Cornelis, hugo.cornelis@gmail.com
13: //'
14: //' coding ............ Hugo Cornelis, hugo.cornelis@gmail.com
15: //'
16: //////////////////////////////////////////////////////////////////////////////
17:
18:
19: #include <stdio.h>
20: #include <stdlib.h>
21: #include <string.h>
22:
23: #include "heccer/compartment.h"
24: #include "heccer/heccer.h"
25:
26:
27: /// **************************************************************************
28: ///
29: /// SHORT: HeccerCompileCompartments()
30: ///
31: /// ARGS.:
32: ///
33: /// pheccer...: a heccer.
34: ///
35: /// RTN..: int
36: ///
37: /// success of operation.
38: ///
39: /// DESCR: Compile the intermediary of the compartments to byte code.
40: ///
41: /// **************************************************************************
42:
43: int HeccerCompartmentCompile(struct Heccer *pheccer)
44: {
45: //- set default result : ok
46:
47: int iResult = TRUE;
48:
49: //v time step
50:
51: double dt;
52:
53: //- for backward euler integration
54:
55: if (pheccer->ho.iOptions & HECCER_OPTION_BACKWARD_EULER)
56: {
57: //- remember to do full time step
58:
59: dt = pheccer->dStep;
60: }
61:
62: //- else : crank-nicholson
63:
64: else
65: {
66: //- halve the time step
67:
68: dt = pheccer->dStep / 2.0;
69: }
70:
71: //v make minimum degree indices accessible
72:
73: int *piChildren = pheccer->indexers.md.piChildren;
74:
75: int **ppiChildren = pheccer->indexers.md.ppiChildren;
76:
77: int *piForward = pheccer->indexers.md.piForward;
78:
79: int *piBackward = pheccer->indexers.md.piBackward;
80:
81: //- create storage for diagterms
82:
83: pheccer->vm.iDiagonals = pheccer->inter.iCompartments;
84:
85: double *pdDiagonals = (double *)calloc(pheccer->vm.iDiagonals, sizeof(double));
86:
87: pheccer->vm.pdDiagonals = pdDiagonals;
88:
89: //- loop over all compartments (in flow order)
90:
91: int i;
92:
93: for (i = 0 ; i < pheccer->inter.iCompartments ; i++)
94: {
95: //- get intermediary number for the current compartment
96:
97: int iIntermediary = piBackward[i];
98:
99: //- get compartment parameters
100:
101: double dCm = pheccer->inter.pcomp[iIntermediary].dCm;
102:
103: double dRm = pheccer->inter.pcomp[iIntermediary].dRm;
104:
105: //- compute diagonal term : time constant
106:
107: pdDiagonals[i] = 1.0 + dt / (dCm * dRm);
108:
109: #ifdef DEBUG_OPERATIONS_DIAGONALS
110:
111: fprintf(stdout, "pdDiagonals[%i] : %e\t= 1.0 + %e / (%e * %e)\n", i, pdDiagonals[i], dt, dCm, dRm);
112:
113: #endif
114:
115: //- get compartment's parent number
116:
117: int iParent = pheccer->inter.pcomp[iIntermediary].iParent;
118:
119: //- get number of children of current compartment
120:
121: int iChildren = piChildren[iIntermediary];
122:
123: //- loop over all children
124:
125: int iChild;
126:
127: for (iChild = 0 ; iChild < iChildren ; iChild++)
128: {
129: //- get index number for child of current compartment
130:
131: int iModelNumber = ppiChildren[iIntermediary][iChild];
132:
133: //- compute reverse contribution (genesis raxial)
134:
135: //t check old hsolve code, this computation seems
136: //t contradictory to me, mixes up child and parent
137: //t properties. ...
138: //t ... ok, done, perhaps this is a bug in hsolve ?
139: //t have to check the computations carefully again.
140:
141: double dReverse = -dt / (pheccer->inter.pcomp[iIntermediary].dCm * pheccer->inter.pcomp[iModelNumber].dRa);
142:
143: pdDiagonals[i] -= dReverse;
144:
145: #ifdef DEBUG_OPERATIONS_DIAGONALS
146:
147: fprintf(stdout, "pdDiagonals[%i] : %e (- %e)\n", i, pdDiagonals[i], dReverse);
148:
149: #endif
150:
151: //- compute verse contribution, comes from child (genesis axial)
152:
153: double dVerse = -dt / (pheccer->inter.pcomp[iModelNumber].dCm * pheccer->inter.pcomp[iModelNumber].dRa);
154:
155: //- convert child compartment number to schedule number
156:
157: int iScheduleNumber = piForward[iModelNumber];
158:
159: pdDiagonals[iScheduleNumber] -= dVerse;
160:
161: #ifdef DEBUG_OPERATIONS_DIAGONALS
162:
163: fprintf(stdout,"pdDiagonals[%i] : %e (- %e)\n", iScheduleNumber, pdDiagonals[i], dVerse);
164:
165: #endif
166: }
167: }
168:
169: #define SETCOP1(array,current,operator) ((array) ? (array)[(current)++] = (operator) : (current)++)
170:
171: #define SETCOP2(array,current,operator,operand) ((array) ? ({(array)[(current)++] = (operator); (array)[(current)++] = (operand); }) : ({ (current) += 2; }))
172:
173: #define SETAXRES(array,current,value) ((array) ? ((array)[(current)++] = (value)) : (current)++)
174:
175: //- first count, next compile the following block
176:
177: int iCountCompile;
178:
179: //- first counting by setting array to NULL
180:
181: int *piCops = NULL;
182:
183: double *pdAxres = NULL;
184:
185: for (iCountCompile = 0 ; iCountCompile < 2 ; iCountCompile++)
186: {
187: //- counters always start at zero
188:
189: int iCops = 0;
190:
191: int iAxres = 0;
192:
193: //v for leaves in the mid of the array, immediate skip to diagonal
194:
195: int iSkipped = 0;
196:
197: //- loop over all compartments : start forward elimination
198:
199: int i;
200:
201: for (i = 0 ; i < pheccer->inter.iCompartments ; i++)
202: {
203: //- get intermediary number for the current compartment
204:
205: int iIntermediary = piBackward[i];
206:
207: //- get compartment capacity
208:
209: double dCm = pheccer->inter.pcomp[iIntermediary].dCm;
210:
211: //- get compartment's parent number
212:
213: int iParent = pheccer->inter.pcomp[iIntermediary].iParent;
214:
215: //- get number of children of current compartment
216:
217: int iChildren = piChildren[iIntermediary];
218:
219: //- if skipped previous compartment
220:
221: if (iSkipped)
222: {
223: //- can only occur once in a row, so next sets diagonal
224:
225: //t easy consistency check
226:
227: iSkipped = 0;
228: }
229:
230: //- previous compartment has set diagonal
231:
232: else
233: {
234: //- if current compartment has children
235:
236: if (iChildren)
237: {
238: //- skip first two compartments, other always set the diagonal
239:
240: if (i > 1)
241: {
242: SETCOP1(piCops, iCops, HECCER_COP_SET_DIAGONAL);
243: }
244: }
245:
246: //- if no children for the current compartment
247:
248: else
249: {
250: //- starting from the second compartment
251:
252: if (i > 1)
253: {
254: SETCOP1(piCops, iCops, HECCER_COP_NEXT_ROW);
255:
256: //- remember skipped to the diagonal for the next compartment
257:
258: iSkipped = 1;
259: }
260: }
261: }
262:
263: //- loop over all children
264:
265: int j;
266:
267: for (j = 0; j < iChildren; j++)
268: {
269: //- get index number for child of current compartment
270:
271: int iModelNumber = ppiChildren[iIntermediary][j];
272:
273: //- convert to schedule number
274:
275: int iScheduleNumber = piForward[iModelNumber];
276:
277: SETCOP2(piCops, iCops, HECCER_COP_FORWARD_ELIMINATION, 2 * iScheduleNumber);
278:
279: //- compute reverse contribution (genesis raxial)
280:
281: //t check old hsolve code, this computation seems
282: //t contradictory to me, mixes up child and parent
283: //t properties. ...
284: //t ... ok, done, perhaps this is a bug in hsolve ?
285: //t have to check the computations carefully again.
286:
287: double dReverse = -dt / (pheccer->inter.pcomp[iIntermediary].dCm * pheccer->inter.pcomp[iModelNumber].dRa);
288:
289: SETAXRES(pdAxres, iAxres, dReverse);
290:
291: //- compute verse contribution, comes from child (genesis axial)
292:
293: double dVerse = -dt / (pheccer->inter.pcomp[iModelNumber].dCm * pheccer->inter.pcomp[iModelNumber].dRa);
294:
295: SETAXRES(pdAxres, iAxres, dVerse);
296: }
297: }
298:
299: //- finishing forward elimination
300:
301: SETCOP1(piCops, iCops, HECCER_COP_FINISH);
302:
303: //- start backward substitution
304:
305: for (i = pheccer->inter.iCompartments - 2 ; i >= 0 ; i--)
306: {
307: //- get intermediary number for the current compartment
308:
309: int iIntermediary = piBackward[i];
310:
311: //- get compartment capacity
312:
313: double dCm = pheccer->inter.pcomp[iIntermediary].dCm;
314:
315: //- get compartment's parent number
316:
317: int iParent = pheccer->inter.pcomp[iIntermediary].iParent;
318:
319: //- get number of children of current compartment
320:
321: int iChildren = piChildren[iIntermediary];
322:
323: //- convert to schedule number
324:
325: int iScheduleNumber = piForward[iParent];
326:
327: SETCOP2(piCops, iCops, HECCER_COP_BACKWARD_SUBSTITUTION, 2 * iScheduleNumber);
328:
329: SETCOP1(piCops, iCops, HECCER_COP_FINISH_ROW);
330:
331: //- set axial resistance value for child
332:
333: double dRa = pheccer->inter.pcomp[iIntermediary].dRa;
334:
335: SETAXRES(pdAxres, iAxres, -dt / (dCm * dRa));
336: }
337:
338: //- finish backward substitution
339:
340: SETCOP1(piCops, iCops, HECCER_COP_FINISH);
341:
342: //- if we were counting during the previous loop
343:
344: if (!pheccer->vm.piCops)
345: {
346: //- prepare for compilation : allocate ->piCops and ->pdAxres, set counters
347:
348: pheccer->vm.piCops
349: = (int *)calloc(iCops, sizeof(int));
350:
351: piCops = pheccer->vm.piCops;
352:
353: pheccer->vm.iCops = iCops;
354:
355: pheccer->vm.pdAxres
356: = (double *)calloc(iAxres, sizeof(double));
357:
358: pdAxres = pheccer->vm.pdAxres;
359:
360: pheccer->vm.iAxres = iAxres;
361: }
362: }
363:
364: pheccer->vm.iResults = 2 * pheccer->inter.iCompartments;
365:
366: pheccer->vm.pdResults = (double *)calloc(pheccer->vm.iResults, sizeof(double));
367:
368: pheccer->vm.iVms = pheccer->inter.iCompartments;
369:
370: pheccer->vm.pdVms = (double *)calloc(pheccer->vm.iVms, sizeof(double));
371:
372: //! the diagonal terms are also needed for solving the current
373: //! equations, but I guess it is ok to deal with that in the
374: //! mechanism compiler. Check needed.
375:
376: //- return result
377:
378: return(iResult);
379: }
380:
381:
382: /// **************************************************************************
383: ///
384: /// SHORT: HeccerCompartmentDump()
385: ///
386: /// ARGS.:
387: ///
388: /// pcomp......: a compartment.
389: /// pfile......: stdio file.
390: /// iSelection.: selection to dump.
391: ///
392: /// RTN..: int
393: ///
394: /// success of operation.
395: ///
396: /// DESCR: Dump compartment parameters to the given stream.
397: ///
398: /// **************************************************************************
399:
400: int
401: HeccerCompartmentDump
402: (struct Compartment *pcomp, FILE *pfile, int iSelection)
403: {
404: //- set default result
405:
406: int iResult = TRUE;
407:
408: //- administrative overhead
409:
410: //! this makes testing quite hard, needs careful thought
411:
412: #ifdef lkjlkslkjasdf
413:
414: HECCER_SOURCE_NEUROSPACES
415:
416: fprintf(pfile, "Compartment (mc.iSerial, mc.iType) : (%i, %i)\n", pcomp->mc.iSerial, pcomp->mc.iType);
417:
418: #else
419:
420: fprintf(pfile, "Compartment (mc.iType) : (%i)\n", pcomp->mc.iType);
421:
422: #endif
423:
424: //- index of parent compartment, -1 for none
425:
426: if (iSelection & HECCER_DUMP_INTERMEDIARY_STRUCTURE)
427: {
428: fprintf(pfile, "Compartment (iParent) : (%i)\n", pcomp->iParent);
429: }
430:
431: /* //m number of mechanisms */
432:
433: /* if (iSelection & HECCER_DUMP_INTERMEDIARY_STRUCTURE) */
434: /* { */
435: /* fprintf(pfile, "Compartment (iMechanisms) : (%i)\n", pcomp->iMechanisms); */
436: /* } */
437:
438: //m descriptive values, alphabetical order
439:
440: if (iSelection & HECCER_DUMP_INTERMEDIARY_COMPARTMENTS_PARAMETERS)
441: {
442: fprintf(pfile, "Compartment (dCm) : (%g)\n", pcomp->dCm);
443:
444: fprintf(pfile, "Compartment (dEm) : (%g)\n", pcomp->dEm);
445:
446: fprintf(pfile, "Compartment (dInitVm) : (%g)\n", pcomp->dInitVm);
447:
448: fprintf(pfile, "Compartment (dInject) : (%g)\n", pcomp->dInject);
449:
450: fprintf(pfile, "Compartment (dRa) : (%g)\n", pcomp->dRa);
451:
452: fprintf(pfile, "Compartment (dRm) : (%g)\n", pcomp->dRm);
453: }
454:
455: //- return result
456:
457: return(iResult);
458: }
459:
460:
461: /// **************************************************************************
462: ///
463: /// SHORT: HeccerCompartmentInitiate()
464: ///
465: /// ARGS.:
466: ///
467: /// pheccer...: a heccer.
468: ///
469: /// RTN..: int
470: ///
471: /// success of operation.
472: ///
473: /// DESCR: Fill the compartment arrays with initial values.
474: ///
475: /// **************************************************************************
476:
477: int HeccerCompartmentInitiate(struct Heccer *pheccer)
478: {
479: //- set default result
480:
481: int iResult = TRUE;
482:
483: //- loop over all compartments
484:
485: int i;
486:
487: for (i = 0 ; i < pheccer->inter.iCompartments ; i++)
488: {
489: //- get schedule number
490:
491: int iSchedule = pheccer->indexers.md.piForward[i];
492:
493: //- fill in initial membrane potential
494:
495: //! if I use initial vm of zero, I get NaN for vm ?
496:
497: pheccer->vm.pdVms[iSchedule] = pheccer->inter.pcomp[i].dInitVm;
498:
499: //t fill in precomputed values for channels
500:
501: }
502:
503: //- return result
504:
505: return(iResult);
506: }
507:
508:
509: /// **************************************************************************
510: ///
511: /// SHORT: HeccerCompartmentSolveCN()
512: ///
513: /// ARGS.:
514: ///
515: /// pheccer...: a heccer.
516: ///
517: /// RTN..: int
518: ///
519: /// success of operation.
520: ///
521: /// DESCR: Perform the compartment operators once.
522: ///
523: /// **************************************************************************
524:
525: int HeccerCompartmentSolveCN(struct Heccer *pheccer)
526: {
527: int *piCops = &pheccer->vm.piCops[0];
528: int iCop;
529:
530: double *pdAxres = &pheccer->vm.pdAxres[0];
531: double *pdResults = &pheccer->vm.pdResults[0];
532: double *pdResult = &pdResults[2];
533: double *pdVms = &pheccer->vm.pdVms[pheccer->inter.iCompartments];
534: double dDiagonal = pdResult[1];
535: double dResult = pdResult[0];
536:
537: //- start forward elimination at row 1
538:
539: int iDone = FALSE;
540:
541: while (!iDone)
542: {
543: iCop = piCops[0];
544: piCops++;
545:
546: if (iCop == HECCER_COP_FORWARD_ELIMINATION)
547: {
548: //! pdResults[piCops[0] + 1] contains other diagonal term
549:
550: double d = pdAxres[0] / pdResults[piCops[0] + 1];
551: dDiagonal -= pdAxres[1] * d;
552: pdAxres += 2;
553: dResult -= pdResults[piCops[0]] * d;
554: piCops++;
555: }
556: else if (iCop == HECCER_COP_SET_DIAGONAL)
557: {
558: //- store right side
559:
560: pdResult[0] = dResult;
561:
562: //- store left side
563:
564: pdResult[1] = dDiagonal;
565:
566: //- next row
567:
568: pdResult += 2;
569: dResult = pdResult[0];
570: dDiagonal = pdResult[1];
571: }
572: else if (iCop == HECCER_COP_NEXT_ROW)
573: {
574: //- store right side
575:
576: pdResult[0] = dResult;
577:
578: //- store left side
579:
580: pdResult[1] = dDiagonal;
581:
582: //- skip zero coefficients
583:
584: pdResult += 2;
585:
586: //- next row
587:
588: pdResult += 2;
589: dResult = pdResult[0];
590: dDiagonal = pdResult[1];
591: }
592: else
593: {
594: iDone = TRUE;
595: }
596: }
597:
598: //- last row done separately
599:
600: //t should try to incorporate this in the next loop
601: //t I have the feeling the loop is coded in the wrong way.
602:
603: /* SOLVE_ROW(dResult, dResult / dDiagonal); */
604:
605: dResult = dResult / dDiagonal;
606:
607: //- explicit forwards step for last row
608:
609: /* COMPUTE_VM(pdVms, dResult); */
610:
611: pdVms--;
612: pdVms[0] = dResult + dResult - pdVms[0];
613:
614: /* REGISTER_ROW_RIGHT(pdResult[0], dResult); */
615:
616: pdResult[0] = dResult;
617: pdResult -= 2;
618:
619: //- start new right side
620:
621: /* INITIATE_ROW(dResult, pdResult[0]); */
622:
623: dResult = pdResult[0];
624:
625: //- backwards elimination
626:
627: iDone = FALSE;
628:
629: while (!iDone)
630: {
631: //- load the operation
632:
633: iCop = piCops[0];
634: piCops++;
635:
636: if (iCop == HECCER_COP_BACKWARD_SUBSTITUTION)
637: {
638: //- fill in the off diagonal contribution
639:
640: /* OFF_DIAGONAL_CONTRIBUTION(dResult, pdAxres[0], pdResults, piCops[0]); */
641:
642: dResult -= pdAxres[0] * pdResults[piCops[0]];
643: pdAxres++;
644: piCops++;
645: }
646: else if (iCop == HECCER_COP_FINISH_ROW)
647: {
648: //- solve the matrix row
649:
650: //! pdResult[1] contains the diagonal
651:
652: /* SOLVE_ROW(dResult, dResult / pdResult[1]); */
653:
654: dResult = dResult / pdResult[1];
655:
656: //- explicit forward step to get the membrane potential
657:
658: /* COMPUTE_VM(pdVms, dResult); */
659:
660: pdVms--;
661: pdVms[0] = dResult + dResult - pdVms[0];
662:
663: /* REGISTER_ROW_RIGHT(pdResult[0], dResult); */
664:
665: pdResult[0] = dResult;
666: pdResult -= 2;
667:
668: //- start new right side
669:
670: /* INITIATE_ROW(dResult, pdResult[0]); */
671:
672: dResult = pdResult[0];
673: }
674: else
675: {
676: iDone = TRUE;
677: }
678: }
679:
680: //- return success
681:
682: return(TRUE);
683: }
684:
685:
686:
Generated by Xrefactory version 2.0.14 on Thu Jul 24 22:41:20 2008