R/denpro.R

Defines functions addnode allokoi.new allokoi alloroot alpha.complex blokitus2 blokitus boundbox branchmap ccentebag ccentedya ccente cenone cfrekv change cintemul cinte colo2scem coloallo colobary.merge colobary.nodes colobary colobary.roots colors2data colorsofdata2 colorsofdata3 colorsofdata.adagrid2 colorsofdata.adagrid.new colorsofdata.adagrid colorsofdata.new colorsofdata colorsofdata.tail complex.rips cumu cutmut cutvalue cvolumbag cvolumdya cvolum declevdya declevgen declevnew declev decombag decomdya decom den2reg dend2parent dendat2lst depth2com depth digit dist.func dotouchgen dotouch drawgene drawhist draw.kern draw.levset draw.pcf epane epan etais etaisrec eva.clayton eva.cop6 eva.copula eva.gauss eva.hat eval.func.1D eval.func.dD eva.lognormal evanor eva.plackett eva.prod eva.skewgauss eva.student eva.t excmas.bylevel excmas exmap explo.compa findbnodes findbranch.pare findbranch findend findleafs findneighbor fs.calc.parti fs.calc graph.matrix.level graph.matrix hgrid histo1d histo2data histo intersec.edges intersec intersec.simpces2 intersec.simpces is.inside is.inside.simp.bary is.inside.simp.long is.inside.simp joincongen joinconne joingene kereva kergrid kernesti.dens lambda2emass leafsfirst.adagrid leafsfirst.bondary leafsfirst.complex leafsfirst.complex.volu leafsfirst.delaunay leafsfirst.intpol leafsfirst.intpol.volu leafsfirst.lst leafsfirst.new leafsfirst.nn leafsfirst leafsfirst.shape leafsfirst.tail leafsfirst.visu leikkaa levord liketree listchange locofmax lst2xy lstseq.kern makehis makeinfo makeparent massat massone maxnodenum modecent modegraph modetestgauss modetest modetestydin montecarlo.ball montecarlo.complex moodilkm mtest multitree negapart nn.indit nn.likeset nn.radit nnt omaind omamax omamin omaord2 omaord onko onsetissa paraclus paracoor.dens paracoor pcf.boundary pcf.func pcf.histo pcf.kernC pcf.kern pcf.kern.vech pcf.matrix perspec.dyna pituus plotbary plotbary.slide plotbranchmap plot.complex plotdata plotdelineator plotexmap plot.histdata plot.histo plotinfo plot.kernscale plotmodet plotprof plottext plottree plottwin plotvecs plotvolu2d plotvolu.new plotvolu point.eval posipart pp.plot preprocess prof2vecs profgene profhist profkernC profkernCRC profkern profkernR proftree proftreeR prunemodes pruneprof qq.plot quanti rota.seq rotation2d rotation3d rotation scaletable shape2d shapetree siborder.new siborder siborToModor sim.1d2modal sim.claw sim.cross sim.data sim.fox sim.fssk simmix1d simmix sim.mulmodII sim.mulmod sim.nested sim.peaks sim.penta4d sim.tetra3d simukern slicing sphere.map sphere.para stseq support tailfunc taillevel tail.plot.dens tail.plot til1 til2 til touchi.boundary touchi.dela touchi touchi.simp touchi.tail touch touchstep.boundary touchstep.complex touchstep.delaunay touchstep touchstep.tail toucrec travel.tree treedisc.ada treedisc.intpol treedisc tree.segme vectomatch volball vols.complex volsimplex voltriangle weightsit

Documented in blokitus branchmap colors2data colorsofdata draw.levset draw.pcf etais evanor excmas exmap findbnodes fs.calc.parti hgrid kernesti.dens leafsfirst leafsfirst.adagrid locofmax lstseq.kern modecent modegraph moodilkm paracoor pcf.func pcf.kern plotbary plotbranchmap plotexmap plotmodet plottree plotvolu plotvolu2d profgene profhist profkern proftree prunemodes scaletable shape2d sim.data slicing stseq treedisc tree.segme

addnode<-function(inde,curre,curdep,left,right,parent,low,upp,N,numnode){
#(inde,curre,curdep,left,right,deplink,low,upp,enofatdep,N,numnode){
#
#inde is d-vector: index (gridpoint) to be added
#curre is pointer to vectors left,right,...
#
d<-length(inde)
apu<-depth2com(curdep,N)
curdir<-apu$direc
depatd<-apu$depind
depit<-log(N,base=2)
#depit[d]<-depit[d]+1
#
while (curdir<=(d-1)){
    ind<-inde[curdir]
    while (depatd<=depit[curdir]){
        mid<-(low[curre]+upp[curre])/2
        if (ind<=mid){
           left[curre]<-numnode+1
           parent[numnode+1]<-curre
           low[numnode+1]<-low[curre]
           upp[numnode+1]<-floor(mid)
        }
        else{
           right[curre]<-numnode+1
           parent[numnode+1]<-curre
           low[numnode+1]<-ceiling(mid)
           upp[numnode+1]<-upp[curre]
        }
        numnode<-numnode+1
        curre<-numnode
        depatd<-depatd+1
        curdep<-curdep+1
#        deplink[endofatdep[curdep]]<-numnode
#        deplink[numnode]<-0
#        endofatdep[curdep]<-numnode
    }
    #
    # Last node of this dimension (first node of next dimension)
    #
    curdir<-curdir+1
    ind<-inde[curdir]
    low[curre]<-1
    upp[curre]<-N[curdir]
    mid<-(low[curre]+upp[curre])/2
    if (ind<=mid){
           left[curre]<-numnode+1
           parent[numnode+1]<-curre
           low[numnode+1]<-low[curre]
           upp[numnode+1]<-floor(mid)
    }
    else{
           right[curre]<-numnode+1
           parent[numnode+1]<-curre
           low[numnode+1]<-ceiling(mid)
           upp[numnode+1]<-upp[curre]
    }
    depatd<-2
    numnode<-numnode+1
    curre<-numnode
    curdep<-curdep+1
#    deplink[endofatdep[curdep]]<-numnode
#    deplink[curre]<-0
#    endofatdep[curdep]<-numnode
}
#
# Last dimension 
#
ind<-inde[curdir]
while (depatd<=depit[curdir]){
        mid<-(low[curre]+upp[curre])/2
        if (ind<=mid){
           left[curre]<-numnode+1
           parent[numnode+1]<-curre
           low[numnode+1]<-low[curre]
           upp[numnode+1]<-floor(mid)
        }
        else{
           right[curre]<-numnode+1
           parent[numnode+1]<-curre
           low[numnode+1]<-ceiling(mid)
           upp[numnode+1]<-upp[curre]
        }
        numnode<-numnode+1
        curre<-numnode
        depatd<-depatd+1
        curdep<-curdep+1
#        deplink[endofatdep[curdep]]<-numnode
#        deplink[curre]<-0
#        endofatdep[curdep]<-numnode
}
#
# Last node of last dimension
#
#left[curre]<-0
#right[curre]<-0
#
#return(list(numnode=numnode,left=left,right=right,deplink=deplink,low=low,
#upp=upp,endofatdep=endofatdep))
return(list(numnode=numnode,left=left,right=right,parent=parent,low=low,
upp=upp,nodeloc=numnode))
}










allokoi.new<-function(cur,vecs,lst,left,right,sibord)
{
# allocates space for all children of "cur"

# Calculate the number of childs and sum of volumes of childs
now<-left[cur]
childnum<-1
childvolume<-lst$volume[now]
while (right[now]>0){
  now<-right[now]
  childnum<-childnum+1
  childvolume<-childvolume+lst$volume[now]
}
 
gaplen<-(lst$volume[cur]-childvolume)/(childnum+1)

if (childnum==1){
   now<-left[cur]
   xbeg<-gaplen+vecs[cur,1]
   xend<-xbeg+lst$volume[now]
   ycoo<-lst$level[now]
   vecs[now,]<-c(xbeg,xend,ycoo)
}
else{
  siblinks<-matrix(0,childnum,1)  #make siblinks in right order
  now<-left[cur] 
  sior<-sibord[now]
  siblinks[sior]<-now
  while (right[now]>0){
    now<-right[now]
    sior<-sibord[now]
    siblinks[sior]<-now
  }
  xend<-vecs[cur,1]      #initialize xend 
  for (i in 1:childnum){
     now<-siblinks[i]
     xbeg<-gaplen+xend
     xend<-xbeg+lst$volume[now]
     ycoo<-lst$level[now]
     vecs[now,]<-c(xbeg,xend,ycoo)
  }
} 


return(vecs)
}


allokoi<-function(vecs,cur,child,sibling,sibord,levels,volumes)
{
#Finds coordinates of a node
#sibord,levels,volumes are nodenum-vector

# Calculate the number of childs and sum of volumes of childs
now<-child[cur]
childnum<-1
childvolume<-volumes[now]
while (sibling[now]>0){
  now<-sibling[now]
  childnum<-childnum+1
  childvolume<-childvolume+volumes[now]
}
 
gaplen<-(volumes[cur]-childvolume)/(childnum+1)

if (childnum==1){
   now<-child[cur]
   xbeg<-gaplen+vecs[cur,1]
   xend<-xbeg+volumes[now]
   ycoo<-levels[now]
   vecs[now,]<-c(xbeg,ycoo,xend,ycoo)
}
else{
  siblinks<-matrix(0,childnum,1)  #make siblinks in right order
  now<-child[cur] 
  sior<-sibord[now]
  siblinks[sior]<-now
  while (sibling[now]>0){
    now<-sibling[now]
    sior<-sibord[now]
    siblinks[sior]<-now
  }
  xend<-vecs[cur,1]      #initialize xend 
  for (i in 1:childnum){
     now<-siblinks[i]
     xbeg<-gaplen+xend
     xend<-xbeg+volumes[now]
     ycoo<-levels[now]
     vecs[now,]<-c(xbeg,ycoo,xend,ycoo)
  }
} 
return(vecs)
}



alloroot<-function(vecs,roots,sibord,levels,volumes)
{
rootnum<-length(roots)

# Calculate sum of volumes of roots
rootsvolume<-0
for (i in 1:rootnum){
  now<-roots[i]
  rootsvolume<-rootsvolume+volumes[now]
}

basis<-rootsvolume+rootsvolume/4
 
gaplen<-(basis-rootsvolume)/(rootnum+1)

rootlinks<-matrix(0,rootnum,1)  #make links in right order

if (rootnum==1) rootlinks[1]<-roots[1]  #1
else{
for (i in 1:rootnum){
  now<-roots[i]
  roor<-sibord[now]
  rootlinks[roor]<-now
}
}
xbeg<-0
xend<-0
for (i in 1:rootnum){
  now<-rootlinks[i]
  xbeg<-gaplen+xend
  xend<-xbeg+volumes[now]
  ycoo<-levels[now]
  vecs[now,]<-c(xbeg,ycoo,xend,ycoo)
}
return(vecs)
}
alpha.complex<-function(complex,dendat,alpha)
{
M<-dim(complex)[1]
n<-dim(dendat)[1]
d<-dim(dendat)[2]  # d<-dim(complex)[2]-1  

acomplex<-matrix(0,M,d+1)
lkm<-0
for (m in 1:M){
    simindex<-complex[m,]
    simplex<-dendat[simindex,]

    tulos<-0
    i<-1
    while ((i<=d) && (tulos==0)){
       v1<-simplex[i,]
       j<-i+1
       while ((j<=(d+1)) && (tulos==0)){
         v2<-simplex[j,]
         etais2<-sum((v1-v2)^2)
         if (etais2>alpha^2) tulos<-1
         j<-j+1
       }
       i<-i+1
    }
    if (tulos==0){ 
       lkm<-lkm+1
       acomplex[lkm,]<-complex[m,]
    }
}
acomplex<-acomplex[1:lkm,]

return(acomplex)
}

blokitus2<-function(obj,blokki){
#
sar<-length(obj[1,]) #sarakkeiden maara 
riv<-length(obj[,1]) #rivien maara 
#
uusobj<-matrix(0,riv,sar+blokki)
uusobj[,1:sar]<-obj
#
return(uusobj)
}
blokitus<-function(obj,blokki){
#
if (dim(t(obj))[1]==1) k<-1 else k<-length(obj[,1]) #rivien maara 
if (k==1){
  len<-length(obj)
  uusobj<-matrix(0,len+blokki,1)
  uusobj[1:len]<-obj
}
else{
  lev<-length(obj[1,])
  uusobj<-matrix(0,k+blokki,lev)
  uusobj[1:k,]<-obj
}
return(uusobj)
}
boundbox<-function(rec1,rec2)
{
# rec:s are 2*d-vectors

d<-length(rec1)/2
rec<-matrix(0,2*d,1)

for (i in 1:d){
    rec[2*i-1]<-min(rec1[2*i-1],rec2[2*i-1])
    rec[2*i]<-max(rec1[2*i],rec2[2*i])
}

return(rec)
}

branchmap<-function(estiseq,hseq=NULL,levnum=80,paletti=NULL,rootpaletti=NULL,
type="jump")
{
#type= "smooth", "jump", "diffe" 

if (is.null(paletti))
 paletti<-c("red","blue","green",
 "orange","navy","darkgreen",
 "orchid","aquamarine","turquoise",
 "pink","violet","magenta","chocolate","cyan",
 colors()[50:100])
if (is.null(rootpaletti)) rootpaletti<-colors()[102:110]

lstseq<-estiseq$lstseq
if (is.null(hseq))
   if (!is.null(estiseq$type)){
       if (estiseq$type=="bagghisto") hseq<--estiseq$hseq
       if (estiseq$type=="carthisto")  hseq<--estiseq$leaf
       if (estiseq$type=="kernel")  hseq<-estiseq$hseq
   }
   else hseq<-estiseq$hseq
hnum<-length(hseq)

if (hseq[1]>hseq[2]){  
    hseq<-hseq[seq(hnum,1)]
    apuseq<-list(lstseq[[hnum]])
    i<-2
    while (i <= hnum){
         apuseq<-c(apuseq,list(lstseq[[hnum-i+1]]))
         i<-i+1 
   }
   lstseq<-apuseq
}

maxlevel<-0
i<-1
while (i<=hnum){
   lst<-lstseq[[i]]
   maxlevel<-max(max(lst$level),maxlevel)
   i<-i+1
}

levstep<-maxlevel/(levnum-1)
level<-seq(0,maxlevel,levstep)
z<-matrix(0,length(level)+1,length(hseq)+1)
#col<-matrix("white",length(level),length(hseq))
colot<-matrix("white",(length(level))*(length(hseq)),1)

i<-1
while (i<=hnum){
    lst<-lstseq[[i]]  #[[hnum-i+1]]

    if ((type=="smooth") || (type=="diffe")){
          eb<-excmas.bylevel(lst,length(level)+1)
          if (type=="smooth")  z[,i]<-eb$levexc
          else z[,i]<-eb$diffe
    }
    mut<-multitree(lst$parent)
    ex<-excmas(lst)
    fb<-findbranch.pare(lst$parent)
    if (is.null(fb)) branchnum<-0 else branchnum<-length(fb)

    if (branchnum==0) toplevel<-max(lst$level) else toplevel<-min(lst$level[fb])
    # toplevel is the level of the next branch
    rootnum<-length(mut$roots)
    rootstep<-toplevel/rootnum
    ordroots<-order(ex[mut$roots])  #order(lst$level[mut$roots])

    exmassa<-0     
    k<-1
    while (k<=rootnum){
        ind<-mut$roots[ordroots[k]]
        exmassa<-exmassa+ex[ind]  
        k<-k+1
    }

    leveka<-1
    levend<-max(leveka,min(round(levnum*toplevel/maxlevel),levnum))
    curleveka<-leveka
    k<-1
    while (k<=rootnum){
        ind<-mut$roots[ordroots[k]]
        curexma<-ex[ind] 
curlevend<-max(curleveka,min(round(curleveka+(levend-leveka)*curexma/exmassa),levend))
        if (type=="jump") z[curleveka:curlevend,i]<-exmassa
        ##col[curleveka:curlevend,i]<-rootpaletti[k]
        aa<-(i-1)*(levnum)+min(curleveka,levnum)
        bb<-(i-1)*(levnum)+min(curlevend,levnum)
        colot[aa:bb]<-rootpaletti[k]
        curleveka<-curlevend+1
        k<-k+1
    }

    curlevel<-toplevel   # curlevel is the level of the previous branching
    ordbranches<-order(lst$level[fb])
    k<-1
    while (k<=branchnum){

        branchind<-ordbranches[k]
        branch<-fb[branchind]

        if (k==branchnum) toplevel<-max(lst$level) 
        else{
            nextbranch<-fb[ordbranches[k+1]]
            toplevel<-lst$level[nextbranch]
        }
        childnum<-2
        children<-c(mut$child[branch],mut$sibling[mut$child[branch]])
        ordchild<-order(ex[children])  #order(lst$level[children])

        exmassa<-0     
        l<-1
        while (l<=childnum){
             ind<-children[ordchild[l]]
             exmassa<-exmassa+ex[ind]  
             l<-l+1
        }

        leveka<-curlevend+1
        levend<-max(leveka,min(leveka+round(levnum*(toplevel-curlevel)/maxlevel),levnum))
        curleveka<-leveka
        l<-1
        while (l<=childnum){
           ind<-children[ordchild[l]]
           curexma<-ex[ind]
  curlevend<-max(curleveka,min(round(curleveka+(levend-leveka)*curexma/exmassa),levend)) 
           if (type=="jump") z[curleveka:curlevend,i]<-exmassa
           ##col[curleveka:curlevend,i]<-paletti[l]
           aa<-(i-1)*(levnum)+min(curleveka,levnum)
           bb<-(i-1)*(levnum)+min(curlevend,levnum)
           colot[aa:bb]<-paletti[l]
           curleveka<-curlevend+1
           l<-l+1
        }
        curlevel<-toplevel
        k<-k+1
    }
    i<-i+1
}

z[,dim(z)[2]]<-z[,dim(z)[2]-1]
z[dim(z)[1],]<-z[dim(z)[1]-1,]
hseq[length(hseq)+1]<-hseq[length(hseq)]+hseq[length(hseq)]-hseq[length(hseq)-1]
level[length(level)+1]<-level[length(level)]+level[length(level)]-level[length(level)-1]

z<-z/max(z)

# add one column to the matrix: a new first column

lisa<-1
zapu<-matrix(0,dim(z)[1],dim(z)[2]+lisa)
zapu[,1:lisa]<-0
zapu[,(lisa+1):(lisa+dim(z)[2])]<-z

yapu<-matrix(0,length(hseq)+lisa,1)
ystep<-hseq[2]-hseq[1]
yapu[lisa:1]<-seq(hseq[1]-ystep,hseq[1]-ystep*lisa,-ystep)
yapu[(lisa+1):(length(hseq)+lisa)]<-hseq

# add colors to the end 
levelo<-lisa*(dim(zapu)[1]-1)
colapu<-matrix("",length(colot)+levelo,1)
colapu[1:length(colot)]<-colot
colapu[(length(colot)+1):length(colapu)]<-colot[(length(colot)-levelo+1):length(colot)]

return(list(level=level,h=yapu,z=zapu,col=colapu))
}





