/*---------------------------------------------
USEFUL PROCEDURES
---------------------------------------------*/
func max(){
if ($1 > $2){
return $1
} else {
return $2
}
}
objectvar factGmaxSynMorphVect
factGmaxSynMorphVect = new Vector(100,1.0) //set/reset
factGmaxSynMorphVect.x(0)= 340
factGmaxSynMorphVect.x(1)= 180
factGmaxSynMorphVect.x(2)= 160
factGmaxSynMorphVect.x(3)= 220
factGmaxSynMorphVect.x(4)= 200
factGmaxSynMorphVect.x(5)= 200
factGmaxSynMorphVect.x(6)= 200
factGmaxSynMorphVect.x(7)= 200
factGmaxSynMorphVect.x(8)= 240
factGmaxSynMorphVect.x(9)= 240
factGmaxSynMorphVect.x(10)= 180
factGmaxSynMorphVect.x(11)= 160
factGmaxSynMorphVect.x(12)= 200
factGmaxSynMorphVect.x(13)= 180
factGmaxSynMorphVect.x(14)= 200
factGmaxSynMorphVect.x(15)= 220
factGmaxSynMorphVect.x(16)= 200
strdef str,strFile
proc checkFile(){
if (strcmp(strFile,"./morphologies/geoc70963.hoc")==0) {iMorph = 0
}else if (strcmp(strFile,"./morphologies/geoc81462.hoc")==0) {iMorph = 1
}else if (strcmp(strFile,"./morphologies/geoc70863.hoc")==0) {iMorph = 2
}else if (strcmp(strFile,"./morphologies/geoc20466.hoc")==0) {iMorph = 3
}else if (strcmp(strFile,"./morphologies/geocd2351.hoc")==0) {iMorph = 4
}else if (strcmp(strFile,"./morphologies/geoc8076.hoc")==0) {iMorph = 5
}else if (strcmp(strFile,"./morphologies/geoc9236.hoc")==0) {iMorph = 6
}else if (strcmp(strFile,"./morphologies/geoc91665.hoc")==0) {iMorph = 7
}else if (strcmp(strFile,"./morphologies/geoc80761.hoc")==0) {iMorph = 8
}else if (strcmp(strFile,"./morphologies/geocd0351.hoc")==0) {iMorph = 9
}else if (strcmp(strFile,"./morphologies/geoc30465.hoc")==0) {iMorph = 10
}else if (strcmp(strFile,"./morphologies/geocd1152.hoc")==0) {iMorph = 11
}else if (strcmp(strFile,"./morphologies/geoc73162.hoc")==0) {iMorph = 12
}else if (strcmp(strFile,"./morphologies/geoc72965.hoc")==0) {iMorph = 13
}else if (strcmp(strFile,"./morphologies/geoc91662.hoc")==0) {iMorph = 14
}else if (strcmp(strFile,"./morphologies/geoc62564mod.hoc")==0) {iMorph = 15
}else if (strcmp(strFile,"./morphologies/geoc73166.hoc")==0) {iMorph = 16
}else {iMorph = 99}
maxStim = factGmaxSynMorphVect.x(iMorph)
if (iMorph == 99) maxStimOnOff = 0
}
objref f
proc readingFile(){
f = new File()
f.chooser("", "Reading Morphologic Model", "", "Accept", "", ".\\morphologies")
if (f.chooser()) {
firstInit = 0
f.getname(strFile)
checkFile()
}else firstInit = 1
}
proc AleFunctionalSetting(){
if (optPrint==1) print "AleFunctionalSetting Setting"
forall delete_section()
xopen(strFile) // read geometry file
firstInit = 0
LENGTHMERGINGserialMETHOD = -1
LENGTHMERGINGparallelMETHOD = 2
DIAMMERGINGserialMETHOD = 0
DIAMMERGINGparallelMETHOD = 0
RESCALEMOD = 0
CUTTINGFACT = 5
PRESERVINGsomaMETHOD = 0
PRESERVINGMETHOD = 0
RaMERGINGMETHOD=1
SURFFACT_PARTIAL_TOTAL = 1
SURFFACT_PARTIAL_METHOD=2
withLocation = 1
Ra0 = 150
RaAx = 50
g_pas0 = 1/28000 //Siemens/cm^2 1/30000
cm0 = 1
forall {insert pas e_pas=Vrest g_pas=g_pas0 Ra=Ra0 cm=cm0}
forsec "axon[0]" {insert pas e_pas=Vrest g_pas = g_pas0 Ra=RaAx cm=cm0}
geom_nseg()
totOr=0
forall {totOr=totOr+nseg}
if(optPrint) print "totSegIni: ", totOr
}
/*---------------------------------------------
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,0) //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 cellAnalysis(){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,L)
secCm.set(sec,cm)
secdiam.set(sec,diam)
mrgL.set(sec,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 if (issection("axon[0]")){
if(optPrint) print "pippo: ",sec
for(x) if(x>0) {
secSurf.set(sec,secSurf.x(sec)+area(x))
secri.set(sec, secri.x(sec)+ri(x))
}
axonVector.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 reduction() 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 reduction(){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
if (axonVector.size()>0){
newClusterList.append(axonVector)
mergingAxon(newClusterList.o(1))
iniClusterList = 2
}
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)
}
}
// 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
clusterListApp = $o1
nCluster = clusterListApp.count
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)}
sh.show(0)
}
// data to compute path from each section to the soma
objectvar eqrvp
objectvar eqslsoma
objectvar eqslroot
{create eqSection[100]}
factEqRa = 1//0.298
RaMERGINGMETHOD = 0
proc createEqCell(){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)!=secApp){
orMinpathri.set(i,secpathri.x(parsec.x(secApp)))
}else{orMinpathri.set(i,0)}
orMaxpathri.set(i,secpathri.x(secApp))
}
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)}
eqSection[i] {L = eqL.x(i) diam = eqdiam.x(i)
if($2==0){RaEq = factEqRa*eqRa.x(i)}\
else if($2==1){RaEq = factEqRa*eqri2.x(i)*PI*(diam/2)^2*100/L}\
else if($2==2){RaEq = factEqRa*(orMaxpathri.x(i)-orMinpathri.x(i))*PI*(diam/2)^2*100/L}
/*psection()*/ }
if (i>=0){
eqSection[i] {insert pas e_pas=Vrest \
g_pas=factRm*surfFactRmCm.x(i)*g_pas0 Ra=RaEq cm=factCm*cm0*surfFactRmCm.x(i)}
}else if(i==0 || i==1){
eqSection[i] {insert pas e_pas=Vrest \
g_pas=surfFactRmCm.x(i)*g_pas0 Ra=RaEq cm=surfFactRmCm.x(i)*cm0}
}
}
totEq=0
forsec "eqSection.*" {
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))}
appLun = eqSection[0].L
{access eqSection[0]}
{pt3dclear()}
{pt3dadd(300,0,0,eqSection[0].diam)}
{pt3dadd(300,appLun,0,eqSection[0].diam)}
appLun = eqSection[1].L
{access eqSection[1]}
{pt3dclear()}
{pt3dadd(300,0,0,eqSection[1].diam)}
{pt3dadd(300+appLun,0,0,eqSection[1].diam)}
appLun = eqSection[2].L
{access eqSection[2]}
{pt3dclear()}
{pt3dadd(300,0,0,eqSection[2].diam)}
{pt3dadd(300,appLun,0,eqSection[2].diam)}
vbox0 = new VBox()
vbox1 = new VBox()
strdef strApp
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)
}
if(optPrint1) {
vbox0.intercept(1)
sh = new Shape()
vbox0.intercept(0)
sprint(strApp, "Plotting %s", str)
vbox0.map(strApp, 0, 630, 300, 180)
app = ""
coloringClusterList(clusterListApp)
sh.action("app = secname() secApp=-1 forall {secApp=secApp+1 if (issection(app)) break} print secApp,app coloringClusterList(clusterListApp)")
sh.show(0)
}
}
objref ff,ffMorph
objectvar synSecVect, synPosVect, synGmaxVect
objectvar synSecMorphVect, synPosMorphVect, synGmaxMorphVect
proc onfile(){
nCluster = $o1.count
ff = new File("./outEq.txt")
ff.wopen("./outEq.txt")
ffMorph = new File("./outMorph.txt")
ffMorph.wopen("./outMorph.txt")
ff.printf("%s\n",strFile)
ffMorph.printf("%s\n",strFile)
ff.printf("%d\n",nCluster)
ind = 0
forsec "eqSection.*" {
ff.printf("%d\t%d\t%-8.8g\t%-8.8g\t%-8.8g\t%-8.8g\t%-8.8g\n",ind,nseg,L,Ra,diam,cm,g_pas)
ind=ind+1
}
ff.printf("%d\t%d\t%d\n",noise0,poissonSeed,positionSeed)
ffMorph.printf("%d\t%d\t%d\n",noise0,poissonSeed,positionSeed)
ff.printf("%d\n",synNum)
ffMorph.printf("%d\n",synNum)
//print synNum
if (synNum>0){
for(i = 0; i < synNum; i = i + 1){
ff.printf("%d\t%-8.8g\t%-8.8g\n",synSecVect.x(i),synPosVect.x(i),synGmaxVect.x(i))
ffMorph.printf("%d\t%-8.8g\t%-8.8g\n",synSecMorphVect.x(i),synPosMorphVect.x(i),synGmaxMorphVect.x(i))
}
}
for (i = 0; i < nCluster; i = i + 1){
ff.printf("%-8.8g\n",surfFact.x(i))
}
ffMorph.printf("%-6.6g\n",freq)
ff.close("./outEq.txt")
ffMorph.close("./outMorph.txt")
}
/*---------------------------------------------
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(100,1.0) //set/reset
