Method calculateRatio

src/GenomeCopyNumber.cpp:1312–1602

}

int GenomeCopyNumber::calculateRatio( GenomeCopyNumber & controlCopyNumber, int degree, bool intercept) {


    int maximalNumberOfIterations = 300;
    int maximalNumberOfCopies = ploidy_*2;
    int realNumberOfIterations = maximalNumberOfIterations;

    int successfulFit = 0;

	if(isRatioLogged_) {
        intercept=1; degree=1;//because it is loglogscale
        vector <float> y; //y ~ a0x+a1
        vector <float> x;

        //fill x and y:
        vector<ChrCopyNumber>::iterator it;
        for ( it=chrCopyNumber_.begin() ; it != chrCopyNumber_.end(); it++ ) {
            if (! (it->getChromosome().find("X")!=string::npos || it->getChromosome().find("Y")!=string::npos)) {
                vector <float> controlcounts = controlCopyNumber.getChrCopyNumber(it->getChromosome()).getValues() ;
                //check that everything is all right:
                if (int(controlcounts.size())!=it->getLength()) {
                    cerr << "Possible Error: calculateMedianAround ()\n";
                }
                for (int i = 0; i< it->getLength(); i++) {
                    if (it->getValueAt(i)>0 && controlcounts[i]>0) {
                        x.push_back(log(controlcounts[i]));
                        y.push_back(log(it->getValueAt(i)));
                    }
                }
                controlcounts.clear();
            }
        }
        cout << "Initial guess for polynomial:\n";

        int nvars = degree; //1 if the fit is linear
		ap::real_2d_array xy;
		int npoints = x.size();
		xy.setlength(npoints,nvars+1);
		int pos = 0;
		for (int i = 0; i <(int)x.size(); i++) {
				xy(pos,degree) = y[i];
				xy(pos,degree-1) = x[i];;
				for (int j = degree-2; j>=0; j--) {
					xy(pos,j)=xy(pos,j+1)*x[i];
				}
				pos++;
		}
		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 linear regression, code: " << info <<"\n";
		}