/*---------------------------------------------------------------------
Written by Alessandro Limongiello
University "Federico II" of Naples, 2012-2013
-----------------------------------------------------------------------*/
/*---------------------------------------------
USEFUL PROCEDURES
---------------------------------------------*/
func max(){
if ($1 > $2){
return $1
} else {
return $2
}
}
/*---------------------------------------------
Strahler Numbers Computation
---------------------------------------------*/
// Useful variables
objectvar secNameList0
objectvar secaddress0 // section address
objectvar paraddress0 // parentsection address
objectvar strahlerNumbersVect // Vector for Strahler Numbers
objectvar stimulatedGroupVect // Vector for Stimulated Group
objref sRef
objectvar childrenSecListVector // a List of Vector for children of each section
objectvar parSecVect // parent for a section
strdef str,strFile
func strahlerNumbersFunc(){ local sizeVect,ii,sN,egualBool localobj vectApp
sizeVect = childrenSecListVector.o($1).size()
if (sizeVect==0){
strahlerNumbersVect.set($1,1)
return 1}
vectApp = new Vector(sizeVect,0)
egualBool = 1
for ii = 0,sizeVect-1{
sN = strahlerNumbersFunc(childrenSecListVector.o($1).x(ii))
if (ii>0 && vectApp.contains(sN)==0) egualBool=0
vectApp.set(ii,sN)
}
if (egualBool==1){
strahlerNumbersVect.set($1,vectApp.x(0)+1)
return vectApp.x(0)+1
}
strahlerNumbersVect.set($1,vectApp.max())
return vectApp.max()
}
proc strahlerNumbersComputation(){local sec,sec1,sec2,ii
stoprun = 1
if (optPrint==1) print "Strahler Numbers Computation"
if (firstInit ==1){
forall delete_section()
readingFile()
xopen(strFile) // read geometry file
firstInit == 0
}else{xopen(strFile)}
{NSEC=0 forall if (!issection("eqSection.*")){NSEC=NSEC+1}}
secNameList0=new List()
strahlerNumbersVect=new Vector(NSEC,0) //set/reset
stimulatedGroupVect=new Vector(NSEC,0) //set/reset
secaddress0=new Vector(NSEC,0) //set/reset
paraddress0=new Vector(NSEC,0) //set/reset
childrenSecListVector = new List()
parSecVect = new Vector(NSEC,0) //set/reset
sec = 0
forall if (!issection("eqSection.*")){
secNameList0.append(new String(secname()))
secaddress0.set(sec, this_section())
paraddress0.set(sec, parent_section(0))
childrenSecListVector.append(new Vector())
sec = sec+1
}
for sec1=1,NSEC-1 {
for sec2=0,NSEC-1 {
if (secaddress0.x(sec2)==paraddress0.x(sec1)) {
parSecVect.set(sec1,sec2)
childrenSecListVector.o(sec2).append(sec1)
}
}
}
sNroot = strahlerNumbersFunc(0)
if (optPrint){
for sec1=0,100 /*NSEC-1*/{
sizeVect = childrenSecListVector.o(sec1).size()
if (sizeVect==0){
print "(sec,secname,sNum): ", sec1,secNameList0.o(sec1).s,"Leaf",strahlerNumbersVect.x(sec1)
}else{
printf("\nsec,secname,sNum,children): (%d,%s,%d", sec1,secNameList0.o(sec1).s,strahlerNumbersVect.x(sec1))
for ii=0,sizeVect-1 printf(",%d",childrenSecListVector.o(sec1).x(ii))
printf("\n")
}
}
}
}
objref vbox00
objectvar sh00
vbox00 = new VBox()
proc coloringStrahler(){local sec
vbox00 = new VBox()
vbox00.intercept(1)
sh00 = new Shape()
vbox00.intercept(0)
vbox00.map("coloringStrahler", 200, 430, 300, 250)
sec = 0
NumSmooth = 0
NumSpiny = 0
sNorder1 = 0
sNorder2 = 0
sNorder3 = 0
sNorder4 = 0
sNorder5 = 0
sNorder6 = 0
sNorder7 = 0
sNorderMajor7 = 0
forall if (!issection("eqSection.*")){
if (strahlerNumbersVect.x(sec)==1) {NumSpiny = NumSpiny+1 sNorder1 = sNorder1+1 sh00.color(2)
}else if (strahlerNumbersVect.x(sec)==2) {NumSpiny = NumSpiny+1 sNorder2 = sNorder2+1 sh00.color(3)
}else if (strahlerNumbersVect.x(sec)==3) {NumSpiny = NumSpiny+1 sNorder3 = sNorder3+1 sh00.color(4)
}else if (strahlerNumbersVect.x(sec)==4) {NumSmooth = NumSmooth+1 sNorder4 = sNorder4+1 sh00.color(5)
}else if (strahlerNumbersVect.x(sec)==5) {NumSmooth = NumSmooth+1 sNorder5 = sNorder5+1 sh00.color(6)
}else if (strahlerNumbersVect.x(sec)==6) {NumSmooth = NumSmooth+1 sNorder6 = sNorder6+1 sh00.color(7)
}else if (strahlerNumbersVect.x(sec)==7) {NumSmooth = NumSmooth+1 sNorder7 = sNorder7+1 sh00.color(8)
}else if (strahlerNumbersVect.x(sec)>7) {NumSmooth = NumSmooth+1 sNorderMajor7 = sNorderMajor7+1 sh00.color(9)}
sec=sec+1
}
if (optPrint) print "(NumSmooth,NumSpiny,sNorder1,sNorder2,sNorder3,sNorder4,sNorder5,sNorder6,sNorder7,sNorderMajor7): ",\
NumSmooth,NumSpiny,sNorder1,sNorder2,sNorder3,sNorder4,sNorder5,sNorder6,sNorder7,sNorderMajor7
sh00.show(0)
}
objref vbox01
objectvar sh01
vbox01 = new VBox()
objectvar stimulatedGroup1Vect // Vector for Stimulated Group
proc coloringStimulatedClusters(){local sec
vbox01 = new VBox()
vbox01.intercept(1)
sh01 = new Shape()
vbox01.intercept(0)
vbox01.map("coloringStimulatedClusters", 0, 130, 270, 250)
sec = 0
forall if (!issection("eqSection.*")){
if (stimulatedGroup1Vect.x(sec)==1) {sh01.color(2)
}else if (stimulatedGroup1Vect.x(sec)==2) {sh01.color(3)
}else if (stimulatedGroup1Vect.x(sec)==3) {sh01.color(4)
}else if (stimulatedGroup1Vect.x(sec)==4) {sh01.color(5)
}else if (stimulatedGroup1Vect.x(sec)==5) {sh01.color(6)
}else if (stimulatedGroup1Vect.x(sec)==6) {sh01.color(7)
}else if (stimulatedGroup1Vect.x(sec)==7) {sh01.color(8)
}else if (stimulatedGroup1Vect.x(sec)==8) {sh01.color(9)
}else if (stimulatedGroup1Vect.x(sec)==9) {sh01.color(11)
}else if (stimulatedGroup1Vect.x(sec)==10) {sh01.color(12)
}else if (stimulatedGroup1Vect.x(sec)==11) {sh01.color(13)
}else if (stimulatedGroup1Vect.x(sec)==12) {sh01.color(14)
}else if (stimulatedGroup1Vect.x(sec)==13) {sh01.color(15)
}else if (stimulatedGroup1Vect.x(sec)==14) {sh01.color(16)
}else if (stimulatedGroup1Vect.x(sec)==15) {sh01.color(17)
}else if (stimulatedGroup1Vect.x(sec)==16) {sh01.color(18)
}else if (stimulatedGroup1Vect.x(sec)>16) {sh01.color(19)}
sec=sec+1
}
sh01.show(0)
}
indexStimulatedCluster = 1
objref vbox02
objectvar sh02
vbox02 = new VBox()
proc coloringSingleStimulatedCluster(){local nCluster,sec,col,i
vbox02 = new VBox()
vbox02.intercept(1)
sh02 = new Shape()
vbox02.intercept(0)
vbox02.map("coloringSingleStimulatedCluster", 300, 230, 500, 450)
sec = 0
//colNum = indexStimulatedCluster+1
forall if (!issection("eqSection.*")){
if (stimulatedGroup1Vect.x(sec)==indexStimulatedCluster) {sh02.color(2)}
sec=sec+1
}
sh02.show(0)
}
/*---------------------------------------------
Subroutine for setting
---------------------------------------------*/
objectvar factGmaxSynMorphVect
factGmaxSynMorphVect = new Vector(100,1.0) //set/reset
objref fIn
strdef directory
proc readingFile(){
fIn = new File()
fIn.chooser("", "Reading Morphologic Model", "", "Accept", "", ".\\morphologies")
if (fIn.chooser()) {
forall delete_section()
firstInit = 0
fIn.getname(strFile)
sprint(directory,"%s",fIn.dir())
}else firstInit = 1
}
proc PurkinjeSetting(){local sec
vbox00.unmap()
if (optPrint==1) print "AlePurkinjeFunctional Setting"
forall delete_section()
xopen(strFile) // read geometry file
firstInit = 0
LENGTHMERGINGserialMETHOD = -1
LENGTHMERGINGparallelMETHOD = 2
DIAMMERGINGserialMETHOD = 0
DIAMMERGINGparallelMETHOD = 0
RESCALEMOD = 0
CUTTINGFACT = 100 //5
PRESERVINGsomaMETHOD = 0
// PRESERVINGMETHOD = 0 //0
// RaMERGINGMETHOD=1
// SURFFACT_PARTIAL_TOTAL = 1
// SURFFACT_PARTIAL_METHOD=2
withLocation = 1
strahlerNumbersComputation()
sec = 0
forall if (!issection("eqSection.*")){
if (issection("soma.*")){
cm = cmSoma
Ra = Ra0
insert Leak gl_Leak = gl_LeakSoma el_Leak = el_Leak0
}else if (strahlerNumbersVect.x(sec)>3){
cm= cmSmoothSec
Ra = Ra0
insert Leak gl_Leak = gl_LeakSmoothSec el_Leak = el_Leak0
}else if (strahlerNumbersVect.x(sec)<=3){
cm = cmSpinySect
Ra = Ra0
insert Leak gl_Leak = gl_LeakSpinySect el_Leak = el_Leak0
}
sec = sec+1
}
geom_nseg()
totOr=0
forall {totOr=totOr+nseg}
if(optPrint) print "totSegIni: ", totOr
//coloringStrahler()
}
/*---------------------------------------------
VARIABLE DEFINITION
---------------------------------------------*/
objectvar order // number of sections between each section and soma
objectvar parsec // parentsection id for each section
// section data input
objectvar secri
objectvar secpathri
objectvar secpathri0
objectvar secpathImp
objectvar secpathImp0
objectvar secL
objectvar secRa
objectvar secCm
objectvar secdiam
objectvar secSurf
//objectvar secApicalBasal
// section data during the merging
objectvar mrgL
objectvar mrgdiam
objectvar mrgdiamSer
objectvar mrgdiamPar
objectvar mrgri
objectvar mrgri2
objectvar mrgSurf
// equivalent section data after the merging
objectvar eqL
objectvar eqdiam
objectvar eqdiamSer
objectvar eqdiamPar
objectvar eqRa
objectvar surfFactVect
objectvar eqri2
objectvar eqsecpathri
objectvar eqsecpathri0
objectvar eqsecpathImp
objectvar eqsecpathImp0
objectvar orMinpathri
objectvar orMaxpathri
//objectvar eqri
objectvar eqSurf
objectvar orSurf
objectvar orSurf1
objectvar orSurf2
objectvar surfFact
objectvar surfFactRmCm
// cluster data
objectvar clusterParent // link between each cluster and its parent
objectvar clusterParentPos // position of each cluster as regards its parent
objectvar clusterList // List of clusters
objectvar newClusterList // List of clusters after merging Y_group_sections
objectvar clusterMarker // Marker to identify each section belonging to its cluster
objectvar somaVector // Vector of sections belonging to the soma
//objectvar axonVector // Vector of sections belonging to the axon
/*---------------------------------------------
INIT PROCEDURE
---------------------------------------------*/
proc initDataStruct(){
order = new Vector(NSEC,0) //set/reset
parsec = new Vector(NSEC,0) //set/reset
secri = new Vector(NSEC,0) //set/reset
secpathri = new Vector(NSEC,0.0) //set/reset
secpathri0 = new Vector(NSEC,0.0) //set/reset
secpathImp = new Vector(NSEC,0.0) //set/reset
secpathImp0 = new Vector(NSEC,0.0) //set/reset
secL = new Vector(NSEC,0) //set/reset
secRa = new Vector(NSEC,0) //set/reset
secCm = new Vector(NSEC,0) //set/reset
secdiam = new Vector(NSEC,0) //set/reset
secSurf = new Vector(NSEC,0) //set/reset
// secApicalBasal = new Vector(NSEC,0) //set/reset if 1=Apical elseif 0=Basal
mrgL = new Vector(NSEC,0) //set/reset
mrgdiam = new Vector(NSEC,0) //set/reset
mrgdiamSer = new Vector(NSEC,0) //set/reset
mrgdiamPar = new Vector(NSEC,0) //set/reset
mrgri = new Vector(NSEC,0) //set/reset
mrgri2 = new Vector(NSEC,0) //set/reset
mrgSurf = new Vector(NSEC,0) //set/reset
clusterList = new List()
newClusterList = new List()
clusterMarker = new Vector(NSEC,-1) //set/reset
somaVector = new Vector(0,0)
// axonVector = new Vector(0,0)
}
/*---------------------------------------------
CELL ANALYSIS
---------------------------------------------*/
// Useful variables
objectvar secNameList
objectvar secName
secNameList=new List()
objectvar secaddress // section address
objectvar paraddress // parentsection address
// data to compute path from each section to the soma
objectvar rvp
objectvar slsoma
objectvar slroot
proc cellPurkAnalysis(){local sec,sec1,sec2, x
strdef sname
{NSEC=0 forall if (!issection("eqSection.*")){NSEC=NSEC+1}}
initDataStruct()
secaddress=new Vector(NSEC,0) //set/reset
paraddress=new Vector(NSEC,0) //set/reset
sec = 0
forall if (!issection("eqSection.*")){
sname =secname() secName = new String(sname) secNameList.append(secName)
secaddress.set(sec, this_section())
paraddress.set(sec, parent_section(0))
secRa.set(sec,Ra)
secL.set(sec,abs(L))
secCm.set(sec,cm)
secdiam.set(sec,diam)
mrgL.set(sec,abs(L))
// if (y3d(1)>=2) secApicalBasal.set(sec,1)
if (issection("soma.*")){
for(x) if(x>0) {
secSurf.set(sec,secSurf.x(sec)+area(x))
secri.set(sec, secri.x(sec)+ri(x))
}
//mrgdiam.set(sec,aaa/PI/L)
//mrgri.set(sec,secri.x(0))
//print sec,totSomaArea,mrgL.x(sec),mrgdiam.x(sec),secri.x(sec)
somaVector.append(sec)
} else {
{slsoma=new SectionList() {rvp=new RangeVarPlot("v")}
{soma rvp.begin(.5)} rvp.end(.5) rvp.list(slsoma)}
{slroot=new SectionList() {rvp=new RangeVarPlot("v")}
{ss=0 forsec slsoma{ {if(ss==1) rvp.begin(.5)} ss=ss+1}}
rvp.end(.5) rvp.list(slroot)}
forsec slroot {order.set(sec, order.x(sec)+1) /*print sec,order.x(sec) psection()*/}
for(x) if(x>0) {
secSurf.set(sec,secSurf.x(sec)+area(x))
secri.set(sec, secri.x(sec)+ri(x))
/*print sec,secri.x(sec)*/
}
if(secri.x(sec)>9999) secri.set(sec, -9999)
nsecPath = 0
forsec slroot { nsecPath=nsecPath+1
for(x) {if(x>0){secpathri.set(sec,secpathri.x(sec)+ri(x))}}
}
nsecPath1 = 0
forsec slroot {
nsecPath1=nsecPath1+1
for(x) {if(x>0 && nsecPath1<nsecPath) secpathri0.set(sec,secpathri0.x(sec)+ri(x))}
}
if (PRESERVINGMETHOD==0){ //PRESERVINGMETHOD = 0 by Destexhe
mrgdiamSer.set(sec,secdiam.x(sec))
mrgdiamPar.set(sec,secdiam.x(sec)^2)
//print mrgdiamPar.x(sec)
mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI))
}else if (PRESERVINGMETHOD==1){
mrgdiam.set(sec,secSurf.x(sec)/secL.x(sec)/PI)
}
mrgri.set(sec,secri.x(sec))
mrgri2.set(sec,secri.x(sec))
mrgSurf.set(sec, secSurf.x(sec))
}
sec=sec+1}
for sec1=1,NSEC-1 {
for sec2=0,NSEC-1 {
if(secaddress.x(sec2)==paraddress.x(sec1)) {parsec.set(sec1,sec2)}
}
}
}
/*---------------------------------------------
CLUSTERISING PROCEDURES
---------------------------------------------*/
xopen("clusterisingMethods.hoc")
/*---------------------------------------------
MERGING PROCEDURES
---------------------------------------------*/
xopen("mergingMethods.hoc")
/*---------------------------------------------
REDUCTION PROCEDURES
---------------------------------------------*/
// Useful variables
objectvar searchAvailable
objectvar mergeAvailable
// seaching the next Y_group_sections / output: {secA,secB,secP}
proc getNextY(){local ii, nsect, idA, idB, sec //input: Vector of id sections
// seaching the last available section
secA=0
nsect = $o1.size()-1
idA = 0
for ii=0,nsect {
if(searchAvailable.x(ii)){
sec = $o1.x(ii)
if(order.x(sec)>=order.x(secA)) {secA=sec idA = ii}
}
}
// seaching the brother section of secA
secB=0
idB = 0
for ii=0,nsect {
if(searchAvailable.x(ii)){
sec = $o1.x(ii)
if(parsec.x(sec)==parsec.x(secA)&&sec!=secA) {secB=sec idB = ii}
}
}
// seaching the parent section of secA
secP=0
if ($o1.contains(parsec.x(secA))) secP=parsec.x(secA)
// marking unavailable sections already visited
searchAvailable.set(idA,0)
if (secB!=0) searchAvailable.set(idB,0)
if (secA!=0){
if (secP==0){
mergeAvailable.set(idA,1)
if (secB!=0) mergeAvailable.set(idB,1)
}
}
}
// Useful proc used in PurkReduction() after each call of getNextY() and before each call of merge2() or mergeY()
// input: {secA,secB,secP}
// output: {AmRa,AmL,Amdiam,Amri,BmRa,BmL,Bmdiam,Bmri,PmRa,PmL,Pmdiam,Pmri}
proc shortnms(){
{AmRa=0 AmL=0 Amdiam=0 Amri=0 Amsurf=0 AmdiamSer=0 AmdiamPar=0}
{BmRa=0 BmL=0 Bmdiam=0 Bmri=0 Bmsurf=0 BmdiamSer=0 BmdiamPar=0}
{PmRa=0 PmL=0 Pmdiam=0 Pmri=0 Pmsurf=0 PmdiamSer=0 PmdiamPar=0}
if(secA>0) {
AmRa=secRa.x(secA)
AmL=mrgL.x(secA)
Amdiam=mrgdiam.x(secA)
AmdiamSer=mrgdiamSer.x(secA)
AmdiamPar=mrgdiamPar.x(secA)
Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100)
Amsurf = AmL*Amdiam*PI
}
if(secB>0) {
BmRa=secRa.x(secB)
BmL=mrgL.x(secB)
Bmdiam=mrgdiam.x(secB)
BmdiamSer=mrgdiamSer.x(secB)
BmdiamPar=mrgdiamPar.x(secB)
Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100)
Bmsurf = BmL*Bmdiam*PI
}
if(secP>0) {
PmRa=secRa.x(secP)
PmL=mrgL.x(secP)
Pmdiam=mrgdiam.x(secP)
PmdiamSer=mrgdiamSer.x(secP)
PmdiamPar=mrgdiamPar.x(secP)
Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100)
Pmsurf = PmL*Pmdiam*PI
}
if(RESCALEMOD==1){
if(secA>0 && secB>0){
if(AmL>BmL){
{BmL=AmL mrgL.set(secB,BmL)}
if (PRESERVINGMETHOD==0){ //PRESERVINGMETHOD = 0 by Destexhe
Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100))
BmdiamSer=Bmdiam
BmdiamPar=Bmdiam
}else if (PRESERVINGMETHOD==1){
Bmdiam=Bmsurf/BmL/PI }
mrgdiam.set(secB,Bmdiam)
{Bmsurf = BmL*Bmdiam*PI}
}else if(BmL>AmL){
{AmL=BmL mrgL.set(secA,AmL)}
if (PRESERVINGMETHOD==0){ //PRESERVINGMETHOD = 0 by Destexhe
Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100))
AmdiamSer=Amdiam
AmdiamPar=Amdiam
}else if (PRESERVINGMETHOD==1){
Amdiam=Amsurf/AmL/PI }
mrgdiam.set(secA,Amdiam)
{Amsurf = AmL*Amdiam*PI}
}
}
}
}
// merging procedure of the section secA with its parent secP
proc mergingParentSec(){
mergingSerialMethod(AmRa,PmRa,AmL,PmL,Amdiam,Pmdiam,Amri,Pmri,Amsurf,Pmsurf,AmdiamSer,PmdiamSer,AmdiamPar,PmdiamPar)
mrgL.set(secP, newL)
mrgri.set(secP, newri)
mrgri2.set(secP, Amri+Pmri)
mrgSurf.set(secP, max(Amsurf,Pmsurf))
mrgdiam.set(secP, newdiam)
mrgdiamSer.set(secP, newdiamSer)
mrgdiamPar.set(secP, newdiamPar)
//print "newdiamPar",newdiamPar,"secP",secNameList.o(secP).s
}
// merging procedure of Y_group_sections
proc mergingYSec(){local AmL2BmL,BmL2AmL
AmL2BmL=AmL/BmL
BmL2AmL=BmL/AmL
if(AmL2BmL>=CUTTINGFACT){mergingSerialMethod(AmRa,PmRa,AmL,PmL,Amdiam,Pmdiam,Amri,Pmri,Amsurf,Pmsurf,AmdiamSer,PmdiamSer,AmdiamPar,PmdiamPar)}
if(BmL2AmL>=CUTTINGFACT){mergingSerialMethod(BmRa,PmRa,BmL,PmL,Bmdiam,Pmdiam,Bmri,Pmri,Bmsurf,Pmsurf,BmdiamSer,PmdiamSer,BmdiamPar,PmdiamPar)}
if(AmL2BmL<CUTTINGFACT && BmL2AmL<CUTTINGFACT){
mergingYMethod(AmRa,BmRa,PmRa,AmL,BmL,PmL,Amdiam,Bmdiam,Pmdiam,Amri,Bmri,Pmri,AmdiamSer,BmdiamSer,PmdiamSer,AmdiamPar,BmdiamPar,PmdiamPar)}
mrgL.set(secP, newL)
mrgri.set(secP, newri)
mrgri2.set(secP, newri2)
mrgSurf.set(secP, newSurf)
mrgdiam.set(secP, newdiam)
mrgdiamSer.set(secP, newdiamSer)
mrgdiamPar.set(secP, newdiamPar)
//print "newdiamPar da Y",newdiamPar,"secP",secNameList.o(secP).s
}
// Useful variables
objectvar clusterVectApp
objectvar mergeAvailableApp
// procedure to create a new Cluster List after merging Y_group_sections
proc newClusters(){local nsect, sec
{clusterVectApp = $o1 mergeAvailableApp = $o2}
nsect = clusterVectApp.size()
clusterVect = new Vector(0,0)
for ii=0,nsect-1 {
if(mergeAvailableApp.x(ii)) {/*print clusterVectApp.x(ii)*/ clusterVect.append(clusterVectApp.x(ii))}
}
newClusterList.append(clusterVect)
}
objectvar appVect
proc PurkReduction(){local i, ii, nsect
eqL=new Vector(clusterList.count,0) //set/reset
eqdiam=new Vector(clusterList.count,0) //set/reset
eqdiamSer=new Vector(clusterList.count,0) //set/reset
eqdiamPar=new Vector(clusterList.count,0) //set/reset
//eqRa=new Vector(clusterList.count,0) //set/reset
eqri2=new Vector(clusterList.count,0) //set/reset
eqSurf=new Vector(clusterList.count,0) //set/reset
orSurf=new Vector(clusterList.count,0) //set/reset
orSurf1=new Vector(clusterList.count,0) //set/reset
orSurf2=new Vector(clusterList.count,0) //set/reset
surfFact=new Vector(clusterList.count,1) //set/reset
surfFactRmCm=new Vector(clusterList.count,1) //set/reset
//eqri=new Vector(clusterList.count,0) //set/reset
newClusterList.append(somaVector)
mergingSoma(newClusterList.o(0))
iniClusterList = 1
for i=iniClusterList,clusterList.count-1 {
nsect = clusterList.o(i).size()
searchAvailable = new Vector(nsect,1)
mergeAvailable = new Vector(nsect,0)
getNextY(clusterList.o(i))
shortnms()
while(secA>0) {
if (secB==0 && secP!=0) mergingParentSec()
if (secB!=0 && secP!=0) mergingYSec()
//print secNameList.o(secA).s,secA,secNameList.o(secB).s,secB,secNameList.o(secP).s,secP
getNextY(clusterList.o(i))
shortnms()
}
newClusters(clusterList.o(i),mergeAvailable)
mergingAllMethod(newClusterList.o(i),i,nsect)
totSurf = 0
appVect = new Vector(0,0)
for ii=0,nsect-1 {
appVect.append(secSurf.x(clusterList.o(i).x(ii)))
totSurf = totSurf + secSurf.x(clusterList.o(i).x(ii))
}
//print totSurf
if (SURFFACT_PARTIAL_METHOD==0 || SURFFACT_PARTIAL_METHOD==3){
orSurf.set(i,appVect.sum())
surfFact.set(i,surfFactVect.x(i)*orSurf.x(i)/eqSurf.x(i))
surfFactRmCm.set(i,surfFactVect.x(i)*orSurf.x(i)/eqSurf.x(i))
}else if (SURFFACT_PARTIAL_METHOD==1 || SURFFACT_PARTIAL_METHOD==4){
orSurf.set(i,orSurf1.x(i))
surfFact.set(i,surfFactVect.x(i)*orSurf.x(i)/eqSurf.x(i))
surfFactRmCm.set(i,surfFactVect.x(i)*orSurf.x(i)/eqSurf.x(i))
//print surfFact.x(i),surfFactRmCm.x(i)
}
}
}
/*---------------------------------------------
OUTPUT DATA:
PRINT CLUSTERLIST
CREATING EQUIVALENT CELL
---------------------------------------------*/
// Useful variables
objectvar clusterListApp
proc printClusterList(){local nCluster,secApp,i,j,ii
clusterListApp = $o1
nCluster = clusterListApp.count
for i=0,nCluster-1{
print "_____________________"
print "cluster n. ", i
for j=0,clusterListApp.o(i).size()-1{
secApp = clusterListApp.o(i).x(j)
print secApp, secNameList.o(secApp).s, secL.x(secApp),secdiam.x(secApp),mrgdiamPar.x(secApp), secSurf.x(secApp),secpathri0.x(secApp),secpathri.x(secApp)
}
}
print "Section\tLength (um)\tDiameter (um)\tOld Area (um^2)\tNew Area (um^2)\tArea Fact"
for ii=0,nCluster-1 print ii, "\t", eqL.x(ii),"\t", eqdiam.x(ii),"\t",orSurf.x(ii),"\t",eqSurf.x(ii),"\t",surfFactRmCm.x(ii)
print "eqRa\t orMinpathri\t orMaxpathri\t eqri2\t eqriEffettiva\t eqpathri0\t eqpathri"
for ii=0,nCluster-1 {
eqSection[ii] {
if(RaMERGINGMETHOD==0){RaEq = factEqRa*eqRa.x(ii)}\
else if(RaMERGINGMETHOD==1){RaEq = factEqRa*eqri2.x(ii)*PI*(diam/2)^2*100/L}\
else if(RaMERGINGMETHOD==2){RaEq = factEqRa*(orMaxpathri.x(ii)-orMinpathri.x(ii))*PI*(diam/2)^2*100/L}
}
print RaEq,"\t",orMinpathri.x(ii),"\t",orMaxpathri.x(ii),"\t",eqri2.x(ii),"\t",(eqRa.x(ii)*eqL.x(ii))/(PI*((eqdiam.x(ii)^2)/4)*100),"\t",eqsecpathri0.x(ii),"\t",eqsecpathri.x(ii)
}
}
// data to compute path from each section to the soma
objectvar eqrvp
objectvar eqslsoma
objectvar eqslroot
{create eqSection[2000]}
factEqRa = 1
objref tgAlph
objref signDx
objref signDy
proc createPurkEqCell(){local nCluster,appLun,sec,col,i,j,jj,secApp
clusterListApp = $o1
nCluster = clusterListApp.count
eqsecpathri = new Vector(nCluster,0.0) //set/reset
eqsecpathri0 = new Vector(nCluster,0.0) //set/reset
orMinpathri = new Vector(nCluster,10^6) //set/reset
orMaxpathri = new Vector(nCluster,0.0) //set/reset
eqRa = new Vector(nCluster,0.0) //set/reset
TotOrSurf = orSurf.sum(1,orSurf.size-1)
TotEqSurf = eqSurf.sum(1,eqSurf.size-1)
for (i = 0; i < nCluster; i = i + 1){
for jj=0,clusterListApp.o(i).size()-1{
secApp = clusterListApp.o(i).x(jj)
eqRa.set(i,eqRa.x(i)+secRa.x(secApp))
if (secpathri.x(secApp)>orMaxpathri.x(i)){
orMaxpathri.set(i,secpathri.x(secApp))
}
if (secpathri0.x(secApp)<orMinpathri.x(i)){
orMinpathri.set(i,secpathri0.x(secApp))
}
}
if (clusterListApp.o(i).size()==1){
if (parsec.x(secApp)!=0){ //Mod Ale: era secApp invece di 0
orMinpathri.set(i,secpathri.x(parsec.x(secApp)))
}else{orMinpathri.set(i,0)}
orMaxpathri.set(i,secpathri.x(secApp))
}
appp = orMinpathri.x(i)
if (orMinpathri.x(i)>orMaxpathri.x(i)){
orMinpathri.x(i) = orMaxpathri.x(i)
orMaxpathri.x(i) = appp
}
if(optPrint) print "orMinpathri di ",i,"è: ", orMinpathri.x(i)
if(optPrint) print "orMaxpathri di ",i,"è: ", orMaxpathri.x(i)
eqRa.set(i,eqRa.x(i)/clusterListApp.o(i).size())
for (j = 1; j < nCluster; j = j + 1){
if (clusterParent.x(j)==i) {eqSection[i] connect eqSection[j](0), clusterParentPos.x(j)}
}
aaa = 0
if (SURFFACT_PARTIAL_TOTAL==0 && i>0) {surfFact.set(i,TotOrSurf/TotEqSurf) surfFactRmCm.set(i,TotOrSurf/TotEqSurf)}
}
for (j = 0; j < eqSomaVector.size(); j = j + 1){
jj = eqSomaVector.x(j)
eqSection[jj] {insert Leak el_Leak = el_Leak0 gl_Leak = factRm*surfFactRmCm.x(jj)*gl_LeakSoma el_Leak = el_Leak0
L = eqL.x(jj) diam = eqdiam.x(jj)
cm=factCm*cmSoma*surfFactRmCm.x(jj)
if($2==0){Ra = factEqRa*eqRa.x(jj)}\
else if($2==1){Ra = factEqRa*eqri2.x(jj)*PI*(diam/2)^2*100/L}\
else if($2==2){Ra = factEqRa*(orMaxpathri.x(jj)-orMinpathri.x(jj))*PI*(diam/2)^2*100/L}
}
}
for (j = 0; j < eqTrunkVector.size(); j = j + 1){
jj = eqTrunkVector.x(j)
eqSection[jj] {insert Leak el_Leak = el_Leak0 gl_Leak = factRm*surfFactRmCm.x(jj)*gl_LeakSmoothSec el_Leak = el_Leak0
L = eqL.x(jj) diam = eqdiam.x(jj)
cm=factCm*cmSmoothSec*surfFactRmCm.x(jj)
if($2==0){Ra = factEqRa*eqRa.x(jj)}\
else if($2==1){Ra = factEqRa*eqri2.x(jj)*PI*(diam/2)^2*100/L}\
else if($2==2){Ra = factEqRa*(orMaxpathri.x(jj)-orMinpathri.x(jj))*PI*(diam/2)^2*100/L}
}
}
for (j = 0; j < eqSmoothVector.size(); j = j + 1){
jj = eqSmoothVector.x(j)
eqSection[jj] {insert Leak el_Leak = el_Leak0 gl_Leak = factRm*surfFactRmCm.x(jj)*gl_LeakSmoothSec el_Leak = el_Leak0
L = eqL.x(jj) diam = eqdiam.x(jj)
cm=factCm*cmSmoothSec*surfFactRmCm.x(jj)
if($2==0){Ra = factEqRa*eqRa.x(jj)}\
else if($2==1){Ra = factEqRa*eqri2.x(jj)*PI*(diam/2)^2*100/L}\
else if($2==2){Ra = factEqRa*(orMaxpathri.x(jj)-orMinpathri.x(jj))*PI*(diam/2)^2*100/L}
}
}
for (j = 0; j < eqSpinyVector.size(); j = j + 1){
jj = eqSpinyVector.x(j)
eqSection[jj] {insert Leak el_Leak = el_Leak0 gl_Leak = factRm*surfFactRmCm.x(jj)*gl_LeakSpinySect el_Leak = el_Leak0
L = eqL.x(jj) diam = eqdiam.x(jj)
cm=factCm*cmSpinySect*surfFactRmCm.x(jj)
if($2==0){Ra = factEqRa*eqRa.x(jj)}\
else if($2==1){Ra = factEqRa*eqri2.x(jj)*PI*(diam/2)^2*100/L}\
else if($2==2){Ra = factEqRa*(orMaxpathri.x(jj)-orMinpathri.x(jj))*PI*(diam/2)^2*100/L}
}
}
totEq=0
for (i = 0; i < nCluster; i = i + 1){
eqSection[i] {
nseg = int((L/(d_lambda*lambda_f(freqEq))+0.9)/2)*2 + 1
totEq=totEq+nseg
}
}
if(optPrint) print "totSegEq: ", totEq
sec = 1
forsec "eqSection.*" if (!issection("eqSection[0]")){
{eqslsoma=new SectionList() eqrvp=new RangeVarPlot("v")
{eqSection[0] eqrvp.begin(.5)} eqrvp.end(.5) eqrvp.list(eqslsoma) }
{eqslroot=new SectionList() eqrvp=new RangeVarPlot("v")
{ss=0 forsec eqslsoma{ {if(ss==1) eqrvp.begin(.5)} ss=ss+1}}
eqrvp.end(.5) eqrvp.list(eqslroot)}
nsecPath = 0
forsec eqslroot {nsecPath=nsecPath+1 for(x) if(x>0) eqsecpathri.set(sec,eqsecpathri.x(sec)+ri(x))}
nsecPath1 = 0
forsec eqslroot {
nsecPath1=nsecPath1+1
for(x) if(x>0 && nsecPath1<nsecPath) eqsecpathri0.set(sec,eqsecpathri0.x(sec)+ri(x))
}
sec = sec+1
}
forsec "eqSection[0]" {for(x) if(x>0) eqsecpathri.set(0,eqsecpathri.x(0)+ri(x))}
tgAlph = new Vector(0,0)
signDx = new Vector(0,0)
signDy = new Vector(0,0)
access soma[0]
distance()
xx0 = x3d(1)
yy0 = y3d(1)
yyPrec = yy0
iii = 0
forsec TrunkSectList{
if (distance(0)==0){yyPrec = yy0}
if (x3d(1)-xx0==0){
tgAlph.append(1e10)
}else {tgAlph.append(abs((y3d(1)-yyPrec)/(x3d(1)-xx0)))}
if (x3d(1)-xx0>=0){signDx.append(1)}else{signDx.append(-1)}
if (y3d(1)-yyPrec>=0){signDy.append(1)}else{signDy.append(-1)}
//print tgAlph.x(iii)
iii =iii+1
yyPrec = yyPrec+y3d(1)-y3d(0)
}
access soma[0]
if (x3d(1)-x3d(0)==0){
tgAlph0 = 1e10
}else {tgAlph0 = abs((y3d(1)-y3d(0))/(x3d(1)-x3d(0)))}
appLun = eqSection[0].L
access eqSection[0]
pt3dclear()
DeltaY = appLun*sin(atan(tgAlph0))
DeltaX = appLun*cos(atan(tgAlph0))
pt3dadd(300,0,0,eqSection[0].diam)
pt3dadd(300+DeltaX,DeltaY,0,eqSection[0].diam)
for (j = 0; j < 1; j = j + 1){ //eqTrunkVector.size()
jj = eqTrunkVector.x(j)
jjParent = clusterParent.x(jj)
eqSection[jjParent]{xxParent = x3d(1) yyParent = y3d(1)}
eqSection[jj] {
pt3dclear()
DeltaY = L*sin(atan(tgAlph.x(j)))
DeltaX = L*cos(atan(tgAlph.x(j)))/10
pt3dadd(xxParent,yyParent,0,diam)
pt3dadd(xxParent+signDx.x(j)*DeltaX,yyParent+signDy.x(j)*DeltaY,0,diam)
}
}
}
// Useful variables
objectvar sh
objectvar shPlot
objref vbox0, vbox1
strdef app
objref strobj,rr,gg,bb
strobj = new StringFunctions()
proc coloringClusterList(){local nCluster,sec,col,i,j
clusterListApp = $o1
nCluster = clusterListApp.count
vbox0 = new VBox()
vbox1 = new VBox()
strdef strApp
vbox0.intercept(1)
sh = new Shape() //PlotShape
app = ""
sh.action("app = secname() secApp=-1 forall {secApp=secApp+1 if (issection(app)) break} print secApp,app coloringClusterList(clusterListApp)")
vbox0.intercept(0)
sprint(strApp, "Plotting %s", str)
vbox0.map(strApp, 0, 530, 270, 250)
sec = 0
forall if (!issection("eqSection.*")){
col = 1
for (i = 0; i < nCluster; i = i + 1){
//print "pippo",i,col,col-int(col/10)
if (col-int(col/10)*10==2 || col-int(col/10)*10==0) col = col-int(col/10)*10+1
//print "pluto",i,col
if (clusterListApp.o(i).contains(sec)&&!issection(app)) sh.color(col)
col = col + 1
}
sec=sec+1
}
col = 1
forsec "eqSection.*" {if (col-int(col/10)*10==2 || col-int(col/10)*10==0) {col = col-int(col/10)*10+1} sh.color(col) col=col+1}
length = strobj.len(app)
if (length>0) forsec "eqSection.*" if(issection(app)){sh.color(2)}
if(optPrint2) {
vbox1.intercept(1)
shPlot = new PlotShape()
shPlot.scale(-70, 0) //-70,30,-50
shPlot.exec_menu("Show Diam")
shPlot.exec_menu("Shape Plot")
vbox1.intercept(0)
vbox1.map(strApp, 400, 295, 500, 420)
}
sh.show(0)
}
/*---------------------------------------------
ACCURACY FUNCTION
---------------------------------------------*/
accu = -1
fact2 = 0.35
TmaxNeg = 10
TmaxPos = 10
// Useful variables
objectvar vecTTSpikeOrApp, vecTTSpikeEqApp
func accuracy(){local i, j, nn
tstop = t
vecTTSpikeOrApp = $o1
vecTTSpikeEqApp = $o2
TP = 0
TN = 0
FP = 0
FN = 0
if(vecTTSpikeOrApp.size()==0) vecTTSpikeOrApp.append(tstop)
if(vecTTSpikeEqApp.size()==0) vecTTSpikeEqApp.append(tstop)
for i = 0,vecTTSpikeOrApp.size()-1{
if (i == 0){tprec = 0
}else{tprec = vecTTSpikeOrApp.x(i - 1)}
if (i == vecTTSpikeOrApp.size()-1){tsucc = tstop
}else{tsucc = vecTTSpikeOrApp.x(i + 1)}
if (fact2*(vecTTSpikeOrApp.x(i) - tprec) <= TmaxPos){DeltaT1 = fact2*(vecTTSpikeOrApp.x(i) - tprec)
}else{DeltaT1 = TmaxPos}
if (i == 0){initNeg = 0
}else{initNeg = tprec + DeltaT1}
fineNeg = vecTTSpikeOrApp.x(i) - DeltaT1
initPos = fineNeg
if (fact2*(tsucc - vecTTSpikeOrApp.x(i)) <= TmaxPos){DeltaT2 = fact2*(tsucc - vecTTSpikeOrApp.x(i))
}else{DeltaT2 = TmaxPos}
finePos = vecTTSpikeOrApp.x(i) + DeltaT2
for (nn = initNeg; nn <= fineNeg; nn = nn + TmaxNeg) {
inizioInter = nn
if (nn + TmaxNeg <= fineNeg) {fineInter = nn + TmaxNeg
}else(fineInter = fineNeg)
nEq = 0
for j = 0,vecTTSpikeEqApp.size()-1{
if (vecTTSpikeEqApp.x(j) > inizioInter && vecTTSpikeEqApp.x(j) < fineInter) nEq = nEq+1
}
if (nEq == 0) {TN = TN + 1
}else{FP = FP + nEq}
}
nEq = 0
for j = 0,vecTTSpikeEqApp.size()-1{
if (vecTTSpikeEqApp.x(j) >= initPos && vecTTSpikeEqApp.x(j) <= finePos) nEq = nEq+1
}
if (nEq == 0) { FN = FN + 1
}else{TP = TP + 1 FP = FP + nEq - 1}
}
initNeg = finePos
fineNeg = tstop
for (nn = initNeg; nn <= fineNeg; nn = nn + TmaxNeg) {
inizioInter = nn
if (nn + TmaxNeg <= fineNeg) {fineInter = nn + TmaxNeg
}else(fineInter = fineNeg)
nEq = 0
for j = 0,vecTTSpikeEqApp.size()-1{
if (vecTTSpikeEqApp.x(j) > inizioInter && vecTTSpikeEqApp.x(j) < fineInter) nEq = nEq+1
}
if (nEq == 0) {TN = TN + 1
}else{FP = FP + nEq}
}
if(TP + TN + FP + FN != 0){ accu = (TP + TN)/(TP + TN + FP + FN) }
return accu
}
surfFactVect = new Vector(2000,1.0) //set/reset
