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 <math.h> 20: #include <stdio.h> 21: #include <stdlib.h> 22: #include <string.h> 23: 24: #include "heccer/heccer.h" 25: #include "heccer/table.h" 26: 27: 28: int 29: static 30: HeccerChannelPersistentSteadyStateDualTauTabulate 31: (struct ChannelPersistentSteadyStateDualTau *pcpsdt, struct Heccer *pheccer); 32: 33: int 34: static 35: HeccerChannelPersistentSteadyStateTauTabulate 36: (struct ChannelPersistentSteadyStateTau *pcpst, struct Heccer *pheccer); 37: 38: int 39: static 40: HeccerChannelSteadyStateSteppedTauTabulate 41: (struct ChannelSteadyStateSteppedTau *pcsst, struct Heccer *pheccer); 42: 43: int 44: static 45: HeccerDiscretizeConcentrationGate 46: (struct Heccer *pheccer, struct ConcentrationActivator *pac); 47: 48: int 49: static 50: HeccerDiscretizeGateConcept 51: (struct Heccer *pheccer, struct GateConcept *pgc); 52: 53: static int 54: HeccerTableDump 55: (struct HeccerTabulatedGate *phtg, int iIndex, FILE *pfile, int iSelection); 56: 57: static 58: int 59: HeccerTabulatedGateCompareParameters 60: (struct HeccerTabulatedGate *phtg, void *pv, size_t iSize); 61: 62: static 63: int 64: HeccerTabulatedGateLookup 65: (struct Heccer *pheccer, void *pv, size_t iSize); 66: 67: static 68: int 69: HeccerTabulatedGateNew 70: (struct Heccer *pheccer, 71: void *pvParameters, 72: size_t iSize, 73: double dStart, 74: double dEnd, 75: int iEntries); 76: 77: int 78: static 79: HeccerTabulateSpringMass(struct Heccer *pheccer, struct ChannelSpringMass *pcsm); 80: 81: static 82: int 83: HeccerTabulatedSpringMassCompareParameters 84: (struct HeccerTabulatedSpringMass *phtsm, void *pv, size_t iSize); 85: 86: static 87: int 88: HeccerTabulatedSpringMassLookup 89: (struct Heccer *pheccer, void *pv, size_t iSize); 90: 91: static 92: int 93: HeccerTabulatedSpringMassNew 94: (struct Heccer *pheccer, void *pvParameters, size_t iSize); 95: 96: static 97: int 98: HeccerTabulatedGateStore 99: (struct Heccer *pheccer, 100: struct HeccerTabulatedGate *phtgNew); 101: 102: 103: /// ************************************************************************** 104: /// 105: /// SHORT: HeccerConcentrationGateTabulate() 106: /// 107: /// ARGS.: 108: /// 109: /// pca.....: a concentration gate concept. 110: /// pheccer.: a heccer. 111: /// 112: /// RTN..: int 113: /// 114: /// success of operation. 115: /// 116: /// phtg....: a filled heccer tabulated gate. 117: /// 118: /// DESCR: Fill the tables with a discretization of concentration gate 119: /// kinetics. 120: /// 121: /// ************************************************************************** 122: 123: int 124: HeccerConcentrationGateTabulate 125: (struct ConcentrationActivator *pca, struct Heccer *pheccer) 126: { 127: //- set default result : ok 128: 129: int iResult = TRUE; 130: 131: //- get access to the tabulated gate structure 132: 133: int iIndex = pca->iTable; 134: 135: struct HeccerTabulatedGate *phtg = &pheccer->tgt.phtg[iIndex]; 136: 137: //- get step size 138: 139: double dStep = phtg->hi.dStep; 140: 141: //- loop over all entries in the table 142: 143: int i; 144: double dx; 145: 146: for (dx = phtg->hi.dStart, i = 0 ; i <= phtg->iEntries ; i++, dx += dStep) 147: { 148: //- compute steady state 149: 150: double dEquilibrium = 1 / (1 + (pca->parameters.dBasalLevel / dx)); 151: 152: //- fill in A and B table 153: 154: //t perhaps the table names should be redefined ? 155: 156: //! time step normalization done elsewhere 157: 158: phtg->pdA[i] = dEquilibrium / pca->parameters.dTau; 159: 160: phtg->pdB[i] = 1 / pca->parameters.dTau; 161: } 162: 163: //- return result 164: 165: return(iResult); 166: } 167: 168: 169: /// ************************************************************************** 170: /// 171: /// SHORT: HeccerDiscretizeConcentrationGate() 172: /// 173: /// ARGS.: 174: /// 175: /// pheccer...: a heccer. 176: /// pca.......: concentration gate concept description. 177: /// 178: /// RTN..: int : success of operation. 179: /// 180: /// DESCR: Discretize the given concentration gate concept. 181: /// 182: /// ************************************************************************** 183: 184: int 185: static 186: HeccerDiscretizeConcentrationGate 187: (struct Heccer *pheccer, struct ConcentrationActivator *pca) 188: { 189: //- set default result : ok 190: 191: int iResult = TRUE; 192: 193: //- if already registered 194: 195: if (pca->iTable != -1) 196: { 197: return(TRUE); 198: } 199: 200: //- allocate structures 201: 202: double dStart = pheccer->ho.dConcentrationGateStart; 203: double dEnd = pheccer->ho.dConcentrationGateEnd; 204: int iEntries = pheccer->ho.iIntervalEntries; 205: 206: //- lookup the table parameters ... 207: 208: int i = HeccerTabulatedGateLookup(pheccer, &pca->parameters, sizeof(pca->parameters)); 209: 210: if (i == -1) 211: { 212: //- ... or create a new table for these parameters 213: 214: i = HeccerTabulatedGateNew(pheccer, &pca->parameters, sizeof(pca->parameters), dStart, dEnd, iEntries); 215: 216: if (i == -1) 217: { 218: return(FALSE); 219: } 220: 221: //- register the index 222: 223: pca->iTable = i; 224: 225: //- fill the table with the discretized gate kinetics 226: 227: iResult = iResult && HeccerConcentrationGateTabulate(pca, pheccer); 228: } 229: else 230: { 231: pca->iTable = i; 232: } 233: 234: //- return result 235: 236: return(iResult); 237: } 238: 239: 240: /// ************************************************************************** 241: /// 242: /// SHORT: HeccerDiscretizeGateConcept() 243: /// 244: /// ARGS.: 245: /// 246: /// pheccer...: a heccer. 247: /// pgc.......: gate concept description. 248: /// 249: /// RTN..: int : success of operation. 250: /// 251: /// DESCR: Discretize the given gate concept. 252: /// 253: /// Note that for a single gate, this function sets up the A 254: /// and B kinetic discretization. 255: /// 256: /// ************************************************************************** 257: 258: int 259: static 260: HeccerDiscretizeGateConcept 261: (struct Heccer *pheccer, struct GateConcept *pgc) 262: { 263: //- set default result : ok 264: 265: int iResult = TRUE; 266: 267: //- if already registered 268: 269: if (pgc->iTable != -1) 270: { 271: return(TRUE); 272: } 273: 274: //- if table is hardcoded 275: 276: if (pgc->htg.pdA 277: || pgc->htg.pdB) 278: { 279: if (!pgc->htg.pdA 280: || !pgc->htg.pdB) 281: { 282: HeccerError 283: (pheccer, 284: NULL, 285: "HeccerDiscretizeGateConcept(): if a gate has a hardcoded table for one kinetic, it must have hardcoded tables for all kinetics."); 286: 287: return(FALSE); 288: } 289: 290: //t first should do the lookup 291: 292: //- store the table as is 293: 294: int i = HeccerTabulatedGateStore(pheccer, &pgc->htg); 295: 296: if (i == -1) 297: { 298: return(FALSE); 299: } 300: 301: //- register the index 302: 303: pgc->iTable = i; 304: 305: //- return success 306: 307: return(TRUE); 308: } 309: 310: //- allocate structures 311: 312: double dStart = pheccer->ho.dIntervalStart; 313: double dEnd = pheccer->ho.dIntervalEnd; 314: int iEntries = pheccer->ho.iIntervalEntries; 315: 316: //- lookup the table parameters ... 317: 318: int i = HeccerTabulatedGateLookup(pheccer, &pgc->parameters, sizeof(pgc->parameters)); 319: 320: if (i == -1) 321: { 322: //- ... or create a new table for these parameters 323: 324: i = HeccerTabulatedGateNew(pheccer, &pgc->parameters, sizeof(pgc->parameters), dStart, dEnd, iEntries); 325: 326: if (i == -1) 327: { 328: return(FALSE); 329: } 330: 331: //- register the index 332: 333: pgc->iTable = i; 334: 335: //- fill the table with the discretized gate kinetics 336: 337: iResult = iResult && HeccerGateConceptTabulate(pgc, pheccer); 338: } 339: else 340: { 341: pgc->iTable = i; 342: } 343: 344: //- return result 345: 346: return(iResult); 347: } 348: 349: 350: /// ************************************************************************** 351: /// 352: /// SHORT: HeccerGateConceptTabulate() 353: /// 354: /// ARGS.: 355: /// 356: /// pgc.....: a heccer gate concept. 357: /// pheccer.: a heccer. 358: /// 359: /// RTN..: int 360: /// 361: /// success of operation. 362: /// 363: /// phtg....: a filled heccer tabulated gate. 364: /// 365: /// DESCR: Fill the tables with a discretization of the gate kinetics. 366: /// 367: /// ************************************************************************** 368: 369: int 370: HeccerGateConceptTabulate 371: (struct GateConcept *pgc, struct Heccer *pheccer) 372: { 373: //- set default result : ok 374: 375: int iResult = TRUE; 376: 377: //- get access to the tabulated gate structure 378: 379: int iIndex = pgc->iTable; 380: 381: struct HeccerTabulatedGate *phtg = &pheccer->tgt.phtg[iIndex]; 382: 383: //- get step size 384: 385: double dStep = phtg->hi.dStep; 386: 387: //- loop over all entries in the table 388: 389: int i; 390: double dx; 391: 392: for (dx = phtg->hi.dStart, i = 0 ; i <= phtg->iEntries ; i++, dx += dStep) 393: { 394: //v alpha carry over to beta 395: 396: double dCarryOver; 397: 398: if (phtg->pdA) 399: { 400: double dMultiplier = pgc->parameters.gkA.dHHScale; 401: double dMembraneDependence = pgc->parameters.gkA.dHHMult; 402: double dMembraneDependenceOffset = pgc->parameters.gkA.dHHOffsetM; 403: int iNominator = pgc->parameters.gkA.iHHFactorFlag; 404: double dDeNominatorOffset = pgc->parameters.gkA.dHHAdd; 405: double dMembraneOffset = pgc->parameters.gkA.dHHOffsetE; 406: double dTauDenormalizer = pgc->parameters.gkA.dHHTau; 407: 408: //t check the MCAD MMGLT macro to see how it deals with 409: //t relative errors. The current implementation is magnitude 410: //t dependent, and obviously completely add hoc. 411: 412: if (fabs(dTauDenormalizer) < 1e-17) 413: { 414: dCarryOver = 0.0; 415: phtg->pdA[i] = 0.0; 416: } 417: else 418: { 419: double dDeNominator = dDeNominatorOffset + exp((dx + dMembraneOffset) / dTauDenormalizer); 420: 421: if (fabs(dDeNominator) < 1e-17) 422: { 423: phtg->pdA[i] = phtg->pdA[i - 1]; 424: } 425: else 426: { 427: if (iNominator == 1) 428: { 429: phtg->pdA[i] = (dMultiplier + (dMembraneDependenceOffset - dMembraneDependence * dx)) * dDeNominator; 430: } 431: else if (iNominator == -1) 432: { 433: phtg->pdA[i] = (dMultiplier + (dMembraneDependenceOffset - dMembraneDependence * dx)) / dDeNominator; 434: } 435: else 436: { 437: iResult = FALSE; 438: } 439: } 440: 441: //- remember alpha, such that we can add it to beta below 442: 443: dCarryOver = phtg->pdA[i]; 444: } 445: } 446: 447: if (phtg->pdA && phtg->pdB) 448: { 449: double dMultiplier = pgc->parameters.gkB.dHHScale; 450: double dMembraneDependence = pgc->parameters.gkB.dHHMult; 451: double dMembraneDependenceOffset = pgc->parameters.gkA.dHHOffsetM; 452: int iNominator = pgc->parameters.gkB.iHHFactorFlag; 453: double dDeNominatorOffset = pgc->parameters.gkB.dHHAdd; 454: double dMembraneOffset = pgc->parameters.gkB.dHHOffsetE; 455: double dTauDenormalizer = pgc->parameters.gkB.dHHTau; 456: 457: if (fabs(dTauDenormalizer) < 1e-17) 458: { 459: phtg->pdB[i] = 0.0; 460: } 461: else 462: { 463: double dDeNominator = dDeNominatorOffset + exp((dx + dMembraneOffset) / dTauDenormalizer); 464: 465: if (fabs(dDeNominator) < 1e-17) 466: { 467: phtg->pdB[i] = phtg->pdB[i - 1]; 468: } 469: else 470: { 471: if (iNominator == 1) 472: { 473: phtg->pdB[i] = (dMultiplier + dMembraneDependenceOffset - (dMembraneDependence * dx)) * dDeNominator; 474: } 475: else if (iNominator == -1) 476: { 477: phtg->pdB[i] = (dMultiplier + dMembraneDependenceOffset - (dMembraneDependence * dx)) / dDeNominator; 478: } 479: else 480: { 481: iResult = FALSE; 482: } 483: } 484: } 485: 486: //- add alpha 487: 488: if (iResult == TRUE) 489: { 490: phtg->pdB[i] += dCarryOver; 491: } 492: } 493: 494: } 495: 496: //- return result 497: 498: return(iResult); 499: } 500: 501: 502: /// ************************************************************************** 503: /// 504: /// SHORT: HeccerChannelPersistentSteadyStateDualTauTabulate() 505: /// 506: /// ARGS.: 507: /// 508: /// pcpsdt..: a channel with steady state and two time constants. 509: /// pheccer.: a heccer. 510: /// 511: /// RTN..: int 512: /// 513: /// success of operation. 514: /// 515: /// DESCR: Fill the tables with a discretization of the gate kinetics. 516: /// 517: /// First tables with 149 entries are created using the commonly 518: /// defined formulas for steady state and time constant (see 519: /// EDS1994). Next the full sized tables are created using the 520: /// small tables, by using b-spline interpolation. 521: /// 522: /// This is the way it was done in the Purkinje cell model, but is 523: /// only here for reasons of legacy. This function should 524: /// normally not be used. 525: /// 526: /// ************************************************************************** 527: 528: int 529: static 530: HeccerChannelPersistentSteadyStateDualTauTabulate 531: (struct ChannelPersistentSteadyStateDualTau *pcpsdt, struct Heccer *pheccer) 532: { 533: //- set default result : ok 534: 535: int iResult = TRUE; 536: 537: //- if already registered 538: 539: if (pcpsdt->iFirstTable != -1) 540: { 541: return(TRUE); 542: } 543: 544: //- allocate structures 545: 546: double dStart = pheccer->ho.dIntervalStart; 547: double dEnd = pheccer->ho.dIntervalEnd; 548: int iEntries = pheccer->ho.iIntervalEntries; 549: 550: //- lookup the table parameters ... 551: 552: int i = HeccerTabulatedGateLookup(pheccer, &pcpsdt->parameters1, sizeof(pcpsdt->parameters1)); 553: 554: if (i == -1) 555: { 556: //- ... or create a new table for these parameters 557: 558: i = HeccerTabulatedGateNew(pheccer, &pcpsdt->parameters1, sizeof(pcpsdt->parameters1), dStart, dEnd, iEntries); 559: 560: if (i == -1) 561: { 562: return(FALSE); 563: } 564: 565: //- register the index 566: 567: pcpsdt->iFirstTable = i; 568: 569: //- get access to the tabulated gate structures 570: 571: int iIndex1 = pcpsdt->iFirstTable; 572: 573: struct HeccerTabulatedGate *phtg1 = &pheccer->tgt.phtg[iIndex1]; 574: 575: //- allocate the small tables 576: 577: int iSmallTableSize = pheccer->ho.iSmallTableSize; 578: 579: //20; 580: 581: //- get step size 582: 583: double dSmallStep = (phtg1->hi.dEnd - phtg1->hi.dStart) / (iSmallTableSize); 584: 585: double *pd1 = (double *)calloc(iSmallTableSize + 1, sizeof(double)); 586: double *pd2 = (double *)calloc(iSmallTableSize + 1, sizeof(double)); 587: 588: //- loop over all entries in the first small table 589: 590: double dx; 591: 592: for (dx = phtg1->hi.dStart, i = 0 ; i <= iSmallTableSize ; i++, dx += dSmallStep) 593: { 594: double dA = pcpsdt->parameters1.dSteadyState; 595: 596: double dC 597: = (pcpsdt->parameters1.tau.dMultiplier 598: / (pcpsdt->parameters1.tau.dDeNominatorOffset 599: + exp((dx + pcpsdt->parameters1.tau.dMembraneOffset) 600: / pcpsdt->parameters1.tau.dTauDenormalizer))); 601: 602: double dA1 = dA; 603: double dB1 = dC; 604: 605: //t check the MCAD MMGLT macro to see how it deals with 606: //t relative errors. The current implementation is magnitude 607: //t dependent, and obviously completely add hoc. 608: 609: if (fabs(dA1) < 1e-17) 610: { 611: if (dA1 < 0.0) 612: { 613: dA1 = -1e-17; 614: } 615: else 616: { 617: dA1 = 1e-17; 618: } 619: } 620: 621: pd1[i] = dB1 / dA1; 622: pd2[i] = 1.0 / dA1; 623: } 624: 625: //- interpolate the tables 626: 627: double *ppdSources[] 628: = { pd1, pd2, NULL, }; 629: double *ppdDestinations[] 630: = { phtg1->pdA, phtg1->pdB, NULL, }; 631: 632: iResult = iResult && HeccerTableInterpolate(ppdSources, ppdDestinations, iSmallTableSize, iEntries); 633: 634: //- free allocated memory 635: 636: free(pd1); 637: free(pd2); 638: } 639: else 640: { 641: pcpsdt->iFirstTable = i; 642: } 643: 644: //- lookup the table parameters ... 645: 646: int j = HeccerTabulatedGateLookup(pheccer, &pcpsdt->parameters2, sizeof(pcpsdt->parameters2)); 647: 648: if (j == -1) 649: { 650: //- ... or create a new table for these parameters 651: 652: j = HeccerTabulatedGateNew(pheccer, &pcpsdt->parameters2, sizeof(pcpsdt->parameters2), dStart, dEnd, iEntries); 653: 654: if (j == -1) 655: { 656: return(FALSE); 657: } 658: 659: //- register the index 660: 661: pcpsdt->iSecondTable = j; 662: 663: int iIndex2 = pcpsdt->iSecondTable; 664: 665: struct HeccerTabulatedGate *phtg2 = &pheccer->tgt.phtg[iIndex2]; 666: 667: //- allocate the small tables 668: 669: int iSmallTableSize = pheccer->ho.iSmallTableSize; 670: 671: //20; 672: 673: //- get step size 674: 675: double dSmallStep = (phtg2->hi.dEnd - phtg2->hi.dStart) / (iSmallTableSize); 676: 677: double *pd3 = (double *)calloc(iSmallTableSize + 1, sizeof(double)); 678: double *pd4 = (double *)calloc(iSmallTableSize + 1, sizeof(double)); 679: 680: //- loop over all entries in the first small table 681: 682: double dx; 683: 684: for (dx = phtg2->hi.dStart, i = 0 ; i <= iSmallTableSize ; i++, dx += dSmallStep) 685: { 686: double dB = pcpsdt->parameters2.dSteadyState; 687: 688: double dD 689: = (pcpsdt->parameters2.tau.dMultiplier 690: / (pcpsdt->parameters2.tau.dDeNominatorOffset 691: + exp((dx + pcpsdt->parameters2.tau.dMembraneOffset) 692: / pcpsdt->parameters2.tau.dTauDenormalizer))); 693: 694: double dA2 = dB; 695: double dB2 = dD; 696: 697: //t check the MCAD MMGLT macro to see how it deals with 698: //t relative errors. The current implementation is magnitude 699: //t dependent, and obviously completely add hoc. 700: 701: if (fabs(dA2) < 1e-17) 702: { 703: if (dA2 < 0.0) 704: { 705: dA2 = -1e-17; 706: } 707: else 708: { 709: dA2 = 1e-17; 710: } 711: } 712: 713: pd3[i] = dB2 / dA2; 714: pd4[i] = 1.0 / dA2; 715: } 716: 717: //- interpolate the tables 718: 719: double *ppdSources[] 720: = { pd3, pd4, NULL, }; 721: double *ppdDestinations[] 722: = { phtg2->pdA, phtg2->pdB, NULL, }; 723: 724: iResult = iResult && HeccerTableInterpolate(ppdSources, ppdDestinations, iSmallTableSize, iEntries); 725: 726: //- free allocated memory 727: 728: free(pd3); 729: free(pd4); 730: } 731: else 732: { 733: pcpsdt->iSecondTable = j; 734: } 735: 736: //- return result 737: 738: return(iResult); 739: } 740: 741: 742: /// ************************************************************************** 743: /// 744: /// SHORT: HeccerChannelPersistentSteadyStateTauTabulate() 745: /// 746: /// ARGS.: 747: /// 748: /// pcpst...: a heccer channel with steady state and time constant. 749: /// pheccer.: a heccer. 750: /// 751: /// RTN..: int 752: /// 753: /// success of operation. 754: /// 755: /// DESCR: Fill the tables with a discretization of the gate kinetics. 756: /// 757: /// First tables with 149 entries are created using the commonly 758: /// defined formulas for steady state and time constant (see 759: /// EDS1994). Next the full sized tables are created using the 760: /// small tables, by using b-spline interpolation. 761: /// 762: /// This is the way it was done in the Purkinje cell model, but is 763: /// only here for reasons of legacy. This function should 764: /// normally not be used. 765: /// 766: /// ************************************************************************** 767: 768: int 769: static 770: HeccerChannelPersistentSteadyStateTauTabulate 771: (struct ChannelPersistentSteadyStateTau *pcpst, struct Heccer *pheccer) 772: { 773: //- set default result : ok 774: 775: int iResult = TRUE; 776: 777: //- if already registered 778: 779: if (pcpst->iTable != -1) 780: { 781: return(TRUE); 782: } 783: 784: //- allocate structures 785: 786: double dStart = pheccer->ho.dIntervalStart; 787: double dEnd = pheccer->ho.dIntervalEnd; 788: int iEntries = pheccer->ho.iIntervalEntries; 789: 790: //- lookup the table parameters ... 791: 792: int i = HeccerTabulatedGateLookup(pheccer, &pcpst->parameters, sizeof(pcpst->parameters)); 793: 794: if (i == -1) 795: { 796: //- ... or create a new table for these parameters 797: 798: i = HeccerTabulatedGateNew(pheccer, &pcpst->parameters, sizeof(pcpst->parameters), dStart, dEnd, iEntries); 799: 800: if (i == -1) 801: { 802: return(FALSE); 803: } 804: 805: //- register the index 806: 807: pcpst->iTable = i; 808: 809: //- get access to the tabulated gate structures 810: 811: int iIndex = pcpst->iTable; 812: 813: struct HeccerTabulatedGate *phtg = &pheccer->tgt.phtg[iIndex]; 814: 815: //- allocate the small tables 816: 817: int iSmallTableSize = pheccer->ho.iSmallTableSize; 818: 819: //- get step size 820: 821: double dSmallStep = (phtg->hi.dEnd - phtg->hi.dStart) / (iSmallTableSize); 822: 823: double *pd1 = (double *)calloc(iSmallTableSize + 1, sizeof(double)); 824: double *pd2 = (double *)calloc(iSmallTableSize + 1, sizeof(double)); 825: 826: //- loop over all entries in the first small table 827: 828: double dx; 829: 830: for (dx = phtg->hi.dStart, i = 0 ; i <= iSmallTableSize ; i++, dx += dSmallStep) 831: { 832: double dA 833: = (pcpst->parameters.ss.dNominator 834: / ( pcpst->parameters.ss.dMultiplier1 835: * ((exp((dx + pcpst->parameters.ss.dMembraneOffset1) / pcpst->parameters.ss.dTauDenormalizer1))) 836: + pcpst->parameters.ss.dMultiplier2 837: * (exp ((dx + pcpst->parameters.ss.dMembraneOffset2) / pcpst->parameters.ss.dTauDenormalizer2)))); 838: 839: double dB 840: = (pcpst->parameters.tc.dNominator 841: / (pcpst->parameters.tc.dDeNominatorOffset 842: + (exp((dx + pcpst->parameters.tc.dMembraneOffset) / pcpst->parameters.tc.dTauDenormalizer)))); 843: 844: double dA3 = dA; 845: 846: double dB3 = dB; 847: 848: //t check the MCAD MMGLT macro to see how it deals with 849: //t relative errors. The current implementation is magnitude 850: //t dependent, and obviously completely add hoc. 851: 852: if (fabs(dA3) < 1e-17) 853: { 854: if (dA3 < 0.0) 855: { 856: dA3 = -1e-17; 857: } 858: else 859: { 860: dA3 = 1e-17; 861: } 862: } 863: 864: pd1[i] = dB3 / dA3; 865: pd2[i] = 1.0 / dA3; 866: } 867: 868: //- interpolate the tables 869: 870: double *ppdSources[] 871: = { pd1, pd2, NULL, }; 872: double *ppdDestinations[] 873: = { phtg->pdA, phtg->pdB, NULL, }; 874: 875: iResult = iResult && HeccerTableInterpolate(ppdSources, ppdDestinations, iSmallTableSize, iEntries); 876: 877: //- free allocated memory 878: 879: free(pd1); 880: free(pd2); 881: } 882: else 883: { 884: pcpst->iTable = i; 885: } 886: 887: //- return result 888: 889: return(iResult); 890: } 891: 892: 893: /// ************************************************************************** 894: /// 895: /// SHORT: HeccerTabulateAny() 896: /// 897: /// ARGS.: 898: /// 899: /// pheccer.: a heccer. 900: /// iType...: type of intermediary. 901: /// pv......: a heccer intermediary with table parameters. 902: /// 903: /// RTN..: int 904: /// 905: /// success of operation. 906: /// 907: /// DESCR: Create the table for a channel. 908: /// 909: /// ************************************************************************** 910: 911: int 912: HeccerTabulateAny 913: (struct Heccer *pheccer, void *pv, int iType) 914: { 915: //- set default result: failure 916: 917: int iResult = FALSE; 918: 919: //- tabulate depending on tables 920: 921: switch (iType) 922: { 923: case MATH_TYPE_Concentration: 924: { 925: struct ConcentrationActivator *pca = (struct ConcentrationActivator *)pv; 926: 927: iResult = HeccerDiscretizeConcentrationGate(pheccer, pca); 928: 929: break; 930: } 931: case MATH_TYPE_GateConcept: 932: { 933: struct GateConcept *pgc = (struct GateConcept *)pv; 934: 935: iResult = HeccerDiscretizeGateConcept(pheccer, pgc); 936: 937: break; 938: } 939: case MATH_TYPE_ChannelPersistentSteadyStateDualTau: 940: { 941: struct ChannelPersistentSteadyStateDualTau *pcpsdt = (struct ChannelPersistentSteadyStateDualTau *)pv; 942: 943: iResult = HeccerChannelPersistentSteadyStateDualTauTabulate(pcpsdt, pheccer); 944: 945: break; 946: } 947: case MATH_TYPE_ChannelPersistentSteadyStateTau: 948: { 949: struct ChannelPersistentSteadyStateTau *pcpst = (struct ChannelPersistentSteadyStateTau *)pv; 950: 951: iResult = HeccerChannelPersistentSteadyStateTauTabulate(pcpst, pheccer); 952: 953: break; 954: } 955: case MATH_TYPE_ChannelSpringMass: 956: { 957: struct ChannelSpringMass *pcsm = (struct ChannelSpringMass *)pv; 958: 959: iResult = HeccerTabulateSpringMass(pheccer, pcsm); 960: 961: break; 962: } 963: case MATH_TYPE_ChannelSteadyStateSteppedTau: 964: { 965: struct ChannelSteadyStateSteppedTau *pcsst = (struct ChannelSteadyStateSteppedTau *)pv; 966: 967: iResult = HeccerChannelSteadyStateSteppedTauTabulate(pcsst, pheccer); 968: 969: break; 970: } 971: } 972: 973: //- return result 974: 975: return(iResult); 976: } 977: 978: 979: /// ************************************************************************** 980: /// 981: /// SHORT: HeccerTabulateSpringMass() 982: /// 983: /// ARGS.: 984: /// 985: /// pheccer.: a heccer. 986: /// pcsm....: a heccer channel with steady state and time constant. 987: /// 988: /// RTN..: int 989: /// 990: /// success of operation. 991: /// 992: /// DESCR: Create the table for a spring mass channel. 993: /// 994: /// ************************************************************************** 995: 996: int 997: static 998: HeccerTabulateSpringMass(struct Heccer *pheccer, struct ChannelSpringMass *pcsm) 999: { 1000: //- set default result : ok 1001: 1002: int iResult = TRUE; 1003: 1004: //- if already registered 1005: 1006: if (pcsm->iTable != -1) 1007: { 1008: return(TRUE); 1009: } 1010: 1011: //- lookup the table parameters ... 1012: 1013: int i = HeccerTabulatedSpringMassLookup(pheccer, &pcsm->parameters, sizeof(pcsm->parameters)); 1014: 1015: if (i == -1) 1016: { 1017: //- ... or create a new table for these parameters 1018: 1019: i = HeccerTabulatedSpringMassNew(pheccer, &pcsm->parameters, sizeof(pcsm->parameters)); 1020: 1021: if (i == -1) 1022: { 1023: return(FALSE); 1024: } 1025: 1026: //- register the index 1027: 1028: pcsm->iTable = i; 1029: 1030: //- get access to the tabulated spring mass structure 1031: 1032: struct HeccerTabulatedSpringMass *phtsm = &pheccer->tsmt.phtsm[i]; 1033: 1034: //m compute time constants derived coefficients 1035: 1036: //! we normalize the first coefficient over the time step, 1037: //! propagates through all equations. 1038: 1039: phtsm->dX1 = (pcsm->parameters.dTau1 * (1.0 - exp(- pheccer->dStep / pcsm->parameters.dTau1))) / pheccer->dStep; 1040: phtsm->dX2 = exp(- pheccer->dStep / pcsm->parameters.dTau1); 1041: 1042: phtsm->dY1 = pcsm->parameters.dTau2 * (1.0 - exp(- pheccer->dStep / pcsm->parameters.dTau2)); 1043: phtsm->dY2 = exp(- pheccer->dStep / pcsm->parameters.dTau2); 1044: 1045: } 1046: else 1047: { 1048: pcsm->iTable = i; 1049: } 1050: 1051: //- return result 1052: 1053: return(iResult); 1054: } 1055: 1056: 1057: /// ************************************************************************** 1058: /// 1059: /// SHORT: HeccerChannelSteadyStateSteppedTauTabulate() 1060: /// 1061: /// ARGS.: 1062: /// 1063: /// pheccer.: a heccer. 1064: /// pcsst...: a heccer channel with steady state and time constant. 1065: /// 1066: /// RTN..: int 1067: /// 1068: /// success of operation. 1069: /// 1070: /// DESCR: Fill the tables with a discretization of the gate kinetics. 1071: /// 1072: /// First tables with 149 entries are created using the commonly 1073: /// defined formulas for steady state and time constant (see 1074: /// EDS1994). Next the full sized tables are created using the 1075: /// small tables, by using b-spline interpolation. 1076: /// 1077: /// This is the way it was done in the Purkinje cell model, but is 1078: /// only here for reasons of legacy. This function should 1079: /// normally not be used. 1080: /// 1081: /// ************************************************************************** 1082: 1083: int 1084: static 1085: HeccerChannelSteadyStateSteppedTauTabulate 1086: (struct ChannelSteadyStateSteppedTau *pcsst, struct Heccer *pheccer) 1087: { 1088: //- set default result : ok 1089: 1090: int iResult = TRUE; 1091: 1092: //- if already registered 1093: 1094: if (pcsst->iFirstTable != -1) 1095: { 1096: return(TRUE); 1097: } 1098: 1099: //- allocate structures 1100: 1101: double dStart = pheccer->ho.dIntervalStart; 1102: double dEnd = pheccer->ho.dIntervalEnd; 1103: int iEntries = pheccer->ho.iIntervalEntries; 1104: 1105: //- lookup the table parameters ... 1106: 1107: int i = HeccerTabulatedGateLookup(pheccer, &pcsst->ss_parameters, sizeof(pcsst->ss_parameters)); 1108: 1109: if (i == -1) 1110: { 1111: //- ... or create a new table for these parameters 1112: 1113: i = HeccerTabulatedGateNew(pheccer, &pcsst->ss_parameters, sizeof(pcsst->ss_parameters), dStart, dEnd, iEntries); 1114: 1115: if (i == -1) 1116: { 1117: return(FALSE); 1118: } 1119: 1120: //- register the index 1121: 1122: pcsst->iFirstTable = i; 1123: 1124: //- get access to the tabulated gate structures 1125: 1126: int iIndex1 = pcsst->iFirstTable; 1127: 1128: struct HeccerTabulatedGate *phtg1 = &pheccer->tgt.phtg[iIndex1]; 1129: 1130: //- allocate the small tables 1131: 1132: int iSmallTableSize = pheccer->ho.iSmallTableSize; 1133: 1134: //- get step size 1135: 1136: double dSmallStep = (phtg1->hi.dEnd - phtg1->hi.dStart) / (iSmallTableSize); 1137: 1138: double *pd1 = (double *)calloc(iSmallTableSize + 1, sizeof(double)); 1139: double *pd2 = (double *)calloc(iSmallTableSize + 1, sizeof(double)); 1140: 1141: //- loop over all entries in the first small table 1142: 1143: double dx; 1144: 1145: for (dx = phtg1->hi.dStart, i = 0 ; i <= iSmallTableSize ; i++, dx += dSmallStep) 1146: { 1147: double dA 1148: = (pcsst->ss_parameters.first.a.dMultiplier 1149: * (dx + pcsst->ss_parameters.first.a.dMembraneDependenceOffset) 1150: / ((exp((dx + pcsst->ss_parameters.first.a.dMembraneOffset) 1151: / pcsst->ss_parameters.first.a.dTauDenormalizer)) 1152: + pcsst->ss_parameters.first.a.dDeNominatorOffset)); 1153: 1154: double dB 1155: = (pcsst->ss_parameters.first.b.dMultiplier 1156: * (exp( - (dx + pcsst->ss_parameters.first.b.dMembraneDependenceOffset) 1157: / pcsst->ss_parameters.first.b.dTauDenormalizer))); 1158: 1159: double dC 1160: = (pcsst->ss_parameters.second.a.dMultiplier 1161: * (dx + pcsst->ss_parameters.second.a.dMembraneDependenceOffset) 1162: / ((exp((dx + pcsst->ss_parameters.second.a.dMembraneOffset) 1163: / pcsst->ss_parameters.second.a.dTauDenormalizer)) 1164: + pcsst->ss_parameters.second.a.dDeNominatorOffset)); 1165: 1166: double dD 1167: = (pcsst->ss_parameters.second.b.dMultiplier 1168: * (exp( - (dx + pcsst->ss_parameters.second.b.dMembraneDependenceOffset) 1169: / pcsst->ss_parameters.second.b.dTauDenormalizer))); 1170: 1171: double dA4 = (1.0 / (dA + dB)); 1172: 1173: double dB4 = (dC / (dC + dD)); 1174: 1175: //t check the MCAD MMGLT macro to see how it deals with 1176: //t relative errors. The current implementation is magnitude 1177: //t dependent, and obviously completely add hoc. 1178: 1179: if (fabs(dA4) < 1e-17) 1180: { 1181: if (dA4 < 0.0) 1182: { 1183: dA4 = -1e-17; 1184: } 1185: else 1186: { 1187: dA4 = 1e-17; 1188: } 1189: } 1190: 1191: pd1[i] = dB4 / dA4; 1192: pd2[i] = 1.0 / dA4; 1193: } 1194: 1195: //- interpolate the tables 1196: 1197: double *ppdSources[] 1198: = { pd1, pd2, NULL, }; 1199: double *ppdDestinations[] 1200: = { phtg1->pdA, phtg1->pdB, NULL, }; 1201: 1202: iResult = iResult && HeccerTableInterpolate(ppdSources, ppdDestinations, iSmallTableSize, iEntries); 1203: 1204: //- free allocated memory 1205: 1206: free(pd1); 1207: free(pd2); 1208: } 1209: else 1210: { 1211: pcsst->iFirstTable = i; 1212: } 1213: 1214: //- lookup the table parameters ... 1215: 1216: int j = HeccerTabulatedGateLookup(pheccer, &pcsst->tc_parameters, sizeof(pcsst->tc_parameters)); 1217: 1218: if (j == -1) 1219: { 1220: //- ... or create a new table for these parameters 1221: 1222: j = HeccerTabulatedGateNew(pheccer, &pcsst->tc_parameters, sizeof(pcsst->tc_parameters), dStart, dEnd, iEntries); 1223: 1224: if (j == -1) 1225: { 1226: return(FALSE); 1227: } 1228: 1229: //- register the index 1230: 1231: pcsst->iSecondTable = j; 1232: 1233: //- get access to the tabulated gate structures 1234: 1235: int iIndex2 = pcsst->iSecondTable; 1236: 1237: struct HeccerTabulatedGate *phtg2 = &pheccer->tgt.phtg[iIndex2]; 1238: 1239: //- allocate the small tables 1240: 1241: int iSmallTableSize = pheccer->ho.iSmallTableSize; 1242: 1243: //- get step size 1244: 1245: double dSmallStep = (phtg2->hi.dEnd - phtg2->hi.dStart) / (iSmallTableSize); 1246: 1247: double *pd3 = (double *)calloc(iSmallTableSize + 1, sizeof(double)); 1248: double *pd4 = (double *)calloc(iSmallTableSize + 1, sizeof(double)); 1249: 1250: //- loop over all entries in the first small table 1251: 1252: double dx; 1253: 1254: for (dx = phtg2->hi.dStart, i = 0 ; i <= iSmallTableSize ; i++, dx += dSmallStep) 1255: { 1256: double dA; 1257: 1258: if (dx < pcsst->tc_parameters.a.dThreshold) 1259: { 1260: dA = pcsst->tc_parameters.a.dLowTarget; 1261: } 1262: else 1263: { 1264: dA = pcsst->tc_parameters.a.dHighTarget; 1265: } 1266: 1267: double dY 1268: = (pcsst->tc_parameters.b.dDeNominatorOffset 1269: + (exp((dx + pcsst->tc_parameters.b.dMembraneOffset) 1270: / pcsst->tc_parameters.b.dTauDenormalizer))); 1271: 1272: double dB = (1.0 / dY); 1273: 1274: //t check the MCAD MMGLT macro to see how it deals with 1275: //t relative errors. The current implementation is magnitude 1276: //t dependent, and obviously completely add hoc. 1277: 1278: if (fabs(dA) < 1e-17) 1279: { 1280: if (dA < 0.0) 1281: { 1282: dA = -1e-17; 1283: } 1284: else 1285: { 1286: dA = 1e-17; 1287: } 1288: } 1289: 1290: pd3[i] = dB / dA; 1291: pd4[i] = 1.0 / dA; 1292: } 1293: 1294: //- interpolate the tables 1295: 1296: double *ppdSources[] 1297: = { pd3, pd4, NULL, }; 1298: double *ppdDestinations[] 1299: = { phtg2->pdA, phtg2->pdB, NULL, }; 1300: 1301: iResult = iResult && HeccerTableInterpolate(ppdSources, ppdDestinations, iSmallTableSize, iEntries); 1302: 1303: //- free allocated memory 1304: 1305: free(pd3); 1306: free(pd4); 1307: } 1308: else 1309: { 1310: pcsst->iSecondTable = j; 1311: } 1312: 1313: //- return result 1314: 1315: return(iResult); 1316: } 1317: 1318: 1319: /// ************************************************************************** 1320: /// 1321: /// SHORT: HeccerTableDump() 1322: /// 1323: /// ARGS.: 1324: /// 1325: /// phtg.......: a tabulated gate. 1326: /// iIndex.....: index of this gate. 1327: /// pfile......: stdio file. 1328: /// iSelection.: selection to dump. 1329: /// 1330: /// RTN..: int 1331: /// 1332: /// success of operation. 1333: /// 1334: /// DESCR: Dump tables to given stream, respecting given selection. 1335: /// 1336: /// ************************************************************************** 1337: 1338: static int 1339: HeccerTableDump 1340: (struct HeccerTabulatedGate *phtg, int iIndex, FILE *pfile, int iSelection) 1341: { 1342: //- set default result 1343: 1344: int iResult = TRUE; 1345: 1346: if (iSelection & HECCER_DUMP_TABLE_GATE_SUMMARY) 1347: { 1348: fprintf(pfile, "Tabulated gate %i, interval (dStart) : (%g)\n", iIndex, phtg->hi.dStart); 1349: fprintf(pfile, "Tabulated gate %i, interval (dEnd) : (%g)\n", iIndex, phtg->hi.dEnd); 1350: fprintf(pfile, "Tabulated gate %i, interval (dStep) : (%g)\n", iIndex, phtg->hi.dStep); 1351: } 1352: 1353: if (iSelection & HECCER_DUMP_TABLE_GATE_SUMMARY) 1354: { 1355: fprintf(pfile, "Tabulated gate %i, interpolation (iShape) : (%i)\n", iIndex, phtg->htao.iShape); 1356: } 1357: 1358: if (iSelection & HECCER_DUMP_TABLE_GATE_SUMMARY) 1359: { 1360: fprintf(pfile, "Tabulated gate %i, (iEntries) : (%i)\n", iIndex, phtg->iEntries); 1361: } 1362: 1363: if (iSelection & HECCER_DUMP_TABLE_GATE_TABLES) 1364: { 1365: int i; 1366: 1367: for (i = 0 ; i < phtg->iEntries ; i++) 1368: { 1369: fprintf(pfile, "Tabulated gate %i, A (iEntry %i) : (%g)\n", iIndex, i, phtg->pdA[i]); 1370: } 1371: 1372: for (i = 0 ; i < phtg->iEntries ; i++) 1373: { 1374: fprintf(pfile, "Tabulated gate %i, B (iEntry %i) : (%g)\n", iIndex, i, phtg->pdB[i]); 1375: } 1376: } 1377: 1378: //- return result 1379: 1380: return(iResult); 1381: } 1382: 1383: 1384: /// ************************************************************************** 1385: /// 1386: /// SHORT: HeccerTableInterpolate() 1387: /// 1388: /// ARGS.: 1389: /// 1390: /// ppdSources.......: source tables for interpolation. 1391: /// ppdDestinations..: destination tables for interpolation. 1392: /// iSourceSize......: size of source tables. 1393: /// iDestinationSize.: size of destination tables. 1394: /// 1395: /// RTN..: int 1396: /// 1397: /// success of operation. 1398: /// 1399: /// DESCR: Interpolate given tables, uses bezier curves. 1400: /// 1401: /// ppdSources and ppdDestinations must be NULL terminated. 1402: /// ppdDestinations must be preallocated, presumably using 1403: /// HeccerTabulatedGateNew(). 1404: /// 1405: /// ************************************************************************** 1406: 1407: int 1408: HeccerTableInterpolate 1409: (double *ppdSources[], 1410: double *ppdDestinations[], 1411: int iSourceSize, 1412: int iDestinationSize) 1413: { 1414: //- set default result : ok 1415: 1416: int iResult = TRUE; 1417: 1418: //- loop over all gates 1419: 1420: int iGate; 1421: 1422: for (iGate = 0 ; ppdSources[iGate] ; iGate++) 1423: { 1424: //! modified and optimized for heccer from 1425: //! http://people.scs.fsu.edu/~burkardt/cpp_src/spline/spline.C 1426: 1427: double dRangeStep = (double)iSourceSize / (double)iDestinationSize; 1428: double dActual; 1429: 1430: double dNSA = 1.0 / 6.0; 1431: double dNSB = 2.0 / 3.0; 1432: 1433: //- loop over the entries of this gate 1434: 1435: double *pdSource = ppdSources[iGate]; 1436: double *pdDestination = ppdDestinations[iGate]; 1437: 1438: int iSource; 1439: int iDestination; 1440: 1441: //- fill the destination till the first element that comes from the source 1442: 1443: for (dActual = 0.0, iDestination = 0, iSource = 0 ; iSource <= 1 ; dActual += dRangeStep, iSource = (int)dActual, iDestination++) 1444: { 1445: pdDestination[iDestination] 1446: = ((dActual - (double)iSource) 1447: * (pdSource[iSource + 1] - pdSource[iSource]) 1448: + pdSource[iSource]); 1449: } 1450: 1451: //- do the interpolation 1452: 1453: for( ; iSource <= iSourceSize - 2 ; dActual += dRangeStep, iSource = (int)dActual, iDestination++) 1454: { 1455: double dSource = dActual - (double)iSource; 1456: double dSource2 = dSource / 2.0; 1457: double dSource4 = dSource * dSource; 1458: double dSource42 = dSource4 / 2.0; 1459: double dSource8 = dSource4 * dSource; 1460: double dSource82 = dSource8 / 2.0; 1461: 1462: double dWeight1 = - dNSA * dSource8 + dSource42 - dSource2 + dNSA; 1463: double dWeight2 = dSource82 - dSource4 + dNSB; 1464: double dWeight3 = - dSource82 + dSource42 + dSource2 + dNSA; 1465: double dWeight4 = dNSA * dSource8; 1466: 1467: pdDestination[iDestination] 1468: = (dWeight1 * pdSource[iSource - 1] 1469: + dWeight2 * pdSource[iSource] 1470: + dWeight3 * pdSource[iSource + 1] + dWeight4 * pdSource[iSource + 2]); 1471: } 1472: 1473: //- fill the destination till the end 1474: 1475: for (; iDestination <= iDestinationSize ; dActual += dRangeStep, iSource = (int)dActual, iDestination++) 1476: { 1477: if (iSource < iSourceSize) 1478: { 1479: pdDestination[iDestination] 1480: = ((dActual - (double)iSource) 1481: * (pdSource[iSource + 1] - pdSource[iSource]) + pdSource[iSource]); 1482: } 1483: else 1484: { 1485: pdDestination[iDestination] = pdSource[iSourceSize]; 1486: } 1487: } 1488: } 1489: 1490: //- return result 1491: 1492: return(iResult); 1493: } 1494: 1495: 1496: /// ************************************************************************** 1497: /// 1498: /// SHORT: HeccerTablesDump() 1499: /// 1500: /// ARGS.: 1501: /// 1502: /// ptgt.......: tabulated gate table. 1503: /// pfile......: stdio file. 1504: /// iSelection.: selection to dump. 1505: /// 1506: /// RTN..: int 1507: /// 1508: /// success of operation. 1509: /// 1510: /// DESCR: Dump intermediary functions. 1511: /// 1512: /// ************************************************************************** 1513: 1514: int 1515: HeccerTablesDump 1516: (struct TabulatedGateTable *ptgt, FILE *pfile, int iSelection) 1517: { 1518: //- set default result : ok 1519: 1520: int iResult = TRUE; 1521: 1522: //- number of tables 1523: 1524: fprintf(pfile, "Tables (iTabulatedGateCount) : (%i)\n", ptgt->iTabulatedGateCount); 1525: 1526: //- loop over tables 1527: 1528: int i; 1529: 1530: for (i = 0 ; i < ptgt->iTabulatedGateCount ; i++) 1531: { 1532: //- dump compartment 1533: 1534: iResult = iResult && HeccerTableDump(&ptgt->phtg[i], i, pfile, iSelection); 1535: } 1536: 1537: //- return result 1538: 1539: return(iResult); 1540: } 1541: 1542: 1543: /// ************************************************************************** 1544: /// 1545: /// SHORT: HeccerTabulatedGateCompareParameters() 1546: /// 1547: /// ARGS.: 1548: /// 1549: /// phtg..: an initialized gate table. 1550: /// pv....: parameters to use for comparison. 1551: /// iSize.: size of parameter block. 1552: /// 1553: /// RTN..: int 1554: /// 1555: /// See memcmp() manual. 1556: /// 1557: /// DESCR: 1558: /// 1559: /// Compare parameters for a candidate gate with the parameters 1560: /// of an existing table. 1561: /// 1562: /// ************************************************************************** 1563: 1564: static 1565: int 1566: HeccerTabulatedGateCompareParameters 1567: (struct HeccerTabulatedGate *phtg, void *pv, size_t iSize) 1568: { 1569: //- set default result : match 1570: 1571: int iResult = 0; 1572: 1573: //- set result : compare memory regions, using smallest size 1574: 1575: iSize = iSize < phtg->iSizeParameters ? iSize : phtg->iSizeParameters; 1576: 1577: iResult = memcmp(phtg->pvParameters, pv, iSize); 1578: 1579: //- return result 1580: 1581: return(iResult); 1582: } 1583: 1584: 1585: /// ************************************************************************** 1586: /// 1587: /// SHORT: HeccerTabulatedGateLookup() 1588: /// 1589: /// ARGS.: 1590: /// 1591: /// pheccer...: a heccer. 1592: /// pv........: table parameter block. 1593: /// iSize.....: size of parameter block. 1594: /// 1595: /// RTN..: int 1596: /// 1597: /// tabulated gate index, -1 for failure. 1598: /// 1599: /// DESCR: Lookup an existing table. 1600: /// 1601: /// ************************************************************************** 1602: 1603: static 1604: int 1605: HeccerTabulatedGateLookup 1606: (struct Heccer *pheccer, void *pv, size_t iSize) 1607: { 1608: //- set default result : not found 1609: 1610: int iResult = -1; 1611: 1612: //! a protection for the case where you accidently forget to dereference 1613: 1614: if (iSize < 10) 1615: { 1616: ((int *)0)[0] = 0; 1617: } 1618: 1619: //- loop over all tables 1620: 1621: int i; 1622: 1623: for (i = 0 ; i < pheccer->tgt.iTabulatedGateCount; i++) 1624: { 1625: //- get pointer to current table 1626: 1627: struct HeccerTabulatedGate *phtg = &pheccer->tgt.phtg[i]; 1628: 1629: //- if match 1630: 1631: if (HeccerTabulatedGateCompareParameters(phtg, pv, iSize) == 0) 1632: { 1633: //- set table index 1634: 1635: iResult = i; 1636: 1637: //- break searching loop 1638: 1639: break; 1640: } 1641: } 1642: 1643: //- return result 1644: 1645: return(iResult); 1646: } 1647: 1648: 1649: /// ************************************************************************** 1650: /// 1651: /// SHORT: HeccerTabulatedGateNew() 1652: /// 1653: /// ARGS.: 1654: /// 1655: /// pheccer...: a heccer. 1656: /// dStart....: start value for table. 1657: /// dEnd......: end value for table. 1658: /// iEntries..: number of entries in the table. 1659: /// 1660: /// RTN..: int 1661: /// 1662: /// tabulated gate index, -1 for failure. 1663: /// 1664: /// DESCR: Allocate a new table. 1665: /// 1666: /// ************************************************************************** 1667: 1668: static 1669: int 1670: HeccerTabulatedGateNew 1671: (struct Heccer *pheccer, 1672: void *pvParameters, 1673: size_t iSize, 1674: double dStart, 1675: double dEnd, 1676: int iEntries) 1677: { 1678: struct HeccerTabulatedGate *phtg = NULL; 1679: 1680: if (pheccer->tgt.iTabulatedGateCount >= HECCER_TABULATED_GATES_MAX) 1681: { 1682: return(-1); 1683: } 1684: 1685: #define HECCER_STATIC_TABULATED_GATES 1686: #ifdef HECCER_STATIC_TABULATED_GATES 1687: 1688: //- set result : from pool 1689: 1690: phtg = &pheccer->tgt.phtg[pheccer->tgt.iTabulatedGateCount]; 1691: 1692: #else 1693: 1694: //- set result : allocate a new tabulated gate 1695: 1696: phtg = calloc(1, sizeof(*phtg)); 1697: 1698: #endif 1699: 1700: //- increment registry count 1701: 1702: pheccer->tgt.iTabulatedGateCount++; 1703: 1704: //- initialize tabulation parameters 1705: 1706: phtg->iSizeParameters = iSize; 1707: phtg->pvParameters = pvParameters; 1708: 1709: //- initialize interval discretizer 1710: 1711: phtg->hi.dStart = dStart; 1712: phtg->hi.dEnd = dEnd; 1713: phtg->hi.dStep = (phtg->hi.dEnd - phtg->hi.dStart) / iEntries; 1714: 1715: //- initialize discrete function entries 1716: 1717: phtg->iEntries = iEntries; 1718: phtg->pdA = calloc(phtg->iEntries + 1, sizeof(*phtg->pdA)); 1719: phtg->pdB = calloc(phtg->iEntries + 1, sizeof(*phtg->pdB)); 1720: 1721: //- initialize the tao 1722: 1723: phtg->htao.iShape = 0; 1724: 1725: //- return result 1726: 1727: return(pheccer->tgt.iTabulatedGateCount - 1); 1728: } 1729: 1730: 1731: /// ************************************************************************** 1732: /// 1733: /// SHORT: HeccerTabulatedSpringMassCompareParameters() 1734: /// 1735: /// ARGS.: 1736: /// 1737: /// phtsm.: an initialized gate table. 1738: /// pv....: parameters to use for comparison. 1739: /// iSize.: size of parameter block. 1740: /// 1741: /// RTN..: int 1742: /// 1743: /// See memcmp() manual. 1744: /// 1745: /// DESCR: 1746: /// 1747: /// Compare parameters for a candidate gate with the parameters 1748: /// of an existing table. 1749: /// 1750: /// ************************************************************************** 1751: 1752: static 1753: int 1754: HeccerTabulatedSpringMassCompareParameters 1755: (struct HeccerTabulatedSpringMass *phtsm, void *pv, size_t iSize) 1756: { 1757: //- set default result : match 1758: 1759: int iResult = 0; 1760: 1761: //- set result : compare memory regions, using smallest size 1762: 1763: iSize = iSize < phtsm->iSizeParameters ? iSize : phtsm->iSizeParameters; 1764: 1765: iResult = memcmp(phtsm->pvParameters, pv, iSize); 1766: 1767: //- return result 1768: 1769: return(iResult); 1770: } 1771: 1772: 1773: /// ************************************************************************** 1774: /// 1775: /// SHORT: HeccerTabulatedSpringMassLookup() 1776: /// 1777: /// ARGS.: 1778: /// 1779: /// pheccer...: a heccer. 1780: /// pv........: table parameter block. 1781: /// iSize.....: size of parameter block. 1782: /// 1783: /// RTN..: int 1784: /// 1785: /// tabulated spring mass index, -1 for failure. 1786: /// 1787: /// DESCR: Lookup an existing table. 1788: /// 1789: /// ************************************************************************** 1790: 1791: static 1792: int 1793: HeccerTabulatedSpringMassLookup 1794: (struct Heccer *pheccer, void *pv, size_t iSize) 1795: { 1796: //- set default result : not found 1797: 1798: int iResult = -1; 1799: 1800: //! a protection for the case where you accidently forget to dereference 1801: 1802: if (iSize < 10) 1803: { 1804: ((int *)0)[0] = 0; 1805: } 1806: 1807: //- loop over all tables 1808: 1809: int i; 1810: 1811: for (i = 0 ; i < pheccer->tsmt.iTabulatedSpringMassCount; i++) 1812: { 1813: //- get pointer to current table 1814: 1815: struct HeccerTabulatedSpringMass *phtsm = &pheccer->tsmt.phtsm[i]; 1816: 1817: //- if match 1818: 1819: if (HeccerTabulatedSpringMassCompareParameters(phtsm, pv, iSize) == 0) 1820: { 1821: //- set table index 1822: 1823: iResult = i; 1824: 1825: //- break searching loop 1826: 1827: break; 1828: } 1829: } 1830: 1831: //- return result 1832: 1833: return(iResult); 1834: } 1835: 1836: 1837: /// ************************************************************************** 1838: /// 1839: /// SHORT: HeccerTabulatedSpringMassNew() 1840: /// 1841: /// ARGS.: 1842: /// 1843: /// pheccer...: a heccer. 1844: /// 1845: /// RTN..: int 1846: /// 1847: /// tabulated spring mass index, -1 for failure. 1848: /// 1849: /// DESCR: Allocate a new table. 1850: /// 1851: /// ************************************************************************** 1852: 1853: static 1854: int 1855: HeccerTabulatedSpringMassNew 1856: (struct Heccer *pheccer, void *pvParameters, size_t iSize) 1857: { 1858: struct HeccerTabulatedSpringMass *phtsm = NULL; 1859: 1860: if (pheccer->tsmt.iTabulatedSpringMassCount >= HECCER_TABULATED_SPRINGMASS_MAX) 1861: { 1862: return(-1); 1863: } 1864: 1865: //- set result : from pool 1866: 1867: phtsm = &pheccer->tsmt.phtsm[pheccer->tsmt.iTabulatedSpringMassCount]; 1868: 1869: //- increment registry count 1870: 1871: pheccer->tsmt.iTabulatedSpringMassCount++; 1872: 1873: //- initialize tabulation parameters 1874: 1875: phtsm->iSizeParameters = iSize; 1876: phtsm->pvParameters = pvParameters; 1877: 1878: //- return result 1879: 1880: return(pheccer->tsmt.iTabulatedSpringMassCount - 1); 1881: } 1882: 1883: 1884: /// ************************************************************************** 1885: /// 1886: /// SHORT: HeccerTabulatedGateStore() 1887: /// 1888: /// ARGS.: 1889: /// 1890: /// pheccer...: a heccer. 1891: /// phtgNew...: tabulated gate. 1892: /// 1893: /// RTN..: int 1894: /// 1895: /// tabulated gate index, -1 for failure. 1896: /// 1897: /// DESCR: Store a new table. 1898: /// 1899: /// Note that the parameters that identify the table must be 1900: /// initialized correctly before this function is called. 1901: /// 1902: /// ************************************************************************** 1903: 1904: static 1905: int 1906: HeccerTabulatedGateStore 1907: (struct Heccer *pheccer, 1908: struct HeccerTabulatedGate *phtgNew) 1909: { 1910: struct HeccerTabulatedGate *phtg = NULL; 1911: 1912: if (pheccer->tgt.iTabulatedGateCount >= HECCER_TABULATED_GATES_MAX) 1913: { 1914: return(-1); 1915: } 1916: 1917: #define HECCER_STATIC_TABULATED_GATES 1918: #ifdef HECCER_STATIC_TABULATED_GATES 1919: 1920: //- set result : from pool 1921: 1922: phtg = &pheccer->tgt.phtg[pheccer->tgt.iTabulatedGateCount]; 1923: 1924: #else 1925: 1926: //- set result : allocate a new tabulated gate 1927: 1928: phtg = calloc(1, sizeof(*phtg)); 1929: 1930: #endif 1931: 1932: //- increment registry count 1933: 1934: pheccer->tgt.iTabulatedGateCount++; 1935: 1936: //- we simply copy all the parameters 1937: 1938: *phtg = *phtgNew; 1939: 1940: //- return result 1941: 1942: return(pheccer->tgt.iTabulatedGateCount - 1); 1943: } 1944: 1945: 1946: