ソース

# 依存パッケージ
library(kinship)
library(MCMCpack)
library(gtools)
library(sets)
library(paramlink)

# 家系情報からkinshipパッケージのpedigreeオブジェクトを作る
# p 
# 第1カラム:ID
# 第2カラム:母
# 第3カラム:父
# 第4カラム:性別 0:男 1:女
# 第5カラム:タイプ (1はジェノタイプ提供者、2は探されている人、3は情報なしの人)
MakePedigreeFromFamilyInfo<-function(p){
	ns<-length(p[,1])
	affected<-status<-rep(1,ns)
	affected[which(p[,5]==2)]<-0
	affected[which(p[,5]==3)]<-0
	status[which(p[,5]==1)]<-0
	status[which(p[,5]==2)]<-0
	ptemp<-pedigree(id=p[,1],dadid=p[,3],momid=p[,2],sex=p[,4],affected=affected,status=status)
	if(sum(ptemp$affected)==0)ptemp$affected<-affected
	ptemp
}

# NL個のローカスに、リピート数アレルとその確率を適当に与える

MakeAlleleProb<-function(NL=15){

	Alleles<-NULL
	Probs<-NULL
	library(MCMCpack)
	for(i in 1:NL){
		na<-sample(2:5,1)
		st<-sample(10:20,1)
		tmp<-NULL  
		for(j in 1:na){
			a<-paste("",st+j,sep="")
			tmp<-c(tmp,a)
		}
		Alleles[[i]]<-tmp
		Probs[[i]]<-c(rdirichlet(1,rep(1,na)))
		#print(sum(Probs[[i]]))
	}
	list(Alleles=Alleles,Probs=Probs)
}
# ランダムにジェノタイプを発生しよう
RandomGenotype<-function(N,n,Alleles,Probs){
	ns<-N
	g<-array(0,c(ns,2,n))
	noparents<-1:ns
	for(i in noparents){
		for(j in 1:n){
			g[i,,j]<-sample(Alleles[[j]],2,replace=TRUE,Probs[[j]])
		}
	}
	g
}
# 家系図を見ながら、ジェノタイプを割り当てよう
# missing =TRUEの場合は、サンプル提供者のみジェノタイプあり
RandomGenotypeFamily<-function(d,n,Alleles,Probs,missing=FALSE){
	ns<-length(d[,1])
	g<-array(0,c(ns,2,n))
	noparents<-which(d[,2]==0,d[,3]==0)
	for(i in noparents){
		for(j in 1:n){
			g[i,,j]<-sample(Alleles[[j]],2,replace=TRUE,Probs[[j]])
		}
	}
	given<-rep(0,ns)
	given[noparents]<-1
	loop<-TRUE
	if(sum(given)==ns)loop<-FALSE
	while(loop){
		yet<-which(given==0)
		for(i in yet){
			prt1<-d[i,2]
			prt2<-d[i,3]
			if(given[prt1]+given[prt2]==2){
				for(j in 1:n){
					g[i,1,j]<-sample(g[prt1,,j],1)
					g[i,2,j]<-sample(g[prt2,,j],1)
				}
				given[i]<-1
			}
		}
		if(sum(given)==ns)loop<-FALSE
	}
	if(missing){
		Type<-p[,5]
		G2<-array(0,c(length(p[,1]),2,NL))
		
		G2[which(Type==1),,]<-g[which(Type==1),,]
		g<-G2
	}
	g
}
# 調べたい個人が何人か、いる
# そのうち、1人ずつだけをサンプルプールに対応づけるか
# 2,3,...人ずつ対応付けるかを考慮する必要がある
# そのパターンを決める

MakeSearchPattern<-function(v){
	library(gtools)
	nt<-length(v)
	x<-0
	X<-NULL
	for(i in 1:nt){
		tmp<-combinations(nt,i)
		tmp2<-tmp
		tmp2[1:length(tmp2)]<-c(v[tmp])
		#X[[i]]<-matrix(tmp2,nrow=length(tmp[,1]))
		X[[i]]<-tmp2
		x<-sum(x,length(as.matrix(X[[i]])[,1]))
	}
	list(X=X,x=x)
}

# ある候補者が複数から指名されていたら、テストしない
MakeTestPattern<-function(candidates){
	gr<-expand.grid(candidates)
	n<-length(gr[,1])
	ok<-rep(1,n)
	for(i in 1:n){
		tmp<-unlist(gr[i,])
		tmp2<-outer(tmp,tmp,"-")
		diag(tmp2)<-1
		if(prod(tmp2)==0){
			ok[i]<-0
		}
	}
	as.matrix(gr[which(ok==1),])
}

OffspringGenotypeProb<-function(p1,p2){
	V<-p1%*%t(p2)
	U<-V+t(V)
	U[lower.tri(U)]<-0
	diag(U)<-diag(U)/2
	U
}

LogLikelihoodTrio<-function(p1,p2,c){
	V<-p1%*%t(p2)
	U<-V+t(V)
	U[lower.tri(U)]<-0
	diag(U)<-diag(U)/2
}


CalcLogLikelihoodFamily<-function(p,g,Alleles,Probs){
	#print(g)
	nl<-length(g[1,1,])
	ret<-rep(0,nl)
	
	for(i in 1:length(g[1,1,])){
		tmp<-CalcLogLikelihoodFamilyLocus(p,g[,,i],Alleles[[i]],Probs[[i]])
		ret[i]<-tmp$loglikesum
	}
	list(loglikelist=ret,loglikesum=sum(ret))
}
LimitDiplotypeZ<-function(p,g,A){
	unknown="0"
	ns<-length(p[,1])
	na<-length(A)
	#D ディプロタイプ
	#H ハプロタイプ
	#0 不明 1必ず持つ -1決して持たない
	D<-matrix(0,ns,na)
	H<-array(0,c(ns,2,na))
		# 確定ディプロタイプで代入
		# 確定ディプロタイプが示すアレルを持つことは確定し
		# 持たないアレルを持たないことも確定する
		# 持たないアレルが伝達することは決してない
		for(i in 1:ns){
			if(g[i,1]!=unknown){
				tmp<-rep(0,na)
				tmp[which(A==g[i,1])]<-1
				tmp[which(A==g[i,2])]<-1
				D[i,which(tmp==1)]<-1
				D[i,which(tmp==0)]<-(-1)
				H[i,1,which(tmp==0)]<-(-1)
				H[i,2,which(tmp==0)]<-(-1)
			}
		}

	# ループ
	loop<-TRUE
	while(loop){
		tmpH<-H
		tmpD<-D
		# 親子トリオ
		for(i in 1:ns){
			# 由来親の持っていないアレルは受け取れない
			# 由来親から受け取れないアレルの数がアレル総数-1に等しいとき、
			# そのアレルは受け取るアレルとして確定する
			if(p[i,2]!=0){# mother known
				H[i,1,which(D[p[i,2],]==-1)]<-(-1)
				if(length(which(D[p[i,2],]>=0))==1){
					#print(i)
					#print(H[i,1,])
					H[i,1,which(D[p[i,2],]==1)]<-1
					H[i,1,which(D[p[i,2],]!=1)]<-(-1)
					#print(H[i,1,])
				}
			}
			if(length(which(H[i,1,]==-1))==(na-1))H[i,1,which(H[i,1,]>(-1))]<-1
			if(length(which(H[i,1,]==1))==1)H[i,1,which(H[i,1,]!=1)]<--1
			D[i,which(H[i,1,]==1)]<-1
			if(p[i,3]!=0){ # father known
				H[i,2,which(D[p[i,3],]==-1)]<-(-1)
				if(length(which(D[p[i,3],]>=0))==1){
					H[i,2,which(D[p[i,3],]==1)]<-1
					H[i,2,which(D[p[i,3],]!=1)]<-(-1)
				}
			}
			if(length(which(H[i,2,]==-1))==(na-1))H[i,2,which(H[i,2,]>(-1))]<-1
			if(length(which(H[i,2,]==1))==1)H[i,2,which(H[i,2,]!=1)]<--1
			D[i,which(H[i,2,]==1)]<-1

			# 子の必ず持つアレルのうち、父(母)由来でないアレルは、必ず母(父)からの伝達である
			H[i,1,which(((D[i,]==1) * (H[i,2,]==-1))==1)]<-1
			H[i,2,which(((D[i,]==1) * (H[i,1,]==-1))==1)]<-1

			D[p[i,2],which(H[i,1,]==1)]<-1
			D[p[i,3],which(H[i,2,]==1)]<-1
			if(length(which(D[i,]==1))==2)D[i,which(D[i,]!=1)]<-(-1)
			
			
			# 両親から伝達される可能性のないアレルはディプロタイプに持ちえない
			for(i in 1:ns){
				D[i,which((H[i,1,]==-1) * (H[i,2,]==-1)==1)]<-(-1)
				if(length(which(D[i,]==-1))>=(na-1))D[i,which(D[i,]>(-1))]<-1
			}
		}
		if(isTRUE(all.equal(D,tmpD)) & isTRUE(all.equal(H,tmpH)))loop<-FALSE
		
	}
	dset<-dsetA<-list()
	dtypeset<-dtypesetA<-list()
	hset1<-hset2<-hset1A<-hset2A<-list()
	for(i in 1:length(D[,1])){
		tmpd1<-as.set(which(D[i,]==1))
		tmpd0<-as.set(which(D[i,]==0))
		tmpd1A<-as.set(A[which(D[i,]==1)])
		tmpd0A<-as.set(A[which(D[i,]==0)])
		if(set_is_empty(tmpd1)){
			#dset[[i]]<-tmpd0*tmpd0
			dset[[i]]<-set_combn(tmpd0,1) | set_combn(tmpd0,2)
			#dsetA[[i]]<-tmpd0A*tmpd0A
			dsetA[[i]]<-set_combn(tmpd0A,1) | set_combn(tmpd0A,2)

		}else{
			if(set_is_empty(tmpd0)){
				dset[[i]]<-set(tmpd1)
				dsetA[[i]]<-set(tmpd1A)
			}else{
				#dset[[i]]<-tmpd1*tmpd0
				dset[[i]]<-dsetA[[i]]<-set()
				tmpd01<-tmpd0 | tmpd1
				for(j in tmpd01){
					dset[[i]]<-dset[[i]] | set(tmpd1 | j)

				}
				tmpd01A<-tmpd0A | tmpd1A
				for(j in tmpd01A){
					dsetA[[i]]<-dsetA[[i]] | set(tmpd1A | j)
				}
				#dset[[i]]<-as.set(sapply(set_combn(tmpd0,1),"|",tmpd1))
				#dsetA[[i]]<-tmpd1A*tmpd0A
				#dsetA[[i]]<-as.set(sapply(set_combn(tmpd0A,1),"|",tmpd1A))
			}
			
		}
		#dset[[i]]<-as.set(which(D[i,]>=0))
		hset1[[i]]<-as.set(which(H[i,1,]>=0))
		hset2[[i]]<-as.set(which(H[i,2,]>=0))
		dtypeset[[i]]<-as.set(hset1[[i]]*hset2[[i]])
		#dsetA[[i]]<-as.set(A[which(D[i,]>=0)])
		hset1A[[i]]<-as.set(A[which(H[i,1,]>=0)])
		hset2A[[i]]<-as.set(A[which(H[i,2,]>=0)])
		dtypesetA[[i]]<-as.set(hset1A[[i]]*hset2A[[i]])
	}
	# 子が持つアレルを1本も持たないディプロタイプを親は持たない
	for(i in 1:length(D[,1])){
		tmpdset<-set()

		for(j in dset[[i]]){
			offsOK<-TRUE
			oyajudge<-FALSE

			oyanashi<-TRUE
			for(k in 1:length(p[,1])){
				if(p[k,2] == i | p[k,3] == i){
				offsOKperchild<-FALSE

					oyajudge<-TRUE
					oyanashi<-FALSE
						for(l in dset[[k]]){
							if(!set_is_empty(j & l))offsOKperchild<-TRUE
							#print(j)
							#print(l)
							#print(offsOKperchild)
							#print("---")
						}
						if(!offsOKperchild)offsOK<-FALSE

				}

			}
			hapsOK<-FALSE
			if(!(set_is_empty(j & hset1[[i]])) & !(set_is_empty(j & hset2[[i]])) &
			length(j & (hset1[[i]] | hset2[[i]]) )==length(j) )hapsOK<-TRUE
			if(offsOK & hapsOK){
				#oyajudge=TRUE
				tmpdset<-tmpdset | set(j)
			}
			if(oyanashi & hapsOK){
				oyajudge=TRUE
				tmpdset<-tmpdset | set(j)
			}
		}
		#if(oyajudge){
			dset[[i]]<-tmpdset
		#}
		
	}
	for(i in 1:length(D[,1])){
		tmpdset<-set()

		for(j in dsetA[[i]]){
			offsOK<-TRUE
			oyajudge<-FALSE

			oyanashi<-TRUE
			for(k in 1:length(p[,1])){
				if(p[k,2] == i | p[k,3] == i){
				offsOKperchild<-FALSE

					oyajudge<-TRUE
					oyanashi<-FALSE
						for(l in dsetA[[k]]){
							if(!set_is_empty(j & l))offsOKperchild<-TRUE
							#print(j)
							#print(l)
							#print(offsOKperchild)
							#print("---")
						}
						if(!offsOKperchild)offsOK<-FALSE

				}

			}
			hapsOK<-FALSE
			if(!(set_is_empty(j & hset1A[[i]])) & !(set_is_empty(j & hset2A[[i]])) &
			length(j & (hset1A[[i]] | hset2A[[i]]) )==length(j) )hapsOK<-TRUE
			if(offsOK & hapsOK){
				#oyajudge=TRUE
				tmpdset<-tmpdset | set(j)
			}
			if(oyanashi & hapsOK){
				oyajudge=TRUE
				tmpdset<-tmpdset | set(j)
			}
		}
		#if(oyajudge){
			dsetA[[i]]<-tmpdset
		#}
	}
	dlist<-dlistA<-list()
	
	#for(i in 1:length(dtypeset)){
	for(i in 1:length(dset)){
		cnt<-1
		dlist[[i]]<-list()
		#for(j in dtypeset[[i]]){
		for(j in dset[[i]]){
			dlist[[i]][[cnt]]<-as.set(j)
			cnt<-cnt+1
		}
	}
	#for(i in 1:length(dtypesetA)){
	for(i in 1:length(dsetA)){
		cnt<-1
		dlistA[[i]]<-list()
		#for(j in dtypesetA[[i]]){
		for(j in dsetA[[i]]){
			dlistA[[i]][[cnt]]<-as.set(j)
			cnt<-cnt+1
		}
	}
	list(D=D,H=H,dlist=dlist,dlistA=dlistA,dset=dset,dsetA=dsetA,dtypeset=dtypeset,dtypesetA=dtypesetA,
	hset1=hset1,hset2=hset2,hset1A=hset1A,hset2A=hset2A)
}
TestAllPatterns<-function(p,G,Gpool,candidates,Alleles,Probs){
	ptemp<-MakePedigreeFromFamilyInfo(p)

	tobeSearched<-which(p[,5]==2)
	SearchPattern<-MakeSearchPattern(tobeSearched)
	ret<-NULL
	ret2<-NULL
	ret3<-NULL
	counter<-1
	for(i in 1:length(SearchPattern$X)){
		spx<-as.matrix(SearchPattern$X[[i]])
		
		for(j in 1:length(spx[,1])){
			
			ToTestCurrent<-spx[j,]
			tmpG2<-G
			#print(tmpG2[,,1])
			tobeEvaluated<-rep(0,length(p[,1]))
			tobeEvaluated[which(p[,5]==1)]<-1
			tobeEvaluated[which(p[,5]==3)]<-1
			tobeEvaluated[ToTestCurrent]<-1
			tmpp<-p[which(tobeEvaluated==1),]
			tmptmpG2<-tmpG2[which(tobeEvaluated==1),,]
			if(length(tmpG2[1,1,])==1){
				tmptmpG2<-array(0,c(length(which(tobeEvaluated==1)),2,1))
				tmptmpG2[,,1]<-tmpG2[which(tobeEvaluated==1),,]
			}

			#print(tmptmpG2)

			#print("---")
			#print(tmpG2)
			#print(tobeEvaluated)
			#print(tmpp)
			#print(tmptmpG2)
			#ret[[counter]]<-CalcLogLikelihoodFamily(tmpp,tmptmpG2,Alleles,Probs)
			ret[[counter]]<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs)
			ret2[[counter]]<-list(ToTestCurrent)
			ret3[[counter]]<-list()
			#print(ret[[counter]])
		
			counter<-counter+1
			
			numtested<-length(spx[j,])
			tmpcandidates<-NULL

			for(k in 1:numtested){
				selectedToBeTested<-which(tobeSearched==spx[j,k])
				
				tmpcandidates[[k]]<-candidates[[selectedToBeTested]]
			}
			tmpTestPattern<-MakeTestPattern(tmpcandidates)
			counter<-counter+1
			
			numtested<-length(spx[j,])
			tmpcandidates<-NULL

			for(k in 1:numtested){
				selectedToBeTested<-which(tobeSearched==spx[j,k])
				
				tmpcandidates[[k]]<-candidates[[selectedToBeTested]]
			}
			tmpTestPattern<-MakeTestPattern(tmpcandidates)
			for(k in 1:length(tmpTestPattern[,1])){
				#print("to be tested IDs in pedigree")
				#print(ToTestCurrent)
				#print("to be tested sampleIDs in pool")
				#print(tmpTestPattern[k,])

				tmpG2<-G
				tmpG2[ToTestCurrent,,]<-Gpool[tmpTestPattern[k,],,]
				#print(tmpG2)
				tobeEvaluated<-rep(0,length(p[,1]))
				tobeEvaluated[which(p[,5]==1)]<-1
				tobeEvaluated[which(p[,5]==3)]<-1
				tobeEvaluated[ToTestCurrent]<-1
				tmpp<-p[which(tobeEvaluated==1),]
				tmptmpG2<-tmpG2[which(tobeEvaluated==1),,]
				if(length(tmpG2[1,1,])==1){
					tmptmpG2<-array(0,c(length(which(tobeEvaluated==1)),2,1))
					tmptmpG2[,,1]<-tmpG2[which(tobeEvaluated==1),,]
				}

				#print("---")
				#print(tobeEvaluated)
				#print(tmpp)
				#print(tmptmpG2)
				#ret[[counter]]<-CalcLogLikelihoodFamily(tmpp,tmptmpG2,Alleles,Probs)
				ret[[counter]]<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs)
				
				#ret[[counter]]<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs)
				ret2[[counter]]<-list(ToTestCurrent)
				ret3[[counter]]<-list(tmpTestPattern[k,])

				#print(ret[[counter]])
				counter<-counter+1
			}
		}
	}
	LRmatrix<-matrix(0,counter,counter)
	for(i in 1:(counter-1)){
		for(j in 1:(counter-1)){
			LRmatrix[i,j]<-ret[[i]][[2]][1]-ret[[j]][[2]][1]
		}
	}
	
	list(LogLike=ret,LRmatrix=LRmatrix,SearchedID=ret2,CandidateID=ret3,familyinfo=p,pedigree=ptemp,tobeSearched=tobeSearched,candidates=candidates,offeredGenotype=G,poolGenotype=Gpool,Alleles=Alleles,Probs=Probs)
}
counterplus<-function(x,t){
	x[1]<-x[1]+1
	for(i in 1:(length(x))){
		if(x[i]==t[i]){
			x[i]<-0
			if(i<length(x)){
				x[i+1]<-x[i+1]+1
			}
			
		}else{
			return(x)
		}
	}
	return(x)
}
counterplus2<-function(x,t){
	x[1]<-x[1]+1
	for(i in 1:(length(x))){
		if(x[i]==(t[i]+1)){
			x[i]<-1
			if(i<length(x)){
				x[i+1]<-x[i+1]+1
			}
			
		}else{
			return(x)
		}
	}
	return(x)
}

subnucs <- function(ped) {# subfunction in linkdat() in package "paramlink"
        parents <- unique(ped[, 2:3])
        parents = parents[-match(0, parents[, 1]), , drop = FALSE]
        list1 <- lapply(nrow(parents):1, function(i) {
            par = parents[i, ]
            c(fa = par[[1]], mo = par[[2]], offs = as.vector(ped[, 
                1])[which(ped[, 2] == par[[1]] & ped[, 3] == 
                par[[2]], useNames = FALSE)])
        })
        res = list()
        i = 1
        k = 1
        while (length(list1) > 1) {
            if (i > length(list1)) {
                warning("Loop detected, likelihood calculations will not work.")
                return(NULL)
            }
            link = ((sub <- list1[[i]]) %in% unlist(list1[-i]))
            if (sum(link) == 1) {
                res[[k]] <- c(pivot = sub[[which(link)]], sub)
                list1 <- list1[-i]
                k <- k + 1
                i <- 1
            }
            else i <- i + 1
        }
        res[[k]] <- c(pivot = 0, list1[[1]])
        res
    }
CalcLogLikelihoodFamilyLocus2<-function(p,g,A,P){
	#print(g)
	unknown<-"0"
	Unfixed<-which(g[,1]==unknown)
	Nunfixed<-length(Unfixed) # 不明ジェノタイプ人数
	ns<-length(p[,1])
	na<-length(A)
	if(Nunfixed==0){# 全員のジェノタイプが指定されている場合
	#print("all fixed")
		tmpret<-rep(1,ns)
		currentg<-g
		gvector<-array(0,c(ns,na,2))
		ret<-0
		#given<-rep(0,ns)
		for(i in 1:ns){
			#tmp<-0
			#if(g[i,1]!=unknown){
				gvector[i,which(A==currentg[i,1]),1]<-1
				#tmp<-tmp+1
			#}
			#if(g[i,2]!=unknown){
				gvector[i,which(A==currentg[i,2]),2]<-1
				#tmp<-tmp+1
			#}
			#if(tmp==2)given[i]<-1
		}
		#print(gvector)
		# すべての親子トリオについて計算
		trioloop<-TRUE
		triocnt<-1
		while(trioloop){
			
			#print(i)
			#P1<-P2<-P
			P1<-gvector[triocnt,,1]
			P2<-gvector[triocnt,,2]
			tmp1<-p[triocnt,2]
			tmp2<-p[triocnt,3]
			if(tmp1!=0){
				P1<-(gvector[tmp1,,1]+gvector[tmp1,,2])/2
			}
			if(tmp2!=0){
				P2<-(gvector[tmp2,,1]+gvector[tmp2,,2])/2
			}
			
			FromParents<-OffspringGenotypeProb(P1,P2)
			tmpallele1<-which(A==currentg[triocnt,1])
			tmpallele2<-which(A==currentg[triocnt,2])
			lessa<-min(tmpallele1,tmpallele2)[1]
			morea<-max(tmpallele1,tmpallele2)[1]
			tmptmp<-FromParents[lessa,morea]
			#Offspring<-OffspringGenotypeProb(gvector[triocnt,,1],gvector[triocnt,,2])
			#ret[i]<-log(sum(FromParents*Offspring))
			#tmptmp<-sum(FromParents*Offspring)
			tmpret[triocnt]<-tmptmp
			if(tmptmp==0)trioloop<-FALSE

			#print(tmptmp)
			triocnt<-triocnt+1
			if(triocnt==(ns+1))trioloop<-FALSE
		}
		#tmprob<-prod(tmpret)*noParentProb
		tmprob<-prod(tmpret)
		#print(tmpret)
		#print(tmprob)
		ret<-ret+tmprob

	}else{# ジェノタイプ不明者が1人以上含まれる場合
		
	# すべての未確定ジェノタイプについて場合分けして計算し、その和をとる
	# 集団ジェノタイプ頻度の影響を受けるのは親が与えられていない人のみ
	# それ以外の個人で不明ジェノタイプの確率は「組み込まれる」
	Ng<-na*(na+1)/2
	#print(Unfixed)
	#print(Nunfixed)
	#print(Ng)
	#limitedDiplotype<-LimitDiplotype(p,g,A)
	#Gvec<-list()
	#for(i in 1:length(limitedDiplotype[,1,1])){
	#	Gvec[[i]]<-DiplotypeListFromMat(limitedDiplotype[i,,])
	#}

	diplotype<-matrix(0,Ng,2)
	diplotypeProb<-rep(0,Ng)
	cnt<-1
	for(i in 1:na){
		for(j in 1:i){
			diplotype[cnt,1]<-A[i]
			diplotype[cnt,2]<-A[j]
			diplotypeProb[cnt]<-P[i]*P[j]
			if(i!=j)diplotypeProb[cnt]<-2*diplotypeProb[cnt]
			cnt<-cnt+1
		}
	}
	
	counters<-rep(0,Nunfixed)
	nchild<-rep(Ng,Nunfixed)
	#print(counters)
	#print(nchild)
	loop<-TRUE
	ret<-0
	noParentProb<-1
	#print("here")
	while(loop){
		#print(loop)
		#print("again?")
		tmpret<-rep(1,ns)
		currentgenotype<-counters+1
		currentg<-g
		#print(counters)
		#print(nchild)
		counters<-counterplus(counters,nchild)
		#print("post")
		#print(counters)
		#print(sum(counters))
		currentg[Unfixed,]<-diplotype[currentgenotype,]
		#print(currentg)
		gvector<-array(0,c(ns,na,2))
		#given<-rep(0,ns)
		for(i in 1:ns){
			#tmp<-0
			#if(g[i,1]!=unknown){
				gvector[i,which(A==currentg[i,1]),1]<-1
				#tmp<-tmp+1
			#}
			#if(g[i,2]!=unknown){
				gvector[i,which(A==currentg[i,2]),2]<-1
				#tmp<-tmp+1
			#}
			#if(tmp==2)given[i]<-1
		}
		#print(gvector)
		# 被捜索対象者のジェノタイプがないときというのは
		# その人が「集団からの任意の個人」であるときの尤度と比較したいからなので
		# そのように与える
		# ここが難しい!
		#tobegeneral<-rep(0,ns)
		#for(i in 1:ns){
		#	if(g[i,1]==unknown & p[i,5]==2){
		#		gvector[i,,1]<-gvector[i,,2]<-P
		#		print(i)
		#		print(g[i,1])
		#		print(p[i,5])
		#	}
		#}
		# すべての親子トリオについて計算
		trioloop<-TRUE
		triocnt<-1
		while(trioloop){
			
			#print(i)
			P1<-P2<-P
			tmp1<-p[triocnt,2]
			tmp2<-p[triocnt,3]
			if(tmp1!=0){
				P1<-(gvector[tmp1,,1]+gvector[tmp1,,2])/2
			}
			if(tmp2!=0){
				P2<-(gvector[tmp2,,1]+gvector[tmp2,,2])/2
			}
			
			FromParents<-OffspringGenotypeProb(P1,P2)
			#tmpallele1<-which(A==currentg[triocnt,1])
			#tmpallele2<-which(A==currentg[triocnt,2])
			#lessa<-min(tmpallele1,tmpallele2)[1]
			#morea<-max(tmpallele1,tmpallele2)[1]
			#tmptmp<-FromParents[lessa,morea]
			Offspring<-OffspringGenotypeProb(gvector[triocnt,,1],gvector[triocnt,,2])
			#ret[i]<-log(sum(FromParents*Offspring))
			tmptmp<-sum(FromParents*Offspring)
			tmpret[triocnt]<-tmptmp
			if(tmptmp==0)trioloop<-FALSE

			#print(tmptmp)
			triocnt<-triocnt+1
			if(triocnt==(ns+1))trioloop<-FALSE
		}
		#tmprob<-prod(tmpret)*noParentProb
		tmprob<-prod(tmpret)
		#print(tmprob)
		ret<-ret+tmprob
		
		#print(ret)
		
		#print(currentg)
		#print(counters)
		if(sum(counters)==0){
			loop<-FALSE
			#break
		}
		#print(loop)
		#if(!loop){print("YYY")}
		#if(loop){print("bbbbb")}
	}
	}
	
	list(loglikesum=log(ret,10),like=ret)
}
gameteFreq<-function(P){
 (apply(P,1,sum)+apply(P,2,sum))/2
}
OffspringDiplotypeFreq<-function(m,f){# m:mother,f:father
 gameteM<-gameteFreq(m)
 gameteF<-gameteFreq(f)
 outer(gameteM,gameteF,FUN="*")
}
LimitDiplotype<-function(p,g,A){
	unknown="0"
	ns<-length(p[,1])
	na<-length(A)
	a<-array(0,c(ns,na,na))
	for(i in 1:ns){
		if(g[i,1]==unknown){
			a[i,,]<-1/(na^2)
		}else{
			tmpa1<-which(A==g[i,1])
			tmpa2<-which(A==g[i,2])
			a[i,tmpa1,tmpa2]<-a[i,tmpa2,tmpa1]<-1
		}
		#print(a[i,,])
	}
	#print(a)
	for(i in 1:ns){
		tmpm<-p[i,2]
		tmpf<-p[i,3]
		if(tmpm!=0){
			print(i)
			print(tmpm)
			print(tmpf)
			print(a[tmpm,,])
			print(a[tmpf,,])
			#print(OffspringDiplotypeFreq(a[tmpm,,],a[tmpf,,]))
			#print("a")
			#print(a[i,,])
			#print(a[i,,]*OffspringDiplotypeFreq(a[tmpm,,],a[tmpf,,]))
			a[i,,]<-a[i,,]*OffspringDiplotypeFreq(a[tmpm,,],a[tmpf,,])
		}
		a[i,,]<-a[i,,]+t(a[i,,])
	}
	sign(a)
}
MakeVectorFromDiplotypeMat<-function(m){
	tmp<-m
	tmp[,]<-1
	tmp[upper.tri(tmp)]<-0
	m[which(tmp==1,arr.ind=TRUE)]
}

DiplotypeListFromMat<-function(m){
	which(MakeVectorFromDiplotypeMat(m)==1)
}

CalcLogLikelihoodFamily<-function(p,g,Alleles,Probs){
	#print("family")
	#print(g)
	nl<-length(g[1,1,])
	ret<-rep(0,nl)
	
	for(i in 1:length(g[1,1,])){
		tmp<-CalcLogLikelihoodFamilyLocus2(p,g[,,i],Alleles[[i]],Probs[[i]])
		#if(tmp$like==0){
		#	ret[i]<-tmp$loglikesum
		#	return(list(loglikelist=ret,loglikesum=sum(ret)))
		#}
		ret[i]<-tmp$loglikesum
	}
	list(loglikelist=ret,loglikesum=sum(ret))
}
printTAPout<-function(tap){
	n<-length(tap$LogLike)
	
	for(i in 1:n){
		if(tap$LogLike[[i]][[2]][1]!=-Inf){
			tmp<-paste("LogLike",tap$LogLike[[i]][[2]][1],"Prob",10^(tap$LogLike[[i]][[2]][1]),"SearchedID",c(tap$SearchedID[[i]]),"CandidateID",c(tap$CandidateID[[i]]))
			print(tmp)
		}
		
	}
}

printSMPFout<-function(SMMFout){
	n<-length(SMMFout$LogLike)
	plot(SMMFout$pedigree,id=unlist(SMMFout$IndNames),main=SMMFout$FamilyName)


	n<-length(SMMFout$LogLike)
	for(i in 1:n){
		if(SMMFout$LogLike[[i]]!=-Inf){
			tmp<-paste("LogLike",SMMFout$LogLike[[i]],"Prob",10^(SMMFout$LogLike[[i]]),"SearchedID",c(SMMFout$SearchedID[[i]]),"CandidateID",c(SMMFout$CandidateID[[i]]))
			print(tmp)
		}
		
	}
}


TestAllPatterns2<-function(p,G,Gpool,candidates,Alleles,Probs){
	# 家系情報から家系図と核家族を作る
	ptemp<-MakePedigreeFromFamilyInfo(p)
	nucs<-subnucs(p)
	# 被捜索者の捜索パターン
	tobeSearched<-which(p[,5]==2)
	SearchPattern<-MakeSearchPattern(tobeSearched)
	
	# Gpoolの候補者が「ポピュレーションで観察される確率」をあらかじめ算出しよう
	
	LikeFromGenPop<-rep(0,length(Gpool[,1,1]))
	for(i in 1:length(Alleles)){
		tmp<-OffspringGenotypeProb(Probs[[i]],Probs[[i]])
		for(j in 1:length(LikeFromGenPop)){
			tmpallele1<-which(Alleles[[i]]==Gpool[j,1,i])
			tmpallele2<-which(Alleles[[i]]==Gpool[j,2,i])
			tmp2<-log(tmp[min(tmpallele1,tmpallele2)[1],max(tmpallele1,tmpallele2)[1]],10)
			LikeFromGenPop[j]<-LikeFromGenPop[j]+tmp2
		}
	}
	# print(LikeFromGenPop)
	
	
	ret<-NULL
	ret2<-NULL
	ret3<-NULL
	counter<-1
	
	for(i in 1:length(SearchPattern$X)){
		spx<-as.matrix(SearchPattern$X[[i]])
		
		for(j in 1:length(spx[,1])){
			ToTestCurrent<-spx[j,]
			tmpG2<-G
			#BaseLike<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs)$loglikesum
			BaseLike<-CalcLogLikelihoodFamilyX(p,nucs,G,Alleles,Probs)$loglikesum
			#ret[[counter]]<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs)
			#ret2[[counter]]<-list(ToTestCurrent)
			#ret3[[counter]]<-list()
			#counter<-counter+1
			
			numtested<-length(spx[j,])
			tmpcandidates<-NULL

			for(k in 1:numtested){
				selectedToBeTested<-which(tobeSearched==spx[j,k])
				
				tmpcandidates[[k]]<-candidates[[selectedToBeTested]]
			}
			print("done for none")
			tmpTestPattern<-MakeTestPattern(tmpcandidates)
			for(k in 1:length(tmpTestPattern[,1])){
				#print("to be tested IDs in pedigree")
				#print(ToTestCurrent)
				#print("to be tested sampleIDs in pool")
				#print(tmpTestPattern[k,])

				tmpG2<-G
				tmpG2[ToTestCurrent,,]<-Gpool[tmpTestPattern[k,],,]
				LoutTmp<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs)$loglikesum
				#print(LoutTmp)
				ret[[counter]]<-LoutTmp-BaseLike-sum(LikeFromGenPop[tmpTestPattern[k,]])
				ret2[[counter]]<-list(ToTestCurrent)
				ret3[[counter]]<-list(tmpTestPattern[k,])
				counter<-counter+1
			}

		}
	}
	LRmatrix<-matrix(0,counter,counter)
	for(i in 1:(counter-1)){
		for(j in 1:(counter-1)){
			LRmatrix[i,j]<-ret[[i]]-ret[[j]]
		}
	}
	
	list(LogLike=ret,LRmatrix=LRmatrix,SearchedID=ret2,CandidateID=ret3,familyinfo=p,pedigree=ptemp,tobeSearched=tobeSearched,candidates=candidates,offeredGenotype=G,poolGenotype=Gpool,Alleles=Alleles,Probs=Probs)
}

printTAPout2<-function(tap){
	n<-length(tap$LogLike)
	for(i in 1:n){
		if(tap$LogLike[[i]]!=-Inf){
			tmp<-paste("LogLike",tap$LogLike[[i]],"Prob",10^(tap$LogLike[[i]]),"SearchedID",c(tap$SearchedID[[i]]),"CandidateID",c(tap$CandidateID[[i]]))
			print(tmp)
		}
		
	}
}

MakeInfo<-function(LDZout,gnew){
	ret<-LDZout$dlistA
	unknown="0"
	for(i in 1:length(gnew[,1])){
		if(gnew[i,1]!=unknown){
			ret[[i]]<-list(as.set(gnew[i,]))
		}
	}
	ret
}


MakeDiplotypePrior<-function(p,g,A,P,LDZout,info){
	ret<-list()
	for(i in 1:length(p[,1])){
		ret[[i]]<-rep(1,length(LDZout$dlist[[i]]))
		#print(ret[[i]])
		if(length(info[[i]])>1){
			if(p[i,2]=="0"){
				tmph1<-LDZout$hset1A[[i]]
				tmph2<-LDZout$hset2A[[i]]
				
				for(j in 1:length(LDZout$dlist[[i]])){
					tmpp<-1
					#print(LDZout$dlistA[[i]][[j]])
					heterocheck<-TRUE
					for(k in LDZout$dlistA[[i]][[j]]){
						#print(k)
						if(!(k %e% (tmph1 & tmph2)))heterocheck<-FALSE
						tmpp<-tmpp*P[which(A==k)]
					}

					if(length(LDZout$dlistA[[i]][[j]])==2){
						if(heterocheck){
							#print("heterocheck")
					#print(heterocheck)
							tmpp<-tmpp*2	
						}
					}
					
					ret[[i]][j]<-tmpp
				}
				ret[[i]]<-ret[[i]]/sum(ret[[i]])
			}

			
		}else{
			ret[[i]]<-rep(0,length(LDZout$dlist[[i]]))
			for(j in 1:length(LDZout$dlist[[i]])){
				#print(info[[i]])
				#print(LDZout$dlistA[[i]][[j]])
				if(info[[i]][[1]] == LDZout$dlistA[[i]][[j]]){
					ret[[i]][j]=1
				}
			}
		}
	}
	ret
}

CalcProbNucs<-function(p,g,A,P,nucs,LDZout){
	rets<-list()
	for(nn in 1:length(nucs)){
		nuc<-nucs[[nn]]
		dimvector<-nuc[2:length(nuc)]
		#print(nuc)
		#print(dimvector)
		#print("Nmembers")
		#print(length(p[,1]))
		fa<-nuc[2]
		mo<-nuc[3]
		offs<-nuc[4:length(nuc)]
		pivot<-nuc[1]
		if(pivot!=0){
			who<-which(nuc[2:length(nuc)]==pivot)
			dimvector<-nuc[2:length(nuc)]
			pivotLoc<-which(dimvector==pivot)
			dimvector<-c(dimvector[-pivotLoc],pivot)
		}
		#print(dimvector)
		genotypeNumVector<-rep(0,length(dimvector))
		#print(LDZout$dlist)
		for(i in 1:length(dimvector)){
			#genotypeNumVector[i]<-length(LDZout$dset[[dimvector[i]]])
			#genotypeNumVector[i]<-length(LDZout$dtypeset[[dimvector[i]]])
			genotypeNumVector[i]<-length(LDZout$dlist[[dimvector[i]]])
		}
		genotypeNumVector
		ret<-array(1,genotypeNumVector)
		faid<-which(dimvector==fa)
		moid<-which(dimvector==mo)
		ofid<-which(dimvector!=fa & dimvector!=mo)

		counter<-array(1:prod(genotypeNumVector),genotypeNumVector)
		for(i in 1:length(ret)){
			x<-which(counter==i,arr.ind=TRUE)
			#print("----")
			tmpprob<-1
			faset<-LDZout$dlist[[dimvector[faid]]][[x[faid]]]
			moset<-LDZout$dlist[[dimvector[moid]]][[x[moid]]]
			for(j in ofid){
				ofset<-LDZout$dlist[[dimvector[j]]][[x[j]]]
				#print(faset)
				#print(moset)
				#print(ofset)
				tmptmpprob<-LikeTrio(faset,moset,ofset)
				#print(tmptmpprob)
				#print(tmpprob)
				tmpprob<-tmpprob*tmptmpprob
				#print(tmpprob)
				
			}
			ret[x]<-tmpprob
		}
		#print("====")
		#print(ret)
		rets[[nn]]<-list(ret,dimvector)
	}
	rets
}
	
LikeTrio<-function(setf,setm,seto){
	couple<-setf*setm
	cnt<-0
	for(i in couple){
		if(as.set(i) == seto)cnt<-cnt+1
	}
	cnt/length(couple)
}

LimitAlleles<-function(g,A,P){
	tmp<-rep(0,length(A))
	for(i in 1:length(g)){
		tmp[which(A==g[i])]<-1
	}
	A2<-c(A[which(tmp==1)],"pool")
	P2<-c(P[which(tmp==1)])
	P2<-c(P2,1-sum(P2))
	list(A=A2,P=P2)
}
LikeNucWithPrior<-function(cpnout,nucs,DiplotypePrior){
	prob<-list()
	for(nn in 1:length(nucs)){
		nucDimVector<-cpnout[[nn]][[2]]
		tmp<-DiplotypePrior[[nucDimVector[1]]]
		for(i in 2:length(nucDimVector)){
			tmp<-tmp%o%DiplotypePrior[[nucDimVector[i]]]
		}
		#print(dim(cpnout[[nn]][[1]]))
		#print(dim(tmp))
		#print("---")
		prob[[nn]]<-tmp * cpnout[[nn]][[1]]
	}
	prob
}


# pivotで足し合わせ

SumPivot<-function(cpnout,like,info){
	ret<-NULL
	#pivotted<-set()
	#pivotted<-list()
	pivotted<-rep(0,length(info))
	cumulProb<-list()
	for(nn in 1:length(cpnout)){
		tmpdim<-dim(cpnout[[nn]][[1]])
		tmp<-rep(1,tmpdim[1])
		if(pivotted[[cpnout[[nn]][[2]][[1]]]]==1){
			#print("yesagain")
			tmp<-cumulProb[[cpnout[[nn]][[2]][[1]]]]
		}
		for(i in 2:length(tmpdim)){
			tmp2<-rep(1,tmpdim[i])
			if(pivotted[[cpnout[[nn]][[2]][[i]]]]==1){
				#print("yesagain")
				tmp2<-cumulProb[[cpnout[[nn]][[2]][[i]]]]
			}
			tmp<-tmp%o%tmp2
		}
		#print(tmp)
		#print(like[[nn]])
		#print(dim(tmp))
		#print(dim(like[[nn]]))
		tmp3<-like[[nn]]*tmp
		#print(tmp3)
		#ret<-apply(like[[nn]],length(cpnout[[nn]][[2]]),sum)
		ret<-apply(tmp3,length(cpnout[[nn]][[2]]),sum)
		#pivotted<-pivotted | as.set(cpnout[[nn]][[2]][[length(cpnout[[nn]][[2]])]])
		#pitotted[[cpnout[[nn]][[2]][[length(cpnout[[nn]][[2]])]]]]<-1
		pivotted[cpnout[[nn]][[2]][[length(cpnout[[nn]][[2]])]]]<-1
		cumulProb[[cpnout[[nn]][[2]][[length(cpnout[[nn]][[2]])]]]]<-ret
		#print(ret)
		#print(pivotted)
		#print(cumulProb)
		#print("-----")
	}
	prob<-sum(ret)
	list(like=prob,loglike=log(prob,10))
}
CalcLikeZ<-function(p,G,nucs,Alleles,Probs){
	tmpret<-0
	for(na in 1:length(Alleles)){
	
		A<-unlist(Alleles[[na]])
		P<-Probs[[na]]
		g<-G[,,na]
		gpool<-G[,,na]
		#A2P2<-LimitAlleles(g,A,P)
		A2P2<-LimitAlleles(gpool,A,P)
		A2<-A2P2[[1]]
		P2<-A2P2[[2]]
		#LDZout<-LimitDiplotypeZ(p,g,A)

		LDZout2<-LimitDiplotypeZ(p,g,A2)
		MCheck<-TRUE
		for(i in 1:length(LDZout2$dlistA)){
			#print(LDZout2$distA[[i]])
			if(length(LDZout2$dlistA[[i]])==0){
				print("checkFALSE")
				print(p)
				print(g)
				print(LDZout2$dlistA)
				MCheck<-FALSE
			}
		}
		if(!MCheck){
			return(-Inf)
		}
		#cpnout<-CalcProbNucs(p,g,A,P,nucs,LDZout)
		#info<-LDZout$dlistA
		cpnout2<-CalcProbNucs(p,g,A2,P2,nucs,LDZout2)
		info2<-LDZout2$dlistA
		#print(g)
#print(info2)
		DiplotypePrior2<-MakeDiplotypePrior(p,g,A2,P2,LDZout2,info2)
		like2<-LikeNucWithPrior(cpnout2,nucs,DiplotypePrior2)
		tmpret<-tmpret+SumPivot(cpnout2,like2,info2)$loglike

	}
	tmpret
}
# pedigreeごとに
# 家系情報
# 生存・協力者のジェノタイプ情報
# 被捜索者ごとに、身元不明者リスト中の候補者がリストアップされている
SearchMissingsMultiFamily<-function(pedigrees,genotypesFamily,Gpool,candidatesList,FamilyNames,IndNames,Alleles,Probs){
	output<-list()
	# Gpoolの候補者が「ポピュレーションで観察される確率」をあらかじめ算出しよう
	
	LikeFromGenPop<-rep(0,length(Gpool[,1,1]))
	for(i in 1:length(Alleles)){
		tmp<-OffspringGenotypeProb(Probs[[i]],Probs[[i]])
		for(j in 1:length(LikeFromGenPop)){
			tmpallele1<-which(Alleles[[i]]==Gpool[j,1,i])
			tmpallele2<-which(Alleles[[i]]==Gpool[j,2,i])
			tmp2<-log(tmp[min(tmpallele1,tmpallele2)[1],max(tmpallele1,tmpallele2)[1]],10)
			LikeFromGenPop[j]<-LikeFromGenPop[j]+tmp2
		}
	}

	for(ip in 1:length(pedigrees)){
		output[[ip]]<-SearchMissingsPerFamily(pedigrees[[ip]],genotypesFamily[[ip]],
		Gpool,candidatesList[[ip]],FamilyNames[[ip]],IndNames[[ip]],Alleles,Probs,LikeFromGenPop)
	}
	output
}



SearchMissingsPerFamily<-function(p,G,Gpool,candidatesL,FamilyNs,IndNs,Alleles,Probs,LikeFromGenPop){
	# 出力は家系ごと
	
	# 家系情報を遺伝系パッケージ用のオブジェクトにする
	# kinshipパッケージオブジェクト:家系図描図用
	ptemp<-MakePedigreeFromFamilyInfo(p)
	# paramlinkパッケージオブジェクト:メンデリアンチェック用
	tmpPed<-p[,1:5]
	tmpPed<-data.frame(ID=tmpPed[,1],FID=tmpPed[,3],MID=tmpPed[,2],SEX=tmpPed[,4]+1,AFF=tmpPed[,5]-1)
	tmpPed<-linkdat(tmpPed,model=1)

	# 家系図を描いてPDF出力する
	plot(ptemp,id=unlist(IndNs))
	pdffile<-paste(FamilyNs,".pdf")
	#print(IndNs)
	pdf(pdffile,family="Japan1")
	plot(ptemp,id=unlist(IndNs))
	dev.off()
	
	# 核家族処理をする
	nucs<-subnucs(p)

	# 尤度を計算しよう
	searched<-which(p[,5]==2)
	candidates<-list()
	for(i in 1:length(searched)){
		candidates[[i]]<-candidatesL[[searched[[i]]]]
	}
	#print(candidates)
	# BaseLikeの計算(すべての被捜索者に候補者を当てはめない場合)は
	# 1チェックパターンでも「メンデリアンチェックを通ったら」実施する
	# そのためのチェッカー
	BaseMade<-FALSE
	#BaseLike<-CalcLikeZ(p,G,nucs,Alleles,Probs)
	
	SearchPattern<-MakeSearchPattern(searched)
	ret<-NULL
	ret2<-NULL
	ret3<-NULL
	counter<-1

	for(i in 1:length(SearchPattern$X)){
		spx<-as.matrix(SearchPattern$X[[i]])
		
		for(j in 1:length(spx[,1])){
			ToTestCurrent<-spx[j,]
			#tmpG2<-G
			#BaseLike<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs)$loglikesum
			#BaseLike<-CalcLogLikelihoodFamilyX(p,nucs,G,Alleles,Probs)$loglikesum
			#ret[[counter]]<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs)
			#ret2[[counter]]<-list(ToTestCurrent)
			#ret3[[counter]]<-list()
			#counter<-counter+1
			
			#numtested<-length(spx[j,])
			numtested<-length(ToTestCurrent)
			tmpcandidates<-list()
			if(numtested>0){
				for(k in 1:numtested){
					selectedToBeTested<-which(searched==ToTestCurrent[k])
					tmpcandidates[[k]]<-candidates[[selectedToBeTested]]
				}
			}

			#print("done for none")
			tmpTestPattern<-MakeTestPattern(tmpcandidates)
			for(k in 1:length(tmpTestPattern[,1])){
				#print("to be tested IDs in pedigree")
				#print(ToTestCurrent)
				#print("to be tested sampleIDs in pool")
				#print(tmpTestPattern[k,])

				tmpG2<-G
				tmpG2[ToTestCurrent,,]<-Gpool[tmpTestPattern[k,],,]
				#print(tmpG2)
				# mendelian Check
				tttmpG<-tmpG2[,,1]
				for(i in 2:length(tmpG2[1,1,])){
					tttmpG<-cbind(tttmpG,tmpG2[,,i])
				}


				tmpPed2<-setMarkers(tmpPed, m=tttmpG)
				
				mCout<-mendelianCheck(tmpPed2)
				#print("mendelianChecked")
				#print(mCout)
				#print(length(mCout))
				if(length(mCout)!=0){
					ret[[counter]]<--Inf
					ret2[[counter]]<-list(ToTestCurrent)
					ret3[[counter]]<-list(tmpTestPattern[k,])
					counter<-counter+1
				}else{
					if(!BaseMade){
						BaseLike<-CalcLikeZ(p,G,nucs,Alleles,Probs)
						BaseMade<-TRUE
					}
					LoutTmp<-CalcLikeZ(p,tmpG2,nucs,Alleles,Probs)
					#LoutTmp<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs)$loglikesum
					#print(LoutTmp)
					# tmpTestPattern[k,]って変?
					ret[[counter]]<-LoutTmp-BaseLike-sum(LikeFromGenPop[tmpTestPattern[k,]])
					ret2[[counter]]<-list(ToTestCurrent)
					ret3[[counter]]<-list(tmpTestPattern[k,])
					counter<-counter+1
				}
				
			}

		}
	}
	#print(ret)
	list(LogLike=ret,SearchedID=ret2,CandidateID=ret3,pedigree=ptemp,FamilyName=FamilyNs,IndNames=IndNs)
}

CalcLikeForCandidates<-function(p,G,candidates,Gcandidates,Alleles,Probs){
	# Gpoolの候補者が「ポピュレーションで観察される確率」をあらかじめ算出しよう
	
	LikeFromGenPop<-rep(0,length(Gpool2[,1,1]))
	for(i in 1:length(Alleles)){
		tmp<-OffspringGenotypeProb(Probs[[i]],Probs[[i]])
		for(j in 1:length(LikeFromGenPop)){
			tmpallele1<-which(Alleles[[i]]==Gpool2[j,1,i])
			tmpallele2<-which(Alleles[[i]]==Gpool2[j,2,i])
			tmp2<-log(tmp[min(tmpallele1,tmpallele2)[1],max(tmpallele1,tmpallele2)[1]],10)
			LikeFromGenPop[j]<-LikeFromGenPop[j]+tmp2
		}
	}
	
	# BaseLikeの計算(すべての被捜索者に候補者を当てはめない場合)
	nucs<-subnucs(p)
	BaseLike<-CalcLikeZ(p,G2,nucs,Alleles,Probs)
	
	SearchPattern<-MakeSearchPattern(searched)
	ret<-NULL
	ret2<-NULL
	ret3<-NULL
	counter<-1

	#SearchPattern<-2^(as.set(searched))
	for(i in 1:length(SearchPattern$X)){
	#for(i in SearchPattern){
		spx<-as.matrix(SearchPattern$X[[i]])
		
		for(j in 1:length(spx[,1])){
			ToTestCurrent<-spx[j,]
			tmpG2<-G
			#BaseLike<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs)$loglikesum
			#BaseLike<-CalcLogLikelihoodFamilyX(p,nucs,G,Alleles,Probs)$loglikesum
			#ret[[counter]]<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs)
			#ret2[[counter]]<-list(ToTestCurrent)
			#ret3[[counter]]<-list()
			#counter<-counter+1
			
			#numtested<-length(spx[j,])
			numtested<-length(ToTestCurrent)
			tmpcandidates<-list()
			if(numtested>0){
				for(k in 1:numtested){
					selectedToBeTested<-which(searched==ToTestCurrent[k])
					tmpcandidates[[k]]<-candidates[[selectedToBeTested]]
				}
			}

			#print("done for none")
			tmpTestPattern<-MakeTestPattern(tmpcandidates)
			for(k in 1:length(tmpTestPattern[,1])){
				#print("to be tested IDs in pedigree")
				#print(ToTestCurrent)
				#print("to be tested sampleIDs in pool")
				#print(tmpTestPattern[k,])

				tmpG2<-G
				tmpG2[ToTestCurrent,,]<-Gpool2[tmpTestPattern[k,],,]
				#print(tmpG2)
				LoutTmp<-CalcLikeZ(p,tmpG2,nucs,Alleles,Probs)
				#LoutTmp<-CalcLogLikelihoodFamily(p,tmpG2,Alleles,Probs)$loglikesum
				#print(LoutTmp)
				ret[[counter]]<-LoutTmp-BaseLike-sum(LikeFromGenPop[tmpTestPattern[k,]])
				ret2[[counter]]<-list(ToTestCurrent)
				ret3[[counter]]<-list(tmpTestPattern[k,])
				counter<-counter+1
			}

		}
	}
	list(LogLike=ret,SearchedID=ret2,CandidateID=ret3)

}

# AmpFlSTRIdentifilerの15 STR

Identifiler15<-c("D8S1179","D21S11","D7S820","CSF1PO","D3S1358","TH01","D13S317","D16S539","D2S1338","D19S433","vWA","TPROX","D18S51","D5S818","FGA")

AllelesId<-NULL
ProbsIDYoshidaN<-NULL
ProbsIDYoshidaFreq<-NULL
ProbsIDYoshidaFreq2<-NULL
AllelesId[[1]]<-c("7",9,10,11,12,13,14,15,16,17,18)
AllelesId[[2]]<-c("27",28,28.2,29,30,30.2,30.3,31,31.2,32,32.2,33,33.1,33.2,34,34.2)
AllelesId[[3]]<-c("7",6,9,9.1,10,10.1,10.3,11,12,13,14,15)
AllelesId[[4]]<-c("7",8,9,10,11,12,13,14,15,16)
AllelesId[[5]]<-c("12",13,14,15,16,17,18,19,21)
AllelesId[[6]]<-c("5",6,7,8,9,9.3,10)
AllelesId[[7]]<-c("7",8,9,10,11,12,13,14,15)
AllelesId[[8]]<-c("7",8,9,10,11,12,13,14,15)
AllelesId[[9]]<-c("16",17,18,19,20,21,22,23,24,25,26,27,28)
AllelesId[[10]]<-c("9.2",10.2,11,11.2,12,12.2,13,13.2,14,14.2,15,15.2,16,16.2,17.2,18)
AllelesId[[11]]<-c("13",14,15,16,17,18,19,20,21,22)
AllelesId[[12]]<-c("7",8,9,10,11,12,13,14)
AllelesId[[13]]<-c("10",11,12,13,14,15,16,17,17.1,18,19,20,21,22,23,24,25,26,27)
AllelesId[[14]]<-c("7",8,9,10,11,12,13,14,15)
AllelesId[[15]]<-c("17",18,19,20,21,22,22.2,23,23.2,24,24.2,25,25.2,26,27,28)

ProbsIDYoshidaN[[1]]<-c(1,2,358,294,331,607,553,363,173,17,1)
ProbsIDYoshidaN[[2]]<-c(3,114,13,666,910,12,3,279,192,52,321,7,6,109,2,11)
ProbsIDYoshidaN[[3]]<-c(8,342,123,1,591,1,2,887,634,94,16,1)
ProbsIDYoshidaN[[4]]<-c(27,3,126,603,561,1130,186,48,14,2)
ProbsIDYoshidaN[[5]]<-c(5,2,78,1069,827,539,171,8,1)
ProbsIDYoshidaN[[6]]<-c(4,603,720,179,1076,94,24)
ProbsIDYoshidaN[[7]]<-c(3,721,348,310,597,546,138,35,2)
ProbsIDYoshidaN[[8]]<-c(1,5,955,531,505,482,196,22,3)
ProbsIDYoshidaN[[9]]<-c(23,263,429,564,285,40,136,396,291,166,78,23,6)
ProbsIDYoshidaN[[10]]<-c(1,3,10,2,109,14,776,81,944,238,137,310,14,52,8,1)
ProbsIDYoshidaN[[11]]<-c(1,524,72,497,763,609,199,27,7,1)
ProbsIDYoshidaN[[12]]<-c(1,1227,308,83,980,96,3,2)
ProbsIDYoshidaN[[13]]<-c(6,12,129,538,600,454,339,220,1,130,98,59,42,35,20,10,4,2,1)
ProbsIDYoshidaN[[14]]<-c(7,18,232,542,789,635,450,24,3)
ProbsIDYoshidaN[[15]]<-c(9,58,181,240,353,544,5,554,13,424,3,198,5,86,21,6)
for(i in 1:15){
	ProbsIDYoshidaFreq[[i]]<-ProbsIDYoshidaN[[i]]/sum(ProbsIDYoshidaN[[i]])
	ProbsIDYoshidaFreq2[[i]]<-(ProbsIDYoshidaN[[i]]+1)/(sum(ProbsIDYoshidaN[[i]])+1)
	ProbsIDYoshidaFreq2[[i]][which(ProbsIDYoshidaN[[i]]<5)]<-5/sum(ProbsIDYoshidaN[[i]])
}