/*---------------------------------------------------------------------

   Written by Alessandro Limongiello
   University "Federico II" of Naples, 2012-2013

-----------------------------------------------------------------------*/


/*---------------------------------------------
                   CLUSTERISING PROCEDURES
  ---------------------------------------------*/

objectvar clusterVect         // Vector for the i-th cluster, containing the id-sections 
objectvar clusterParentVect   // Vector containing the parent sections for the i-th cluster
objectvar nclusterVect        // Vector containing the number of sections for each cluster
objectvar cutOffVect          // Vector containing the cutOff-values to clusterise

objectvar eqSomaVector
objectvar eqTrunkVector
objectvar eqSmoothVector
objectvar eqSpinyVector
objref sectListApp
objref sectListApp1
objref SomaSectList
objref TrunkSectList
objref TrunkSectListApp
objref SmoothSectList
objref SpinySectList
objref stimulationClusterList 
objref stimulationClusterVect


proc clusterisePurkinjeCell(){local ii,icluster,sec

    eqSomaVector = new Vector(1,0)
    eqTrunkVector = new Vector(0,0)    
    eqSmoothVector = new Vector(1,1)
    eqSpinyVector = new Vector(1,2)
    
    SomaSectList = new SectionList()
    TrunkSectList = new SectionList()
    SmoothSectList = new SectionList()
    SpinySectList = new SectionList()  


	clusterParent = new Vector(3,0)
	clusterParentPos = new Vector(3,0)

    // the first cluster is the soma
    clusterList.append(somaVector)      
            
    nclusterVect = new Vector(3,0)

    for ii=0,somaVector.size()-1{
	  clusterMarker.set(somaVector.x[ii], 0)
	  nclusterVect.x(0)=nclusterVect.x(0)+1        
    }

    sec = 0    
    forall{
       if (issection("soma.*")){
           SomaSectList.append()
       }else if (!issection("soma.*")){            
           if (strahlerNumbersVect.x(sec)>3){
              clusterMarker.set(sec, 1)
        	  nclusterVect.x(1)=nclusterVect.x(1)+1 
        	  SmoothSectList.append()
           }else{
              clusterMarker.set(sec, 2)
        	  nclusterVect.x(2)=nclusterVect.x(2)+1  
              SpinySectList.append()          
           }
       }
       sec = sec+1    
    }
          
    clusterParent.set(1,0)
    clusterParentPos.set(1,1)    
    clusterParent.set(2,1)
    clusterParentPos.set(2,0.5)     
    
  
    for ii=1,2{
        clusterVect = new Vector(nclusterVect.x(ii),0)
        icluster=0
        for sec=1,NSEC-1{
            if (clusterMarker.x(sec)==ii){clusterVect.set(icluster, sec)  icluster=icluster+1}
        }
        
        if (icluster>0) {
            clusterList.append(clusterVect)
        }
    }  
        
    
}

sNumberLimit = 5
sizeLimit = 80
possSpiny = 1

