Method recalculateRatioUsingCG

src/GenomeCopyNumber.cpp:442–559

}

void GenomeCopyNumber::recalculateRatioUsingCG (int degree, bool intercept, float minExpectedGC, float maxExpectedGC) {

	if (degree > MAXDEGREE) {
        cerr << "polynomial degree should be < 10\n";
        exit (-1);
	}
	int maximalNumberOfCopies = ploidy_*2;
	float interval = float (0.01) ;

	//first guess about parameters
	const int npoints = degree+2;


	double around [MAXDEGREE+2];
	for (int i = 0; i<npoints; i++) {
		around[i] = minExpectedGC + (maxExpectedGC-minExpectedGC)/(npoints-1)*i; //0.55-0.35
	}  // for degree = 3 one will get : { 0.35, 0.40, 0.45, 0.5, 0.55 }; for 0.35 and 0.55


	double yValues [MAXDEGREE+2];
	for (int i = 0; i <npoints; i++)
		yValues[i] = calculateMedianRatioAround(interval, float(around[i]));
	int nvars = degree; //fit by cubic polynomial
	ap::real_2d_array xy;
	xy.setlength(npoints,nvars+1);
	for (int i = 0; i <npoints; i++) {
		xy(i,degree) = yValues[i];
		xy(i,degree-1) = around[i];
		for (int j = degree-2; j>=0; j--) {
			xy(i,j)=xy(i,j+1)*around[i];
		}
		/* this is equal to
		xy(i,0) = around[i]*around[i]*around[i];
		xy(i,1) = around[i]*around[i];
		xy(i,2) = around[i];
		xy(i,3) = yValues[i]; */
	}
	linearmodel lm;
	int info;
	lrreport ar;

    if  (intercept)
        lrbuild(xy,npoints,nvars,info,lm,ar);
	else
		lrbuildz(xy,npoints,nvars,info,lm,ar);
	if (info != 1) {
		cerr << "Error in the first linear regression (the first guess about parameters), code: " << info <<"\n";
	}
	ap::real_1d_array v;
    v.setlength(nvars+int(intercept));
	lrunpack(lm,v,nvars);
	//cout << v(0) << "\t" << v(1)<< "\t" <<v(2)<< "\t" <<v(3)<< "\n";
	double a[MAXDEGREE+1];
	for (int i = 0; i <degree; i++) {
			a[i] = v(i);
	} /* this is equal to
	double a0 = v(0);
	double a1 = v(1);