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 <string.h> 20: 21: #include "heccer/heccer.h" 22: #include "heccer/vm.h" 23: 24: 25: struct HeccerCommandInfo 26: { 27: //m operation value 28: 29: int iValue; 30: 31: //m name of command 32: 33: char *pcName; 34: 35: //m length of command (number of operands + 1) 36: 37: int iLength; 38: 39: //m number of secondary operands, expressed in doubles 40: 41: int iSecondaries; 42: 43: //m sizeof secondaries 44: 45: int iSecondariesSize; 46: }; 47: 48: 49: struct HeccerCommandTable 50: { 51: //m initialization status 52: 53: int iStatus; 54: 55: //m number of commands 56: 57: int iCommands; 58: 59: //m -1 means no formatted output 60: 61: int iFormatted; 62: 63: //m command info 64: 65: struct HeccerCommandInfo * phci; 66: }; 67: 68: 69: static struct HeccerCommandInfo phciCops[] = 70: { 71: { HECCER_COP_FORWARD_ELIMINATION, "FORWARD_ELIMINATION", 2 * sizeof(int), }, 72: { HECCER_COP_BACKWARD_SUBSTITUTION, "BACKWARD_SUBSTITUTION", 2 * sizeof(int), }, 73: { HECCER_COP_FINISH_ROW, "FINISH_ROW", 1 * sizeof(int), }, 74: { HECCER_COP_FINISH, "FINISH", 1 * sizeof(int), }, 75: { HECCER_COP_SET_DIAGONAL, "SET_DIAGONAL", 1 * sizeof(int), }, 76: { HECCER_COP_NEXT_ROW, "NEXT_ROW", 1 * sizeof(int), }, 77: { -1, NULL, -1, }, 78: }; 79: 80: 81: static struct HeccerCommandTable hctCops = 82: { 83: 0, 84: sizeof(phciCops) / sizeof(struct HeccerCommandInfo) - 1, 85: 0, 86: phciCops, 87: }; 88: 89: 90: static struct HeccerCommandInfo phciMops[] = 91: { 92: { HECCER_MOP_CALLOUT, "CALLOUT", 1 * sizeof(int), 0, sizeof(struct MatsCallout), }, 93: { HECCER_MOP_COMPARTMENT, "COMPARTMENT", 1 * sizeof(int), 4, sizeof(struct MatsCompartment), }, 94: { HECCER_MOP_CONCEPTGATE, "CONCEPTGATE", sizeof(struct MopsSingleGateConcept), 1, sizeof(struct MatsSingleGateConcept), }, 95: { HECCER_MOP_EVENTGENERATE, "EVENTGENERATE", sizeof(struct MopsEventGenerate), 2, sizeof(struct MatsEventGenerate), }, 96: { HECCER_MOP_EXPONENTIALDECAY, "EXPONENTIALDECAY", sizeof(struct MopsExponentialDecay), 1, sizeof(struct MatsExponentialDecay), }, 97: { HECCER_MOP_FINISH, "FINISH", 1 * sizeof(int), 0, 0, }, 98: { HECCER_MOP_FLUXPOOL, "FLUXPOOL", sizeof(struct MopsFluxPool), 1, sizeof(struct MatsFluxPool), }, 99: { HECCER_MOP_INITIALIZECHANNEL, "INITIALIZECHANNEL", sizeof(struct MopsInitializeChannel), 0, 0, }, 100: { HECCER_MOP_INITIALIZECHANNELEREV, "INITIALIZECHANNELEREV", sizeof(struct MopsInitializeChannelErev), 0, 0, }, 101: { HECCER_MOP_LOADVOLTAGETABLE, "LOADVOLTAGETABLE", sizeof(struct MopsVoltageTableDependency), 0, 0, }, 102: { HECCER_MOP_REGISTERCHANNELCURRENT, "REGISTERCHANNELCURRENT", sizeof(struct MopsRegisterChannelCurrent), 0, 0, }, 103: { HECCER_MOP_UPDATECOMPARTMENTCURRENT, "UPDATECOMPARTMENTCURRENT", sizeof(struct MopsUpdateCompartmentCurrent), 0, 0, }, 104: { HECCER_MOP_AGGREGATECURRENT, "AGGREGATECURRENT", sizeof(struct MopsAggregateCurrent), 0, 0, }, 105: { HECCER_MOP_INTERNALNERNST, "INTERNALNERNST", sizeof(struct MopsInternalNernst), 1, sizeof(struct MatsInternalNernst), }, 106: { HECCER_MOP_SPRINGMASS, "SPRINGMASS", sizeof(struct MopsSpringMass), 3, sizeof(struct MatsSpringMass), }, 107: { HECCER_MOP_STORESINGLECHANNELCURRENT, "STORESINGLECHANNELCURRENT", sizeof(struct MopsStoreSingleChannelCurrent), 2, sizeof(struct MatsStoreSingleChannelCurrent), }, 108: { -1, NULL, -1, }, 109: }; 110: 111: 112: static struct HeccerCommandTable hctMops = 113: { 114: 0, 115: sizeof(phciMops) / sizeof(struct HeccerCommandInfo) - 1, 116: 1, 117: phciMops, 118: }; 119: 120: 121: static int operator_comparator(const void *pv1,const void *pv2) 122: { 123: struct HeccerCommandInfo *phci1 = (struct HeccerCommandInfo *)pv1; 124: struct HeccerCommandInfo *phci2 = (struct HeccerCommandInfo *)pv2; 125: 126: return(phci1->iValue < phci2->iValue 127: ? -1 128: : phci1->iValue > phci2->iValue ? 1 : 0); 129: } 130: 131: 132: static void SortTable(struct HeccerCommandTable *phct) 133: { 134: int iEntries = 0; 135: 136: while (phct->phci[iEntries].pcName != NULL) 137: { 138: iEntries++; 139: } 140: 141: phct->iCommands = iEntries; 142: 143: qsort 144: ((void *)phct->phci, 145: phct->iCommands, 146: sizeof(struct HeccerCommandInfo), 147: operator_comparator); 148: 149: phct->iStatus |= 1; 150: } 151: 152: 153: //f static prototypes 154: 155: static 156: struct HeccerCommandInfo * 157: HeccerCommandInfoLookup 158: (struct HeccerCommandTable *phct, 159: int iOp); 160: 161: static int 162: HeccerVMDumpOperators 163: (char * pcDescription, 164: struct VM *pvm, 165: int piArray[], 166: void *pvOperands, 167: struct HeccerCommandTable *phct, 168: int iStart, 169: int iEnd, 170: FILE *pfile); 171: 172: 173: /// *************************************************************************** 174: /// 175: /// NAME : HeccerCommandInfoLookup() 176: /// 177: /// SHORT : lookup info on given command 178: /// 179: /// PARAMETERS : 180: /// 181: /// phci....: array with info 182: /// iOp.....: operation to lookup 183: /// 184: /// RETURN : struct HeccerCommandInfo * : info on command 185: /// 186: /// DESCRIPTION : lookup info on given command 187: /// 188: /// *************************************************************************** 189: 190: static 191: struct HeccerCommandInfo * 192: HeccerCommandInfoLookup 193: (struct HeccerCommandTable *phct, 194: int iOp) 195: { 196: //- set default result : not found 197: 198: struct HeccerCommandInfo *phciResult = NULL; 199: 200: //- init top and bottom counters 201: 202: int iLower = 0; 203: int iUpper = phct->iCommands - 1; 204: 205: if (!(phct->iStatus & 1)) 206: { 207: SortTable(phct); 208: } 209: 210: iLower = 0; 211: iUpper = phct->iCommands - 1; 212: 213: //! binary search 214: 215: //- search until the range to search in becomes invalid 216: 217: while (iUpper - iLower >= 0) 218: { 219: //- determine the middle of the search range 220: 221: int iMiddle = (iLower + iUpper) / 2; 222: 223: //- get pointer to this command info entry 224: 225: struct HeccerCommandInfo *phci = &phct->phci[iMiddle]; 226: 227: //- set result and break out loop if search value is found here 228: 229: if (iOp == phci->iValue) 230: { 231: phciResult = phci; 232: 233: break; 234: } 235: 236: //- set a new lower or upper limit 237: 238: if (iOp > phci->iValue) 239: { 240: iLower = iMiddle + 1; 241: } 242: else 243: { 244: iUpper = iMiddle - 1; 245: } 246: } 247: 248: //- return result 249: 250: return(phciResult); 251: } 252: 253: 254: /// ************************************************************************** 255: /// 256: /// SHORT: HeccerVMDump() 257: /// 258: /// ARGS.: 259: /// 260: /// pvm........: heccer vm. 261: /// pfile......: stdio file. 262: /// iSelection.: selection to dump. 263: /// 264: /// RTN..: int 265: /// 266: /// success of operation. 267: /// 268: /// DESCR: Dump VM functions. 269: /// 270: /// ************************************************************************** 271: 272: int HeccerVMDump(struct VM *pvm, FILE *pfile, int iSelection) 273: { 274: //- set default result : ok 275: 276: int iResult = TRUE; 277: 278: //- dump compartment data 279: 280: if (iSelection & HECCER_DUMP_VM_COMPARTMENT_DATA) 281: { 282: //! nothing here 283: } 284: 285: //- dump compartment operations 286: 287: if (iSelection & HECCER_DUMP_VM_COMPARTMENT_OPERATIONS) 288: { 289: int iCops = pvm->iCops; 290: 291: int *piCops = pvm->piCops; 292: 293: //! cops are allocated as integers, dump array expects char's, 294: //! so that is why we have to multiply here. 295: 296: HeccerVMDumpOperators("Compartment operations", pvm, &piCops[0], NULL, &hctCops, 0, iCops * sizeof(int), pfile); 297: } 298: 299: //- dump mechanism data 300: 301: if (iSelection & HECCER_DUMP_VM_MECHANISM_DATA) 302: { 303: int iMats = pvm->iMats; 304: 305: void *pvMats = pvm->pvMats; 306: 307: /* HeccerVMDumpData("Mechanism data", pvm, pvMats, NULL, 0, iMats, pfile); */ 308: } 309: 310: //- dump mechanism operations 311: 312: if (iSelection & HECCER_DUMP_VM_MECHANISM_OPERATIONS) 313: { 314: int iMops = pvm->iMops; 315: 316: int *piMops = (int *)pvm->pvMops; 317: 318: HeccerVMDumpOperators("Mechanism operations", pvm, &piMops[0], pvm->pvMats, &hctMops, 0, iMops, pfile); 319: } 320: 321: /* //- dump channel to pool fluxes */ 322: 323: /* if (iSelection & HECCER_DUMP_VM_CHANNEL_POOL_FLUXES) */ 324: /* { */ 325: /* int i; */ 326: 327: /* for (i = 0 ; i < pvm->iFluxes ; i++) */ 328: /* { */ 329: /* fprintf(pfile, "VM Fluxes (pdFluxes[%i]) : (%g)\n", i, pvm->pdFluxes[i]); */ 330: /* } */ 331: /* } */ 332: 333: //- compartment data : diagonals 334: 335: if (iSelection & HECCER_DUMP_VM_COMPARTMENT_MATRIX) 336: { 337: int i; 338: 339: for (i = 0 ; i < pvm->iDiagonals ; i++) 340: { 341: fprintf(pfile, "VM Diagonals (pdDiagonals[%i]) : (%g)\n", i, pvm->pdDiagonals[i]); 342: } 343: } 344: 345: //- compartment data : axial resistances 346: 347: if (iSelection & HECCER_DUMP_VM_COMPARTMENT_MATRIX) 348: { 349: int i; 350: 351: for (i = 0 ; i < pvm->iAxres ; i++) 352: { 353: fprintf(pfile, "VM Axial Resistances (pdAxres[%i]) : (%g)\n", i, pvm->pdAxres[i]); 354: } 355: } 356: 357: //- results : intermediate 358: 359: if (iSelection & HECCER_DUMP_VM_COMPARTMENT_MATRIX) 360: { 361: int i; 362: 363: for (i = 0 ; i < pvm->iResults ; i++) 364: { 365: fprintf(pfile, "VM Axial Resistances (pdResults[%i]) : (%g)\n", i, pvm->pdResults[i]); 366: } 367: } 368: 369: //- results : membrane potentials 370: 371: if (iSelection & HECCER_DUMP_VM_COMPARTMENT_MATRIX) 372: { 373: int i; 374: 375: for (i = 0 ; i < pvm->iVms ; i++) 376: { 377: fprintf(pfile, "VM Membrane Potentials (pdVms[%i]) : (%g)\n", i, pvm->pdVms[i]); 378: } 379: } 380: 381: //- results : aggregators 382: 383: if (iSelection & HECCER_DUMP_VM_AGGREGATORS) 384: { 385: int i; 386: 387: for (i = 0 ; i < pvm->iAggregators ; i++) 388: { 389: fprintf(pfile, "VM Membrane Aggregator results (pdAggregators[%i]) : (%g)\n", i, pvm->pdAggregators[i]); 390: } 391: } 392: 393: //- return result 394: 395: return(iResult); 396: } 397: 398: 399: /// *************************************************************************** 400: /// 401: /// NAME : HeccerVMDumpOperators() 402: /// 403: /// SHORT : print operations in given array according to given info array 404: /// 405: /// PARAMETERS : 406: /// 407: /// pvm.....: heccer vm. 408: /// piArray.: array to print 409: /// iSize...: size of one array entry 410: /// phct....: table with info on array to print 411: /// iStart..: start index 412: /// iEnd....: ending index 413: /// 414: /// RETURN : int 415: /// 416: /// success of operation. 417: /// 418: /// DESCRIPTION : print operations in given array. 419: /// 420: /// NOTE : length of operations is ignored at the moment 421: /// 422: /// *************************************************************************** 423: 424: static 425: int 426: HeccerVMDumpOperators 427: (char * pcDescription, 428: struct VM *pvm, 429: int piOperators[], 430: void *pvOperands, 431: struct HeccerCommandTable *phct, 432: int iStart, 433: int iEnd, 434: FILE *pfile) 435: { 436: //- set default result : ok 437: 438: int iResult = TRUE; 439: 440: //- give the banner 441: 442: fprintf(pfile, "%s\n", pcDescription); 443: 444: fprintf(pfile, "-----\n"); 445: 446: //- initial line count is zero 447: 448: int iLineCount = 0; 449: 450: //- loop from start to end 451: 452: //! dump array expects char's but goes over the array as if they 453: //! are int's. 454: 455: //t fix char's int's mismatch somehow. 456: 457: int i; 458: 459: for (i = iStart; i < iEnd; ) 460: { 461: char pcOutput1[500]; 462: char pcOutput2[500]; 463: 464: //- get current command 465: 466: int iCommand = piOperators[i / sizeof(int)]; 467: 468: //- lookup info for current operand 469: 470: struct HeccerCommandInfo *phciCurrent 471: = HeccerCommandInfoLookup(phct, iCommand); 472: 473: int iCommandLength = phciCurrent->iLength; 474: 475: //- print counter 476: 477: sprintf(pcOutput1, "%5.5i ::", iLineCount); 478: 479: /* //- print numerical info */ 480: 481: /* fprintf(pcOutput2, " %i", iCommand); */ 482: 483: /* strcat(pcOutput1, pcOutput2); */ 484: 485: //- if found 486: 487: if (phciCurrent) 488: { 489: //- print name of operator 490: 491: sprintf(pcOutput2, " %s", phciCurrent->pcName); 492: 493: strcat(pcOutput1, pcOutput2); 494: 495: //- if operand length is valid 496: 497: if (iCommandLength >= 2) 498: { 499: if (phct->iFormatted) 500: { 501: //t so only for five doubles ... 502: 503: char pc[100] = ""; 504: 505: { 506: void *pv = (void *)&piOperators[i / sizeof(int) + 1]; 507: 508: if (phciCurrent->iValue == HECCER_MOP_INITIALIZECHANNEL) 509: { 510: struct MopsInitializeChannel *pmops 511: = (struct MopsInitializeChannel *)&piOperators[i / sizeof(int)]; 512: 513: sprintf(pc, " %g %g", pmops->dReversalPotential, pmops->dMaximalConductance); 514: } 515: else if (phciCurrent->iValue == HECCER_MOP_EXPONENTIALDECAY) 516: { 517: struct MopsExponentialDecay *pmops 518: = (struct MopsExponentialDecay *)&piOperators[i / sizeof(int)]; 519: 520: sprintf(pc, " %g %g %g\n\t\t\t", pmops->dBeta, pmops->dSteadyState, pmops->dTau); 521: 522: uMC2Mat *puExternal = &pmops->puExternal[0]; 523: 524: int k; 525: 526: char pc2[100]; 527: 528: for (k = 0 ; k < EXPONENTIALDECAY_CONTRIBUTORS ; k++) 529: { 530: if (puExternal[k].pdValue) 531: { 532: sprintf(pc2, " (%g)", *puExternal[k].pdValue); 533: } 534: else 535: { 536: sprintf(pc2, " (nil)"); 537: } 538: 539: strcat(pc, pc2); 540: } 541: } 542: else if (phciCurrent->iValue == HECCER_MOP_CONCEPTGATE) 543: { 544: struct MopsSingleGateConcept *pmops 545: = (struct MopsSingleGateConcept *)&piOperators[i / sizeof(int)]; 546: 547: if (pmops->uState.pdValue) 548: { 549: sprintf(pc, " %i %i (%g)", pmops->iTableIndex, pmops->iPower, *pmops->uState.pdValue); 550: } 551: else 552: { 553: sprintf(pc, " %i %i (nil)", pmops->iTableIndex, pmops->iPower); 554: } 555: } 556: else if (phciCurrent->iValue == HECCER_MOP_INTERNALNERNST) 557: { 558: struct MopsInternalNernst *pmops 559: = (struct MopsInternalNernst *)&piOperators[i / sizeof(int)]; 560: 561: sprintf(pc, " %g %g", pmops->dConstant, pmops->dExternal); 562: 563: int k; 564: 565: char pc2[100]; 566: 567: for (k = 0 ; k < 1 ; k++) 568: { 569: if (pmops->uInternal.pdValue) 570: { 571: sprintf(pc2, " (%g)", *pmops->uInternal.pdValue); 572: } 573: else 574: { 575: sprintf(pc2, " (nil)"); 576: } 577: 578: strcat(pc, pc2); 579: } 580: } 581: else if (phciCurrent->iValue == HECCER_MOP_INITIALIZECHANNELEREV) 582: { 583: struct MopsInitializeChannelErev *pmops 584: = (struct MopsInitializeChannelErev *)&piOperators[i / sizeof(int)]; 585: 586: if (pmops->uReversalPotential.pdValue) 587: { 588: sprintf(pc, " (%g)", *pmops->uReversalPotential.pdValue); 589: } 590: else 591: { 592: sprintf(pc, "(nil)"); 593: } 594: 595: char pc2[100]; 596: 597: sprintf(pc2, " %g", pmops->dMaximalConductance); 598: 599: strcat(pc, pc2); 600: } 601: else if (phciCurrent->iValue == HECCER_MOP_SPRINGMASS) 602: { 603: struct MopsSpringMass *pmops 604: = (struct MopsSpringMass *)&piOperators[i / sizeof(int)]; 605: 606: if (pmops->pdEvents 607: && pmops->pdEvents[0]) 608: { 609: sprintf(pc, " %i %g %i %i %g", pmops->iEvent, pmops->pdEvents[pmops->iEvent], pmops->iDiscreteTarget, pmops->iTable, pmops->dFrequency); 610: } 611: else 612: { 613: sprintf(pc, " %i (nil) %i %i %g", pmops->iEvent, pmops->iDiscreteTarget, pmops->iTable, pmops->dFrequency); 614: } 615: } 616: else if (phciCurrent->iValue == HECCER_MOP_EVENTGENERATE) 617: { 618: struct MopsEventGenerate *pmops 619: = (struct MopsEventGenerate *)&piOperators[i / sizeof(int)]; 620: 621: if (pmops->uSource.pdValue == (double *)-1) 622: { 623: sprintf(pc, " (%s) %g %g %i", "dVm", pmops->dThreshold, pmops->dRefractoryReset, pmops->iTable); 624: } 625: else if (pmops->uSource.pdValue) 626: { 627: sprintf(pc, " (%g) %g %g %i", *pmops->uSource.pdValue, pmops->dThreshold, pmops->dRefractoryReset, pmops->iTable); 628: } 629: else 630: { 631: sprintf(pc, " (nil) %g %g %i", pmops->dThreshold, pmops->dRefractoryReset, pmops->iTable); 632: } 633: } 634: else if (phciCurrent->iValue == HECCER_MOP_AGGREGATECURRENT) 635: { 636: struct MopsAggregateCurrent *pmops 637: = (struct MopsAggregateCurrent *)&piOperators[i / sizeof(int)]; 638: 639: sprintf(pc, " %i (%g)", pmops->iIndex, pvm->pdAggregators[pmops->iIndex]); 640: } 641: } 642: 643: sprintf(pcOutput2, "%s", pc); 644: 645: strcat(pcOutput1, pcOutput2); 646: } 647: else 648: { 649: //- print operands numerically 650: 651: int j; 652: 653: for (j = sizeof(int) ; j < iCommandLength ; j += sizeof(int)) 654: { 655: int iOperand; 656: 657: //- get current operand 658: 659: iOperand = piOperators[(i + j) / sizeof(int)]; 660: 661: sprintf(pcOutput2, " %4i", iOperand); 662: 663: strcat(pcOutput1, pcOutput2); 664: } 665: } 666: 667: //- if a secondary operands array is given 668: 669: if (pvOperands) 670: { 671: if (phciCurrent->iSecondaries) 672: { 673: sprintf(pcOutput2, "\t\t\t\t\t\t\t"); 674: 675: strcat(pcOutput1, pcOutput2); 676: 677: double *pdOperands = (double *)pvOperands; 678: 679: int iSecondary; 680: 681: for (iSecondary = 0; iSecondary < phciCurrent->iSecondaries; iSecondary++) 682: { 683: sprintf(pcOutput2, " %g", pdOperands[iSecondary]); 684: 685: strcat(pcOutput1, pcOutput2); 686: } 687: } 688: 689: pvOperands = (void *)&((char *)pvOperands)[phciCurrent->iSecondariesSize]; 690: } 691: } 692: } 693: else 694: { 695: fprintf(pfile, " Error: %4i command not found in command table", iCommand); 696: } 697: 698: //- terminate line 699: 700: fprintf(pfile, "%s\n", pcOutput1); 701: 702: //- increment line count 703: 704: iLineCount++; 705: 706: //- add size of operands 707: 708: if (phciCurrent && iCommandLength >= 2) 709: { 710: i += iCommandLength; 711: } 712: else 713: { 714: i += sizeof(int); 715: } 716: } 717: 718: //- return result 719: 720: return(iResult); 721: } 722: 723: 724: