2009-12-08

高次元分割表

R 分割分割表

昨日、N個のものをk種類に分けるわけ方について、各種類に上限数があるときに、何通りの分け方があるかについて書いた。
少し、変える。
N個のものがあって、第1の因子では、あるわけ方をされ、第2の因子では、別の分け方をされ、。。。と高次に分けられているときの、分け方について考える。まずは単純にするために、それぞれの因子では、2種類に分けられるものとする。すると、k因子があるときにはその因子の持ち方が $2^k$ 通り存在する。これらの持ち方には、上限があるが、この上限は、k個の因子の2種類の数に依存する。このような条件で、さて、何通りあるか・・・

Rのコマンド
- 任意次元2^k分割。総数Ｎにおいて、第１因子で(n1,N-n1) に分かれていたとする。その後、第２因子で(n2,N-n2),第３因子で(n3,N-n3),..,第i因子で(ni,N-ni)に分けたときに、1,...,k要素が(0,...,0)(0,..,0,1),(0,,,1,0),...,(1,1,...,1)の2^iパターンがあるが、その数の割り振り方がいくつあるかを算出。重いです
- 全サンプル数が15で、因子数が３つ、第1因子が(10,5),第2因子が（3,12)、第3因子が(8,7)であるときに、次のようにする。depthは２から変えないこと。silent=FALSEにすると、場合数だけでなく、分配パターンも表示される
- PatternsHighDim(nlist=c(10,5),k=c(10,30,8),depth=2,silent=TRUE)
Rのソース

PatternsHighDim<-function(nlist=c(10,5),k=c(10,3,8),depth=2,silent=TRUE){
 #print("nlist,depth,val,count")
 #print(nlist)
 #print(depth)
 #print(val)
 #print(count)
 #print("------")
 if(depth>length(k)){
  if(!silent){
    print(nlist)
  }
  return(list(count=1))
 }
 #write.table(res,file="TEXTFILE1.txt",quote=F,sep="\t",row.names=F,col.names=F)
 #sink("TEXTFILE1.txt")
 lenNlist<-length(nlist)

 MinMaxList<-LimitRange(nlist,k[depth])
 #newNlist<-rep(0,lenNlist*2)
 N<-sum(nlist)
 patterncounter<-0
 plusQ<-InitPlusMinMax(MinMaxList$MinList,MinMaxList$MaxList,k[depth])
 loop=plusQ$judge
 while(loop!=0){
  #print("plusQplus")
  #print(plusQ$plus)
  if(plusQ$judge==1){
   newNlist<-UpDateNlist(nlist,plusQ$plus)
   
   #print(newNlist)
   #newval<-lExp(newNlist,LogFact)
   #print("newval")
   #print(newval)
   if(FALSE){
    #return(newval)
    #val<-val+newval
    #patterncounter<-patterncounter+1
   }else{
    ans<-PatternsHighDim(newNlist,k,depth+1,silent)
    #val<-val+ans$val
    patterncounter<-patterncounter+ans$count

   }
  }
  plusQ<-UpdatePlusMinMax(plusQ$plus,MinMaxList$MinList,MinMaxList$MaxList,k[depth])
  loop<-plusQ$judge

 }
 return(list(count=patterncounter))
}
Patterns<-function(nlist=c(10,20,6,8),k=20,silent=TRUE){
 #write.table(res,file="TEXTFILE1.txt",quote=F,sep="\t",row.names=F,col.names=F)
 #sink("TEXTFILE1.txt")
 lenNlist<-length(nlist)

 MinMaxList<-LimitRange(nlist,k)
 #newNlist<-rep(0,lenNlist*2)
 N<-sum(nlist)
 patterncounter<-0
 plusQ<-InitPlusMinMax(MinMaxList$MinList,MinMaxList$MaxList,k)
 loop=plusQ$judge
 while(loop!=0){
  #newNlist<-UpDateNlist(nlist,plusQ$plus)
   #print(plusQ$plus)
  if(plusQ$judge==1){
   if(!silent){
    print(plusQ$plus)
   }
  
   patterncounter<-patterncounter+1
  }
  

  plusQ<-UpdatePlusMinMax(plusQ$plus,MinMaxList$MinList,MinMaxList$MaxList,k)
  loop<-plusQ$judge
  
 }
 #sink()
 return(list(NumPatterns=patterncounter))
 
}
InitPlusMinMax<-function(MinList,MaxList,k){
 plus<-MinList
 plus[length(plus)]<-k-sum(plus)+plus[length(plus)]
 if(plus[length(plus)]<=MaxList[length(plus)] && plus[length(plus)]>=MinList[length(plus)]){
  return(list(plus=plus,judge=1))
 }else{
  ans<-UpdatePlusMinMax(plus,MinList,MaxList,k)
  return(list(plus=ans$plus,judge=ans$judge))
 }
}
UpdatePlusMinMax<-function(plus,MinList,MaxList,k){
 #print("update")
 plusLen<-length(plus)
 retplus<-plus
 success<-1
 while(success==1){
  retplus[1]<-retplus[1]+1
  retplus[plusLen]<-retplus[plusLen]-1
  for(i in 1:(plusLen-1)){
   
   if(retplus[i]<=MaxList[i]){
    #if(retplus[plusLen]<MinList[plusLen]){
    # loop<-0
    #}else if(retplus[plusLen]>MaxList[plusLen]){
    # retplus[i]<-retplus[i]+retplus[plusLen]-MaxList[plusLen]
    # retplus[plusLen]<-MaxList[plusLen]
    # if(retplus[i]<=MaxList[i] && retplus[i]>=MinList[i]){
    #  return(list(plus=retplus,judge=1))
    # }else{
    #  loop<-0
    # }
    #}else 
    if(retplus[plusLen]<=MaxList[plusLen] && retplus[plusLen]>=MinList[plusLen]){
     return(list(plus=retplus,judge=1))
    }else{
     return(list(plus=retplus,judge=2))
    }
    #retplus[i]<-retplus[i]+1
    #retplus[plusLen]<-retplus[plusLen]-1
   }
   retplus[plusLen]<-retplus[plusLen]+retplus[i]-MinList[i]-1
   retplus[i]<-MinList[i]
   retplus[i+1]<-retplus[i+1]+1
   if(i==(plusLen-1)){
    #if(retplus[i+1]>MaxList[i+1]){
     return(list(plus=retplus,judge=0))
    #}
   }
  }
 }
 
 
 #return(list(plus=retplus,judge=0))
}
LimitRangeKeta<-function(nlist,k,keta,loopMax=10){
 MinList<-rep(0,keta)
 MaxList<-nlist[1:keta]
 dif<-1
 loop<-0
 while(dif!=0 && loop<loopMax){
  loop<-loop+1
  tmpMinList<-UpdateMinList(MinList,MaxList,k)
  tmpMaxList<-UpdateMaxList(MinList,MaxList,k)
  dif<-sum(abs(tmpMinList-MinList)+abs(tmpMaxList-MaxList))
  MinList<-tmpMinList
  MaxList<-tmpMaxList
 }
 return(list(MinList=MinList,MaxList=MaxList))

}
LimitRange<-function(nlist,k,loopMax=10){
 MinList<-rep(0,length(nlist))
 MaxList<-nlist
 LEN<-length(nlist)
 dif<-1
 loop<-0
 while(dif!=0 && loop<loopMax){
  #print("#")
  loop<-loop+1
  #print(loop)
  tmpMinList<-UpdateMinList(MinList,MaxList,k)
  tmpMaxList<-UpdateMaxList(MinList,MaxList,k)
  #print(tmpMinList)
  #print(tmpMaxList)
  dif<-sum(abs(tmpMinList-MinList)+abs(tmpMaxList-MaxList))
  MinList<-tmpMinList
  MaxList<-tmpMaxList
 }
 return(list(MinList=MinList,MaxList=MaxList))
}



UpdateMinList<-function(MinList,MaxList,k){
 list<-MinList
 sumMax<-sum(MaxList)
 for(i in 1:(length(list))){
  if(MinList[i]<k-sumMax+MaxList[i]){
   list[i]<-k-sumMax+MaxList[i]
  }
 }
 return(list)
}
UpdateMaxList<-function(MinList,MaxList,k){
 list<-MaxList
 sumMin<-sum(MinList)
 for(i in 1:(length(list))){
  if(MaxList[i]>k-sumMin+MinList[i]){
   list[i]<-k-sumMin+MinList[i]
  }
 }
 return(list)
}

InitPlus<-function(nlist,k){
 retplus<-rep(0,length(nlist))
 sofar<-k
 for(i in 1:length(nlist)){
  if(sofar>=nlist[i]){
   retplus[i]<-nlist[i]
   sofar<-sofar-nlist[i]
  }else{
   retplus[i]<-sofar
   return(retplus)
  }
 }
 return(retplus)
}
UpdatePlus<-function(plus,nlist){
 plusLen<-length(plus)
 totalplus<-sum(plus)
 plus[1]<-plus[1]+1
# if(plusLen>2){
  for(i in 1:plusLen){
   #print(plus)
   if(plus[i]>nlist[i]){
    #plus[plusLen]<-plus[plusLen]+plus[i]-1
    plus[i]<-0
    if(i<plusLen){
     plus[i+1]<-plus[i+1]+1
    }else{
     plus[i]<-plus[i]+1
     return(plus)
    }
   }else{
    return(plus)
   }
  }
 #}
 return(plus)
}
UpDateNlist<-function(nlist,plus){
 newNlist<-rep(0,length(nlist)*2)
 for(i in 1:length(nlist)){
  newNlist[i*2-1]<-plus[i]
  newNlist[i*2]<-nlist[i]-plus[i]
 }
 return(newNlist)
}
lExp<-function(nlist,LogFact){
 ret<-0
 for(i in 1:length(nlist)){
  ret<-ret+LogFact[nlist[i]+1]
 }
 return(ret)
}
lExpHighDim<-function(nlist,LogFact,k,depth){
 ret<-0
 N<-sum(nlist)
 sortlist<-sort(nlist)
 for(i in 1:length(nlist)){
  ret<-ret-LogFact[sortlist[i]+1]
 }
 for(i in 1:depth){
  ret<-ret+LogFact[k[i]+1]+LogFact[N-k[i]+1]
 }
 ret<-ret-(depth-1)*LogFact[N+1]
 return(ret)
}

MakeLogFact<-function(n){
 ret<-rep(0,n+1)
 ret[1]<-0
 for(i in 2:(n+1)){
  ret[i]<-ret[i-1]+log(i-1)
 }
 return(ret)
}