ccentebag<-function(component,AtomlistAtom,AtomlistNext,low,upp,volume,
step,suppo)
{
d<-dim(low)[2]

componum<-length(component)
center<-matrix(0,componum,d)

for (i in 1:componum){
   curcente<-matrix(0,d,1)
   pointer<-component[i]
   while (pointer>0){
        atompointer<-AtomlistAtom[pointer]
        
        newcente<-matrix(0,d,1)
        for (j in 1:d){
            # calculate 1st volume of d-1 dimensional rectangle where
            # we have removed j:th dimension

            vol<-1
            k<-1
            while (k<=d){
               if (k!=j){
                  vol<-vol*(upp[atompointer,k]-low[atompointer,k])*step[k]
               }
               k<-k+1
            }

            ala<-suppo[2*j-1]+step[j]*low[atompointer,j]
            yla<-suppo[2*j-1]+step[j]*upp[atompointer,j]
            newcente[j]<-vol*(yla^2-ala^2)/2
        }

        curcente<-curcente+newcente
        pointer<-AtomlistNext[pointer]
   }
   center[i,]<-curcente/volume[i]
}
return(t(center))
}

ccentedya<-function(volofatom,component,AtomlistNext,AtomlistAtom,
volume,minim,h,delta,index,d){
#
componum<-length(component)
center<-matrix(0,componum,d)
#
for (i in 1:componum){
   curcente<-0
   pointer<-component[i]
   while (pointer>0){
        atompointer<-AtomlistAtom[pointer]
        inde<-index[atompointer,]
        newcente<-minim-h+delta*inde
        curcente<-curcente+newcente
        pointer<-AtomlistNext[pointer]
   }
   center[i,]<-volofatom*curcente/volume[i]
}
return(t(center))
}
ccente<-function(levels,items,mass){
#Calculates centers from a collection of level sets.
#center is 1st moment didided by volume.
#
#levels is tasolkm*N-matrix of 1:s and 0:s
#items is N*(2*d)-matrix
#mass is tasolkm-vector
#
#returns N*d-matrix of 1st moments.
#
N<-length(levels[,1])
d<-length(items[1,])/2
res<-matrix(0,N,d)
if (dim(t(levels))[1]==1) tasolkm<-1 else tasolkm<-length(levels[,1]) 
for (i in 1:tasolkm){
  lev2<-change(levels[i,])
  m<-length(lev2)
  vol<-matrix(0,d,1)
  for (j in 1:m){
    ind<-lev2[j]
    rec<-items[ind,]
    vol<-vol+cenone(rec)
  }
  res[i,]<-vol/mass[i]
}
return(t(res))
}

 

cenone<-function(rec){
#Calculates the 1st moment of a rectangle.
#
#rec is (2*d)-vector, represents rectangle in d-space
#Returns a d-vector.
#
d<-length(rec)/2
res<-matrix(0,d,1)
for (j in 1:d){
  apurec<-rec      #apurec such that is volume is equal to
  apurec[2*j-1]<-0 #volume of d-1 dimensional rectangle where
  apurec[2*j]<-1   #we have removed j:th dimension
  vajmas<-massone(apurec) 
  res[j]<-vajmas*(rec[2*j]^2-rec[2*j-1]^2)/2  
}
return(res) 
}

cfrekv<-function(levels,arvot){
#laskee tasojoukon osien frekvenssit
#arvo on reaaliluku
#kumu on kork*n-matriisi, n saraketta, kuvaa kork kpl:tta tasojoukon osia
#1 jos vastaava data-matriisin rivin indikoima pallo kuuluu tasojouon osaan
#muodostetaan matriisi, jonka 1. sarakkeessa "arvo", 
#2. sarakkeessa kunkin tasojoukon osan frekvenssi  
#ts. laskettu kuinka monesta pallosta tasojoukko on yhdistetty
#
tasolkm<-length(levels[,1])     #levels:n rivien maara
frek<-matrix(0,tasolkm,1)
a<-1
while (a<=tasolkm){
   frek[a]<-sum(levels[a,]*arvot)
   a<-a+1 
}
return(t(frek))
}







change<-function(levset){
#
#
len<-length(levset)
m<-sum(levset)
rindeksit<-matrix(0,m,1)
j<-1
for (i in 1:len){
    if (levset[i]==1){
       rindeksit[j]<-i
       j<-j+1
    }
}
return(rindeksit)
}

cintemul<-function(roots,child,sibling,volume,level){
#
#integrate function, over the level of roots, in the region of roots
#
itemnum<-length(child)
rootnum<-length(roots)
inte<-0
for (i in 1:rootnum){
    pino<-matrix(0,itemnum,1)
    valpino<-matrix(0,itemnum,1)  #level of parent
    pino[1]<-roots[i]
    valpino[1]<-0
    sibling[roots[i]]<-0
    #    
    pinin<-1
    while (pinin>0){
        cur<-pino[pinin]      #take from stack
        valcur<-valpino[pinin] 
        pinin<-pinin-1
        #
        if (level[cur]>0){
           inte<-inte+(level[cur]-max(valcur,0))*volume[cur]
        }
        #
        if (sibling[cur]>0){
              pinin<-pinin+1
              pino[pinin]<-sibling[cur]
              valpino[pinin]<-valcur
        }
        while (child[cur]>0){    #go to left and put right nodes to stack
              valcur<-level[cur]
              cur<-child[cur]
              #
              if (level[cur]>0){
                 inte<-inte+(level[cur]-max(valcur,0))*volume[cur]
              }
              #
              if (sibling[cur]>0){  #if cur has siblings
                 pinin<-pinin+1
                 pino[pinin]<-sibling[cur]
                 valpino[pinin]<-valcur
             }
        }
    }
}
#
return(inte)
}

cinte<-function(values,volumes,parents){
#Calculates the integral of a piecewise continuous function.
#
len<-length(values)
int<-0
for (i in len:1){
  par<-parents[i]
  if (par==0) valpar<-0 
  else valpar<-values[par]
  int<-int+volumes[i]*(values[i]-valpar)
}
return(int)
}

colo2scem<-function(sp,mt,ca)
{
#sp result of scemprof
#mt result of modegraph
#ca result of coloallo
#origlis translates h-values from sp terminology to mt terminology

len<-length(sp$bigdepths)
mtlen<-length(ca)
colors<-matrix("black",len,1)

for (i in 1:len){
  label<-sp$mlabel[i]       #label for mode
  if (label>0){ 
      smoot<-sp$smoot[i]    #smoothing paramter value/leafnum
      # we find the corresponding slot from "ca" where
      # label corresponds and smoothing parameter value corresponds
      run<-1
      koesmoot<-mt$ycoor[run]
      koelabel<-mt$mlabel[run] 
      while (((koesmoot!=smoot) || (koelabel!=label)) && (run<=mtlen)){
         run<-run+1
         koesmoot<-mt$ycoor[run]
         koelabel<-mt$mlabel[run] 
      }
      # we have found the slot
      colors[i]<-ca[run]
  }
}

return(colors)
}








coloallo<-function(mt,paletti=NULL)
{
# fast allocation of colors (matching of modes)
# mt is mode tree
# paletti gives a list of colors

if (is.null(paletti))
paletti<-c("red","blue","green","turquoise","orange","navy",
"darkgreen","orchid",colors()[50:100])

d<-dim(mt$xcoor)[2]

snum<-0
for (i in 1:length(mt$mlabel)){
  if (mt$mlabel[i]==1) snum<-snum+1
}

xcoor<-mt$xcoor
ycoor<-mt$ycoor
mlabel<-mt$mlabel
lenni<-length(ycoor)

colot<-matrix("",lenni,1)

# find the locations for the information for each h

low<-matrix(0,snum,1)
upp<-matrix(0,snum,1)
low[1]<-1
glob<-2
while ((glob<=lenni) && (mlabel[glob]!=1)){
       glob<-glob+1
}
upp[1]<-glob-1
# now glob is at the start of new block
i<-2
while (i<=snum){
   low[i]<-glob
   glob<-glob+1
   while ((glob<=lenni) && (mlabel[glob]!=1)){
       glob<-glob+1
   }
   upp[i]<-glob-1
   i<-i+1
}

# first we allocate colors for the largest h

run<-1  #low[1]
while (run<=upp[1]){
   colot[run]<-paletti[run]
   run<-run+1
}

firstnewcolo<-run

i<-2

while (i<=snum){
   prenum<-upp[i-1]-low[i-1]+1
   curnum<-upp[i]-low[i]+1

   smallernum<-min(prenum,curnum)
   greaternum<-max(prenum,curnum)

   if (prenum==smallernum){   
        bases<-i
        compa<-i-1
   }
   else{
        bases<-i-1
        compa<-i
   }

      dista<-matrix(NA,smallernum,greaternum)
      for (ap in low[bases]:upp[bases]){
         for (be in low[compa]:upp[compa]){
           if (d==1){
               curcenter<-xcoor[ap]
               precenter<-xcoor[be]
           }
           else{
               curcenter<-xcoor[ap,]
               precenter<-xcoor[be,]
           }
           dista[be-low[compa]+1,ap-low[bases]+1]<-etais(curcenter,precenter)
         }
      }

      match<-matrix(0,smallernum,1)  #for each mode the best match
      findtie<-TRUE

      # find the best match for all and check whether there are ties
      match<-matrix(0,smallernum,1)
      for (bm in 1:smallernum){
          minimi<-min(dista[bm,],na.rm=TRUE)
          match[bm]<-which(minimi==dista[bm,])[1]
      }
      findtie<-FALSE
      bm<-1
      while ((bm<=smallernum) && (findtie==FALSE)){
         koe<-match[bm]
         bm2<-bm+1
         while (bm2<=smallernum){
            if (koe==match[bm2]){
                  findtie<-TRUE
            }
            bm2<-bm2+1
         }
         bm<-bm+1
      }
    
      onkayty<-FALSE

      while (findtie){

      onkayty<-TRUE
      tiematch<-matrix(0,smallernum,1)
      
      # find the best match for all
      bestmatch<-matrix(0,smallernum,1)
      for (bm in 1:smallernum){
          allna<-TRUE
          am<-1
          while ((am<=greaternum) && (allna)){
             if (!is.na(dista[bm,am])) allna<-FALSE
             am<-am+1
          }
          if (!(allna)){
             minimi<-min(dista[bm,],na.rm=TRUE)
             bestmatch[bm]<-which(minimi==dista[bm,])[1]
          }
          else bestmatch[bm]<-match[bm]
      }

      # find the first tie
      findtie<-FALSE

      tieset<-matrix(0,smallernum,1)
      bm<-1
      while ((bm<=smallernum) && (findtie==FALSE)){
         koe<-bestmatch[bm]
         bm2<-bm+1
         while (bm2<=smallernum){
            if (koe==bestmatch[bm2]){
                  findtie<-TRUE
                  tieset[bm]<-1
                  tieset[bm2]<-1
            }
            bm2<-bm2+1
         }
         bm<-bm+1
      }

      # solve the first tie
      if (findtie==TRUE){
         numofties<-sum(tieset)
         kavelija<-0
         tiepointer<-matrix(0,numofties,1) 
         # find the second best
         secondbest<-matrix(0,smallernum,1)
         for (bm in 1:smallernum){
            if (tieset[bm]==1){
               redudista<-dista[bm,]
               redudista[bestmatch[bm]]<-NA
               minimi<-min(redudista,na.rm=TRUE)
               secondbest[bm]<-which(minimi==redudista)[1]

               kavelija<-kavelija+1
               tiepointer[kavelija]<-bm
            }
         }
         # try different combinations       
         # try all subsets of size 2 from the set of ties
         numofsubsets<-choose(numofties,2)
            #gamma(numofties+1)/gamma(numofties-2+1)
         valuelist<-matrix(0,numofsubsets,1)
         vinnerlist<-matrix(0,numofsubsets,1)
         matchlist<-matrix(0,numofsubsets,1)
         runneri<-1
         eka<-1
         while (eka<=numofties){
            ekapo<-tiepointer[eka]
            toka<-eka+1
            while (toka<=numofties){
               tokapo<-tiepointer[toka]
               # try combinations for this subset (there are 2)
               # 1st combination
               fvinner<-ekapo
               fvinnermatch<-bestmatch[fvinner]
               floser<-tokapo
               flosermatch<-secondbest[floser]
               fvalue<-dista[fvinner,fvinnermatch]+dista[floser,flosermatch]
                # 2nd combination
               svinner<-tokapo
               svinnermatch<-bestmatch[svinner]
               sloser<-ekapo
               slosermatch<-secondbest[sloser]
               svalue<-dista[svinner,svinnermatch]+dista[sloser,slosermatch]
               # tournament
               if (fvalue<svalue){
                   valuelist[runneri]<-fvalue
                   vinnerlist[runneri]<-fvinner
                   matchlist[runneri]<-fvinnermatch
               }
               else{ 
                   valuelist[runneri]<-svalue
                   vinnerlist[runneri]<-svinner
                   matchlist[runneri]<-svinnermatch
               }
               runneri<-runneri+1 
               # 
               toka<-toka+1
            }
            eka<-eka+1
         }
         minimi<-min(valuelist,na.rm=TRUE)
         bestsub<-which(minimi==valuelist)[1]
         vinnerson<-vinnerlist[bestsub]
         matcherson<-matchlist[bestsub]

         tiematch[vinnerson]<-matcherson
         dista[vinnerson,]<-NA
         dista[,matcherson]<-NA

      }

      }  #while (findtie)

      if (onkayty){  #there was one tie
          
          for (sepo in 1:smallernum){
               if (tiematch[sepo]!=0) match[sepo]<-tiematch[sepo]
               else match[sepo]<-bestmatch[sepo]
          }
      }

      # finally allocate colors
      run<-1
      while (run<=smallernum){
          
          if (prenum==smallernum){
             xind<-run
             yind<-match[xind]
          }
          else{
             yind<-run
             xind<-match[yind]
          }

          colot[low[i]+yind-1]<-colot[low[i-1]+xind-1]    
          run<-run+1
      }
                    
      if (prenum<greaternum){

        run<-low[bases]
        while (run<=upp[bases]){
            if (colot[run]==""){
               colot[run]<-paletti[firstnewcolo]
               firstnewcolo<-firstnewcolo+1   
            }
            run<-run+1
        }

     }

     i<-i+1
}

return(colot)
}


















colobary.merge<-function(parent,level,colothre=min(level),paletti=NULL)
{
mt<-multitree(parent) #roots<-mt$roots child<-mt$child sibling<-mt$sibling

itemnum<-length(mt$child)
rootnum<-length(mt$roots)
if (is.null(paletti))
 paletti<-c("red","blue","green",
 "orange","navy","darkgreen",
 "orchid","aquamarine","turquoise",
 "pink","violet","magenta","chocolate","cyan",
 colors()[50:657],colors()[50:657])
hep<-1

colot<-matrix("",itemnum,1)
col<-colobary(parent,paletti)

for (i in 1:rootnum){
   curroot<-mt$roots[i]
   colot[curroot]<-col[curroot]  #"grey"  #paletti[hep]
   hep<-hep+1
   if (mt$child[curroot]>0){
      pino<-matrix(0,itemnum,1)
      pino[1]<-mt$child[curroot]
      pinin<-1
      while (pinin>0){
          cur<-pino[pinin]      #take from stack
          pinin<-pinin-1
          #if (level[mt$child[cur]]>colothre)
          #if (level[parent[cur]]>colothre)
          if (level[cur]>colothre)
              colot[cur]<-colot[parent[cur]]
          else{ 
                colot[cur]<-col[cur]  #"grey"  #paletti[hep]
                hep<-hep+1
          } 
          # put to the stack 
          if (mt$sibling[cur]>0){
                 pinin<-pinin+1
                 pino[pinin]<-mt$sibling[cur]
          }
          # go to left and put right nodes to the stack
          while (mt$child[cur]>0){   
             cur<-mt$child[cur]
             #if (level[mt$child[cur]]>colothre)
             #if (level[parent[cur]]>colothre)
             if (level[cur]>colothre) 
                       colot[cur]<-colot[parent[cur]]
             else{ 
                    colot[cur]<-col[cur]  #"grey"  #paletti[hep]
                    hep<-hep+1
             } 
             if (mt$sibling[cur]>0){
                 pinin<-pinin+1
                 pino[pinin]<-mt$sibling[cur]
             }
           }
       }#while (pinin>0)
   }
}       

ind<-(level<colothre)
colot[ind]<-"grey"
                    
return(colot)
}


colobary.nodes<-function(parent,nodes,paletti=NULL)
{
mt<-multitree(parent) #roots<-mt$roots child<-mt$child sibling<-mt$sibling

itemnum<-length(mt$child)
rootnum<-length(mt$roots)
if (is.null(paletti))
 paletti<-c("red","blue","green",
 "orange","navy","darkgreen",
 "orchid","aquamarine","turquoise",
 "pink","violet","magenta","chocolate","cyan",
 colors()[50:657],colors()[50:657])

colot<-matrix("",itemnum,1)
col<-colobary(parent,paletti)

nodenum<-length(parent)
allnodes<-matrix(0,nodenum,1)
#allnodes[1:length(nodes)]<-nodes
counter<-0  #length(nodes)
for (i in 1:length(nodes)){
    node<-nodes[i]
    tt<-travel.tree(parent,node)
    allnodes[(counter+1):(counter+length(tt))]<-tt
    counter<-counter+length(tt)
}
allnodes<-allnodes[1:counter]

for (i in 1:rootnum){
   curroot<-mt$roots[i]
   colot[curroot]<-col[curroot]  #"grey"  #paletti[hep]
   if (mt$child[curroot]>0){
      pino<-matrix(0,itemnum,1)
      pino[1]<-mt$child[curroot]
      pinin<-1
      while (pinin>0){
          cur<-pino[pinin]      #take from stack
          pinin<-pinin-1
          if (sum(cur==allnodes)>0)
              colot[cur]<-colot[parent[cur]]
          else{ 
                colot[cur]<-col[cur]  #"grey"  #paletti[hep]
          } 
          # put to the stack 
          if (mt$sibling[cur]>0){
                 pinin<-pinin+1
                 pino[pinin]<-mt$sibling[cur]
          }
          # go to left and put right nodes to the stack
          while (mt$child[cur]>0){   
             cur<-mt$child[cur]
             if (sum(cur==allnodes)>0)
                       colot[cur]<-colot[parent[cur]]
             else{ 
                    colot[cur]<-col[cur]  #"grey"  #paletti[hep]
             } 
             if (mt$sibling[cur]>0){
                 pinin<-pinin+1
                 pino[pinin]<-mt$sibling[cur]
             }
           }
       }#while (pinin>0)
   }
}       

ind<-setdiff(seq(1:itemnum),allnodes)
colot[ind]<-"grey"
                    
return(colot)
}


colobary<-function(parent,paletti,roots=NULL,
modecolo=NULL,modepointer=NULL #,segtype="char"
)
{
nodenum<-length(parent)
#if (segtype=="char") colot<-matrix("",nodenum,1) 
#else 
colot<-matrix(0,nodenum,1)

fb<-findbranch(parent)$indicator
modloc<-moodilkm(parent)$modloc
#if (repretype=="B"){
#   fb<-findbranchB(parent,roots)$indicator
#   modloc<-moodilkmB(parent)$modloc
#}

moodilkm<-length(modloc)
palerun<-0

# first allocate colors for modes
if (is.null(modecolo)){
   i<-1
   while (i<=moodilkm){
       cur<-modloc[i]
       palerun<-palerun+1
       colot[cur]<-paletti[palerun]
       i<-i+1
   }
}
else{
   # remove modecolo:s from paletti
   indu<-0
   for (pp in 1:length(paletti)) 
       for (ppp in 1:length(modecolo))
          if (paletti[pp]==modecolo[ppp]){ 
                 indu<-indu+1
                 paletti[pp]<-colors()[100+indu] 
          }
   
   i<-1
   while (i<=moodilkm){
       cur<-modepointer[i]
       colot[cur]<-modecolo[i]
       i<-i+1
   } 
}

# then allocate for others
i<-1
while (i<=moodilkm){
 
  cur<-modloc[i]
  while (parent[cur]>0){

     child<-parent[cur]

     if ((fb[cur]==1) && (colot[child]==0)){ #cur is a result of a branch 
           palerun<-palerun+1
           colot[child]<-paletti[palerun]
     }      
     else if (colot[child]==0) colot[child]<-colot[cur]

     cur<-child
  }
  i<-i+1
}

return(colot)
}    


colobary.roots<-function(parent,level,colothre=min(level),paletti=NULL)
{
mt<-multitree(parent) #roots<-mt$roots child<-mt$child sibling<-mt$sibling

itemnum<-length(mt$child)
colot<-matrix("",itemnum,1)
rootnum<-length(mt$roots)
if (is.null(paletti))
 paletti<-c("red","blue","green",
 "orange","navy","darkgreen",
 "orchid","aquamarine","turquoise",
 "pink","violet","magenta","chocolate","cyan",
 colors()[50:657],colors()[50:657])
hep<-1

for (i in 1:rootnum){
   curroot<-mt$roots[i]
   colot[curroot]<-paletti[hep]
   hep<-hep+1
   if (mt$child[curroot]>0){
      pino<-matrix(0,itemnum,1)
      pino[1]<-mt$child[curroot]
      pinin<-1
      while (pinin>0){
          cur<-pino[pinin]      #take from stack
          pinin<-pinin-1
          if ((mt$sibling[cur]==0)||(level[parent[cur]]<colothre)) 
              colot[cur]<-colot[parent[cur]]
          else{ 
                colot[cur]<-paletti[hep]
                hep<-hep+1
          } 
          # put to the stack 
          if (mt$sibling[cur]>0){
                 pinin<-pinin+1
                 pino[pinin]<-mt$sibling[cur]
          }
          # go to left and put right nodes to the stack
          while (mt$child[cur]>0){   
             cur<-mt$child[cur]
             if ((mt$sibling[cur]==0)||(level[parent[cur]]<colothre)) 
                       colot[cur]<-colot[parent[cur]]
             else{ 
                    colot[cur]<-paletti[hep]
                    hep<-hep+1
             } 
             if (mt$sibling[cur]>0){
                 pinin<-pinin+1
                 pino[pinin]<-mt$sibling[cur]
             }
           }
       }#while (pinin>0)
   }
}                           
return(colot)
}




 colors2data<-function(dendat,pcf,lst,paletti=NULL,clusterlevel=NULL,nodes=NULL,
type="regular")
{
if (type=="regular")
return( colorsofdata(dendat,pcf,lst,paletti=paletti,clusterlevel=clusterlevel,nodes=nodes) )
else
return( colorsofdata.adagrid(dendat,pcf,lst,paletti=paletti,clusterlevel=clusterlevel,nodes=nodes) )

}

colorsofdata2<-function(dendat,pcf,lst,paletti=NULL,
clusterlevel=NULL,nodes=NULL)
{
# links from dendat to rec to node to color
# "lst$infopointer" gives links from nodes to recs

n<-dim(dendat)[1]
d<-dim(dendat)[2]
rnum<-length(pcf$value)

step<-matrix(0,d,1)
for (i in 1:d) step[i]<-(pcf$support[2*i]-pcf$support[2*i-1])/pcf$N[i]
    
if (is.null(paletti))
paletti<-c("red","blue","green",
    "orange","navy","darkgreen",
    "orchid","aquamarine","turquoise",
    "pink","violet","magenta","chocolate","cyan",
    colors()[50:657],colors()[50:657])

# links from node to color
if ((is.null(clusterlevel))&&(is.null(nodes))) col<-colobary(lst$parent,paletti)
if (!is.null(clusterlevel))
col<-colobary.merge(lst$parent,lst$level,colothre=clusterlevel,paletti)
if (!is.null(nodes))
col<-colobary.nodes(lst$parent,nodes,paletti)

# links from rec to node (invert the links in infopointer)
nodefinder<-matrix(0,rnum,1)
for (i in 1:rnum) nodefinder[lst$infopointer[i]]<-i

# find links from dendat to rec
den2pcf<-matrix(0,n,1)
pcf2den<-matrix(0,rnum,1)
value<-matrix(0,n,1)
for (i in 1:n){
    j<-1
    while (j<=rnum){
         inside<-TRUE
         coordi<-1
         while ((inside) && (coordi<=d)){
             ala<-pcf$down[j,coordi]
             yla<-pcf$high[j,coordi]
             ala<-pcf$support[2*coordi-1]+ala*step[coordi]
             yla<-pcf$support[2*coordi-1]+yla*step[coordi]
             if ((dendat[i,coordi]<ala) || (dendat[i,coordi]>yla)) 
                         inside<-FALSE
             coordi<-coordi+1
         }
         if (inside){
            den2pcf[i]<-j
            pcf2den[j]<-i
            value[i]<-pcf$value[j]
         }
         j<-j+1
    }
}

datcol<-matrix("white",n,1)
for (i in 1:n){
    eka<-den2pcf[i]
    if (eka>0) tok<-nodefinder[eka]
    if ((eka>0)&&(tok>0)) datcol[i]<-col[tok]
}

or<-order(value,decreasing=FALSE)
return(list(datacolo=datcol,ord=or))
}



colorsofdata3<-function(dendat,pcf,lst,paletti=NULL, clusterlevel=NULL,nodes=NULL){
# this version written made by Sauli Herrala
# links from dendat to rec to node to color
# "lst$infopointer" gives links from nodes to recs

  n<-dim(dendat)[1]
  d<-dim(dendat)[2]
  rnum<-length(pcf$value)
  
  i <- 1:d
  step <-(pcf$support[2*i]-pcf$support[2*i-1])/pcf$N[i]
	  
  if (is.null(paletti))
  paletti<-c("red","blue","green",
      "orange","navy","darkgreen",
      "orchid","aquamarine","turquoise",
      "pink","violet","magenta","chocolate","cyan",
      colors()[50:657],colors()[50:657])
  
  # links from node to color
  if ((is.null(clusterlevel))&&(is.null(nodes))) col<-colobary(lst$parent,paletti)
  if (!is.null(clusterlevel)) col<-colobary.merge(lst$parent,lst$level,colothre=clusterlevel,paletti)
  if (!is.null(nodes)) col<-colobary.nodes(lst$parent,nodes,paletti)
  # links from rec to node (invert the links in infopointer)
 
  nodefinder<-matrix(0,rnum,1)
  for (i in 1:rnum) nodefinder[lst$infopointer[i]]<-i
  
  den2pcf<-matrix(0,n,1)
  pcf2den<-matrix(0,rnum,1)
  value<-matrix(0,n,1)
  ala <- pcf$down
  yla <- pcf$high
  alaTesti <- t(ala) * step + pcf$support[2*1:ncol(ala) -1]
  ylaTesti <- t(yla) * step + pcf$support[2*1:ncol(ala) -1]
  
  for (i in 1:n){
    bol <- (c(dendat[i, ]) < alaTesti) | (unlist(dendat[i,]) >  ylaTesti)
	j <- which.min(colSums(bol))
	den2pcf[i] <- j
    pcf2den[j] <- i
    value[i] <- pcf$value[j]	
  } 	
  
  datcol<-matrix("white",n,1)
  tok <- 0
  for (i in 1:n){
      eka<-den2pcf[i]
      if (eka>0) tok<-nodefinder[eka]
      if (tok>0) datcol[i]<-col[tok]
  }
  
  or<-order(value,decreasing=FALSE)
  return(list(datacolo=datcol,ord=or))
}


colorsofdata.adagrid2<-function(dendat,pcf,lst,paletti=NULL,
clusterlevel=NULL,nodes=NULL)
{
# links from dendat to rec to node to color
# "lst$infopointer" gives links from nodes to recs

n<-dim(dendat)[1]
d<-dim(dendat)[2]
rnum<-length(pcf$value)

if (is.null(paletti))
paletti<-c("red","blue","green",
    "orange","navy","darkgreen",
    "orchid","aquamarine","turquoise",
    "pink","violet","magenta","chocolate","cyan",
    colors()[50:657],colors()[50:657])

# links from node to color
if ((is.null(clusterlevel))&&(is.null(nodes))) col<-colobary(lst$parent,paletti)
if (!is.null(clusterlevel))
col<-colobary.merge(lst$parent,lst$level,colothre=clusterlevel,paletti)
if (!is.null(nodes))
col<-colobary.nodes(lst$parent,nodes,paletti)

# links from rec to node (invert the links in infopointer)
nodefinder<-matrix(0,rnum,1)
for (i in 1:rnum) nodefinder[lst$infopointer[i]]<-i

# find links from dendat to rec
den2pcf<-matrix(0,n,1)
pcf2den<-matrix(0,rnum,1)
value<-matrix(0,n,1)
for (i in 1:n){
    j<-1
    while (j<=rnum){
         inside<-TRUE
         coordi<-1
         while ((inside) && (coordi<=d)){
             ala<-pcf$down[j,coordi]
             yla<-pcf$high[j,coordi]
             ala<-pcf$grid[ala,coordi]   #pcf$support[2*coordi-1]+ala*step[coordi]
             yla<-pcf$grid[yla,coordi]   #pcf$support[2*coordi-1]+yla*step[coordi]
             if ((dendat[i,coordi]<ala) || (dendat[i,coordi]>yla)) 
                         inside<-FALSE
             coordi<-coordi+1
         }
         if (inside){
            den2pcf[i]<-j
            pcf2den[j]<-i
            value[i]<-pcf$value[j]
         }
         j<-j+1
    }
}

datcol<-matrix("white",n,1)
for (i in 1:n){
    eka<-den2pcf[i]
    if (eka>0) tok<-nodefinder[eka]
    if ((eka>0)&&(tok>0)) datcol[i]<-col[tok]
}

or<-order(value,decreasing=FALSE)
return(list(datacolo=datcol,ord=or))
}



colorsofdata.adagrid.new<-function(dendat, pcf, lst, paletti = NULL, clusterlevel=NULL, 
nodes=NULL){
  # links from dendat to rec to node to color
  # "lst$infopointer" gives links from nodes to recs

  n<-dim(dendat)[1]
  d<-dim(dendat)[2]
  rnum<-length(pcf$value)
  
  	  
  if (is.null(paletti)){
    paletti<-c("red","blue","green","orange","navy","darkgreen",
      "orchid","aquamarine","turquoise", "pink","violet","magenta",
      "chocolate","cyan", colors()[50:657],colors()[50:657])
  }
  # links from node to color
  if ((is.null(clusterlevel))&&(is.null(nodes))) col<-colobary(lst$parent,paletti)
  if (!is.null(clusterlevel)) col<-colobary.merge(lst$parent,lst$level,colothre=clusterlevel,paletti)
  if (!is.null(nodes)) col<-colobary.nodes(lst$parent,nodes,paletti)

  # links from rec to node (invert the links in infopointer)
  nodefinder<-matrix(0,rnum,1)
  for (i in 1:rnum) nodefinder[lst$infopointer[i]]<-i

  # find links from dendat to rec
  den2pcf<-matrix(0, n, 1)
  pcf2den<-matrix(0, rnum, 1)
  value<-matrix(0, n, 1)
  ala <- pcf$down
  yla <- pcf$high
  # a bit complex
  ala <- sapply(1:ncol(ala), function(x, pcf, ala) pcf$grid[ala[, x], x], pcf = pcf, ala = ala)
  yla <- sapply(1:ncol(yla), function(x, pcf, yla) pcf$grid[yla[, x], x], pcf = pcf, yla = yla)
  ala <- t(ala)
  yla <- t(yla)
  
  prc <- proc.time()
  for (i in 1:n){
    bol <- (c(dendat[i, ]) < ala) | (c(dendat[i,]) >  yla)
	j <- which.min(colSums(bol))
	den2pcf[i] <- j
    pcf2den[j] <- i
    value[i] <- pcf$value[j]	
  } 	
  datcol <- matrix("white",n,1)
  for (i in 1:n){
    eka<-den2pcf[i]
    if (eka>0) tok<-nodefinder[eka]
    if ((eka>0)&&(tok>0)) datcol[i]<-col[tok]
  } 
  or<-order(value,decreasing=FALSE)
  return(list(datacolo=datcol,ord=or))
}


colorsofdata.adagrid <- function(dendat, pcf, lst, paletti = NULL, clusterlevel=NULL, nodes=NULL){
  # links from dendat to rec to node to color
  # "lst$infopointer" gives links from nodes to recs
# this version written made by Sauli Herrala

  n<-dim(dendat)[1]
  d<-dim(dendat)[2]
  rnum<-length(pcf$value)
  
  	  
  if (is.null(paletti)){
    paletti<-c("red","blue","green","orange","navy","darkgreen",
      "orchid","aquamarine","turquoise", "pink","violet","magenta",
      "chocolate","cyan", colors()[50:657],colors()[50:657])
  }
  # links from node to color
  if ((is.null(clusterlevel))&&(is.null(nodes))) col<-colobary(lst$parent,paletti)
  if (!is.null(clusterlevel)) col<-colobary.merge(lst$parent,lst$level,colothre=clusterlevel,paletti)
  if (!is.null(nodes)) col<-colobary.nodes(lst$parent,nodes,paletti)

  # links from rec to node (invert the links in infopointer)
  nodefinder<-matrix(0,rnum,1)
  for (i in 1:rnum) nodefinder[lst$infopointer[i]]<-i

  # find links from dendat to rec
  den2pcf<-matrix(0, n, 1)
  pcf2den<-matrix(0, rnum, 1)
  value<-matrix(0, n, 1)
  ala <- pcf$down
  yla <- pcf$high
  # a bit complex
  ala <- sapply(1:ncol(ala), function(x, pcf, ala) pcf$grid[ala[, x], x], pcf = pcf, ala = ala)
  yla <- sapply(1:ncol(yla), function(x, pcf, yla) pcf$grid[yla[, x], x], pcf = pcf, yla = yla)
  ala <- t(ala)
  yla <- t(yla)
  
  prc <- proc.time()
  for (i in 1:n){
    bol <- (c(dendat[i, ]) < ala) | (c(dendat[i,]) >  yla)
	j <- which.min(colSums(bol))
	den2pcf[i] <- j
    pcf2den[j] <- i
    value[i] <- pcf$value[j]	
  } 	
  datcol <- matrix("white",n,1)
  for (i in 1:n){
    eka<-den2pcf[i]
    if (eka>0) tok<-nodefinder[eka]
    if ((eka>0)&&(tok>0)) datcol[i]<-col[tok]
  } 
  or<-order(value,decreasing=FALSE)
  return(list(datacolo=datcol,ord=or))
}

colorsofdata.new<-function(dendat, pcf, lst, paletti=NULL, clusterlevel=NULL, nodes=NULL){
# links from dendat to rec to node to color
# "lst$infopointer" gives links from nodes to recs
  n<-dim(dendat)[1]
  d<-dim(dendat)[2]
  rnum<-length(pcf$value)
  
  i <- 1:d
  step <-(pcf$support[2*i]-pcf$support[2*i-1])/pcf$N[i]
	  
  if (is.null(paletti))
  paletti<-c("red","blue","green",
      "orange","navy","darkgreen",
      "orchid","aquamarine","turquoise",
      "pink","violet","magenta","chocolate","cyan",
      colors()[50:657],colors()[50:657])
  
  # links from node to color
  if ((is.null(clusterlevel))&&(is.null(nodes))) col<-colobary(lst$parent,paletti)
  if (!is.null(clusterlevel)) col<-colobary.merge(lst$parent,lst$level,colothre=clusterlevel,paletti)
  if (!is.null(nodes)) col<-colobary.nodes(lst$parent,nodes,paletti)
  # links from rec to node (invert the links in infopointer)
 
  nodefinder<-matrix(0,rnum,1)
  for (i in 1:rnum) nodefinder[lst$infopointer[i]]<-i
  
  den2pcf<-matrix(0,n,1)
  pcf2den<-matrix(0,rnum,1)
  value<-matrix(0,n,1)
  ala <- pcf$down
  yla <- pcf$high
  alaTesti <- t(ala) * step + pcf$support[2*1:ncol(ala) -1]
  ylaTesti <- t(yla) * step + pcf$support[2*1:ncol(ala) -1]
  
  for (i in 1:n){
    bol <- (c(dendat[i, ]) < alaTesti) | (c(dendat[i,]) >  ylaTesti)
	j <- which.min(colSums(bol))
	den2pcf[i] <- j
    pcf2den[j] <- i
    value[i] <- pcf$value[j]	
  } 	
  
  datcol<-matrix("white",n,1)
  tok <- 0
  for (i in 1:n){
      eka<-den2pcf[i]
      if (eka>0) tok<-nodefinder[eka]
      if (tok>0) datcol[i]<-col[tok]
  }
  
  or<-order(value,decreasing=FALSE)
  return(list(datacolo=datcol,ord=or))
}





colorsofdata<-function(dendat, pcf, lst, paletti=NULL, clusterlevel=NULL, nodes=NULL)
{
# links from dendat to rec to node to color
# "lst$infopointer" gives links from nodes to recs
# this version written made by Sauli Herrala

  n<-dim(dendat)[1]
  d<-dim(dendat)[2]
  rnum<-length(pcf$value)
  
  i <- 1:d
  step <-(pcf$support[2*i]-pcf$support[2*i-1])/pcf$N[i]
	  
  if (is.null(paletti))
  paletti<-c("red","blue","green",
      "orange","navy","darkgreen",
      "orchid","aquamarine","turquoise",
      "pink","violet","magenta","chocolate","cyan",
      colors()[50:657],colors()[50:657])
  
  # links from node to color
  if ((is.null(clusterlevel))&&(is.null(nodes))) col<-colobary(lst$parent,paletti)
  if (!is.null(clusterlevel)) col<-colobary.merge(lst$parent,lst$level,colothre=clusterlevel,paletti)
  if (!is.null(nodes)) col<-colobary.nodes(lst$parent,nodes,paletti)
  # links from rec to node (invert the links in infopointer)
 
  nodefinder<-matrix(0,rnum,1)
  for (i in 1:rnum) nodefinder[lst$infopointer[i]]<-i
  
  den2pcf<-matrix(0,n,1)
  pcf2den<-matrix(0,rnum,1)
  value<-matrix(0,n,1)
  ala <- pcf$down
  yla <- pcf$high
  alaTesti <- t(ala) * step + pcf$support[2*1:ncol(ala) -1]
  ylaTesti <- t(yla) * step + pcf$support[2*1:ncol(ala) -1]
  
  for (i in 1:n){
    bol <- (c(dendat[i, ]) < alaTesti) | (c(dendat[i,]) >  ylaTesti)
	j <- which.min(colSums(bol))
	den2pcf[i] <- j
    pcf2den[j] <- i
    value[i] <- pcf$value[j]	
  } 	
  
  datcol<-matrix("white",n,1)
  tok <- 0
  for (i in 1:n){
      eka<-den2pcf[i]
      if (eka>0) tok<-nodefinder[eka]
      if (tok>0) datcol[i]<-col[tok]
  }
  
  or<-order(value,decreasing=FALSE)
  return(list(datacolo=datcol,ord=or))
}