proc clusterise1PurkinjeCell(){local ii,icluster,sec,sec1,secPar

    eqSomaVector = new Vector(0,0)
    eqTrunkVector = new Vector(0,0)    
    eqSmoothVector = new Vector(0,0)
    eqSpinyVector = new Vector(0,0)
    
    SomaSectList = new SectionList()
    TrunkSectList = new SectionList()
    TrunkSectListApp = new SectionList()
    SmoothSectList = new SectionList()
    SpinySectList = new SectionList()
    stimulationClusterList = new List()

	clusterParent = new Vector(0,0)
	clusterParentPos = new Vector(0,0)
    nclusterVect = new Vector(0,0)

    // the first cluster is the soma
    clusterList.append(somaVector)      
            

    for ii=0,somaVector.size()-1{
	  clusterMarker.set(somaVector.x[ii], 0)        
    }
    nclusterVect.append(somaVector.size())
    
    clusterParent.append(0)
    clusterParentPos.append(0)   
    eqSomaVector.append(0) 
    
    access soma[0]
    distance()
    
    sec = 0   
    cMark = 0 
    forall{
       if (issection("soma.*")){
           SomaSectList.append()
       }else if (!issection("soma.*")){             
           if (strahlerNumbersVect.x(sec)>=sNumberLimit || distance(0)==0){
                secPar = parsec.x(sec)
                //print sec,secPar
                if (clusterMarker.x(secPar)!=-1){
                    if (distance(0)==0){
                       clusterParent.append(clusterMarker.x(0)) 
                       clusterParentPos.append(0.5)
                    }else{
                       clusterParent.append(clusterMarker.x(secPar))
                       clusterParentPos.append(1)
                    }
                    
                    cMark = cMark +1
                    clusterMarker.x(sec)= cMark
                    nclusterVect.append(1)
                    eqTrunkVector.append(cMark)
                    TrunkSectList.append()
                }           
           }
       }
       sec = sec+1    
    }
      
 
    sec2 = 0
    nStimulationClusterList = 0
    forsec TrunkSectList{
    
        sectListApp = new SectionList()
        sectListApp.children()
        sectListApp.remove(SomaSectList)
        sectListApp.remove(TrunkSectList)
        
        sectListApp1 = new SectionList()
        
        if (optPrint==1) print "TrunkSect: ",secname()
        forsec sectListApp{
               if (optPrint==1) print "sectListApp: ",secname()
               sectListApp1.subtree()               
        }
               
        sec1 = 0
        okSmooth = 0
        okSpiny = 0
        stimulationClusterVect = new Vector()
        forall{
               ifsec sectListApp1{
                     if (strahlerNumbersVect.x(sec1)>3){
                          okSmooth = 1
                     }else if (strahlerNumbersVect.x(sec1)<=3){
                          stimulationClusterVect.append(sec1)
                          okSpiny = 1
                     } 
               }
               sec1 = sec1+1
        }
        
        if (nStimulationClusterList==0){
           if (stimulationClusterVect.size()>0){
              stimulationClusterList.append(stimulationClusterVect)
              nStimulationClusterList = 1
           }
        }else if (stimulationClusterVect.size()>0){
           if (stimulationClusterList.o(nStimulationClusterList-1).size() < sizeLimit){
               stimulationClusterList.o(nStimulationClusterList-1).append(stimulationClusterVect)                                                
           }else{
               stimulationClusterList.append(stimulationClusterVect)
               nStimulationClusterList = nStimulationClusterList+1           
           }
        }
        
        currCMark = cMark
        
        if (okSmooth){
           if (okSpiny){
                clusterParent.append(eqTrunkVector.x(sec2))
                clusterParentPos.append(1)   
                nclusterVect.append(0)
                eqSmoothVector.append(cMark+1)
                clusterParent.append(cMark+1)
                clusterParentPos.append(possSpiny)   //0.5 1
                nclusterVect.append(0)
                eqSpinyVector.append(cMark+2)  
                currCMark = cMark+2                    
           }else {
                clusterParent.append(eqTrunkVector.x(sec2))
                clusterParentPos.append(1)   
                nclusterVect.append(0)
                eqSmoothVector.append(cMark+1) 
                currCMark = cMark+1          
           }
        }else if (okSpiny){
                clusterParent.append(eqTrunkVector.x(sec2))
                clusterParentPos.append(1)   
                nclusterVect.append(0)
                eqSpinyVector.append(cMark+1) 
                currCMark = cMark+1
        }
        
        
        sec1 = 0
        //if (sec2<10) {print secname(),": "}
        forall{
            ifsec sectListApp1{
               //if (sec2<10) {print sec1,secNameList.o(sec1).s, cMark}
               if (strahlerNumbersVect.x(sec1)>3){  
                    clusterMarker.x(sec1)= cMark+1
                    nclusterVect.x(cMark+1)=nclusterVect.x(cMark+1)+1
                    SmoothSectList.append()           
               }else if (strahlerNumbersVect.x(sec1)<=3){
                    clusterMarker.x(sec1)= currCMark
                    nclusterVect.x(currCMark)=nclusterVect.x(currCMark)+1
                    SpinySectList.append()          
               }
           }
           sec1 = sec1+1 
        }    
        cMark = currCMark
        sec2 = sec2+1 
    }
    
    for ii=1,nclusterVect.size()-1{
        clusterVect = new Vector(nclusterVect.x(ii),0)
        icluster=0
        for sec=1,NSEC-1{
            if (clusterMarker.x(sec)==ii){clusterVect.set(icluster, sec)  icluster=icluster+1}
        }
        //print ii,nclusterVect.x(ii),icluster,clusterParent.x(ii)
        if (icluster>0) {
            clusterList.append(clusterVect)
        }
    }  
    
    for ii = 0,stimulationClusterList.count-1{
        if (optPrint==1) print "stimulationCluster n.",ii," size: ",stimulationClusterList.o(ii).size()
        for jj = 0,stimulationClusterList.o(ii).size()-1{
               stimulatedGroupVect.x(stimulationClusterList.o(ii).x(jj)) = ii+1
               //print stimulationClusterList.o(ii).x(jj)
        }        
    }
    
    

}

proc stimulationClusterise(){local ii,icluster,sec,sec1,secPar

    stoprun = 1
    if (optPrint==1) print "Stimulation Clusterise"

    if (firstInit ==1){
        forall delete_section()                       
        readingFile() 
        xopen(strFile)		// read geometry file
        firstInit == 0
    }else{xopen(strFile)}
    
    if(firstInit == 0){
        PurkinjeSetting()
        cellPurkAnalysis()
        clusterise1PurkinjeCell()   
        firstInit = -1
    }
    
}



