file:/local_home/local_home/hugo/neurospaces_project/heccer/source/c/snapshots/0/minimumdegree.c        (Mon Jun 16 00:04:16 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 <stdlib.h>
  20: 
  21: #include "heccer/heccer.h"
  22: #include "heccer/minimumdegree.h"
  23: 
  24: 
  25: //f static prototypes
  26: 
  27: static
  28: int
  29: HeccerMDFlowEnumerator
  30: (struct Heccer *pheccer, int iCurrent, int iTarget);
  31: 
  32: static
  33: int
  34: HeccerMDFlowEnumeratorB
  35: (struct Heccer *pheccer, int iCurrent, int iTarget);
  36: 
  37: static
  38: int
  39: HeccerMDFlowEnumeratorL
  40: (struct Heccer *pheccer, int iCurrent, int iTarget);
  41: 
  42: static int HeccerMDFindFlow(struct Heccer *pheccer, int iCompartments);
  43: 
  44: static int HeccerMDStructuralyze(struct Heccer *pheccer, int iCompartments);
  45: 
  46: 
  47: /// **************************************************************************
  48: ///
  49: /// SHORT: HeccerMDFindFlow()
  50: ///
  51: /// ARGS.:
  52: ///
  53: ///     pheccer........: a heccer.
  54: ///     iCompartments..: number of compartments.
  55: ///
  56: /// RTN..: int
  57: ///
  58: ///     success of operation.
  59: ///
  60: /// DESCR: Find natural flow in the morphology.
  61: ///
  62: /// NOTE.:
  63: ///
  64: ///     This is assumed to follow reverse current flow for a passive
  65: ///     morphology, hence the name of the function.
  66: ///
  67: /// **************************************************************************
  68: 
  69: /// driver : select and drive the enumerator
  70: ///
  71: ///     A small variant of hines, 1986.  Recursive function, typically
  72: ///     driven with something like HeccerMDFlowEnumerator(p, a, b),
  73: ///     where the a is the soma index and b is the number of
  74: ///     compartments.
  75: ///
  76: /// NOTE.:
  77: ///
  78: ///     This is a very specific implementation of more general minimum
  79: ///     degree algorithms, and is therefore suboptimal, afatiac.
  80: ///
  81: 
  82: static
  83: int
  84: HeccerMDFlowEnumerator
  85: (struct Heccer *pheccer, int iCurrent, int iTarget)
  86: {
  87:     //- register schedule number for current compartment
  88: 
  89:     pheccer->indexers.md.piForward[iCurrent] = iTarget;
  90: 
  91:     //- register reverse schedule number
  92: 
  93:     pheccer->indexers.md.piBackward[iTarget] = iCurrent;
  94: 
  95:     //- decrement schedule number
  96: 
  97:     iTarget -= 1;
  98: 
  99:     //- if branches first scheduling
 100: 
 101:     //t note : leaves first does not work yet, do not use.
 102: 
 103:     if (1 || pheccer->ho.iOptions & HECCER_OPTION_BRANCHES_FIRST_SCHEDULING)
 104:     {
 105:         //- do a simple depth-first ordering
 106: 
 107:         return(HeccerMDFlowEnumeratorB(pheccer, iCurrent, iTarget));
 108:     }
 109: 
 110:     //- else leaves first scheduling
 111: 
 112:     else
 113:     {
 114:         //- loop over children for current compartment
 115: 
 116:         int i;
 117: 
 118:         for (i = 0 ; i < pheccer->indexers.md.piChildren[iCurrent]; i++)
 119:         {
 120:             //- get current child index
 121: 
 122:             int iIndex = pheccer->indexers.md.ppiChildren[iCurrent][i];
 123: 
 124:             //- only if this child has no children
 125: 
 126:             if (pheccer->indexers.md.piChildren[iIndex] == 0)
 127:             {
 128:                 //- register schedule number for current compartment
 129: 
 130:                 pheccer->indexers.md.piForward[iIndex] = iTarget;
 131: 
 132:                 //- register reverse schedule number
 133: 
 134:                 pheccer->indexers.md.piBackward[iTarget] = iIndex;
 135: 
 136:                 //- decrement schedule number
 137: 
 138:                 iTarget -= 1;
 139:             }
 140:         }
 141: 
 142:         //- loop over children for current compartment
 143: 
 144:         for (i = 0 ; i < pheccer->indexers.md.piChildren[iCurrent]; i++)
 145:         {
 146:             //- get current child index
 147: 
 148:             int iIndex = pheccer->indexers.md.ppiChildren[iCurrent][i];
 149: 
 150:             //- find reverse flow starting at this child
 151: 
 152:             //! depth-first is the essence of branches first.
 153: 
 154:             //! depth-first has two flavours : pre-order, post-order,
 155:             //! they are equivalent for this purpose, ie. for a heccer.
 156: 
 157:             iTarget = HeccerMDFlowEnumeratorL(pheccer, iIndex, iTarget);
 158:         }
 159:     }
 160: 
 161:     //- no compartments : should return zero (recursively)
 162: 
 163:     //! check should be made by the driver of this routine
 164: 
 165:     return(iTarget);
 166: }
 167: 
 168: 
 169: /// worker : recursively schedule all compartments, branches first
 170: 
 171: static
 172: int
 173: HeccerMDFlowEnumeratorB
 174: (struct Heccer *pheccer, int iCurrent, int iTarget)
 175: {
 176:     //- loop over children for current compartment
 177: 
 178:     int i;
 179: 
 180:     for (i = 0 ; i < pheccer->indexers.md.piChildren[iCurrent]; i++)
 181:     {
 182:         //- get current child index
 183: 
 184:         int iIndex = pheccer->indexers.md.ppiChildren[iCurrent][i];
 185: 
 186:         //- find reverse flow starting at this child
 187: 
 188:         //! depth-first is the essence of branches first.
 189: 
 190:         //! depth-first has two flavours : pre-order, post-order,
 191:         //! they are equivalent for this purpose, ie. for a heccer.
 192: 
 193:         iTarget = HeccerMDFlowEnumerator(pheccer, iIndex, iTarget);
 194:     }
 195: 
 196:     //- no compartments : should return zero (recursively)
 197: 
 198:     //! check should be made by the driver of this routine
 199: 
 200:     return(iTarget);
 201: }
 202: 
 203: 
 204: /// worker : recursively schedule all compartments, leaves first
 205: 
 206: static
 207: int
 208: HeccerMDFlowEnumeratorL
 209: (struct Heccer *pheccer, int iCurrent, int iTarget)
 210: {
 211:     //- loop over children for current compartment
 212: 
 213:     int i;
 214: 
 215:     for (i = 0 ; i < pheccer->indexers.md.piChildren[iCurrent]; i++)
 216:     {
 217:         //- get current child index
 218: 
 219:         int iIndex = pheccer->indexers.md.ppiChildren[iCurrent][i];
 220: 
 221:         //- only if this child has no children
 222: 
 223:         if (pheccer->indexers.md.piChildren[iIndex] == 0)
 224:         {
 225:             //- register schedule number for current compartment
 226: 
 227:             pheccer->indexers.md.piForward[iIndex] = iTarget;
 228: 
 229:             //- register reverse schedule number
 230: 
 231:             pheccer->indexers.md.piBackward[iTarget] = iIndex;
 232: 
 233:             //- decrement schedule number
 234: 
 235:             iTarget -= 1;
 236:         }
 237:     }
 238: 
 239:     //- loop over children for current compartment
 240: 
 241:     for (i = 0 ; i < pheccer->indexers.md.piChildren[iCurrent]; i++)
 242:     {
 243:         //- get current child index
 244: 
 245:         int iIndex = pheccer->indexers.md.ppiChildren[iCurrent][i];
 246: 
 247:         //- find reverse flow starting at this child
 248: 
 249:         iTarget = HeccerMDFlowEnumerator(pheccer, iIndex, iTarget);
 250:     }
 251: 
 252:     //- no compartments : should return zero (recursively)
 253: 
 254:     //! check should be made by the driver of this routine
 255: 
 256:     return(iTarget);
 257: }
 258: 
 259: 
 260: /// driver : initialize and drive the enumerator
 261: 
 262: static int HeccerMDFindFlow(struct Heccer *pheccer, int iCompartments)
 263: {
 264:     //- set default result : ok
 265: 
 266:     int iResult = TRUE;
 267: 
 268:     struct MinimumDegree *pmd = &pheccer->indexers.md;
 269: 
 270:     //- allocate unordered to flow
 271: 
 272:     int *piForward = (int *)calloc(iCompartments, sizeof(int));
 273: 
 274:     //- allocate flow to unordered
 275: 
 276:     int *piBackward = (int *)calloc(iCompartments, sizeof(int));
 277: 
 278:     //- set pointers
 279: 
 280:     pmd->piForward = piForward;
 281: 
 282:     pmd->piBackward = piBackward;
 283: 
 284:     //- loop over compartments
 285: 
 286:     int i;
 287: 
 288:     int iStart;
 289: 
 290:     for (i = 0; i < iCompartments; i++)
 291:     {
 292:         //- if compartment has no parent
 293: 
 294:         if (pheccer->inter.pcomp[i].iParent == -1)
 295:         {
 296:             //- this is a candidate to start the reverse flow
 297: 
 298:             //! usually it will be the soma, but ok if it is a dendritic tip
 299: 
 300:             iStart = i;
 301:         }
 302:     }
 303: 
 304:     //- schedule this compartment to be the last one solved
 305: 
 306:     int iTarget = iCompartments - 1;
 307: 
 308:     //- find reverse flow using the candidate as starting point
 309: 
 310:     int iEnd = HeccerMDFlowEnumerator(pheccer, iStart, iTarget);
 311: 
 312:     //- check if we ended with something reasonable (can only be -1, afaict)
 313: 
 314:     if (iEnd != -1)
 315:     {
 316:         //t this can happen if the intermediary structure is wrong,
 317:         //t e.g. out of bound parent index or cycles
 318: 
 319:         //t also happens for multiple segments without a parent
 320: 
 321:         //t add something like HeccerError(number, message, varargs);
 322: 
 323:         //! segv
 324: 
 325:         *(int *)0 = 0;
 326:     }
 327: 
 328:     //- set number of entries
 329: 
 330:     pmd->iEntries = iCompartments;
 331: 
 332:     //- return result
 333: 
 334:     return(piForward && piBackward);
 335: }
 336: 
 337: 
 338: /// **************************************************************************
 339: ///
 340: /// SHORT: HeccerMDStructuralyze()
 341: ///
 342: /// ARGS.:
 343: ///
 344: ///     pheccer........: a heccer.
 345: ///     iCompartments..: number of compartments.
 346: ///
 347: /// RTN..: int
 348: ///
 349: ///     success of operation.
 350: ///
 351: /// DESCR: Build structural indices.
 352: ///
 353: /// **************************************************************************
 354: 
 355: static int HeccerMDStructuralyze(struct Heccer *pheccer, int iCompartments)
 356: {
 357:     //- set default result : ok
 358: 
 359:     int iResult = TRUE;
 360: 
 361:     struct MinimumDegree *pmd = &pheccer->indexers.md;
 362: 
 363:     //- allocate structural analyzers
 364: 
 365: /*     int *piParents = (int *)calloc(iCompartments, sizeof(int)); */
 366: 
 367:     int *piChildren = (int *)calloc(iCompartments, sizeof(int));
 368: 
 369:     int **ppiChildren = (int **)calloc(iCompartments, sizeof(int *));
 370: 
 371:     //- analyze
 372: 
 373:     //- loop over all compartments
 374: 
 375:     int i;
 376: 
 377:     for (i = 0; i < iCompartments; i++)
 378:     {
 379:         //- set current compartment
 380: 
 381:         struct Compartment *pcomp = &pheccer->inter.pcomp[i];
 382: 
 383: /*      //- set default : current compartment has no parent */
 384: 
 385: /*      piParents[i] = -1; */
 386: 
 387:         //- if compartment has parent
 388: 
 389:         if (pcomp->iParent != -1)
 390:         {
 391:             if (pcomp->iParent == i
 392:                 || pcomp->iParent >= iCompartments)
 393:             {
 394:                 //t HeccerError(number, message, varargs);
 395: 
 396:                 HeccerError
 397:                     (pheccer,
 398:                      NULL,
 399:                      "the compartment array does not describe a valid tree structure"
 400:                      " at compartment (%i)\n",
 401:                      i);
 402:             }
 403: 
 404: /*          //- register parent compartment index */
 405: 
 406: /*          piParents[i] = pcomp->iParent; */
 407: 
 408:             //- increment number of children for the parent compartment
 409: 
 410:             piChildren[pcomp->iParent] += 1;
 411:         }
 412:     }
 413: 
 414:     //- build indices for children
 415: 
 416:     //- i loop over all compartments
 417: 
 418:     for (i = 0; i < iCompartments; i++)
 419:     {
 420:         //- clear number of children for this compartment
 421: 
 422:         int iChildren = 0;
 423: 
 424:         //- if compartment has no children
 425: 
 426:         if (piChildren[i] == 0)
 427:         {
 428:             //- it's a terminal branch
 429: 
 430:             ppiChildren[i] = NULL;
 431:         }
 432: 
 433:         //- else
 434: 
 435:         else
 436:         {
 437:             //- allocate for number of children
 438: 
 439:             ppiChildren[i] = (int *)calloc(piChildren[i], sizeof(int));
 440: 
 441:             //- j loop over all compartments
 442: 
 443:             int j;
 444: 
 445:             for (j = 0; j < iCompartments; j++)
 446:             {
 447:                 //- if parent of j compartment is i compartment
 448: 
 449:                 if (pheccer->inter.pcomp[j].iParent == i)
 450:                 {
 451:                     //- remember j as a children of i
 452: 
 453:                     ppiChildren[i][iChildren] = j;
 454: 
 455:                     //- increment number of children
 456: 
 457:                     iChildren++;
 458:                 }
 459:             }
 460:         }
 461:     }
 462:         
 463:     //- set pointers
 464: 
 465: /*     pmd->piParents = piParents; */
 466: 
 467:     pmd->piChildren = piChildren;
 468: 
 469:     pmd->ppiChildren = ppiChildren;
 470: 
 471:     //- return result
 472: 
 473:     return(/* piParents &&  */piChildren && ppiChildren);
 474: }
 475: 
 476: 
 477: /// **************************************************************************
 478: ///
 479: /// SHORT: HeccerMinimumDegree()
 480: ///
 481: /// ARGS.:
 482: ///
 483: ///     pheccer...: a heccer.
 484: ///
 485: /// RTN..: int
 486: ///
 487: ///     success of operation.
 488: ///
 489: /// DESCR: Minimum degree enumeration for compartment matrix.
 490: ///
 491: /// **************************************************************************
 492: 
 493: int HeccerMinimumDegree(struct Heccer *pheccer)
 494: {
 495:     //- set default result : ok
 496: 
 497:     int iResult = TRUE;
 498: 
 499:     if (pheccer->inter.iCompartments == 0)
 500:     {
 501:         return(TRUE);
 502:     }
 503: 
 504:     //- build structural indices
 505: 
 506:     iResult = iResult && HeccerMDStructuralyze(pheccer, pheccer->inter.iCompartments);
 507: 
 508:     //- do minimum degree
 509: 
 510:     //! this is assumed to follow reverse current flow for a passive
 511:     //! morphology, hence the name of the function.
 512: 
 513:     iResult = iResult && HeccerMDFindFlow(pheccer, pheccer->inter.iCompartments);
 514: 
 515: /*     //- sort children according to flow */
 516: 
 517: /*     iResult = iResult && HeccerMDSortChildren(pheccer, pheccer->inter.iCompartments); */
 518: 
 519:     //- return result
 520: 
 521:     return(iResult);
 522: }
 523: 
 524: 
 525: /// **************************************************************************
 526: ///
 527: /// SHORT: HeccerMinimumDegreeDump()
 528: ///
 529: /// ARGS.:
 530: ///
 531: ///     pmd........: minimum degree index.
 532: ///     pfile......: stdio file.
 533: ///     iSelection.: selection to dump.
 534: ///
 535: /// RTN..: int
 536: ///
 537: ///     success of operation.
 538: ///
 539: /// DESCR: Perform the compartment operators once.
 540: ///
 541: /// **************************************************************************
 542: 
 543: int
 544: HeccerMinimumDegreeDump
 545: (struct MinimumDegree *pmd, FILE *pfile, int iSelection)
 546: {
 547:     //- set default result
 548: 
 549:     int iResult = TRUE;
 550: 
 551:     //- index of parent compartment, -1 for none
 552: 
 553:     if (iSelection & HECCER_DUMP_INDEXERS_SUMMARY)
 554:     {
 555:         fprintf(pfile, "MinimumDegree (iEntries) : (%i)\n", pmd->iEntries);
 556:     }
 557: 
 558:     //- structural analyzers
 559: 
 560:     if (iSelection & HECCER_DUMP_INDEXERS_STRUCTURE)
 561:     {
 562:         int i;
 563: 
 564:         for (i = 0 ; i < pmd->iEntries ; i++)
 565:         {
 566:             fprintf(pfile, "MinimumDegree (piChildren[%i]) : (%i)\n", i, pmd->piChildren[i]);
 567: 
 568:             int j;
 569: 
 570:             for (j = 0 ; j < pmd->piChildren[i] ; j++)
 571:             {
 572:                 fprintf(pfile, "MinimumDegree (piChildren[%i][%i]) : (%i)\n", i, j, pmd->ppiChildren[i][j]);
 573:             }
 574:         }
 575:     }
 576: 
 577:     //- unordered to flow
 578: 
 579:     if (iSelection & HECCER_DUMP_INDEXERS_STRUCTURE)
 580:     {
 581:         int i;
 582: 
 583:         for (i = 0 ; i < pmd->iEntries ; i++)
 584:         {
 585:             fprintf(pfile, "MinimumDegree (piForward[%i]) : (%i)\n", i, pmd->piForward[i]);
 586:         }
 587:     }
 588: 
 589:     //- flow to unordered
 590: 
 591:     if (iSelection & HECCER_DUMP_INDEXERS_STRUCTURE)
 592:     {
 593:         int i;
 594: 
 595:         for (i = 0 ; i < pmd->iEntries ; i++)
 596:         {
 597:             fprintf(pfile, "MinimumDegree (piBackward[%i]) : (%i)\n", i, pmd->piBackward[i]);
 598:         }
 599:     }
 600: 
 601:     //- return result
 602: 
 603:     return(iResult);
 604: }
 605: 
 606: 
 607: 








































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