colorsofdata.tail<-function(dendat,lst,paletti=NULL)
{
# links from dendat to node to color
# "lst$infopointer" gives links from nodes to data

n<-dim(dendat)[1]
d<-dim(dendat)[2]
    
if (is.null(paletti))
paletti<-c("red","blue","green",
    "orange","navy","darkgreen",
    "orchid","aquamarine","turquoise",
    "pink","violet","magenta","chocolate","cyan",
    colors()[50:657],colors()[50:657])

# links from node to color
col<-colobary(lst$parent,paletti)

# links from dendat to node (invert the links in infopointer)
nodefinder<-matrix(0,n,1)
for (i in 1:n) nodefinder[lst$infopointer[i]]<-i

datcol<-matrix("white",n,1)
for (i in 1:n){
    tok<-nodefinder[i]
    datcol[i]<-col[tok]
}

return(datacolo=datcol)
}



complex.rips<-function(dendat,rho)
{
n<-dim(dendat)[1]
d<-dim(dendat)[2]

lkm<-0
complex<-matrix(0,1,d+1)

nn<-nn.indit(dendat)
maxk<-n-1
nnr<-nn.radit(dendat,maxk)

for (i in 1:n){
   rcur<-nnr[i,1]
   j<-1
   while ( (rcur<=rho) && (j<=(n-1)) ){
        cind<-nn[i,j]

        if (cind>i){
        # find the connection
        rcur1<-nnr[i,j]
        rcur2<-nnr[cind,1]
        kk<-j+1
        found<-FALSE
        while ( (rcur1<=rho) && (kk<=(n-1)) && (!found) ){
           koe1<-nn[i,kk]
           ll<-1
           while ( (rcur2<=rho) && (ll<=(n-1)) && (!found) ){
                koe2<-nn[cind,ll]
                if (koe1==koe2){
                     found<-TRUE
                     addi<-matrix(c(i,cind,koe1),1,3)
                     if (lkm==0) complex<-matrix(c(i,cind,koe1),1,3)
                     else complex<-rbind(complex,addi)
                     lkm<-lkm+1 
                } 
                ll<-ll+1 
                if (ll<=(n-1)) rcur2<-nnr[cind,ll]     
           }
           kk<-kk+1
           if (kk<=(n-1)) rcur1<-nnr[i,kk]
        }        
        # connection search end
        }

        j<-j+1
        if (j<=(n-1)) rcur<-nnr[i,j]
   }
}

if (lkm==0) complex<-NULL
return(complex)
}

cumu<-function(values,recs,frekv=NULL){
#Finds level sets of a piecewise constant function 
#
#values is recnum-vector
#recs is recnum*(2*d)-matrix
#frekv is recnum-vector
#
#returns list(levels,lsets,recs)
#   levels is levnum-vector,
#   lsets is levnum*atomnum-matrix,
#   atoms is recs but rows in different order
#   frekv is also only ordered differently

jarj<-omaord(values,recs,frekv)
values<-jarj$values
recs<-jarj$recs
frekv<-jarj$frekv
recnum<-length(values) #=length(recs[,1])#numb of lev is in the worst case
levels<-matrix(0,recnum,1)      
lsets<-matrix(1,recnum,recnum)  #same as the number of recs
       #at the beginning we mark everything belonging to level sets
       #next we start removing recs from level sets
levels[1]<-values[1]    #smallest values are first, first row of levels
                       #contains already 1:s 
curval<-values[1]
curlev<-1
for (i in 1:recnum){
  if (values[i]<=curval) lsets[(curlev+1):recnum,i]<-0
    else{
      curlev<-curlev+1      
      curval<-values[i]
      levels[curlev]<-values[i]
      if ((curlev+1)<=recnum) lsets[(curlev+1):recnum,i]<-0
    }
}
levels<-levels[1:curlev]
lsets<-lsets[1:curlev,]
return(list(levels=levels,lsets=lsets,atoms=recs,frekv=frekv))
}
cutmut<-function(mut,level,levels){
#
roots<-mut$roots
child<-mut$child
sibling<-mut$sibling
#
itemnum<-length(child)
rootnum<-length(roots)       
#
newroots<-matrix(0,itemnum,1)
newsibling<-sibling
ind<-0
#
for (i in 1:rootnum){
      curroot<-roots[i]
      pino<-matrix(0,itemnum,1)
      pino[1]<-curroot
      pinin<-1
      while (pinin>0){
          cur<-pino[pinin]      #take from stack
          pinin<-pinin-1
          # 
          # if cur acrosses the level, make cur root
          if (levels[cur]>level){
             ind<-ind+1
             newroots[ind]<-cur     # add to list
             newsibling[cur]<-0     # remove siblings
          }               
          # put to the stack
          if (sibling[cur]>0){
                 pinin<-pinin+1
                 pino[pinin]<-sibling[cur]
          }
          # go to left and put right nodes to the stack
          # if we have not already crossed the level
          while ((child[cur]>0) && (levels[cur]<=level)){
             cur<-child[cur]
             # if cur acrosses the level, make cur root
             if (levels[cur]>level){
                ind<-ind+1
                newroots[ind]<-cur  # add to list
                newsibling[cur]<-0     # remove siblings
             }  
             if (sibling[cur]>0){
                 pinin<-pinin+1
                 pino[pinin]<-sibling[cur]
             }
           }
       }
}                                        
if (ind>0){
    newroots<-newroots[1:ind]
}
else{
    newroots<-NULL
}
return(list(roots=newroots,sibling=newsibling))
}





cutvalue<-function(roots,child,sibling,level,component,
          AtomlistAtom,AtomlistNext,valnum){
#
#from the cutted multitree, form a "newvalue",
#which gives quantized values for the kernel estimate,
#in addition the values are cutted, so that one mode is 
#removed (input is cutted multitree)
#
itemnum<-length(child)
rootnum<-length(roots)
newvalue<-matrix(0,valnum,1)

for (i in 1:rootnum){
    pino<-matrix(0,itemnum,1)
    pino[1]<-roots[i]
    #    
    pinin<-1
    while (pinin>0){
        cur<-pino[pinin]      #take from stack
        pinin<-pinin-1
        #
        node<-cur
        compo<-component[node]
        ato<-compo                          #ato is pointer to "value"
        while (ato>0){
           newvalue[AtomlistAtom[ato]]<-level[node]
           ato<-AtomlistNext[ato]
        }
        #
        if (sibling[cur]>0){
              pinin<-pinin+1
              pino[pinin]<-sibling[cur]
        }
        while (child[cur]>0){    #go to left and put right nodes to stack
              cur<-child[cur]
              #
              node<-cur
              compo<-component[node]
              ato<-compo                    #ato is pointer to "value"
              while (ato>0){
                  newvalue[AtomlistAtom[ato]]<-level[node]
                  ato<-AtomlistNext[ato]
              }
              #
              if (sibling[cur]>0){  #if cur has siblings
                 pinin<-pinin+1
                 pino[pinin]<-sibling[cur]
             }
        }
    }
}
#
return(newvalue)
}


cvolumbag<-function(component,AtomlistAtom,AtomlistNext,low,upp,steppi)
{
d<-dim(low)[2]

componum<-length(component)
volume<-matrix(0,componum,1)

for (i in 1:componum){
   curvolu<-0
   pointer<-component[i]
   while (pointer>0){
        atto<-AtomlistAtom[pointer]

        vol<-1
        for (j in 1:d){
            vol<-vol*(upp[atto,j]-low[atto,j])*steppi[j]
        }

        curvolu<-curvolu+vol
        pointer<-AtomlistNext[pointer]
   }
   volume[i]<-curvolu
}
return(volume)
}
cvolumdya<-function(volofatom,component,AtomlistNext){
#
componum<-length(component)
volume<-matrix(0,componum,1)
#
# it is enough to calculate the number of aoms in each component
for (i in 1:componum){
   numofatoms<-0
   pointer<-component[i]
   while (pointer>0){
        numofatoms<-numofatoms+1
        pointer<-AtomlistNext[pointer]
   }
   volume[i]<-numofatoms*volofatom
}
return(volume)
}
cvolum<-function(levels,items){
#Calculates volumes of set of level sets
#
#levels is tasolkm*N-matrix of 1:s and 0:s
#items is N*(2*d)-matrix
#
#returns N-vector of volumes
#
N<-length(levels[,1])
res<-matrix(0,N,1)
if (dim(t(levels))[1]==1) tasolkm<-1 else tasolkm<-length(levels[,1]) 
for (i in 1:tasolkm){
  lev2<-change(levels[i,])
  m<-length(lev2)
  vol<-0
  for (j in 1:m){
    ind<-lev2[j]
    rec<-items[ind,]
    vol<-vol+massone(rec)
  }
  res[i]<-vol
}
return(t(res))
}


 

declevdya<-function(beg,AtomlistNext,AtomlistAtom,kg,N,nodenumOfDyaker,
terminalnum,d){
#
#beg is pointer to AtomlistAtom
#nodenumOfDyaker is the num of nodes of the _original dyaker_
#terminalnum is the num of terminal nodes of the "current dyaker"
#
#kg=kernel estimate is represented as binary tree with "nodenumOfDyaker" nodes:
#   vectors whose length is "nodenum":
#     -left,right,parent
#     -infopointer: pointer to "value" and "index" (only for terminal nodes)
#   additional data structures  
#     -value, index
#
#to be created:
#-separy, vector with length "nodenumOfDyaker", points to begsSepaBegs
#-begsSepaBegs, begsSepaNext, begsLeftBoun, begsRighBoun
#   vectors of same length as "value" and "index",
#   list of separate sets, index gives starting point for set in atomsSepaNext
#-atomsSepaAtom, atomsSepaNext, atomsLBounNext, atomsRBounNext
#   vector of same length as "value" and "index",
#   list of atoms in separate sets, index gives the atom in value and index
#
#return: 
#begs: list of beginnings of lists
#AtomlistAtoms, AtomlistNext: list of lists of atoms
#
left<-kg$left
right<-kg$right
parent<-kg$parent
index<-kg$index
nodefinder<-kg$nodefinder
#infopointer<-kg$infopointer
#
separy<-matrix(0,nodenumOfDyaker,1)
#
begsSepaNext<-matrix(0,terminalnum,1)
begsSepaBegs<-matrix(0,terminalnum,1)   #pointers to begsSepaAtoms
begsLeftBoun<-matrix(0,terminalnum,1)
begsRighBoun<-matrix(0,terminalnum,1)
#
atomsSepaNext<-matrix(0,terminalnum,1)  #pointers to value,index
atomsSepaAtom<-matrix(0,terminalnum,1)
atomsLBounNext<-matrix(0,terminalnum,1)
atomsRBounNext<-matrix(0,terminalnum,1)
#
nextFloor<-matrix(0,terminalnum,1)
currFloor<-matrix(0,terminalnum,1)
already<-matrix(0,nodenumOfDyaker,1)
#
#############################
# INITIALIZING: "we go through the nodes at depth "depoftree""
# we make currFloor to be one over bottom floor and initialize 
# separy, boundary, atomsSepaAtom, atomsBounAtom
##############################
#
lkm<-0
curlkm<-0
curre<-beg
while(curre>0){
    lkm<-lkm+1
    atom<-AtomlistAtom[curre]
    node<-nodefinder[atom]
    #
    separy[node]<-lkm
    atomsSepaAtom[lkm]<-atom
    #
    exists<-parent[node]
    if (already[exists]==0){
         curlkm<-curlkm+1
         currFloor[curlkm]<-exists 
         already[exists]<-1
    }
    #
    curre<-AtomlistNext[curre]
}     # obs terminalnum=lkm
# initialize the rest
begsSepaBegs[1:terminalnum]<-seq(1,terminalnum)       
begsLeftBoun[1:terminalnum]<-seq(1,terminalnum)
begsRighBoun[1:terminalnum]<-seq(1,terminalnum)
#obs: we need not change 
#begsSepaNext, atomsSepaNext, atomsLBounNext, atomsRBounNext
#since at the beginning set consist only one member:
#pointer is always 0, since we do not have followers 
#
###########################
#
# START the MAIN LOOP
###########################
i<-d
while (i >= 2){
   j<-log(N[i],base=2)   #depth at direction d
   while (j>=1){
        nexlkm<-0
        k<-1
        while (k <= curlkm){
            node<-currFloor[k]  
            # we create simultaneously the upper floor
            exists<-parent[node]
            if (already[exists]==0){
                 nexlkm<-nexlkm+1
                 nextFloor[nexlkm]<-exists
                 already[exists]<-1
            }
####################################################
#            now we join childs            
####################################################
            leftbeg<-left[node]
            rightbeg<-right[node]
            direction<-i
            jg<-joingene(node,leftbeg,rightbeg,separy,
begsSepaNext,begsSepaBegs,begsLeftBoun,begsRighBoun,
atomsSepaNext,atomsSepaAtom,atomsLBounNext,atomsRBounNext,
direction,index)
#
separy<-jg$separy
begsSepaNext<-jg$begsSepaNext
begsSepaBegs<-jg$begsSepaBegs
begsLeftBoun<-jg$begsLeftBoun
begsRighBoun<-jg$begsRighBoun
atomsSepaNext<-jg$atomsSepaNext
atomsSepaAtom<-jg$atomsSepaAtom
atomsLBounNext<-jg$atomsLBounNext
atomsRBounNext<-jg$atomsRBounNext
           k<-k+1
        }
        j<-j-1
        curlkm<-nexlkm
        currFloor<-nextFloor
   }
   # now we move to the next direction, correct boundaries
   begsLeftBoun<-begsSepaBegs
   begsRighBoun<-begsSepaBegs
   #
   atomsLBounNext<-atomsSepaNext
   atomsRBounNext<-atomsSepaNext
   #
   i<-i-1
}
#########################
# ENO OF MAIN LOOP
#########################
#
#
###################
# LAST DIMENSION WILL BE handled, (because this contains root node
##################
   i<-1
   j<-log(N[i],base=2)   #depth at direction d
   while (j>=2){
        nexlkm<-0
        k<-1
        while (k <= curlkm){
            node<-currFloor[k]  
            # we create simultaneously the upper floor
            exists<-parent[node]
            if (already[exists]==0){
                 nexlkm<-nexlkm+1
                 nextFloor[nexlkm]<-exists
                 already[exists]<-1
            }
####################################################
#            now we join childs            
#if (right[parent[node]]==node)  #if node is right child 
####################################################
            leftbeg<-left[node]
            rightbeg<-right[node]
            direction<-1 
jg<-joingene(node,leftbeg,rightbeg,separy,
begsSepaNext,begsSepaBegs,begsLeftBoun,begsRighBoun,
atomsSepaNext,atomsSepaAtom,atomsLBounNext,atomsRBounNext,
direction,index)
#
separy<-jg$separy
begsSepaNext<-jg$begsSepaNext
begsSepaBegs<-jg$begsSepaBegs
begsLeftBoun<-jg$begsLeftBoun
begsRighBoun<-jg$begsRighBoun
#
atomsSepaNext<-jg$atomsSepaNext
atomsSepaAtom<-jg$atomsSepaAtom
atomsLBounNext<-jg$atomsLBounNext
atomsRBounNext<-jg$atomsRBounNext
           k<-k+1
        }
        j<-j-1
        curlkm<-nexlkm
        currFloor<-nextFloor
   }
#########################
# ROOT NODE, we do not anymore update boundaries
#########################
            k<-1
            node<-currFloor[k]  
  while (k <= 1){
####################################################
#            now we join childs            
#if (right[parent[node]]==node)  #if node is right child 
####################################################
            leftbeg<-left[node]
            rightbeg<-right[node]
            if ((leftbeg==0) || (separy[leftbeg]==0)){
                              #if left child does not exist    
                separy[node]<-separy[rightbeg]
            }
            else{   #eka else
                if ((rightbeg==0) || (separy[rightbeg]==0)){  
                              #right child does not exist  
                   separy[node]<-separy[leftbeg]
                } 
                else{   #toka else: both children exist
                    #check whether left boundary of right child is empty
                    Lempty<-TRUE
                    note<-separy[rightbeg]
                    while (note>0){
                        if (begsLeftBoun[note]>0){
                              Lempty<-FALSE
                        }
                        note<-begsSepaNext[note]
                     }
                     #check whether right bound of left child is empty     
                     Rempty<-TRUE
                     note<-separy[leftbeg]
                     while (note>0){
                          if (begsRighBoun[note]>0){
                                 Rempty<-FALSE
                          }
                          note<-begsSepaNext[note]
                     }
                     #check whether one of boundaries is empty
                     if (Lempty || Rempty){
                            #one of boundaries is empty         
############
#concatenating separate parts
#separy[node]<- concatenate separy[leftbeg],separy[rightbeg]
###########
akku<-separy[leftbeg]
while (begsSepaNext[akku]>0){
  akku<-begsSepaNext[akku]
}                           
begsSepaNext[akku]<-separy[rightbeg]
#
separy[node]<-separy[leftbeg]
####################
#end of concatenating, handle next boundaries
###################
                    }
                    else{  #both children exist, both boundaries non-empty  
direction<-i
jc<-joinconne(leftbeg,rightbeg,separy,
begsSepaNext,begsSepaBegs,begsLeftBoun,begsRighBoun,
atomsSepaNext,atomsSepaAtom,atomsLBounNext,atomsRBounNext,
direction,index)  #direction<-i
#
separy<-jc$separy
separy[node]<-jc$totbegSepary 
#
begsSepaNext<-jc$begsSepaNext
begsSepaBegs<-jc$begsSepaBegs
begsLeftBoun<-jc$begsLeftBoun
begsRighBoun<-jc$begsRighBoun
#
atomsSepaNext<-jc$atomsSepaNext
atomsSepaAtom<-jc$atomsSepaAtom
atomsLBounNext<-jc$atomsLBounNext
atomsRBounNext<-jc$atomsRBounNext
                        #
                    }
                } #toka else
            } #eka else
 k<-k+1
}
###########################################
#          end of child joining
###########################################
###################################
# END of ROOT
###################################
#
###########################
###########################
totbegSepary<-separy[node]
return(list(totbegSepary=totbegSepary,
begsSepaNext=begsSepaNext,begsSepaBegs=begsSepaBegs,
atomsSepaNext=atomsSepaNext,atomsSepaAtom=atomsSepaAtom))
}













declevgen<-function(tobehandled,curterminalnum,
left,right,val,vec,infopointer,parent,
low,upp,
d)
# source("~/denpro/R/declevgen.R")
# dl<-declevgen(tobehandled,curterminalnum,left,right,val,vec,
# infopointer,parent,low,upp,d)

{
nodelkm<-length(left)

separy<-matrix(0,nodelkm,1)          #=infopointer?

begsSepaNext<-matrix(0,curterminalnum,1)
begsSepaBegs<-matrix(0,curterminalnum,1)    #pointers to begsSepaAtom
begsLeftBoun<-matrix(0,curterminalnum,1)
begsRighBoun<-matrix(0,curterminalnum,1)

atomsSepaNext<-matrix(0,curterminalnum,1)   #pointers to value,index
atomsSepaAtom<-matrix(0,curterminalnum,1)
atomsLBounNext<-matrix(0,curterminalnum,1)
atomsRBounNext<-matrix(0,curterminalnum,1)

leafloc<-findleafs(left,right)

lkm<-0
node<-nodelkm
while (node>=1){ 
          # root is in position 1
    if ((leafloc[node]==1) && (tobehandled[node]==1)){   #we are in leaf 

          tobehandled[parent[node]]<-1

          lkm<-lkm+1
          separy[node]<-lkm
          atomsSepaAtom[lkm]<-infopointer[node]

          begsSepaBegs[lkm]<-lkm

          #begsLeftBoun[lkm]<-lkm
          #begsRighBoun[lkm]<-lkm

          #obs: we need not change
          #begsSepaNext, atomsSepaNext, atomsLBounNext, atomsRBounNext
          #since at the beginning set consist only one member:
          #pointer is always 0, since we do not have followers

    }
    else if (tobehandled[node]==1){   #not a leaf

        tobehandled[parent[node]]<-1

        leftbeg<-left[node]
        rightbeg<-right[node]

        if ((leftbeg==0) || (separy[leftbeg]==0)){
           #if left child does not exist

           #note that since we consider subsets of the
           #terminal nodes of the original tree, it may happen
           #that leftbeg>0 but left child does not exist
          
           separy[node]<-separy[rightbeg]

           #we need that all lists contain as many members
           #left boundary is empty, but we will make it a list
           #of empty lists

           #note<-separy[node]
           #while (note>0){
           #    begsLeftBoun[note]<-0
           #    note<-begsSepaNext[note]
           #}
           # right boundary stays same as for rightbeg
        }
        else{   # eka else
            if ((rightbeg==0) || (separy[rightbeg]==0)){
            #right child does not exist
    
                   separy[node]<-separy[leftbeg]
    
                   # left boundary stays same as for leftbeg
                   # right boundary is empty
    
                   #note<-separy[node]
                   #while (note>0){
                   #    begsRighBoun[note]<-0
                   #    note<-begsSepaNext[note]
                   #}
            }
            else{   #toka else: both children exist
                    #create boundaries 
                    
                    direktiooni<-vec[node]
                    splittiini<-val[node]     
               
                    #left boundary of right child :
                    #create/check whether empty 
                   
                    Lempty<-TRUE
                    note<-separy[rightbeg]
                    while (note>0){
                        thisnoteempty<-TRUE

                        poiju<-begsSepaBegs[note]
                        prevpoiju<-poiju 
                        while (poiju>0){
                           aatto<-atomsSepaAtom[poiju]
                           if (!(splittiini<low[aatto,direktiooni])){
                               #this atom belongs to boundary
                               if (thisnoteempty==TRUE){
                                   #poiju is the 1st non-empty
                                   begsLeftBoun[note]<-poiju
                               }                
                               Lempty<-FALSE   
                               atomsLBounNext[prevpoiju]<-poiju    
                               atomsLBounNext[poiju]<-0
                               prevpoiju<-poiju

                               thisnoteempty<-FALSE
                           }
                           poiju<-atomsSepaNext[poiju]
                        }
                        if (thisnoteempty) begsLeftBoun[note]<-0
                        note<-begsSepaNext[note]
                     }

                     #right boundary of left child

                     Rempty<-TRUE
                     note<-separy[leftbeg]
                     while (note>0){
                        thisnoteempty<-TRUE
                       
                        poiju<-begsSepaBegs[note]
                        prevpoiju<-poiju 
                        while (poiju>0){
                           aatto<-atomsSepaAtom[poiju]
                           if (!(splittiini>upp[aatto,direktiooni])){
                               #this atom belongs to boundary
                               if (thisnoteempty==TRUE){
                                   #poiju is the 1st non-empty
                                   begsRighBoun[note]<-poiju
                               }                
                               Rempty<-FALSE   
                               atomsRBounNext[prevpoiju]<-poiju    
                               atomsRBounNext[poiju]<-0
                               prevpoiju<-poiju 
                           
                               thisnoteempty<-FALSE
                           }
                           poiju<-atomsSepaNext[poiju]
                        }
                        if (thisnoteempty) begsRighBoun[note]<-0
                        note<-begsSepaNext[note]
                     }

                     #check whether one of boundaries is empty

                     if (Lempty || Rempty){
                        #one of boundaries is empty
                        ############
                        #concatenating separate parts
                        ############
                        akku<-separy[leftbeg]
                        #begsRighBoun[akku]<-0 
                        # right boundaries of sets in left child are empty
                        # begsLeftBoun[akku] does not change
                        
                        #while (begsSepaNext[akku]>0){
                        #    akku<-begsSepaNext[akku]
                        #    begsRighBoun[akku]<-0
                        #}
                        begsSepaNext[akku]<-separy[rightbeg] 
                        #concatenate list of separate sets

                        separy[node]<-separy[leftbeg]
                        akku<-separy[rightbeg]
                          
                        #begsLeftBoun[akku]<-0 
                        #left boundaries of sets in right child are empty

                        #while (begsSepaNext[akku]>0){
                        #   akku<-begsSepaNext[akku]
                        #   begsLeftBoun[akku]<-0
                        #}
                        ############# 
                        #end of concatenating
                        #############
                     }
                     else{  #both children exist, both boundaries non-empty
                     direction<-vec[node]
                     jc<-joincongen(leftbeg,rightbeg,separy,
                     begsSepaNext,begsSepaBegs,begsLeftBoun,begsRighBoun,
                     atomsSepaNext,atomsSepaAtom,atomsLBounNext,atomsRBounNext,
                     direction,low,upp)   
                      
                     separy<-jc$separy
                     separy[node]<-jc$totbegSepary

                     begsSepaNext<-jc$begsSepaNext
                     begsSepaBegs<-jc$begsSepaBegs
                     #begsLeftBoun<-jc$begsLeftBoun
                     #begsRighBoun<-jc$begsRighBoun

                     atomsSepaNext<-jc$atomsSepaNext
                     atomsSepaAtom<-jc$atomsSepaAtom
                     #atomsLBounNext<-jc$atomsLBounNext
                     #atomsRBounNext<-jc$atomsRBounNext
                     }
            } #toka else
        } #eka else
    }  #else not a leaf
    node<-node-1
}

return(list(
totbegSepary=separy[1],
begsSepaNext=begsSepaNext,
begsSepaBegs=begsSepaBegs,
atomsSepaNext=atomsSepaNext,
atomsSepaAtom=atomsSepaAtom
))
}













declevnew<-function(rindeksit,linkit,n){
#Splits level set to disconnected subsets
#
#rindeksit is m-vector, links to observations, levset is union of m atoms
#linkit is n*n-matrix, n is the total number of atoms,
#  describes which atoms touch each other
#
#Returns sublevnum*n-matrix, describes disconnected parts of levset  
#
m<-length(rindeksit)
tulos<-matrix(0,m,n)      #in the worst case we split to m parts
#                         #blokit !!!!!!!!!!!!!!!!
merkatut<-matrix(0,m,1) #laitetaan 1 jos on jo sijoitettu johonkin tasojouk.
pino<-matrix(0,m,1) #pinoon laitetaan aina jos koskettaa, max kosketuksia m
pinind<-0           #pinossa viitataan rindeksin elementteihin
curleima<-1
i<-1   #i ja j viittavat rindeksit-vektoriin, jonka alkiot viittavat atomeihin
while (i<=m){
  if (merkatut[i]==0){  #jos ei viela merkattu niin 
    pinind<-pinind+1    #pannaan pinoon
    pino[pinind]<-i    
    while (pinind>0){
      curviite<-pino[pinind]  #otetaan pinosta viite rindeksit-vektoriin
                              #jossa puolestaan viitteet itse palloihin
      curpallo<-rindeksit[curviite]  #haetaan viite palloon
      pinind<-pinind-1  
      tulos[curleima,curpallo]<-1    #laitetaan pallo ko tasojoukkoon
#      merkatut[curviite]<-1          #merkataan kaytetyksi
      j<-1
      while (j<=m){          #pannnaan linkeista pinoon
        ehdokas<-rindeksit[j]         #kaydaan ko tasojoukon atomit lapi
        touch<-(linkit[curpallo,ehdokas]==1)
        if ((touch) && (merkatut[j]==0)){
          pinind<-pinind+1      
          pino[pinind]<-j  
          merkatut[j]<-1
        }
        j<-j+1
      }
    }
    curleima<-curleima+1   #uusi leima
  }
  i<-i+1
} 
tullos<-matrix(0,(curleima-1),n)
tullos[1:(curleima-1),]<-tulos[1:(curleima-1),] #poistetaan ylimaaraiset
return(tullos)
}


















declev<-function(rindeksit,linkit,n){
#Splits level set to disconnected subsets
#
#rindeksit is m-vector, links to observations, levset is union of m atoms
#linkit is n*n-matrix, n is the total number of atoms,
#  describes which atoms touch each other
#
#Returns sublevnum*n-matrix, describes disconnected parts of levset  
#
m<-length(rindeksit)
tulos<-matrix(0,m,n)      #in the worst case we split to m parts
#                         #blokit !!!!!!!!!!!!!!!!
merkatut<-matrix(0,m,1) #laitetaan 1 jos on jo sijoitettu johonkin tasojouk.
pino<-matrix(0,m,1) #pinoon laitetaan aina jos koskettaa, max kosketuksia m
pinind<-0           #pinossa viitataan rindeksin elementteihin
curleima<-1
i<-1   #i ja j viittavat rindeksit-vektoriin, jonka alkiot viittavat atomeihin
while (i<=m){
  if (merkatut[i]==0){  #jos ei viela merkattu niin 
    pinind<-pinind+1    #pannaan pinoon
    pino[pinind]<-i    
    while (pinind>0){
      curviite<-pino[pinind]  #otetaan pinosta viite rindeksit-vektoriin
                              #jossa puolestaan viitteet itse palloihin
      curpallo<-rindeksit[curviite]  #haetaan viite palloon
      pinind<-pinind-1  
      tulos[curleima,curpallo]<-1    #laitetaan pallo ko tasojoukkoon
#      merkatut[curviite]<-1          #merkataan kaytetyksi
      j<-1
      while (j<=m){          #pannnaan linkeista pinoon
        curkoske<-linkit[curpallo,]   #ne joihin curpallo koskettaa
        ehdokas<-rindeksit[j]         #kaydaan ko tasojoukon atomit lapi
        touch<-onko(curkoske,ehdokas) #onko ehdokas rivilla "curkoske"
            #touch<-(linkit[curpallo,rindeksit[j]]==1)
        if ((touch) && (merkatut[j]==0)){
          pinind<-pinind+1      
          pino[pinind]<-j  
          merkatut[j]<-1
        }
        j<-j+1
      }
    }
    curleima<-curleima+1   #uusi leima
  }
  i<-i+1
} 
tullos<-matrix(0,(curleima-1),n)
tullos[1:(curleima-1),]<-tulos[1:(curleima-1),] #poistetaan ylimaaraiset
return(tullos)
}


















decombag<-function(numofall,levseq,
left,right,val,vec,infopointer,parent,
nodenumOfTree,terminalnum,
value,low,upp,nodefinder,
d)
{
#Makes density tree
#returns list(parent,level,component)
#  -component is pointer to AtomlistAtom, AtomlistNext, ei tarvita

#apu2<-0

levnum<-length(levseq)

AtomlistNext<-matrix(0,numofall,1)
AtomlistAtom<-matrix(0,numofall,1) #point to value,..: values in 1,...,atomnum

componum<-numofall
parentLST<-matrix(0,componum,1)
level<-matrix(0,componum,1)
component<-matrix(0,componum,1)    

pinoComponent<-matrix(0,componum,1)   #pointer to component, level,...
pinoTaso<-matrix(0,componum,1)        #ordinal of level (pointer to levseq)

# Initilize the lists

AtomlistAtom[1:terminalnum]<-seq(1,terminalnum)
AtomlistNext[1:(terminalnum-1)]<-seq(2,terminalnum)
AtomlistNext[terminalnum]<-0
listEnd<-terminalnum
beg<-1

# Let us divide the lowest level set to disconnected parts

begi<-1
tobehandled<-matrix(0,nodenumOfTree,1)
atto<-AtomlistAtom[begi]
while (begi>0){
   if (value[atto]>0){
      node<-nodefinder[atto]
      tobehandled[node]<-1
   }    
   begi<-AtomlistNext[begi]
   atto<-AtomlistAtom[begi]
}
curterminalnum<-terminalnum

dld<-declevgen(tobehandled,curterminalnum,
left,right,val,vec,infopointer,parent,
low,upp,
d)
 
totbegSepary<-dld$totbegSepary
begsSepaNext<-dld$begsSepaNext
begsSepaBegs<-dld$begsSepaBegs
atomsSepaNext<-dld$atomsSepaNext
atomsSepaAtom<-dld$atomsSepaAtom

lc<-listchange(AtomlistAtom,AtomlistNext,totbegSepary,
begsSepaNext,begsSepaBegs,atomsSepaNext,atomsSepaAtom,
curterminalnum,beg)
begs<-lc$begs
AtomlistNext<-lc$AtomlistNext
AtomlistAtom<-lc$AtomlistAtom

koko<-length(begs)
# Talletetaan osat 
component[1:koko]<-begs
level[1:koko]<-levseq[1]       #arvo toistuu 
efek<-koko                     #kirjataan uusien osien lkm  ????? jos vain yksi
# Laitetaan kaikki osat pinoon
pinoComponent[1:koko]<-seq(1,koko)      #1,2,...,koko
pinoTaso[1:koko]<-1     #kaikki osat kuuluvat alimpaan tasojoukkoon
pinind<-koko            #indeksi pinoon

if (levnum>1){  while (pinind>=1){
  # Take from stack
  ind<-pinoComponent[pinind]      #indeksi tasoon
  levind<-pinoTaso[pinind]        #ko tason korkeus
  pinind<-pinind-1                #otettu pinosta
  partlevsetbeg<-component[ind]  
  higlev<-levseq[levind+1]
  
  # Make intersection with the curr. component and higher lev.set
  PrevlistEnd<-listEnd
  addnum<-0      #num of atoms in the intersection
  removenum<-0   #num of atoms which have to be removed to get intersec.
  runner<-partlevsetbeg
  origiListEnd<-listEnd
  while ((runner>0) && (runner<=origiListEnd)){
      atom<-AtomlistAtom[runner]
      arvo<-value[atom]
      if (arvo>=higlev){
          listEnd<-listEnd+1    
          AtomlistAtom[listEnd]<-atom
          AtomlistNext[listEnd]<-listEnd+1 
          addnum<-addnum+1     
      }
      else{           
          removenum<-removenum+1
      }                                
      runner<-AtomlistNext[runner] 
  }
  AtomlistNext[listEnd]<-0      # we have to correct the end to terminate
  if (addnum>0){
      AtomlistNext[PrevlistEnd]<-0
      beghigher<-PrevlistEnd+1 
  }
  if (removenum==0){  #jos leikkaus ei muuta, niin tasoj sailyy samana  
      level[ind]<-levseq[levind+1]  #remove lower part
          #component and parentLST stay same, it is enough to change level
      if (levind+1<levnum){ #jos ei olla korkeimmalla tasolla,laita pinoon
        pinoComponent[pinind+1]<-ind     
        pinoTaso[pinind+1]<-levind+1 #tasojouk taso on levind+1 
        pinind<-pinind+1       
      }
  }
  else if (addnum>0){     #leikkaus ei tyhja
      beg<-beghigher

      begi<-beghigher
      tobehandled<-matrix(0,nodenumOfTree,1)
      while (begi>0){
          atto<-AtomlistAtom[begi]
          node<-nodefinder[atto]
          tobehandled[node]<-1
          begi<-AtomlistNext[begi]
      }
      curterminalnum<-addnum

      dld<-declevgen(tobehandled,curterminalnum,
      left,right,val,vec,infopointer,parent,
      low,upp,
      d)
 
      totbegSepary<-dld$totbegSepary
      begsSepaNext<-dld$begsSepaNext
      begsSepaBegs<-dld$begsSepaBegs
      atomsSepaNext<-dld$atomsSepaNext
      atomsSepaAtom<-dld$atomsSepaAtom

#      apu2<-apu2+1
#if (apu2==1) an<-atomsSepaNext     

      lc<-listchange(AtomlistAtom,AtomlistNext,totbegSepary,
      begsSepaNext,begsSepaBegs,atomsSepaNext,atomsSepaAtom,
      curterminalnum,beg)
      begs<-lc$begs
      AtomlistNext<-lc$AtomlistNext
      AtomlistAtom<-lc$AtomlistAtom

      koko<-length(begs)    #jos vain yksi ?????????
      #
      # paivitetaan kumu tulokseen
      level[(efek+1):(efek+koko)]<-levseq[levind+1]   #arvo toistuu  
      parentLST[(efek+1):(efek+koko)]<-ind
      component[(efek+1):(efek+koko)]<-begs
      efek<-efek+koko
      if (levind+1<levnum){ #jos ei olla korkeimmalla tasolla,laita pinoon
        pinoComponent[(pinind+1):(pinind+koko)]<-seq(efek-koko+1,efek)
        pinoTaso[(pinind+1):(pinind+koko)]<-levind+1  #tasjouk tas on levind+1 
        pinind<-pinind+koko       
      }
  }
}} 
level<-t(level[1:efek])
parentLST<-t(parentLST[1:efek])
component<-t(component[1:efek])
#
return(list(level=level,parent=parentLST,
component=component,AtomlistAtom=AtomlistAtom,AtomlistNext=AtomlistNext))
}









decomdya<-function(numofall,atomnum,levseq,kg,N,nodenumOfDyaker){
#Makes density tree
#
#returns list(parent,level,component)
#  -component is pointer to AtomlistAtom, AtomlistNext, ei tarvita
#
d<-length(N)
levnum<-length(levseq)
#
AtomlistNext<-matrix(0,numofall,1)
AtomlistAtom<-matrix(0,numofall,1) #point to value,..: values in 1,...,atomnum
#
componum<-numofall
parent<-matrix(0,componum,1)
level<-matrix(0,componum,1)
component<-matrix(0,componum,1)    
#
pinoComponent<-matrix(0,componum,1)   #pointer to component, level,...
pinoTaso<-matrix(0,componum,1)        #ordinal of level (pointer to levseq)
#
# Initilize the lists
#
AtomlistAtom[1:atomnum]<-seq(1,atomnum)
AtomlistNext[1:(atomnum-1)]<-seq(2,atomnum)
AtomlistNext[atomnum]<-0
listEnd<-atomnum
#
# Let us divide the lowest level set to disconnected parts
#
beg<-1
terminalnum<-atomnum
dld<-declevdya(beg,AtomlistNext,AtomlistAtom,kg,N,nodenumOfDyaker,
terminalnum,d)  #terminalnum<-atomnum 
totbegSepary<-dld$totbegSepary
begsSepaNext<-dld$begsSepaNext
begsSepaBegs<-dld$begsSepaBegs
atomsSepaNext<-dld$atomsSepaNext
atomsSepaAtom<-dld$atomsSepaAtom
#
lc<-listchange(AtomlistAtom,AtomlistNext,totbegSepary,
begsSepaNext,begsSepaBegs,atomsSepaNext,atomsSepaAtom,
terminalnum,beg)
begs<-lc$begs
AtomlistNext<-lc$AtomlistNext
AtomlistAtom<-lc$AtomlistAtom
#
koko<-length(begs)
# Talletetaan osat 
component[1:koko]<-begs
level[1:koko]<-levseq[1]       #arvo toistuu 
efek<-koko                     #kirjataan uusien osien lkm  ????? jos vain yksi
# Laitetaan kaikki osat pinoon
pinoComponent[1:koko]<-seq(1,koko)      #1,2,...,koko
pinoTaso[1:koko]<-1     #kaikki osat kuuluvat alimpaan tasojoukkoon
pinind<-koko            #indeksi pinoon
# 
if (levnum>1){  while (pinind>=1){
  # Take from stack
  ind<-pinoComponent[pinind]      #indeksi tasoon
  levind<-pinoTaso[pinind]        #ko tason korkeus
  pinind<-pinind-1                #otettu pinosta
  partlevsetbeg<-component[ind]  
  higlev<-levseq[levind+1]
  #
  # Make intersection with the curr. component and higher lev.set
  PrevlistEnd<-listEnd
  addnum<-0      #num of atoms in the intersection
  removenum<-0   #num of atoms which have to be removed to get intersec.
  runner<-partlevsetbeg
  origiListEnd<-listEnd
  value<-kg$value
  while ((runner>0) && (runner<=origiListEnd)){
      atom<-AtomlistAtom[runner]
      arvo<-value[atom]
      if (arvo>=higlev){
          listEnd<-listEnd+1    
          AtomlistAtom[listEnd]<-atom
          AtomlistNext[listEnd]<-listEnd+1 
          addnum<-addnum+1     
      }
      else{           
          removenum<-removenum+1
      }                                
      runner<-AtomlistNext[runner] 
  }
  AtomlistNext[listEnd]<-0      # we have to correct the end to terminate
  if (addnum>0){
      AtomlistNext[PrevlistEnd]<-0
      beghigher<-PrevlistEnd+1 
  }
  if (removenum==0){  #jos leikkaus ei muuta, niin tasoj sailyy samana  
      level[ind]<-levseq[levind+1]  #remove lower part
          #component and parent stay same, it is enough to change level
      if (levind+1<levnum){ #jos ei olla korkeimmalla tasolla,laita pinoon
        pinoComponent[pinind+1]<-ind     
        pinoTaso[pinind+1]<-levind+1 #tasojouk taso on levind+1 
        pinind<-pinind+1       
      }
  }
  else if (addnum>0){     #leikkaus ei tyhja
      beg<-beghigher
      terminalnum<-addnum
      dld<-declevdya(beg,AtomlistNext,AtomlistAtom,kg,N,nodenumOfDyaker,
                     terminalnum,d) 
totbegSepary<-dld$totbegSepary
begsSepaNext<-dld$begsSepaNext
begsSepaBegs<-dld$begsSepaBegs
atomsSepaNext<-dld$atomsSepaNext
atomsSepaAtom<-dld$atomsSepaAtom
#
lc<-listchange(AtomlistAtom,AtomlistNext,totbegSepary,
begsSepaNext,begsSepaBegs,atomsSepaNext,atomsSepaAtom,
terminalnum,beg)
begs<-lc$begs
AtomlistNext<-lc$AtomlistNext
AtomlistAtom<-lc$AtomlistAtom
#
      koko<-length(begs)    #jos vain yksi ?????????
      #
      # paivitetaan kumu tulokseen
      level[(efek+1):(efek+koko)]<-levseq[levind+1]   #arvo toistuu  
      parent[(efek+1):(efek+koko)]<-ind
      component[(efek+1):(efek+koko)]<-begs
      efek<-efek+koko
      if (levind+1<levnum){ #jos ei olla korkeimmalla tasolla,laita pinoon
        pinoComponent[(pinind+1):(pinind+koko)]<-seq(efek-koko+1,efek)
        pinoTaso[(pinind+1):(pinind+koko)]<-levind+1  #tasjouk tas on levind+1 
        pinind<-pinind+koko       
      }
  }
}} 
level<-t(level[1:efek])
parent<-t(parent[1:efek])
component<-t(component[1:efek])
#
return(list(level=level,parent=parent,
component=component,AtomlistAtom=AtomlistAtom,AtomlistNext=AtomlistNext))
}









decom<-function(lsets,levels,links,alkublokki,blokki){
#Makes density tree,  edellyttaa jarjestyksen ???????
#
#lsets is levnum*atomnum-matrix
#levels is levnum-vector
#links is  atomnum*maxtouchnum-matrix, pointers to atoms,
#  for each atom we indicate which atoms it touches
#
#lsets matriisiin lisataan riveja, silla aikaisemmat rivit jakautuvat
#erillisiin osiinsa, lisataan levels:n vastaavat alkiot 
#
#Returns list(lsets,levels,parents)
# parents and levels are newlevnum-vectors
# lsets is newlevnum*atomnum
#
if (dim(t(lsets))[1]==1) levnum<-1 else levnum<-length(lsets[,1])
                                              #rows of lsets 
if (levnum==1) atomnum<-length(lsets) else
atomnum<-length(lsets[1,])     #maxalkio on n
newlevnum<-levnum*atomnum       #karkea arvio, jokainen tasojoukko 
#                               #voi periaatteessa jakautua n:aan osaan
newlsets<-matrix(0,alkublokki,atomnum) 
newlevels<-matrix(0,alkublokki,1)
parents<-matrix(0,alkublokki,1)
#
curblokki<-alkublokki
#
pino<-matrix(0,newlevnum,2) #1.col indeksi newlsets:iin, 2.col ind tasoon
a<-1
b<-2
#
if (levnum==1) levset<-lsets else
levset<-lsets[1,]                   #alin tasojoukko
rindeksit<-change(levset)            #change the representation
kumu<-declev(rindeksit,links,atomnum)  #jaetaan alin tasoj. osiin
if (dim(t(kumu))[1]==1) koko<-1 else koko<-length(kumu[,1]) #osien lkm
# Talletetaan osat 
newlsets[1:koko,]<-kumu   
newlevels[1:koko]<-levels[1]     #arvo toistuu 
efek<-koko                       #kirjataan uusien osien lkm
# Laitetaan kaikki osat pinoon
pino[1:koko,a]<-seq(1,koko)      #1,2,...,koko
pino[1:koko,b]<-1   #kaikki osat kuuluvat alimpaan tasojoukkoon
pinind<-koko        #indeksi pinoon
# 
if (levnum>1){  while (pinind>=1){
  ind<-pino[pinind,a]         #indeksi tasoon
  levind<-pino[pinind,b]      #ko tason korkeus
  pinind<-pinind-1            #otettu pinosta
  partlevset<-newlsets[ind,]
  higlevset<-lsets[levind+1,] #huom levind<levnum
  levset<-partlevset*higlevset
  if (sum(partlevset-levset)==0){  
          #jos leikkaus ei muuta, niin tasoj sailyy samana  
      newlevels[ind]<-levels[levind+1]  #poistetaan alempi osa      
      if (levind+1<levnum){ #jos ei olla korkeimmalla tasolla,laita pinoon
        pino[pinind+1,a]<-ind     
        pino[pinind+1,b]<-levind+1 #tasojouk taso on levind+1 
        pinind<-pinind+1       
      }
  }
  else if (sum(levset)>0){   #leikkaus ei tyhja
      rindeksit<-change(levset)
      kumu<-declev(rindeksit,links,atomnum) #jaet tasoj. osa osiin
      if (dim(t(kumu))[1]==1) koko<-1 else koko<-length(kumu[,1]) #osien lkm
      if ((efek+koko)>curblokki){   
        newlsets<-blokitus(newlsets,blokki)
        newlevels<-blokitus(newlevels,blokki)
        parents<-blokitus(parents,blokki)
        curblokki<-curblokki+blokki    
      }
      newlsets[(efek+1):(efek+koko),]<-kumu  #paivitetaan kumu tulokseen
      newlevels[(efek+1):(efek+koko),]<-levels[levind+1]   #arvo toistuu  
      parents[(efek+1):(efek+koko),]<-ind
      efek<-efek+koko
      if (levind+1<levnum){ #jos ei olla korkeimmalla tasolla,laita pinoon
        pino[(pinind+1):(pinind+koko),a]<-seq(efek-koko+1,efek)
        pino[(pinind+1):(pinind+koko),b]<-levind+1 #tasojouk taso on levind+1 
        pinind<-pinind+koko       
      }
  }
}} 
newlevels<-t(newlevels[1:efek])
newlsets<-newlsets[1:efek,]
parents<-t(parents[1:efek])
return(list(lsets=newlsets,levels=newlevels,parents=parents))
}








den2reg<-function(dendat,binlkm,kantaja)
{
#muuntaa tiheysdatan regressiodataksi

#dendat on tiheysdatan sisaltava n*xlkm matriisi
#binlkm on niitten lokeroitten maara joihin yksi muuttuja ositetaan
#otosavaruus ositetaan binlkm^p lokeroon
#kantaja on 2*xlkm vaakavektori, sisaltaa kantajan kullekin 
#muuttujalle, olet etta kaikki dendat:n hav todella sisaltyvat kantajaan.

#palauttaa: listan vast
#vast.dep on saatu*1 vektori, sisaltaa frekvenssit,
#(saatu on erillisten diskretoitujen havaintojen lkm) 
#vast.ind on saatu*xlkm matriisi, sisaltaa niitten lokeroitten 
#keskipisteet joita (diskretoidussa) datassa esiintyy vahintaan yksi, 
#vast.hila on xlkm*3 matriisi, 
#vast.hila:n i:s rivi sis. i:nnelle muuttujalle 
#pisteiden maaran
#ensimmaisen regressiodatan havaintopisteen,
#viimeisen regressiodatan havaintopisteen.

xlkm<-length(dendat[1,]) #dendat:in sarakkeitten lkm on muuttujien lkm
n<-length(dendat[,1])    #dendat:in rivien lkm on havaintojen maara
hila<-matrix(1,xlkm,3)   #hila on xlkm*3 matriisi
binhila<-hila            #binhila on xlkm*3 matriisi
valpit<-matrix(1,xlkm,1)      #hilan valien pituudet
hila[,1]<-binlkm*hila[,1]
binhila[,1]<-binlkm*binhila[,1]
i<-1
while (i<=xlkm){
     binhila[i,2]<-kantaja[i,1]   #min(dendat[,i])-epsi 
     binhila[i,3]<-kantaja[i,2]     #max(dendat[,i])+epsi 
     valpit[i]<-(binhila[i,3]-binhila[i,2])/binlkm
       #binhila:n i:s rivi sis. i:nnelle muuttujalle 
       #lokeroinnin alkupisteen, arvoalueen loppupisteen   
       #valpit: arvoalueen pituus jaettuna lokeroitten lkm:lla
       #eli yhden lokeron leveys. 
     i<-i+1
}
hila[,2]<-binhila[,2]+valpit/2
hila[,3]<-binhila[,3]-valpit/2
       #hila:n i:s rivi sis. i:nnelle muuttujalle 
       #ensimmaisen regressiodatan havaintopisteen,
       #viimeisen regressiodatan havaintopisteen
#if (valpit<=0) stop("in some variable there is no variation")
hiladat<-matrix(1,n,xlkm) #muunnetaan dendat hiladat:iksi
                          #ts. pyoristetaan havainnot 
                          #hiladat sis. diskretoidut havainnot
i<-1
while (i<=n){       #kaydaan lapi aineisto
    #pyoristetaan i:s havainto hilapisteeseen
    j<-1
    while (j<=xlkm){
       alavali<-floor((dendat[i,j]-binhila[j,2])/valpit[j])
           #alavali ilmaisee monennessako lokerossa hav. sijaitsee
       hiladat[i,j]<-binhila[j,2]+alavali*valpit[j]+valpit[j]/2
       j<-j+1
    }
i<-i+1
}
xtulos<-matrix(0,n,xlkm) #periaatteessa mahdollista etta kaikki n
                         #havaintoa ovat eri lokeroissa, siksi
                         #laitetaan xtulos matriisin rivien maaraksi n
ytulos<-matrix(0,n,1)
xtulos[1,]<-hiladat[1,]  #hiladat:in ensimmainen rivi esiintyy ainakin kerran
ytulos[1]<-1             #sen frekvenssi ainakin yksi
saatu<-1                 #toistaiseksi yksi erillinen havainto
i<-1
while (i<n){  #kaydaan lapi aineisto
   i<-i+1     #aloitetaan kakkosesta
   lippu<-0   #apriori kyseessa uusi lajityyppi
   j<-1
   while ((j<=saatu) && (lippu==0)){  #kaydaan lapi keratyt havinnot
       if (all(hiladat[i,]==xtulos[j,])){  #jos on jo saatu
            lippu<-1     #liputetaan etta havaittiin toisto
            jind<-j      #merkataan indeksi frekvenssin paivitysta varten
       }
       j<-j+1       
   }
   if (lippu==1) ytulos[jind]<-ytulos[jind]+1 
            #jos saatiin toisto, paivitetaan frekvenssi 
      else{ 
         saatu<-saatu+1      #jos saatiin uusi, lisataan saatu:un yksi ja
         xtulos[saatu,]<-hiladat[i,]  #merkitaan uusi lajityyppi muistiin
         ytulos[saatu]<-1    #uuden lajityypin frekvenssi on aluksi yksi
      }
}
xtulos<-xtulos[1:saatu,]
ytulos<-ytulos[1:saatu]
ytulos<-t(t(ytulos))
if (xlkm==1) xtulos<-t(t(xtulos))
return(list(dep=ytulos,ind=xtulos,hila=hila))
}

dend2parent<-function(hc,dendat)
{
n<-dim(dendat)[1]
d<-dim(dendat)[2]
nodenum<-length(hc$height)+n
parent<-matrix(0,nodenum,1)
volume<-matrix(0,nodenum,1)
center<-matrix(0,d,nodenum)
level<-matrix(0,nodenum,1)
level[(n+1):nodenum]<-hc$height
pointers<-matrix(0,n,1)      #pointers dendat to tree 

joinnum<-length(hc$height)
lkm<-matrix(0,joinnum,1)
vec1<-matrix(0,1,d)
vec2<-matrix(0,1,d)
v1<-0
v2<-0
vapaa<-1
for (i in 1:joinnum){
   node1<-hc$merge[i,1]
   node2<-hc$merge[i,2]
   if (node1<0){
        parent[vapaa]<-i+n
        for (j in 1:d) vec1[j]<-dendat[-node1,j]
        v1<-1
        center[,vapaa]<-vec1
        volume[vapaa]<-v1
        pointers[-node1]<-vapaa
        lkm1<-1
        vapaa<-vapaa+1
   }
   else{
        parent[node1+n]<-i+n
        vec1<-center[,node1+n]
        v1<-volume[node1+n]
        lkm1<-lkm[node1]
   }
   if (node2<0){
        parent[vapaa]<-i+n
        for (j in 1:d) vec2[j]<-dendat[-node2,j]
        v2<-1
        center[,vapaa]<-vec2
        volume[vapaa]<-v2
        pointers[-node2]<-vapaa
        lkm2<-1
        vapaa<-vapaa+1
   }
   else{
        parent[node2+n]<-i+n
        vec2<-center[,node2+n]
        v2<-volume[node2+n]
        lkm2<-lkm[node2]
   }
   volume[i+n]<-1.1*(v1+v2)
   center[,i+n]<-(lkm1*vec1+lkm2*vec2)/(lkm1+lkm2)
   lkm[i]<-lkm1+lkm2
}

apoin<-matrix(0,n,1)
for (i in 1:n) apoin[i]<-nodenum-pointers[i]+1
apar<-parent[nodenum:1]
apar2<-matrix(0,n,1)
for (i in 1:nodenum) if (apar[i]!=0) apar2[i]<-nodenum-apar[i]+1

return(list(parent=apar2,level=level[nodenum:1],
volume=volume[nodenum:1],center=center[,nodenum:1],pointers=apoin))
}


dendat2lst<-function(dendat,lst,pcf)
{
# compare liketree

rnum<-length(pcf$value)
nodefinder<-matrix(0,rnum,1)
for (i in 1:rnum) nodefinder[lst$infopointer[i]]<-i

n<-dim(dendat)[1]
d<-dim(dendat)[2]
den2lst<-matrix(0,n,1)

step<-matrix(0,d,1)
for (i in 1:d) step[i]<-(pcf$support[2*i]-pcf$support[2*i-1])/pcf$N[i]

# find links from dendat to pcf
den2pcf<-matrix(0,n,1)
pcf2den<-matrix(0,rnum,1)
for (i in 1:n){
    j<-1
    while (j<=rnum){
         inside<-TRUE
         coordi<-1
         while ((inside) && (coordi<=d)){
             ala<-pcf$down[j,coordi]
             yla<-pcf$high[j,coordi]
             ala<-pcf$support[2*coordi-1]+ala*step[coordi]
             yla<-pcf$support[2*coordi-1]+yla*step[coordi]
             if ((dendat[i,coordi]<ala) || (dendat[i,coordi]>yla)) 
                         inside<-FALSE
             coordi<-coordi+1
         }
         if (inside){
                  den2pcf[i]<-j
                  pcf2den[j]<-i
         }
         j<-j+1
    }
}


for (i in 1:n){
    pcfi<-den2pcf[i]
    den2lst[i]<-nodefinder[pcfi]
}

return(den2lst)
}

depth2com<-function(dep,N){
#
d<-length(N)
logn<-log(N,base=2)
cusu<-cumsum(logn)
ind<-1
while ((ind<=d) && ((dep-cusu[ind])>0)){
  ind<-ind+1
}
direc<-min(ind,d)
if (direc==1){
  depind<-dep
}
else{
  depind<-dep-cusu[direc-1]
}
return(list(direc=direc,depind=depind))
}
depth<-function(mt){
#finds the dephts of the nodes
#
#mt is a result from multitree or pruneprof
#
roots<-mt$roots
child<-mt$child
sibling<-mt$sibling
#
itemnum<-length(child)
depths<-matrix(0,itemnum,1)
#
rootnum<-length(roots)
#
for (i in 1:rootnum){
    pino<-matrix(0,itemnum,1)
    pino[1]<-roots[i]  
    pinin<-1
    depths[roots[i]]<-1
    while (pinin>0){
        cur<-pino[pinin]      #take from stack
        pinin<-pinin-1
        if (sibling[cur]>0){    #put right to stack
             pinin<-pinin+1
             pino[pinin]<-sibling[cur]
             depths[sibling[cur]]<-depths[cur]
        }
        while (child[cur]>0){    #go to leaf and put right nodes to stack
             chi<-child[cur]
             depths[chi]<-depths[cur]+1
             if (sibling[chi]>0){ 
                   pinin<-pinin+1
                   pino[pinin]<-sibling[chi]
                   depths[sibling[chi]]<-depths[cur]+1
             }
             cur<-chi
        }
    }
}
return(depths)
}







digit<-function(luku,base){
#Gives representation of luku for system with base
#
#luku is a natural number >=0
#base is d-vector of integers >=2, d>=2, 
#base[d] tarvitaan vain tarkistamaan onko luku rajoissa
#
#Returns d-vector of integers.
#
#example: digit(52,c(10,10)), returns vector (2,5)
#
d<-length(base)
digi<-matrix(0,d,1)
jako<-matrix(0,d,1)
jako[d]<-base[1]
for (i in (d-1):1){
  jako[i]<-base[d-i+1]*jako[i+1]
}
vah<-0
for (i in 1:(d-1)){
  digi[i]<-floor((luku-vah)/jako[i+1]) #if digi[i]>base[i], then ERROR
  vah<-vah+digi[i]*jako[i+1]
}
digi[d]<-luku-vah
# annetaan vastaus kaanteisesti se 2354 annetaan c(4,5,3,2)
# talloin vastaavuus sailyy base:n kanssa 
#apu<-matrix(0,d,1)
#for (i in 1:d){
#  apu[i]<-digi[d-i+1]
#}
apu<-digi[d:1]
return(apu)
}
dist.func<-function(dendat,xepsi=0,yepsi=0,col="black",type="distr",
log="y",cex.axis=1,dendat2=NULL,dendat3=NULL,col2="red",col3="blue",
pch2=20,pch3=20,split=median(dendat),xlim=NULL,xaxt="s",yaxt="s")
{
n<-length(dendat)

if (type=="distr"){

plot(x="",y="",
xlim=c(min(dendat)-xepsi,max(dendat)+xepsi), ylim=c(0-yepsi,1+yepsi),
xlab="",ylab="",cex.axis=cex.axis)
ycur<-0
ordi<-order(dendat)
dendatord<-dendat[ordi]
for(i in 1:(n-1)){
    segments(dendatord[i],ycur,dendatord[i],ycur+1/n,col=col)
    segments(dendatord[i],ycur+1/n,dendatord[i+1],ycur+1/n,col=col)
    ycur<-ycur+1/n
}
segments(dendatord[n],ycur,dendatord[n],ycur+1/n,col=col)
segments(dendatord[n],ycur+1/n,max(dendat)+xepsi,ycur+1/n,col=col)

}
else if ((type=="right.tail") || (type=="left.tail")){

  if (type=="right.tail"){
         redu.ind<-(dendat>split) 
         dendat.redu<-dendat[redu.ind]
  }
  else{
         redu.ind<-(dendat<split)
         dendat.redu<--dendat[redu.ind]
  }
  ordi<-order(dendat.redu)
  dendat.ord<-dendat.redu[ordi]
  nredu<-length(dendat.redu)
  level<-seq(nredu,1)
  if (type=="right.tail")
  plot(dendat.ord,level,log=log,xlab="",ylab="",cex.axis=cex.axis,xlim=xlim)
  else
  plot(-dendat.ord,level,log=log,xlab="",ylab="",cex.axis=cex.axis,xlim=xlim,
  xaxt=xaxt,yaxt=yaxt)

  #ordi<-order(dendat)
  #dendat.ord<-dendat[ordi]
  #medi.ind<-floor(n/2)
  #dendat.redu<-dendat.ord[medi.ind:n]
  #nredu<-length(dendat.redu)
  #level<-seq(nredu,1)
  #plot(dendat.redu,level,log="y",xlab="",ylab="")

  if (!is.null(dendat2)){

     if (type=="right.tail"){
          redu.ind<-(dendat2>split) 
          dendat.redu<-dendat2[redu.ind]
     }
     else{
          redu.ind<-(dendat2<split)
          dendat.redu<--dendat2[redu.ind]
     }
     ordi<-order(dendat.redu)
     dendat.ord<-dendat.redu[ordi]
     nredu<-length(dendat.redu)
     level<-seq(nredu,1)
     if (type=="right.tail")
     matplot(dendat.ord,level,log=log,xlab="",ylab="",cex.axis=cex.axis,
     add=TRUE,col=col2,pch=pch2)
     else
     matplot(-dendat.ord,level,log=log,xlab="",ylab="",cex.axis=cex.axis,
     add=TRUE,col=col2,pch=pch2)

  }

  if (!is.null(dendat3)){

     if (type=="right.tail"){
          redu.ind<-(dendat3>split) 
          dendat.redu<-dendat3[redu.ind]
     }
     else{
          redu.ind<-(dendat3<split)
          dendat.redu<--dendat3[redu.ind]
     }
     ordi<-order(dendat.redu)
     dendat.ord<-dendat.redu[ordi]
     nredu<-length(dendat.redu)
     level<-seq(nredu,1)
      if (type=="right.tail")
     matplot(dendat.ord,level,log=log,xlab="",ylab="",cex.axis=cex.axis,
     add=TRUE,col=col3,pch=pch3)
     else
     matplot(-dendat.ord,level,log=log,xlab="",ylab="",cex.axis=cex.axis,
     add=TRUE,col=col3,pch=pch3)
  }

}

}

dotouchgen<-function(indelow1,indeupp1,indelow2,indeupp2,direction){
#
epsi<-0
d<-length(indelow1)
touch<-TRUE
i<-1
while (i<=d){
     if ((i != direction) &&
         ((indelow1[i]>indeupp2[i]+epsi) || (indeupp1[i]<indelow2[i]-epsi))){
              touch<-FALSE
     }
     i<-i+1
}
return(touch)
}
dotouch<-function(inde1,inde2,direction){
#
d<-length(inde1)
touch<-TRUE
for (i in direction:d){
     if ((inde1[i]>inde2[i]+1) || (inde1[i]<inde2[i]-1)){
           touch<-FALSE
     }
}
return(touch)
}
drawgene<-function(values,recs,plkm=60,ep1=0.5){
#Makes data for drawing a perspective plot.

#plkm on kuvaajan hilan pisteiden lkm
#ep1 makes 0-corona around the support (useful for densities)

#koe<-drawgene(values,recs,plkm=30)
#persp(koe$x,koe$y,koe$z,phi=30,theta=60) 

alkux<-min(recs[,1])-ep1
alkuy<-min(recs[,3])-ep1
loppux<-max(recs[,2])+ep1
loppuy<-max(recs[,4])+ep1
pitx<-(loppux-alkux)/plkm
pity<-(loppuy-alkuy)/plkm
x<-alkux+c(0:plkm)*pitx
y<-alkuy+c(0:plkm)*pity

reclkm<-length(values)
xdim<-length(x)
ydim<-length(y)
arvot<-matrix(0,xdim,ydim)

l<-1
while (l<=reclkm){
   begx<-recs[l,1]
   endx<-recs[l,2]
   begy<-recs[l,3]
   endy<-recs[l,4]

   begxind<-round(plkm*(begx-alkux)/(loppux-alkux))
   endxind<-round(plkm*(endx-alkux)/(loppux-alkux))
   begyind<-round(plkm*(begy-alkuy)/(loppuy-alkuy))
   endyind<-round(plkm*(endy-alkuy)/(loppuy-alkuy))

   arvot[begxind:endxind,begyind:endyind]<-values[l]

   l<-l+1
}

return(list(x=x,y=y,z=arvot))
#persp(x,y,arvot)
}










drawhist<-function(dendat,binlkm,epsi=0,plkm){
#piirtaa 2-ulotteisessa tapauksessa histogramma estimaattorin kuvaajan 
#
#plkm on kuvaajan hilan pisteiden lkm
#
#dendat<-matrix(rnorm(20),10) 
#koe<-drawhist(dendat,binlkm=3,plk=30)
#persp(koe$x,koe$y,koe$z,phi=30,theta=60)
#
hi<-histo(dendat,binlkm,epsi)
recs<-hi$recs
values<-hi$values   
#
ans<-drawgene(values,recs,plkm)
return(list(x=ans$x,y=ans$y,z=ans$z))
}                                      








draw.kern<-function(value,index,N,support,minval=0,dendat=NULL,h=NULL)
{

d<-length(N)

if (d==2){

x<-matrix(0,N[1]+2,1)
y<-matrix(0,N[2]+2,1)
z<-matrix(minval,N[1]+2,N[2]+2)
#col<-matrix("black",dim(z)[1]*dim(z)[2],1)

minim<-matrix(0,d,1)  #minim is d-vector of minimums
maxim<-matrix(0,d,1)
for (i in 1:d){
  minim[i]<-support[2*i-1]
  maxim[i]<-support[2*i]
}
delta<-(maxim-minim)/(N+1)

indenum<-dim(index)[1]

i<-1
while (i<=indenum){
   inde<-index[i,]
   z[1+inde[1],1+inde[2]]<-value[i]
   #col[1+inde[1]+dim(z)[1]*inde[2]]<-ts[i]
   i<-i+1
}

i<-1
while (i<=N[1]){
   x[1+i]<-support[1]+delta[1]*i
   i<-i+1
}

i<-1
while (i<=N[2]){
   y[1+i]<-support[3]+delta[2]*i
   i<-i+1
}

x[1]<-support[1]
x[N[1]+2]<-support[2]
y[1]<-support[3]
y[N[2]+2]<-support[4]

return(list(x=x,y=y,z=z)) #col=col[length(col):1]))

}

else{    #d=1

x<-matrix(0,N+2,1)
y<-matrix(0,N+2,1)

minim<-min(dendat)
maxim<-max(dendat)
hmax<-max(h)
delta<-(maxim-minim+2*hmax)/(N+1)

indenum<-dim(index)[1]

i<-1
while (i<=indenum){

   inde<-index[i]
   point<-minim-hmax+delta*inde
    
   y[1+inde]<-value[i]
   x[1+inde]<-point

   i<-i+1
}
x[1]<-minim-hmax
x[N+2]<-minim-hmax+delta*N+delta

return(list(x=x,y=y))

}

}











draw.levset<-function(pcf,lev=NULL,bary=NULL,propor=0.1,col=NULL,
bound=NULL,dendat=NULL,xaxt="s",yaxt="s",cex.axis=1)
{

if (is.null(lev)) lev<-propor*max(pcf$value)

d<-length(pcf$N)
step<-matrix(0,d,1)
for (i in 1:d) step[i]=(pcf$support[2*i]-pcf$support[2*i-1])/pcf$N[i];

if (is.null(bound)){
  xmin<-pcf$support[1]
  xmax<-pcf$support[2]
  ymin<-pcf$support[3]
  ymax<-pcf$support[4]
}
else{
  xmin<-bound[1]
  xmax<-bound[2]
  ymin<-bound[3]
  ymax<-bound[4]
}

if (is.null(bary))
   plot(xmin,ymin,type="n",xlab="",ylab="",xlim=c(xmin,xmax),ylim=c(ymin,ymax),
   pch=20,col="red",xaxt=xaxt,yaxt=yaxt,cex.axis=cex.axis)
else
  plot(x=bary[1],y=bary[2],
  xlab="",ylab="",xlim=c(xmin,xmax),ylim=c(ymin,ymax),
  pch=20,col="red",xaxt=xaxt,yaxt=yaxt,cex.axis=cex.axis)

lenni<-length(pcf$value)
for (i in 1:lenni){
  if (pcf$value[i]>=lev){

     x1<-pcf$support[1]+step[1]*pcf$down[i,1]
     x2<-pcf$support[1]+step[1]*pcf$high[i,1] 
     y1<-pcf$support[3]+step[2]*pcf$down[i,2]
     y2<-pcf$support[3]+step[2]*pcf$high[i,2] 

     if (is.null(col)) polygon(c(x1,x2,x2,x1),c(y1,y1,y2,y2))  
     else  polygon(c(x1,x2,x2,x1),c(y1,y1,y2,y2),col=col[i],lty="blank")
  }
  i<-i+1
}

if (!is.null(dendat)) points(dendat)

#points(x=bary[1],y=bary[2],pch=20,col="red")

}



draw.pcf<-function(pcf,pnum=rep(32,length(pcf$N)),corona=5,dens=FALSE,minval=0,
drawkern=TRUE)
{
#Makes data for drawing a perspective plot.
#pnum on kuvaajan hilan pisteiden lkm
#corona makes corona around the support (useful for densities)

d<-length(pcf$N)

if (d==2){

#col=matrix("black",length(pcf$value),1)

step<-matrix(0,d,1)
for (i in 1:d){
   step[i]<-(pcf$support[2*i]-pcf$support[2*i-1])/pcf$N[i]
}

if ((drawkern)&&(!is.null(pcf$index))){
   return(draw.kern(pcf$value,pcf$index,pcf$N,pcf$support,minval=minval))
}

else{
     pit<-matrix(0,d,1)
     for (i in 1:d){
         pit[i]<-(pcf$support[2*i]-pcf$support[2*i-1])/pnum[i]
     }
     alkux<-pcf$support[1]+pit[1]/2-corona*pit[1]
     alkuy<-pcf$support[3]+pit[2]/2-corona*pit[2]
     loppux<-pcf$support[2]-pit[1]/2+corona*pit[1]
     loppuy<-pcf$support[4]-pit[2]/2+corona*pit[2]

     pnum2<-pnum+2*corona
     x<-alkux+c(0:(pnum2[1]-1))*pit[1]
     y<-alkuy+c(0:(pnum2[2]-1))*pit[2]

     reclkm<-length(pcf$value)
     xdim<-length(x)
     ydim<-length(y)
     arvot<-matrix(minval,xdim,ydim)

     l<-1
     while (l<=reclkm){
         begx<-pcf$support[1]+step[1]*(pcf$down[l,1])   #recs[l,1]
         endx<-pcf$support[1]+step[1]*(pcf$high[l,1])     #recs[l,2]
         begy<-pcf$support[3]+step[2]*(pcf$down[l,2])   #recs[l,3]
         endy<-pcf$support[3]+step[2]*(pcf$high[l,2])     #recs[l,4]

         begxind<-round(pnum2[1]*(begx-alkux)/(loppux-alkux))
         endxind<-round(pnum2[1]*(endx-alkux)/(loppux-alkux))
         begyind<-round(pnum2[2]*(begy-alkuy)/(loppuy-alkuy))
         endyind<-round(pnum2[2]*(endy-alkuy)/(loppuy-alkuy))

         arvot[begxind:endxind,begyind:endyind]<-pcf$value[l]
         #col[(begxind+(ydim-1)*begxind):(endxind+(ydim-1)*endyind)]<-ts[l]

         l<-l+1
     }

     #return(apu)
     return(list(x=x,y=y,z=arvot))
}  # else

}  # if (d==2)

else{  #d==1

if (dens){ 
  N<-pcf$N+2 
  alku<-1
}
else{
  N<-pcf$N
  alku<-0
}

x<-matrix(0,N,1)
y<-matrix(0,N,1)
step<-(pcf$support[2]-pcf$support[1])/pcf$N
minim<-pcf$support[1]
maxim<-pcf$support[2]
indenum<-length(pcf$value)

i<-1
while (i<=indenum){

   inde<-pcf$high[i]
   point<-minim+step*inde-step/2
    
   y[alku+inde]<-pcf$value[i]
   x[alku+inde]<-point

   i<-i+1
}

if (dens){
  x[1]<-minim-step/2
  x[N]<-maxim+step/2
}

# remove zeros
if (dens){
  numposi<-0
  for (i in 1:length(y)){
    if (y[i]>0){
       numposi<-numposi+1
       y[numposi]<-y[i]
       x[numposi]<-x[i]
    }
  }
  x<-x[1:numposi]
  y<-y[1:numposi]
}

or<-order(x)
xor<-x[or]
yor<-y[or]

return(list(x=xor,y=yor))

}  # else d=1

}











epane<-function(x,h){
#
d<-length(x)
val<-1
for (i in 1:d){
   val<-val*3*(1-(x[i]/h)^2)/2
}
return(val)
}
epan<-function(x)
{
d<-length(x)
val<-1
for (i in 1:d){
   val<-val*3*(1-x[i]^2)/2
}
return(val)
}
etais<-function(x,y){
#laskee euklid etais nelion vektorien x ja y valilla
#
pit<-length(x)
vast<-0
i<-1
while (i<=pit){
  vast<-vast+(x[i]-y[i])^2
  i<-i+1
}
return(vast)
}
etaisrec<-function(point,rec)
{
# calculates the squared diatance of a point to a rectangle

d<-length(rec)/2

res<-0
for (i in 1:d){
   if (point[i]>rec[2*i]) res<-res+(point[i]-rec[2*i])^2
   else if (point[i]<rec[2*i-1]) res<-res+(point[i]-rec[2*i-1])^2
}

return(res)

}



eva.clayton<-function(x,t,marginal="unif",sig=c(1,1),df=1)
# t>0
{
u<-x[1]
v<-x[2]
marg1<-1
marg2<-1

if (marginal=="normal"){
   u<-pnorm(x[1]/sig[1])
   v<-pnorm(x[2]/sig[2])
   marg1<-evanor(x[1]/sig[1])/sig[1]
   marg2<-evanor(x[2]/sig[2])/sig[2]
}
if (marginal=="student"){
   u<-pt(x[1],df)
   v<-pt(x[2],df)
   marg1<-dt(x[1],df)
   marg2<-dt(x[2],df)
}

val<-(1+t)*(u*v)^(-1-t)*(u^(-t)+v^(-t)-1)^(-2-1/t)*marg1*marg2

#if (val<0) val<-0

return(val)
}

eva.cop6<-function(x,t,marginal="unif",sig=c(1,1))
# t in [1,\infty)
{
u<-x[1]
v<-x[2]
marg1<-1
marg2<-1

if (marginal=="normal"){
   u<-pnorm(x[1]/sig[1])
   v<-pnorm(x[2]/sig[2])
   marg1<-evanor(x[1]/sig[2])/sig[1]
   marg2<-evanor(x[2]/sig[2])/sig[2]
}

val<-t*(1-u)^(t-1)*(1-v)^(t-1)*
((1-u)^t+(1-v)^t-(1-u)^t*(1-v)^t)^(1/t-2)*
(-(1/t-1)*(1-(1-u)^t)*(1-(1-v)^t)+
 (1-u)^t+(1-v)^t-(1-u)^t*(1-v)^t)*marg1*marg2

#if (val<0) val<-0

return(val)
}

