Function runEM_linear

src/myFunc.cpp:1395–1461

	return float(a0*x*x*x+a1*x*x+a2*x+a3);
}
float runEM_linear (const vector<float>& x,const vector<float>& y,double & a0,double & a1, int maximalNumberOfIterations,int ploidy, int maximalNumberOfCopies) {

	float rmserror = -1;

	if (x.size() != y.size()) {
		cerr << "Error: object size is different";
		return 0;
	}
	vector <int> cluster (x.size());
	vector <float> res (maximalNumberOfCopies+1);

	int count = 0;

	for (count = 0; count < maximalNumberOfIterations; count++) {

		for (int i = 0; i <(int)x.size(); i++) {
			for (int j = 0; j <= maximalNumberOfCopies; j++)
				res[j] = fabs(j*a0/ploidy*x[i]+ j*a1/ploidy -y[i]);
			cluster[i] = get_min_index(res);
		}
		int npoints = 0;
		for (int i = 0; i <(int)x.size(); i++)
			if (cluster[i] == ploidy)
				npoints++;
		int nvars = 1; //fit by linear function
		ap::real_2d_array xy;
		xy.setlength(npoints,nvars+1);
		int pos = 0;
		for (int i = 0; i <(int)x.size(); i++) {
			if (cluster[i] == ploidy){
				xy(pos,0) = x[i];
				xy(pos,1) = y[i];
				pos++;
			}
		}

		linearmodel lm;
		int info;
		lrreport ar;
		lrbuildz(xy,npoints,nvars,info,lm,ar);
		if (info != 1) {
			cerr << "Error in linear regression, code: " << info <<"\n";
			break;
		}
		rmserror = float(ar.rmserror);
		ap::real_1d_array v;
		v.setlength(nvars+1);
		lrunpack(lm,v,nvars);
		//cout << v(0);  //x
		//cout << v(1);  //x^3
		//cout << v(2);  //x^2
		//cout << v(3);  //intercept

		cout << a0 - v(0) << "\t"; //x
		//cout << a1 - v(1) << "\n"; //intercept
		//if (fabs(a0 - v(0)) + fabs(a1 - v(1)) == 0) {
		if (fabs(a0 - v(0)) == 0) {
			break;