1: //
   2: // Heccer : a compartmental solver that implements efficient Crank-Nicolson
   3: // integration for neuronal models.
   4: //
   5: 
   6: //////////////////////////////////////////////////////////////////////////////
   7: //'
   8: //' Heccer : testbed C implementation
   9: //'
  10: //' Copyright (C) 2006-2008 Hugo Cornelis
  11: //'
  12: //' functional ideas .. Hugo Cornelis, hugo.cornelis@gmail.com
  13: //'
  14: //' coding ............ Hugo Cornelis, hugo.cornelis@gmail.com
  15: //'
  16: //////////////////////////////////////////////////////////////////////////////
  17: 
  18: 
  19: #include <stdio.h>
  20: #include <stdlib.h>
  21: 
  22: #include <neurospaces/pidinstack.h>
  23: #include <neurospaces/treespacetraversal.h>
  24: 
  25: #include "heccer/addressing.h"
  26: #include "heccer/compartment.h"
  27: #include "heccer/intermediary.h"
  28: #include "heccer/neurospaces/heccer.h"
  29: #include "heccer/neurospaces/segments.h"
  30: #include "heccer/service.h"
  31: 
  32: 
  33: static
  34: int
  35: solver_segmentprocessor(struct TreespaceTraversal *ptstr, void *pvUserdata);
  36: 
  37: static int cellsolver_getsegments(struct Heccer *pheccer, struct TranslationService *pts);
  38: 
  39: static int cellsolver_linksegments(struct Heccer *pheccer);
  40: 
  41: static int cellsolver_setupcomp2mech(struct Heccer *pheccer);
  42: 
  43: 
  44: static
  45: int
  46: solver_segmentprocessor(struct TreespaceTraversal *ptstr, void *pvUserdata)
  47: {
  48:     //- set default result : ok
  49: 
  50:     int iResult = TSTR_PROCESSOR_SUCCESS;
  51: 
  52:     //- get actual symbol
  53: 
  54:     struct symtab_HSolveListElement *phsle = (struct symtab_HSolveListElement *)TstrGetActual(ptstr);
  55: 
  56:     //- if segment
  57: 
  58:     if (instanceof_segment(phsle))
  59:     {
  60:         //- get pointer to intermediary
  61: 
  62:         struct Heccer *pheccer = (struct Heccer *)pvUserdata;
  63: 
  64:         struct Intermediary *pinter = &pheccer->inter;
  65: 
  66:         //- register current symbol
  67: 
  68:         int iSegment = pinter->iCompartments;
  69: 
  70:         //- register type
  71: 
  72:         pinter->pcomp[iSegment].mc.iType = MATH_TYPE_Compartment;
  73: 
  74:         //- register identification
  75: 
  76:         //! note: assumes pp define HECCER_SOURCE_NEUROSPACES
  77: 
  78:         int iSerial = PidinStackToSerial(ptstr->ppist);
  79: 
  80:         iSerial = ADDRESSING_NEUROSPACES_2_HECCER(iSerial);
  81: 
  82:         pinter->pcomp[iSegment].mc.iSerial = iSerial;
  83: 
  84: #ifdef HECCER_SOURCE_TYPING
  85: 
  86:         int iModelSourceType
  87:             = SymbolParameterResolveValue(phsle, ptstr->ppist, "SOURCE_TYPE");
  88: 
  89:         pinter->pcomp[iSegment].mc.iModelSourceType = iModelSourceType;
  90: 
  91: #endif
  92: 
  93:         //- register parameters
  94: 
  95:         //t check for error returns, abort traversal if necessary
  96: 
  97:         double dCm
  98:             = SymbolParameterResolveScaledValue(phsle, ptstr->ppist, "CM");
  99: 
 100:         if (dCm == FLT_MAX)
 101:         {
 102:             iResult = TSTR_PROCESSOR_ABORT;
 103:         }
 104: 
 105:         double dEm
 106:             = SymbolParameterResolveValue(phsle, ptstr->ppist, "ELEAK");
 107: 
 108:         if (dEm == FLT_MAX)
 109:         {
 110:             iResult = TSTR_PROCESSOR_ABORT;
 111:         }
 112: 
 113:         double dInitVm
 114:             = SymbolParameterResolveValue(phsle, ptstr->ppist, "Vm_init");
 115: 
 116:         if (dInitVm == FLT_MAX)
 117:         {
 118:             iResult = TSTR_PROCESSOR_ABORT;
 119:         }
 120: 
 121:         double dInject
 122:             = SymbolParameterResolveValue(phsle, ptstr->ppist, "INJECT");
 123: 
 124:         if (dInject == FLT_MAX)
 125:         {
 126:             dInject = 0.0;
 127:         }
 128: 
 129:         double dRa
 130:             = SymbolParameterResolveScaledValue(phsle, ptstr->ppist, "RA");
 131: 
 132:         if (dRa == FLT_MAX)
 133:         {
 134:             iResult = TSTR_PROCESSOR_ABORT;
 135:         }
 136: 
 137:         double dRm
 138:             = SymbolParameterResolveScaledValue(phsle, ptstr->ppist, "RM");
 139: 
 140:         if (dRm == FLT_MAX)
 141:         {
 142:             iResult = TSTR_PROCESSOR_ABORT;
 143:         }
 144: 
 145:         pinter->pcomp[iSegment].dCm = dCm;
 146:         pinter->pcomp[iSegment].dEm = dEm;
 147:         pinter->pcomp[iSegment].dInitVm = dInitVm;
 148:         pinter->pcomp[iSegment].dInject = dInject;
 149:         pinter->pcomp[iSegment].dRa = dRa;
 150:         pinter->pcomp[iSegment].dRm = dRm;
 151: 
 152:         //- register serial of parent
 153: 
 154:         {
 155:             struct symtab_Parameters *pparParent
 156:                 = SymbolFindParameter(phsle, ptstr->ppist, "PARENT");
 157: 
 158:             if (pparParent)
 159:             {
 160:                 //t I can just subtract the cell's segment ID ?
 161: 
 162:                 struct PidinStack *ppistParent
 163:                     = ParameterResolveToPidinStack(pparParent, ptstr->ppist);
 164: 
 165:                 if (ppistParent)
 166:                 {
 167:                     PidinStackLookupTopSymbol(ppistParent);
 168: 
 169:                     int iParent = PidinStackToSerial(ppistParent);
 170: 
 171:                     //t check for error
 172: 
 173:                     if (iParent != INT_MAX)
 174:                     {
 175:                         iParent = ADDRESSING_NEUROSPACES_2_HECCER(iParent);
 176: 
 177:                         pinter->pcomp[iSegment].iParent = iParent;
 178:                     }
 179:                     else
 180:                     {
 181:                         //! parent does not exist
 182: 
 183:                         iResult = TSTR_PROCESSOR_ABORT;
 184:                     }
 185:                 }
 186:                 else
 187:                 {
 188:                     iResult = TSTR_PROCESSOR_ABORT;
 189:                 }
 190:             }
 191:             else
 192:             {
 193:                 pinter->pcomp[iSegment].iParent = -1;
 194:             }
 195:         }
 196: 
 197:         //- increment number of solved segments
 198: 
 199:         pinter->iCompartments++;
 200: 
 201:         if (iResult == TSTR_PROCESSOR_ABORT)
 202:         {
 203:             HeccerError
 204:                 (pheccer,
 205:                  NULL,
 206:                  "compartment array translation failed at compartment (compartment %i, serial %i)\n",
 207:                  pinter->iCompartments,
 208:                  pinter->pcomp[iSegment].mc.iSerial);
 209:         }
 210: 
 211:     }
 212: 
 213:     //- return result
 214: 
 215:     return(iResult);
 216: }
 217: 
 218: 
 219: static int cellsolver_getsegments(struct Heccer *pheccer, struct TranslationService *pts)
 220: {
 221:     //- set default result : ok
 222: 
 223:     int iResult = TRUE;
 224: 
 225:     //- allocate pidin stack pointing to root
 226: 
 227:     struct PidinStack *ppistRoot = pts->ptsd->ppistRoot;
 228: 
 229:     struct symtab_HSolveListElement *phsleRoot = pts->ptsd->phsleRoot;
 230: 
 231:     //- get model context to solve
 232: 
 233:     struct PidinStack *ppistModel
 234:         = SymbolPrincipalSerial2Context(phsleRoot, ppistRoot, pts->ptsd->iModel);
 235: 
 236:     if (ppistModel)
 237:     {
 238:         //- lookup symbol
 239: 
 240:         struct symtab_HSolveListElement *phsleModel
 241:             = PidinStackLookupTopSymbol(ppistModel);
 242: 
 243:         int iSegments = SymbolCountSegments(phsleModel, ppistModel);
 244: 
 245:         if (iSegments == 0)
 246:         {
 247:             char pc[1000];
 248: 
 249:             PidinStackString(ppistModel, pc, sizeof(pc));
 250: 
 251:             HeccerError
 252:                 (pheccer,
 253:                  NULL,
 254:                  "no compartments found inside %s\n",
 255:                  pc);
 256: 
 257:             return(FALSE);
 258:         }
 259: 
 260:         int i;
 261: 
 262:         //- allocate intermediary for segments to solve
 263: 
 264:         struct Intermediary *pinter = &pheccer->inter;
 265: 
 266:         pinter->iCompartments = 0; //iSegments;
 267: 
 268:         pinter->pcomp
 269:             = (struct Compartment *)calloc(iSegments, sizeof(struct Compartment));
 270: 
 271:         //- register solved segments in cell
 272: 
 273:         if (SymbolTraverseSegments
 274:             (phsleModel, ppistModel, solver_segmentprocessor, NULL, pheccer) == -1)
 275:         {
 276:             iResult = FALSE;
 277:         }
 278: 
 279:         //t should use SolverInfoPrincipalSerial2SegmentSerial() for this
 280: 
 281:         //- link the segments together using the parent link
 282: 
 283:         iResult = iResult && cellsolver_linksegments(pheccer);
 284: 
 285:         //- produce a piC2m array: all zeros for now
 286: 
 287:         iResult = iResult && cellsolver_setupcomp2mech(pheccer);
 288:     }
 289:     else
 290:     {
 291:         iResult = FALSE;
 292:     }
 293: 
 294:     //- return result
 295: 
 296:     return(iResult);
 297: }
 298: 
 299: 
 300: static int cellsolver_linksegments(struct Heccer *pheccer)
 301: {
 302:     //- set default result : ok
 303: 
 304:     int iResult = TRUE;
 305: 
 306:     //! this can never fail: multiple parent and multiple root
 307:     //! diagnosis is left to the minimum degree logic.
 308: 
 309:     //- loop over all compartments
 310: 
 311:     struct Intermediary *pinter = &pheccer->inter;
 312: 
 313:     int iCompartment;
 314: 
 315:     for (iCompartment = 0 ; iCompartment < pinter->iCompartments ; iCompartment++)
 316:     {
 317:         //- get current compartment
 318: 
 319:         struct Compartment *pcomp = &pinter->pcomp[iCompartment];
 320: 
 321:         //- get parent serial
 322: 
 323:         int iParent = pcomp->iParent;
 324: 
 325:         if (iParent != -1)
 326:         {
 327:             //- search all compartments
 328: 
 329:             int i;
 330: 
 331:             for (i = 0 ; i < pinter->iCompartments ; i++)
 332:             {
 333:                 //- get current compartment
 334: 
 335:                 struct Compartment *pcompParent = &pinter->pcomp[i];
 336: 
 337:                 //- if is parent
 338: 
 339:                 if (pcompParent->mc.iSerial == iParent)
 340:                 {
 341:                     //- ok, stop searching loop
 342: 
 343:                     break;
 344:                 }
 345:             }
 346: 
 347:             //- set current compartment parent index
 348: 
 349:             if (i < pinter->iCompartments)
 350:             {
 351:                 pcomp->iParent = i;
 352:             }
 353:             else
 354:             {
 355:                 HeccerError
 356:                     (pheccer,
 357:                      NULL,
 358:                      "compartment (compartment %i, serial %i (internal %i)) parent linking failed for parent serial %i (internal %i)\n",
 359:                      iCompartment,
 360:                      ADDRESSING_HECCER_2_NEUROSPACES(pcomp->mc.iSerial),
 361:                      pcomp->mc.iSerial,
 362:                      ADDRESSING_HECCER_2_NEUROSPACES(iParent),
 363:                      iParent);
 364: 
 365:                 iResult = FALSE;
 366:             }
 367:         }
 368:     }
 369: 
 370:     //- return result
 371: 
 372:     return(iResult);
 373: }
 374: 
 375: 
 376: static int cellsolver_setupcomp2mech(struct Heccer *pheccer)
 377: {
 378:     //- allocate
 379: 
 380:     struct Intermediary *pinter = &pheccer->inter;
 381: 
 382:     pinter->piC2m = (int *)calloc(pinter->iCompartments + 1, sizeof(int));
 383: 
 384:     //- sign end of array
 385: 
 386:     pinter->piC2m[pinter->iCompartments] = -1;
 387: 
 388:     //- return result
 389: 
 390:     return(pinter->piC2m != NULL);
 391: }
 392: 
 393: 
 394: int HeccerNeurospacesSegments2Compartments(struct TranslationService *pts)
 395: {
 396:     //- set the service function pointer to cellsolver_getsegments()
 397: 
 398:     pts->segments_inspector = cellsolver_getsegments;
 399: 
 400:     //- return result
 401: 
 402:     return(1);
 403: }
 404: 
 405: 
 406: