bin/IDEAS_2018/IDEAS_APP/tensorHMMbase.cpp
//This function is specific to categorical data
//return maxK
//nh: group size per position
//Q: group size overall
//r: recombination probability
//asz: number of distinct values per position
This function is specific to categorical data
return maxK
nh: group size per position
Q: group size overall
r: recombination probability
asz: number of distinct values per position
This function is specific to categorical data return maxK nh : group size per position Q : group size overall r : recombination probability asz : number of distinct values per position
void tensorHMMbase::_addInd(MYDATA const &mydata, TENSORPARA &tpara, int id, float wb, int st, int ed)
{       int i, j, k, idd = id / mydata.ploidity, idst = mydata.indIndex[idd], ided = mydata.indIndex[idd + 1];
	if(st < 0) st = 0;
	if(ed < 0) ed = mydata.L;
	int const *pp = tpara.pop + id * mydata.L + st, *pp2 = pp + (mydata.ploidity - 1) * mydata.L;
	float const *dp = mydata.data + idst * mydata.ploidity * mydata.L + st, *dp2 = dp + (mydata.ploidity - 1) * mydata.L;
	bool const *rr = tpara.rec + id * mydata.L + st, *rr2 = rr + (mydata.ploidity - 1) * mydata.L + st;
	float *rn = tpara.recombN + st;
	float *nhp = tpara.nh + st * tpara.maxK;
	int K = tpara.maxK;
	int const *tpp = pp, *tpp2 = pp2;
	for(i = st; i < ed; i++, tpp++, tpp2++)
	{       if(*tpp >= K) K = (*tpp) + 1;
		if(*tpp2 >= K) K = (*tpp2) + 1;
	}
	if(K > tpara.maxK)
	{	PARAUNIT pu;
		pu.lapseN = pu.popN = 0;	
		float *newnh = new float[K * mydata.L], *newnhp = newnh; 
		int *aszp = mydata.asz;
		for(i = 0; i < K * mydata.L; i++, newnhp++) *newnhp = 0;
		newnhp = newnh;
		nhp = tpara.nh;
		for(i = 0; i < mydata.L; i++, aszp++)
		{       for(j = 0; j < tpara.maxK; j++, newnhp++, nhp++)
				(*newnhp) = (*nhp);
			newnhp += K - tpara.maxK;
			tpara.param[i].resize(K, pu);
		}
		delete tpara.nh;
		tpara.nh = newnh;
		for(i = tpara.maxK; i < K + 1; i++) tpara.Q[i] = 0;
		tpara.maxK = K;
		nhp = tpara.nh + st * tpara.maxK;
	}
	bool *tlapse = tpara.lapse + idst * mydata.L + st;
	vector<vector<PARAUNIT> >::iterator pit = tpara.param.begin() + st;
	int ksz = mydata.L * mydata.ploidity, kkk = idst * mydata.ploidity, kkksz = mydata.ploidity;
	for(i = st; i < ed; i++, pp++, pp2++, dp++, dp2++, nhp += tpara.maxK, pit++, rr++, rr2++, rn++, tlapse++)
	{ 
		if(*rr) *rn += tpara.idw[id];
		int kk = 0, kks = kkk;
		for(k = 0; k < ided - idst; k++, kk += ksz, kks += kkksz) 
		{	float b = tpara.indWeight[idst + k];
			*(nhp + (*pp)) += b;
			_updatePara_add((*pit)[*pp], (int)(*(dp + kk)), (float)b, kks, (int)(tpara.lapse!=NULL && *(tlapse+kk)));
		}
		if(i == 0 || *rr) tpara.Q[*pp]+=tpara.idw[id];
		if(mydata.ploidity == 2)
		{        
			if(*rr2) *rn += tpara.idw[id];
			kk=0; kks=kkk+1;
			for(k = 0; k < ided - idst; k++, kk+=ksz, kks+=kkksz) 
			{	float b = tpara.indWeight[idst + k];
				*(nhp + (*pp2)) += b;
				_updatePara_add((*pit)[*pp2], (int)(*(dp2 + kk)), (float)b, kks,(int)(tpara.lapse!=NULL && *(tlapse+kk)));
			}
			if(i == 0 || *rr2) tpara.Q[*pp2]+=tpara.idw[id];
		}
	}
}
void tensorHMMbase :: _addInd ( MYDATA const & mydata , TENSORPARA & tpara , int id , float wb , int st , int ed ) { int i , j , k , idd = id / mydata . ploidity , idst = mydata . indIndex [ idd ] , ided = mydata . indIndex [ idd + 1 ] ; if ( st < 0 ) st = 0 ; if ( ed < 0 ) ed = mydata . L ; int const * pp = tpara . pop + id * mydata . L + st , * pp2 = pp + ( mydata . ploidity - 1 ) * mydata . L ; float const * dp = mydata . data + idst * mydata . ploidity * mydata . L + st , * dp2 = dp + ( mydata . ploidity - 1 ) * mydata . L ; bool const * rr = tpara . rec + id * mydata . L + st , * rr2 = rr + ( mydata . ploidity - 1 ) * mydata . L + st ; float * rn = tpara . recombN + st ; float * nhp = tpara . nh + st * tpara . maxK ; int K = tpara . maxK ; int const * tpp = pp , * tpp2 = pp2 ; for ( i = st ; i < ed ; i ++ , tpp ++ , tpp2 ++ ) { if ( * tpp >= K ) K = ( * tpp ) + 1 ; if ( * tpp2 >= K ) K = ( * tpp2 ) + 1 ; } if ( K > tpara . maxK ) { PARAUNIT pu ; pu . lapseN = pu . popN = 0 ; float * newnh = new float [ K * mydata . L ] , * newnhp = newnh ; int * aszp = mydata . asz ; for ( i = 0 ; i < K * mydata . L ; i ++ , newnhp ++ ) * newnhp = 0 ; newnhp = newnh ; nhp = tpara . nh ; for ( i = 0 ; i < mydata . L ; i ++ , aszp ++ ) { for ( j = 0 ; j < tpara . maxK ; j ++ , newnhp ++ , nhp ++ ) ( * newnhp ) = ( * nhp ) ; newnhp += K - tpara . maxK ; tpara . param [ i ] . resize ( K , pu ) ; } delete tpara . nh ; tpara . nh = newnh ; for ( i = tpara . maxK ; i < K + 1 ; i ++ ) tpara . Q [ i ] = 0 ; tpara . maxK = K ; nhp = tpara . nh + st * tpara . maxK ; } bool * tlapse = tpara . lapse + idst * mydata . L + st ; vector < vector < PARAUNIT > > :: iterator pit = tpara . param . begin ( ) + st ; int ksz = mydata . L * mydata . ploidity , kkk = idst * mydata . ploidity , kkksz = mydata . ploidity ; for ( i = st ; i < ed ; i ++ , pp ++ , pp2 ++ , dp ++ , dp2 ++ , nhp += tpara . maxK , pit ++ , rr ++ , rr2 ++ , rn ++ , tlapse ++ ) { if ( * rr ) * rn += tpara . idw [ id ] ; int kk = 0 , kks = kkk ; for ( k = 0 ; k < ided - idst ; k ++ , kk += ksz , kks += kkksz ) { float b = tpara . indWeight [ idst + k ] ; * ( nhp + ( * pp ) ) += b ; _updatePara_add ( ( * pit ) [ * pp ] , ( int ) ( * ( dp + kk ) ) , ( float ) b , kks , ( int ) ( tpara . lapse != NULL && * ( tlapse + kk ) ) ) ; } if ( i == 0 || * rr ) tpara . Q [ * pp ] += tpara . idw [ id ] ; if ( mydata . ploidity == 2 ) { if ( * rr2 ) * rn += tpara . idw [ id ] ; kk = 0 ; kks = kkk + 1 ; for ( k = 0 ; k < ided - idst ; k ++ , kk += ksz , kks += kkksz ) { float b = tpara . indWeight [ idst + k ] ; * ( nhp + ( * pp2 ) ) += b ; _updatePara_add ( ( * pit ) [ * pp2 ] , ( int ) ( * ( dp2 + kk ) ) , ( float ) b , kks , ( int ) ( tpara . lapse != NULL && * ( tlapse + kk ) ) ) ; } if ( i == 0 || * rr2 ) tpara . Q [ * pp2 ] += tpara . idw [ id ] ; } } }
This function is specific to categorical data
return maxK
nh: group size per position
Q: group size overall
r: recombination probability
asz: number of distinct values per position
This function is specific to categorical data return maxK nh : group size per position Q : group size overall r : recombination probability asz : number of distinct values per position
//printf("K=%d maxK=%d\n", K, tpara.maxK);fflush(stdout); //, vector<float>(*aszp, 0)); //	int flipi = 0, curflip = 0; //	int fsz = 0; //	if((int)tpara.flips.size() > 0) fsz = (int)tpara.flips[idd].size(); //      while(flipi < fsz && tpara.flips[idd][flipi] < i) flipi += 2; //	if(flipi < fsz && i == tpara.flips[idd][flipi]) //	{	if(curflip != tpara.flips[idd][flipi + 1]) //		{//	float const *tdp = dp; dp = dp2; dp2 = tdp; //			curflip = tpara.flips[idd][flipi + 1]; //		} //	} //i); //i);