objref stimulationCluster2List 
objref stimulationCluster2Vect

nStimulationClusterList2 = 3
objref stimulationPercVect
stimulationPercVect = new Vector(nStimulationClusterList2+1,1/nStimulationClusterList2)
stimulationPercVect.x(0) = 0
stimulationPercVect.x(1) = 0.5
stimulationPercVect.x(2) = 0.7
stimulationPercVect.x(3) = 1


proc stimulationClusterise2(){local ii,sec,maxDist

    stoprun = 1
    if (optPrint==1) print "Stimulation Clusterise2"
    if (firstInit ==1){strahlerNumbersComputation()}
    
        
    if(firstInit == 0){
        stimulationCluster2List = new List()
        
        access soma[0]
        distance()
                       
        for (ii=0;ii<nStimulationClusterList2;ii=ii+1){
            stimulationCluster2Vect = new Vector()
            stimulationCluster2List.append(stimulationCluster2Vect)
        }
            
        maxDist = 0
        forall if (!issection("eqSection.*")){
               if (distance(1)>maxDist) maxDist = distance(1)
        }
        //print maxDist
                      
            
        sec = 0  
        forall if (!issection("eqSection.*")){
               for (ii=0;ii<nStimulationClusterList2;ii=ii+1){ 
                   if (distance(1)> stimulationPercVect.x(ii)*maxDist && distance(1)<= stimulationPercVect.x(ii+1)*maxDist \
                                    && strahlerNumbersVect.x(sec)<=3 && !issection("soma.*")){
                           stimulationCluster2List.o(ii).append(sec)
                           break
                   }
               }              
               sec = sec + 1
        }
        
        
        
 
        stimulatedGroup1Vect=new Vector(sec,0)		    //set/reset
        
        for ii = 0,stimulationCluster2List.count-1{
            //print "stimulationCluster n.",ii," size: ",stimulationCluster2List.o(ii).size()
            for jj = 0,stimulationCluster2List.o(ii).size()-1{
                   stimulatedGroup1Vect.x(stimulationCluster2List.o(ii).x(jj)) = ii+1
                   //print stimulationCluster2List.o(ii).x(jj)
            }        
        } 
 
        firstInit = -1
    }
    
}


proc findSimAxes(){local mm,x0,y0,x1,y1,dd,alph

     mm = 1
     access soma[0]
     x0 = x3d(0.5)
     y0 = y3d(0.5)
     minDist = 1e10
     mmOk = mm
     for (alph=0; alph<180; alph = alph + 5){
         numPos = 0
         numNeg = 0
         if (alph!=90){
             mm = sin(alph)/cos(alph)
             forall if (!issection("eqSection.*")){
                    x1 = x3d(0.5)
                    y1 = y3d(0.5)
                    dd = (y1-mm*x1-(y0-mm*x0))/(sqrt(1+mm^2))
                    if(dd>=0){numPos = numPos+1
                    }else{numNeg = numNeg+1}
             }
         }
         if (abs(numPos-numNeg) < minDist && abs(numPos-numNeg)!=0){
            minDist = abs(numPos-numNeg)
            mmOk = mm   
         }
         print numPos,numNeg,minDist,mmOk
     }

}


objref stimulationCluster1List 
objref stimulationCluster1Vect

nStimulationClusterList1 = 4

