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: