2008-06-11

等分する(3)〜組み合わせて〜

幾何立体角球

この記事の続き。
n等分したいと同時にm等分したい。 $cos\theta =\frac{-1}{n-1}$ にしつつ、 $cos\theta =\frac{-1}{m-1}$ にするということ。
$(n-1) \times (m-1)$ 次元のうち、 $n-1$ 次元の方で $cos\theta =\frac{-1}{n-1}$ とし、 $m-1$ 次元の方で $cos\theta =\frac{-1}{m-1}$ にする。 $(n-1) \times (m-1)$ 次元のうちの、低次元分割にて、「均等割り」が達成されており、両者は相互に独立している。
これらを連結したときに、「均質」にドッキングするとどうなるか。 $cos\theta=\frac{-1}{(n-1)(m-1)}$ となるらしい。
さらに一般化して、 $\prod n_i$ , $i=1,2,...k$ 空間にて、 $n_i$ 均等分割しつつ、それらを、相互に独立した低次元空間で実現し、「均質」にドッキングするとなると、 $\{1,2,...,k\}$ という集合のべき集合を考えて、その要素(部分集合,subset)、(ちょっと厳密に記載するのはまだるっこしいのだが)、 $cos\theta = \prod_{subset} \frac{-1}{(n_i-1)}$ という関係が、すべての部分集合にて満足するような、そんな分割が存在する(らしい)。もしかすると、 $cos \theta$ は正負の交代などが入り込むかもしれないけれども。。。。

public class TableSpace {

	int na;//No. Axes
	int[] cat;//No. categories per axis
	int df;//Degree of freedom
	int neq;//No. eqdivs
	int[] data;
	int sum;
	double[] exp;
	double chi;
	double[] dif;
	double[][] eq;//coefficient of each equally dividing vectors
	/**
	 * @param args
	 */
	public static void main(String[] args) {
		// TODO 自動生成されたメソッド・スタブ
		int na=3;
		/*
		int[] cat = new int[na];
		for(int i=0;i<cat.length;i++){
			cat[i]=4;
		}
		*/
		int[] cat = {3,4,5};
		
		double[][] testeqdiv=eqdiv(na);
		for(int i=0;i<testeqdiv.length;i++){
			for(int j=0;j<testeqdiv[i].length;j++){
				System.out.print(testeqdiv[i][j]+"\t");
			}
			System.out.println();
		}
		boolean loop=true;
		int[] counter = new int[na];
		while(loop){
			for(int i=0;i<counter.length;i++){
				System.out.print(counter[i]+"\t");
			}
			//System.out.println();
			int index=index(counter,cat);
			System.out.println(index);
			counter=kuriage(counter,cat);
			boolean judge=true;
			for(int i=0;i<counter.length;i++){
				if(counter[i]!=0){
					judge=false;
				}
			}
			if(judge){
				loop=false;
			}
		}
		TableSpace ts=new TableSpace(na,cat);
		for(int i=0;i<ts.eq.length;i++){
			for(int j=0;j<ts.eq[i].length;j++){
				System.out.print(ts.eq[i][j]+"\t");
			}
			System.out.println();
		}
		System.out.println("Eqdiv norm");
		for(int i=0;i<ts.eq.length;i++){
			double norm=0;
			for(int j=0;j<ts.eq[i].length;j++){
			norm+=ts.eq[i][j]*ts.eq[i][j];	
			}
			System.out.println("norm\t"+i+"\t"+norm);
		}
		System.out.println("Sum lattice vector norm");
		for(int i=0;i<ts.eq[0].length;i++){
			double sum=0;
			for(int j=0;j<ts.eq.length;j++){
			sum+=ts.eq[j][i];	
			}
			System.out.println("lattice sum\t"+i+"\t"+sum);
		}
		
		System.out.println("inner product");
		for(int i=0;i<ts.eq.length;i++){
			System.out.print("i=\t"+i+"\t");
			for(int j=0;j<ts.eq.length;j++){
				if(j<=i){
					System.out.print("\t");
				}else{
					double innerProduct=0;
					for(int k=0;k<ts.eq[i].length;k++){
						innerProduct+=ts.eq[i][k]*ts.eq[j][k];
					}
					System.out.print(innerProduct+"\t");
				}
				
			}
			System.out.println();
		}
		

	}
	public TableSpace(int na_,int[] cat_){
		na=na_;
		cat=StatUtilsZ.MiscUtilX.DeepCopyInt1(cat_);
		int[] cat2=new int[cat.length];
		for(int i=0;i<cat2.length;i++){
			cat2[i]=cat[i]-1;
		}
		df=1;
		neq=1;
		for(int i=0;i<cat.length;i++){
			df*=(cat[i]-1);
			neq*=cat[i];
		}
		eq=new double[neq][neq];
		for(int i=0;i<eq.length;i++){
			for(int j=0;j<eq[i].length;j++){
				//eq[i][j]=1;
				eq[i][j]=0;
			}
		}
		
		for(int k=0;k<na;k++){
			int[] counter = new int[na];
			double[][] eqdiv=eqdiv(cat[k]);
			//for(int i=0;i<eqdiv.length;i++){
				//for(int j=0;j<eqdiv[i].length;j++){
					//System.out.print(eqdiv[i][j]+"\t");
				//}
				//System.out.println();
			//}
			boolean loop=true;
			int index=0;
			//System.out.println("-----------");
			//for(int i=0;i<eq.length;i++){
				//for(int j=0;j<eq[i].length;j++){
					//System.out.print(eq[i][j]+"\t");
				//}
				//System.out.println();
			//}
			//System.out.println("-----------");
			while(loop){
				
				
				//System.out.println("index!Pre="+index);
				//for(int i=0;i<eq[index].length;i++){
					//System.out.print(eq[index][i]+"\t");
				//}
				//System.out.println();
				boolean loop2=true;
				int[] counter2 = new int[na];
				int index2=0;
				while(loop2){
					//for(int i=0;i<eq[index].length;i++){
						//for(int i=0;i<eqdiv[0].length;i++){
							//int[] tmpcounter=StatUtilsZ.MiscUtilX.DeepCopyInt1(counter);
							//tmpcounter[k]=i;
							//int tmpindex=index(tmpcounter,cat);
							//int tmpindex=i;
							//int tmpcat=
							//System.out.println("tmpindex="+tmpindex);
							//System.out.println("counter[k]="+counter[k]);
							//System.out.println("i="+i);
							//System.out.println("eq[index][tmpindex]="+eq[index][tmpindex]);
							//eq[index][tmpindex]*=eqdiv[counter[k]][i];
					//System.out.println("index="+index);
					//System.out.println("index2="+index2);
					//System.out.println("counter k ="+counter[k]);
					//System.out.println("counter2 k="+counter2[k]);
					if(counter2[k]<cat[k]-1 && counter2[counter2.length-1]<cat[na-1]-1){
						if(k==0){
							eq[index][index2]=eqdiv[counter[k]][counter2[k]];
									
						}else{
									//if(eq[index][tmpindex]==0){
										//eq[index][tmpindex]=eqdiv[counter[k]][i];
									//}else{
								eq[index][index2]*=eqdiv[counter[k]][counter2[k]];
									//}
						}
					}
					
							
							
							//System.out.println("eq="+eq[index][tmpindex]);
						//}
					counter2=kuriage(counter2,cat);
					boolean judge2=true;
					for(int i=0;i<counter2.length;i++){
						if(counter2[i]!=0){
							judge2=false;
						}
					}
					if(judge2){
						loop2=false;
					}
					index2++;
					if(k==0){
						//loop=false;
					}
				}
				
				//System.out.println("index!Post="+index);
				//for(int i=0;i<eq[index].length;i++){
					//System.out.print(eq[index][i]+"\t");
				//}
				//System.out.println();
				counter=kuriage(counter,cat);
				boolean judge=true;
				for(int i=0;i<counter.length;i++){
					if(counter[i]!=0){
						judge=false;
					}
				}
				if(judge){
					loop=false;
				}
				index++;
				if(k==0){
					//loop=false;
				}
			}
			/*
			for(int i=0;i<eq.length;i++){
				for(int j=0;j<eq[i].length;j++){
					eq[i][j]=1;
				}
			}
			*/
			
		}
		for(int i=0;i<eq.length;i++){
			for(int j=0;j<eq[i].length;j++){
				eq[i][j]/=(Math.pow(2,(double)(na-2)/2));
			}
		}
	}
	void addData(int[] d){
		data=new int[eq.length];
		for(int i=0;i<data.length;i++){
			data[i]=d[i];
			sum+=d[i];
		}
	}
	void chi(){
		for(int i=0;i<data.length;i++){
			chi+=(data[i]-exp[i])*(data[i]-exp[i])/exp[i];
		}
	}
	void calcExp(){
		exp=new double[eq.length];
		int[][] marginal = new int[na][0];
		for(int i=0;i<marginal.length;i++){
			marginal[i]=new int[cat[i]];
			
		}
		boolean loop2=true;
		int[] counter2 = new int[na];
		int index2=0;
		while(loop2){
			for(int i=0;i<marginal.length;i++){
				marginal[i][counter2[i]]+=data[index2];
			}
			
					
					
					//System.out.println("eq="+eq[index][tmpindex]);
				//}
			counter2=kuriage(counter2,cat);
			boolean judge2=true;
			for(int i=0;i<counter2.length;i++){
				if(counter2[i]!=0){
					judge2=false;
				}
			}
			if(judge2){
				loop2=false;
			}
			index2++;
			
		}
		loop2=true;
		counter2 = new int[na];
		index2=0;
		while(loop2){
			//exp[index2]=(double)data[index2];
			exp[index2]=sum;
			for(int i=0;i<marginal.length;i++){
				exp[index2]*=(double)(marginal[i][counter2[i]])/(double)sum;
				
				
			}
			
				
					
					//System.out.println("eq="+eq[index][tmpindex]);
				//}
			counter2=kuriage(counter2,cat);
			boolean judge2=true;
			for(int i=0;i<counter2.length;i++){
				if(counter2[i]!=0){
					judge2=false;
				}
			}
			if(judge2){
				loop2=false;
			}
			index2++;
			
		}
	}
	public static int[] kuriage(int[] a,int[] keta){
		int[] ret = StatUtilsZ.MiscUtilX.DeepCopyInt1(a);
		ret[0]++;
		for(int i=0;i<a.length;i++){
			if(ret[i]==keta[i]){
				ret[i]=0;
				if(i<a.length-1){
					ret[i+1]++;
				}
				
			}
		}
		
		return ret;
	}
	public static int index(int[] counter,int[] cat){
		int kurai=1;
		int index=0;
		for(int i=0;i<counter.length;i++){
			if(i>0){
				kurai*=cat[i-1];
			}
			
			index+=counter[i]*kurai;
		}
		return index;
	}
	public static double[][] eqdiv(int n){
		double[][] ret = new double[n][n-1];
		for(int i=0;i<ret.length;i++){
			if(i<n-1){
				ret[i][i]=Math.sqrt((double)(n)/(double)(n-1))*
				Math.sqrt((double)(n-i-1)/(double)(n-i));
			}
			
			for(int j=0;j<i;j++){
				ret[i][j]=-Math.sqrt((double)(n)/(double)(n-1))*
				Math.sqrt(1/(double)((n-1-j)*(n-j)));
			}
		}
		return ret;
	}

}