file:/local_home/local_home/hugo/neurospaces_project/heccer/source/c/snapshots/0/table.c        (Sun Jul 13 22:29:22 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 <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: 








































Generated by Xrefactory version 2.0.14 on Thu Jul 24 22:41:20 2008