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 
This function is specific to categorical data
return maxK
nh: group size per position
Q: group size overall
r: recombination probability
This function is specific to categorical data return maxK nh : group size per position Q : group size overall r : recombination probability
void tensorHMMbase::_removeInd(MYDATA const &mydata, TENSORPARA &tpara, int id, float wb, int st, int ed, bool clean)
{       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;
	float *rn = tpara.recombN + st;
        float *nhp = tpara.nh + st * tpara.maxK;
        int K = 0;
	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_remove((*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_remove((*pit)[*pp2], (int)(*(dp2 + kk)), (float)b, kks, (int)(tpara.lapse != NULL && *(tlapse+kk)));
			}
			if(i == 0 || (*rr2)) tpara.Q[*pp2]-=tpara.idw[id];
		}
        }
	if(clean)
	{	K = tpara0.popn0;
		nhp = tpara.nh;
		for(i = 0; i < mydata.L; i++, nhp += tpara.maxK)
		{	for(j = tpara.maxK - 1; j >= K; j--) if(*(nhp + j) > NUMPRECISION) break;
			K = j + 1;
		}
		if(K < tpara.maxK)
		{	PARAUNIT pu;
			pu.lapseN = pu.popN = 0;
			float *newnh = new float[K * mydata.L], *newnhp = newnh;
	                nhp = tpara.nh;
        	        for(i = 0; i < mydata.L; i++)
                	{       for(j = 0; j < K; j++, newnhp++, nhp++)
                        	        (*newnhp) = (*nhp);
	                        nhp += tpara.maxK - K;
        	                tpara.param[i].resize(K);
                	}
	                delete tpara.nh;
        	        tpara.nh = newnh;
        	        tpara.maxK = K;
		}
        }
}
void tensorHMMbase :: _removeInd ( MYDATA const & mydata , TENSORPARA & tpara , int id , float wb , int st , int ed , bool clean ) { 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 ; float * rn = tpara . recombN + st ; float * nhp = tpara . nh + st * tpara . maxK ; int K = 0 ; 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_remove ( ( * 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_remove ( ( * pit ) [ * pp2 ] , ( int ) ( * ( dp2 + kk ) ) , ( float ) b , kks , ( int ) ( tpara . lapse != NULL && * ( tlapse + kk ) ) ) ; } if ( i == 0 || ( * rr2 ) ) tpara . Q [ * pp2 ] -= tpara . idw [ id ] ; } } if ( clean ) { K = tpara0 . popn0 ; nhp = tpara . nh ; for ( i = 0 ; i < mydata . L ; i ++ , nhp += tpara . maxK ) { for ( j = tpara . maxK - 1 ; j >= K ; j -- ) if ( * ( nhp + j ) > NUMPRECISION ) break ; K = j + 1 ; } if ( K < tpara . maxK ) { PARAUNIT pu ; pu . lapseN = pu . popN = 0 ; float * newnh = new float [ K * mydata . L ] , * newnhp = newnh ; nhp = tpara . nh ; for ( i = 0 ; i < mydata . L ; i ++ ) { for ( j = 0 ; j < K ; j ++ , newnhp ++ , nhp ++ ) ( * newnhp ) = ( * nhp ) ; nhp += tpara . maxK - K ; tpara . param [ i ] . resize ( K ) ; } delete tpara . nh ; tpara . nh = newnh ; tpara . maxK = K ; } } }
This function is specific to categorical data
return maxK
nh: group size per position
Q: group size overall
r: recombination probability
This function is specific to categorical data return maxK nh : group size per position Q : group size overall r : recombination probability
//	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]; //		} //	}		 //wb;//ploidity; //printf("%d:%d ", id, i);fflush(stdout); //i); //i); //, pu);