eva.copula<-function(x,type="gauss",marginal="unif",sig=rep(1,length(x)),r=0,
t=rep(4,length(x)),g=1)
{
# sig is std of marginals, r is the correlation coeff, 
# t is deg of freedom

d<-length(x)
marg<-matrix(0,d,1)
u<-matrix(0,d,1)

if (marginal=="unif"){
   for (i in 1:d){
      u[i]<-x[i]/sig[i]  #+1/2
      marg[i]<-1/sig[i]
   }
}
if ((marginal=="normal")||(marginal=="gauss")){
   for (i in 1:d){
      u[i]<-pnorm(x[i]/sig[i])
      marg[i]<-evanor(x[i]/sig[i])/sig[i]
   }
}
if (marginal=="student"){
   for (i in 1:d){
      u[i]<-pt(x[i]/sig[i],df=t[i])
      marg[i]<-dt(x[i]/sig[i],df=t[i])/sig[i]
   }
}
if (type=="gauss"){
   d<-2
   x1<-qnorm(u[1],sd=1)
   x2<-qnorm(u[2],sd=1)

#   produ<-dnorm(x1,sd=1)*dnorm(x2,sd=1)
#   nelio<-(x1^2+x2^2-2*r*x1*x2)/(1-r^2)
#   vakio<-(2*pi)^(-d/2) 
#   g<-vakio*(1-r^2)^(-1/2)*exp(-(1/2)*nelio)
#   val<-g/produ*marg[1]*marg[2]

  copuval<-(1-r^2)^(-1/2)*
  exp(-(x1^2+x2^2-2*r*x1*x2)/(2*(1-r^2)))/exp(-(x1^2+x2^2)/2)
  val<-copuval*marg[1]*marg[2]

}
if (type=="gumbel"){
  link<-function(y,g){ return ( (-log(y))^g ) }
  linkinv<-function(y,g){ return ( exp(-y^(1/g)) ) }
  der1<-function(y,g){ return ( -g*(-log(y))^(g-1)/y ) }
  der2<-function(y,g){ return ( g*y^(-2)*(-log(y))^(g-2)*(g-1-log(y)) ) }
  
  linky<-link(u,g)
  a<-sum(linky)
  b<-linkinv(a,g)
  der1b<-der1(b,g)
  der2b<-der2(b,g)
  psi<--der2b*der1b^(-3)
  deriy<-der1(u,g) 
  val<-psi*abs(prod(deriy))*prod(marg)
}
if (type=="frank"){
  link<-function(y,g){ return ( -log((exp(-g*y)-1)/(exp(-g)-1)) ) }
  linkinv<-function(y,g){ return ( -log(1+(exp(-g)-1)/exp(y))/g  ) }
  der1<-function(y,g){ return ( g*exp(-g*y)/(exp(-g*y)-1) ) }
  der2<-function(y,g){ return ( g^2*exp(-g*y)/(exp(-g*y)-1)^2 ) }
  
  linky<-link(u,g)
  a<-sum(linky)
  b<-linkinv(a,g)
  der1b<-der1(b,g)
  der2b<-der2(b,g)
  psi<--der2b*der1b^(-3)
  deriy<-der1(u,g) 
  val<-psi*abs(prod(deriy))*prod(marg)
}
if (type=="clayton"){
  link<-function(y,g){ return ( y^(-g)-1 ) }
  linkinv<-function(y,g){ return ( (y+1)^(-1/g) ) }
  der1<-function(y,g){ return ( -g*y^(-g-1) ) }
  der2<-function(y,g){ return ( g*(g+1)*y^(-g-2) ) }
  
  linky<-link(u,g)
  a<-sum(linky)
  b<-linkinv(a,g)
  der1b<-der1(b,g)
  der2b<-der2(b,g)
  psi<--der2b*der1b^(-3)
  deriy<-der1(u,g) 
  val<-psi*abs(prod(deriy))*prod(marg)
}

return(val)
}






eva.gauss<-function(x,t=1,marginal="unif",sig=c(1,1),r=0,tapa1=TRUE)
{
#  sig is std of marginals

if (marginal=="unif"){
   u<-x[1]/sig[1]+1/2
   v<-x[2]/sig[2]+1/2
   marg1<-1/sig[1]
   marg2<-1/sig[2]
}
if (marginal=="normal"){
   u<-pnorm(x[1]/sig[1])
   v<-pnorm(x[2]/sig[2])
   marg1<-evanor(x[1]/sig[1])/sig[1]
   marg2<-evanor(x[2]/sig[2])/sig[2]
}
if (marginal=="student"){
   u<-pt(x[1]/sig[1],df=t)
   v<-pt(x[2]/sig[2],df=t)
   marg1<-dt(x[1]/sig[1],df=t)/sig[1]
   marg2<-dt(x[2]/sig[2],df=t)/sig[2]
}

d<-2
x1<-qnorm(u,sd=1)
x2<-qnorm(v,sd=1)

if (tapa1){
  produ<-dnorm(x1,sd=1)*dnorm(x2,sd=1)
  nelio<-(x1^2+x2^2-2*r*x1*x2)/(1-r^2)
  vakio<-(2*pi)^(-d/2) 
  g<-vakio*(1-r^2)^(-1/2)*exp(-(1/2)*nelio)
  val<-g/produ*marg1*marg2
}
else{
  # x1,x2 -> copuval
  copuval<-(1-r^2)^(-1/2)*
  exp(-(x1^2+x2^2-2*r*x1*x2)/(2*(1-r^2)))/exp(-(x1^2+x2^2)/2)
  val<-copuval*marg1*marg2
}

return(val)
}

eva.hat<-function(x,a=0.5,b=0.5)
{
# 0<a<1, b<1
# if b<a then marginal is unimodal
# if a^2 < b < a then not star unimodal

d<-length(x)
eta<-sum(x^2)     #vektorin x pituuden nelio
normvakio<-((2*pi)^d*(a^(-d)-b))^(-1)
tulos<-normvakio*(exp(-a^2*eta/2)-b*exp(-eta/2))

return(tulos)
}
eval.func.1D<-function(func,N,support=NULL,g=1,std=1,distr=FALSE,
M=NULL,sig=NULL,p=NULL,a=0.5,b=0.5,d=2)
{
if (func=="gauss"){
   norma<-(2*pi)^(-1/2)
   funni<-function(t){ fu<-exp(-t^2/2); return( norma*fu ) }
}
if (func=="polynomial"){
   support<-c(-std,std)
   norma<-(2*(1-1/(g+1)))^(-1)
   funni<-function(t){ fu<-1-abs(t)^g; return( norma*fu ) }
}
if (func=="student"){
   norma<-gamma((g+1)/2)/((g*pi)^(1/2)*gamma(g/2))
   funni<-function(t){ fu<-(1+t^2/g)^(-(g+1)/2); return( norma*fu ) }
   #y<-dt(x,df=g)
}
if (func=="exponential"){
   norma<-1/2
   funni<-function(t){ fu<-exp(-abs(t)); return( norma*fu ) }
}
if (func=="exponential"){
   norma<-1/2
   funni<-function(t){ fu<-exp(-abs(t)); return( norma*fu ) }
}
if (func=="mixt"){
   funni<-function(t){ 
       mixnum<-length(p)
       val<-0
       for (mi in 1:mixnum){
               evapoint<-(t-M[mi])/sig[mi]
               val<-val+p[mi]*evanor(evapoint)/sig[mi]
        }
        return( val ) 
   }
}
if (func=="hat"){
   normavak<-((2*pi)^d*(a^(-d)-b))^(-1)
   norma<-normavak*(2*pi)^((d-1)/2)
   funni<-function(t){  #(t,a,b,d,...){ 
          fu<-a^(1-d)*exp(-a^2*t^2)-b*exp(-t^2/2); return( norma*fu ) }
}

if (is.null(support)) support<-c(-1,1)

value<-matrix(0,N,1)
step<-(support[2]-support[1])/N
lowsuppo<-support[1]

if (!distr){
   for (i in 1:N){
       inde<-i
       t<-lowsuppo+step*inde-step/2
       value[i]<-funni(t/std)/std
   }
}
else{
   inde<-1
   t<-lowsuppo+step*inde-step/2
   value[1]<-step*funni(t/std)/std
   for (i in 2:N){
       inde<-i
       t<-lowsuppo+step*inde-step/2
       value[i]<-value[i-1]+step*funni(t/std)/std
       #funni(t/std,g=g,a=a,b=b,d=d)/std
   }
}

index<-seq(1:N)
len<-length(index)
down<-matrix(0,len,1)
high<-matrix(0,len,1)
down[,1]<-index-1
high[,1]<-index

res<-list(
value=value,
down=down,high=high,
#down=index-1,high=index,  
support=support,N=N)

return(res)
}

                              

eval.func.dD<-function(func,N,
sig=rep(1,length(N)),support=NULL,theta=NULL,g=1,
M=NULL,p=NULL,mul=3,
t=rep(1,length(N)),marginal="normal",r=0,
mu=NULL,xi=NULL,Omega=NULL,alpha=NULL,df=NULL,a=0.5,b=0.5
)   
# func== "mixt", "epan", "cop1"
{
d<-length(N)
recnum<-prod(N)
value<-matrix(0,recnum,1)
index<-matrix(0,recnum,d)

if (is.null(support)){

   if (func=="mixt"){
     support<-matrix(0,2*d,1)
        for (i in 1:d){
           support[2*i-1]<-min(M[,i]-mul*sig[,i])
           support[2*i]<-max(M[,i]+mul*sig[,i])
        }
   }

   if (func=="epan"){
      if (is.null(sig)) sig<-c(1,1)
      support<-matrix(0,2*d,1)
      for (i in 1:d){
          support[2*i-1]<--sig[i]
          support[2*i]<-sig[i]
      }
   }
}

if ((marginal=="unif")) support<-c(0,sig[1],0,sig[2])
# && (is.null(support))) 
#support<-c(-sig[1]/2,sig[1]/2,-sig[2]/2,sig[2]/2)


lowsuppo<-matrix(0,d,1)
for (i in 1:d) lowsuppo[i]<-support[2*i-1]
step<-matrix(0,d,1)
for (i in 1:d) step[i]<-(support[2*i]-support[2*i-1])/N[i]

numpositive<-0
for (i in 1:recnum){
    inde<-digit(i-1,N)+1
    #if ((inde[1]==0) && (inde[2]==N[2])) inde<-c(0,0)
    point<-lowsuppo+step*inde-step/2

    if (!is.null(theta)){
         rotmat<-matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),2,2)
         point<-rotmat%*%point
    }

    if (func=="prod") val<-eva.prod(point,marginal,g)
    if (func=="skewgauss") val<-eva.skewgauss(point,mu,sig,alpha)
    #if (func=="dmsn") val<-dmsn(point,xi,Omega,alpha)
    if (func=="student") val<-eva.student(point,t,marginal,sig,r,df)
    if (func=="gumbel") val<-eva.copula(point,
        type="gumbel",marginal=marginal,sig=sig,r=r,t=t,g=g)
    if (func=="frank") val<-eva.copula(point,
        type="frank",marginal=marginal,sig=sig,t=t,g=g)
    if (func=="plackett") val<-eva.plackett(point,t,marginal,sig)
    if (func=="clayton2") val<-eva.clayton(point,t,marginal,sig,df)
    if (func=="clayton") val<-eva.copula(point,
        type="clayton",marginal=marginal,sig=sig,r=r,t=t,g=g)
    if (func=="cop6") val<-eva.cop6(point,t,marginal,sig)
    if (func=="epan") val<-epan(point)
    if (func=="gauss") val<-eva.copula(point,
        type="gauss",marginal=marginal,sig=sig,r=r,t=t)
    if (func=="normal") val<-eva.gauss(point,t=t,marginal=marginal,sig=sig,r=r)
    if (func=="mixt"){
        val<-0
        mixnum<-length(p)
        for (mi in 1:mixnum){
            evapoint<-(point-M[mi,])/sig[mi,]
            val<-val+p[mi]*evanor(evapoint)/prod(sig[mi,])
        }
    }
   if (func=="hat") val<-eva.hat(point,a=a,b=b)

    if (val>0){
       numpositive<-numpositive+1
       value[numpositive]<-val
       index[numpositive,]<-inde
    }
}

value<-value[1:numpositive]
index<-index[1:numpositive,]
down<-index-1
high<-index

res<-list(
value=value,index=index,
down=down,high=high,  #step=delta,
support=support,N=N)

return(res)
}                              











eva.lognormal<-function(t)
{

ans<-(2*pi)^(-1/2)*t^(-1)*exp(-(log(t))^2/2)

return(ans)
}

evanor<-function(x){

d<-length(x) 
eta<-sum(x^2)     #vektorin x pituuden nelio
normvakio<-(sqrt(2*pi))^{-d}
tulos<-exp(-eta/2)*normvakio
return(tulos)
}
eva.plackett<-function(x,t,marginal="unif",sig=c(1,1))
# t>=0, t \neq 1
{
u<-x[1]
v<-x[2]
marg1<-1
marg2<-1

if (marginal=="normal"){
   u<-pnorm(x[1]/sig[1])
   v<-pnorm(x[2]/sig[2])
   marg1<-evanor(x[1]/sig[1])/sig[1]
   marg2<-evanor(x[2]/sig[2])/sig[2]
}

val<-t*(1+(t-1)*(u+v-2*u*v))*((1+(t-1)*(u+v))^2-4*t*(t-1)*u*v)^(-3/2)*marg1*marg2
#if (val<0) val<-0

return(val)
}

eva.prod<-function(x,marginal="student",g=1)
{
if (marginal=="student"){
    d<-1
    vakio<-gamma((g+d)/2)/((g*pi)^(d/2)*gamma(g/2))
    y<-vakio*(1+x^2/g)^(-(g+d)/2)
    val<-prod(y)
}
if (marginal=="student.old"){
    d<-1
    vakio<-gamma((g+d)/2)/((g*pi)^(d/2)*gamma(g/2))
    x1<-vakio*(1+x[1]^2/g)^(-(g+d)/2)
    x2<-vakio*(1+x[2]^2/g)^(-(g+d)/2)
    val<-x1*x2
}
if (marginal=="studentR"){
    #x1<-dt(x[1],df=g)
    #x2<-dt(x[2],df=g)
    y<-dt(x,df=g)
    val<-prod(y)  
}
if (marginal=="polyno.old"){
    vakio<-2*(1-1/(g+1))
    y<-vakio*abs(1-x)^g
    val<-prod(y)
}
if (marginal=="polyno"){ 
    vakio<-1/(2*(1-1/(g+1)))
    y<-vakio*abs(1-abs(x)^g)
    val<-prod(y)
}
if (marginal=="double"){
    vakio<-1/2
    y<-exp(-abs(x))
    val<-prod(y)
}
if (marginal=="gauss"){
    vakio<-(2*pi)^(-1/2)
    y<-exp(-x^2/2)
    val<-prod(y)
}

  
return(val)
}

eva.skewgauss<-function(x,mu,sig,alpha)
{

norvak<-prod(sig)^(-1)
point<-(x-mu)/sig
en<-evanor(point)     #dnorm(poi)

point2<-alpha%*%((x-mu)/sig)
pn<-pnorm(point2)

ans<-2*en*pn

return(ans)
}
eva.student<-function(x,t=rep(4,length(x)),
marginal="unif",sig=c(1,1),r=0,df=1)
# t>2 
#  sig is std of marginals
{
if (marginal=="unif"){
   u<-x[1]/sig[1]
   v<-x[2]/sig[2]
   marg1<-1/sig[1]
   marg2<-1/sig[2]
}
if (marginal=="normal"){
   u<-pnorm(x[1]/sig[1])
   v<-pnorm(x[2]/sig[2])
   marg1<-evanor(x[1]/sig[1])/sig[1]
   marg2<-evanor(x[2]/sig[2])/sig[2]
}
if (marginal=="student"){
   u<-pt(x[1]/sig[1],df=t[1])
   v<-pt(x[2]/sig[2],df=t[2])
   marg1<-dt(x[1]/sig[1],df=t[1])/sig[1]
   marg2<-dt(x[2]/sig[2],df=t[2])/sig[2]
}

d<-2
x1<-qt(u,df=df)
x2<-qt(v,df=df)
produ<-dt(x1,df=df)*dt(x2,df=df)

nelio<-(x1^2+x2^2-2*r*x1*x2)/(1-r^2)
vakio<-gamma((df+d)/2)/((df*pi)^(d/2)*gamma(df/2))
ga<-vakio*(1-r^2)^(-1/2)*(1+nelio/df)^(-(df+d)/2)
#ga<-(1-r^2)^(1/2)*(1+(x1^2+x2^2-2*r*x1*x2)/(t*(1-r^2)))^(-(t+d)/2)

val<-ga/produ*marg1*marg2

return(val)
}

eva.t<-function(x,df,mu,Sigma)
{
norma<-gamma((df+1)/2)/((df*pi)^(1/2)*gamma(df/2))
fu<-(1+x^2/df)^(-(df+1)/2)
val<-norma*fu

return(val)
}


excmas.bylevel<-function(lst,levnum)
{
#source("~/denpro/R/excmas.bylevel.R")
#excmas.bylevel(lst,20)

levexc<-matrix(0,levnum,1)

maxlev<-max(lst$level)
step<-maxlev/levnum
nodelkm<-length(lst$parent)

mlkm<-moodilkm(lst$parent)
modloc<-mlkm$modloc    #pointers to modes
lkm<-mlkm$lkm       

added<-matrix(0,nodelkm,1)  #1 if we have visited this node

i<-1
while (i<=lkm){
    node<-modloc[i]
    # calculate curexc
    par<-lst$parent[node]
    if (par==0) valpar<-0 else valpar<-lst$level[par] 
    curexc<-(lst$level[node]-valpar)*lst$volume[node]
    
    nodelevind<-min(max(round(lst$level[node]/step),1),levnum)    
    levexc[1:nodelevind]<-levexc[1:nodelevind]+curexc

    while (lst$parent[node]>0){
         node<-lst$parent[node]
         if (added[node]==0){   
           # calculate curexc
           par<-lst$parent[node]
           if (par==0) valpar<-0 else valpar<-lst$level[par] 
           curexc<-(lst$level[node]-valpar)*lst$volume[node] 
           
           nodelevind<-min(max(round(lst$level[node]/step),1),levnum)    
           levexc[1:nodelevind]<-levexc[1:nodelevind]+curexc

           added[node]<-1 
         }
    }
    i<-i+1
}

levexc<-levexc/levexc[1]

diffe<-matrix(0,length(levexc),1)
for (i in 1:(length(levexc)-1)) diffe[i]<-(levexc[i+1]-levexc[i])/step
diffe[length(diffe)]<-diffe[length(diffe)-1]

return(list(levexc=levexc,diffe=diffe))
}




excmas<-function(lst){
#
parents<-lst$parent
volumes<-lst$volume
levels<-lst$level
#
nodelkm<-length(parents)
excmasses<-matrix(1,nodelkm,1)
#
mlkm<-moodilkm(parents)
modloc<-mlkm$modloc    #pointers to modes
lkm<-mlkm$lkm       
#
added<-matrix(0,nodelkm,1)  #1 if we have visited this node
#
for (i in 1:lkm){
    node<-modloc[i]
    # calculate curexc
    par<-parents[node]
    if (par==0) valpar<-0 else valpar<-levels[par] 
    curexc<-(levels[node]-valpar)*volumes[node]
    #
    excmasses[node]<-curexc
    while (parents[node]>0){
         node<-parents[node]
         if (added[node]==0){   
           # calculate curexc
           par<-parents[node]
           if (par==0) valpar<-0 else valpar<-levels[par] 
           curexc<-curexc+(levels[node]-valpar)*volumes[node] 
           #
           excmasses[node]<-curexc 
           added[node]<-1 
         }
         else{   #add only previous cumulative 
            excmasses[node]<-excmasses[node]+curexc
         }
    }
}
return(t(excmasses))
}




exmap<-function(estiseq,mt,hind=NULL,hseq=NULL,
n=NULL,moteslist=NULL,ylist=NULL)
{
#moteslist is a list of alpha values for every node
#not just for the branch nodes, but it might be nonsense for others

pk<-estiseq$lstseq
if (is.null(hseq)) hseq<-mt$hseq
if (is.null(hind)) hind<-c(1:length(mt$hseq))
slis<-mt$hseq[hind]

d<-dim(pk[[1]]$center)[1]

if (is.null(ylist)) ylist<-c(length(slis):1)

hrun<-1
for (i in 1:length(slis)){
   while (hseq[hrun]!=slis[i]){
      hrun<-hrun+1
   }
   if (i==1) treelist<-list(pk[[hrun]])
   else  treelist=c(treelist,list(pk[[hrun]]))
}

parnum<-length