proc stimulationClusterise1(){local ii,sec,xx0,yy0,xx1,yy1,divx,divy,xx,yy, indd

    stoprun = 1
    if (optPrint==1) print "Stimulation Clusterise1"
    if (firstInit ==1){strahlerNumbersComputation()}
    
    if(firstInit == 0){
        stimulationCluster1List = new List()
        
        xx0 = 1e10   
        yy0 = 1e10
        xx1 = -1e10
        yy1 = -1e10
        
        forall if (!issection("eqSection.*")){
          if (x3d(0)<xx0) xx0 = x3d(0)
          if (x3d(1)<xx0) xx0 = x3d(1)
          if (y3d(0)<yy0) yy0 = y3d(0)
          if (y3d(1)<yy0) yy0 = y3d(1)
          if (x3d(0)>xx1) xx1 = x3d(0)
          if (x3d(1)>xx1) xx1 = x3d(1)
          if (y3d(0)>yy1) yy1 = y3d(0)
          if (y3d(1)>yy1) yy1 = y3d(1)                  
        }

        //print xx0,yy0,xx1,yy1
        divx = 1
        divy = 1
        if (nStimulationClusterList1==1) {divx = 1 divy = 1
        }else if (nStimulationClusterList1==2) {divx = 1 divy = 2
        }else if (nStimulationClusterList1==3) {divx = 1 divy = 3
        }else if (nStimulationClusterList1==4) {divx = 2 divy = 2
        }else if (nStimulationClusterList1==5) {nStimulationClusterList1=4 divx = 2 divy = 2
        }else if (nStimulationClusterList1==6) {divx = 2 divy = 3
        }else if (nStimulationClusterList1==7) {nStimulationClusterList1=6 divx = 2 divy = 3
        }else if (nStimulationClusterList1==8) {divx = 2 divy = 4
        }else if (nStimulationClusterList1==9) {divx = 3 divy = 3
        }else if (nStimulationClusterList1==10) {nStimulationClusterList1=9 divx = 3 divy = 3
        }else if (nStimulationClusterList1==12) {divx = 3 divy = 4
        }else if (nStimulationClusterList1==11) {nStimulationClusterList1=12 divx = 3 divy = 4
        }else {nStimulationClusterList1=12 divx = 3 divy = 4}
        
        for (ii=0;ii<nStimulationClusterList1;ii=ii+1){
            stimulationCluster1Vect = new Vector()
            stimulationCluster1List.append(stimulationCluster1Vect)
        }
            
        sec = 0  
        forall if (!issection("eqSection.*")){
               indd = 0
               inddOk = 0
               for (yy=yy0+(yy1-yy0)/(2*divy);yy<=yy1;yy=yy+(yy1-yy0)/divy){
                   for (xx=xx0+(xx1-xx0)/(2*divx);xx<=xx1;xx=xx+(xx1-xx0)/divx){   
                       if(x3d(0.5)>=xx-(xx1-xx0)/(2*divx) && x3d(0.5)<xx+(xx1-xx0)/(2*divx)\
                       && y3d(0.5)>=yy-(yy1-yy0)/(2*divy) && y3d(0.5)<yy+(yy1-yy0)/(2*divy)){ 
                          inddOk = indd
                       }                      
                       indd = indd + 1     
                   }
               }   
               if (!issection("soma.*") && strahlerNumbersVect.x(sec)<=3){
                  //print secname(),x3d(0.5),y3d(0.5),inddOk
                  stimulationCluster1List.o(inddOk).append(sec)
               }
               
               sec = sec + 1
        }
 
        stimulatedGroup1Vect=new Vector(sec,0)		    //set/reset
        
        for ii = 0,stimulationCluster1List.count-1{
            if (optPrint==1) print "stimulationCluster n.",ii," size: ",stimulationCluster1List.o(ii).size()
            for jj = 0,stimulationCluster1List.o(ii).size()-1{
                   stimulatedGroup1Vect.x(stimulationCluster1List.o(ii).x(jj)) = ii+1
                   //print stimulationCluster1List.o(ii).x(jj)
            }        
        } 
 
        firstInit = -1
    }
    
}


proc clusteriseByDestexhe() { local sec,icluster,ii
     
     
    eqSomaVector = new Vector(1,0)
    eqTrunkVector = new Vector(0,0)    
    eqSmoothVector = new Vector(1,1)
    eqSpinyVector = new Vector(1,2)
    
    SomaSectList = new SectionList()
    TrunkSectList = new SectionList()
    SmoothSectList = new SectionList()
    SpinySectList = new SectionList()     
     
     
	clusterParent = new Vector(3,0)
	clusterParentPos = new Vector(3,0)
                   
    nclusterVect = new Vector(3,0)
    cutOffVect = new Vector(4,0)
    cutOffVect.set(0,0)
    cutOffVect.set(1,CUTOFF1)
    cutOffVect.set(2,CUTOFF2)
    cutOffVect.set(3,1e10)
    
    for sec=0,NSEC-1{
            for ii=1,3{
        		if(secpathri.x(sec)>=cutOffVect.x(ii-1) && secpathri.x(sec)< cutOffVect.x(ii)) {
        		  clusterMarker.set(sec, ii-1)
        		  nclusterVect.x(ii-1)=nclusterVect.x(ii-1)+1
        		}            
            } 
    }  
    for ii=0,2{
        clusterVect = new Vector(nclusterVect.x(ii),0)
        icluster=0
        for sec=1,NSEC-1{
            if (clusterMarker.x(sec)==ii){clusterVect.set(icluster, sec)  icluster=icluster+1}
        }
        if (icluster>0) {
            clusterList.append(clusterVect)    
            clusterParent.set(ii,ii-1)
    	    clusterParentPos.set(ii,1)    
        }  
    }     
         
    sec = 0   
    forall{
           if(secpathri.x(sec)>=cutOffVect.x(0) && secpathri.x(sec)< cutOffVect.x(1)){ SomaSectList.append()
           }else if(secpathri.x(sec)>=cutOffVect.x(1) && secpathri.x(sec)< cutOffVect.x(2)){ SmoothSectList.append()
           }else {SpinySectList.append()}
       sec = sec+1    
    }         
         
         
}