| 1634 | |
| 1635 | |
| 1636 | float runEMlog (const vector<float>& x,const vector<float>& y,double * a, int degree, int & realNumberOfIterations, int maximalNumberOfIterations,int ploidy, int maximalNumberOfCopies, bool intercept, float contamination) { |
| 1637 | |
| 1638 | float rmserror = -1; |
| 1639 | |
| 1640 | if (x.size() != y.size()) { |
| 1641 | cerr << "Error: object size is different"; |
| 1642 | return 0; |
| 1643 | } |
| 1644 | vector <int> cluster (x.size()); |
| 1645 | vector <float> res (maximalNumberOfCopies+1); |
| 1646 | |
| 1647 | int count = 0; |
| 1648 | |
| 1649 | for (count = 0; count < maximalNumberOfIterations; count++) { |
| 1650 | |
| 1651 | for (int i = 0; i <(int)x.size(); i++) { |
| 1652 | for (int j = 0; j <= maximalNumberOfCopies; j++) { |
| 1653 | float conLevel = log(1./10); |
| 1654 | if (j>=0) conLevel = log(float(j)/ploidy); |
| 1655 | res[j] = fabs(polynomial(x[i],a,1,degree)-y[i]+ conLevel); |
| 1656 | } |
| 1657 | cluster[i] = get_min_index(res); |
| 1658 | } |
| 1659 | int npoints = 0; |
| 1660 | for (int i = 0; i <(int)x.size(); i++) |
| 1661 | if (cluster[i] == ploidy) |
| 1662 | npoints++; |
| 1663 | int nvars = degree; //3 if fit by cubic polynomial |
| 1664 | ap::real_2d_array xy; |
| 1665 | xy.setlength(npoints,nvars+1); |
| 1666 | int pos = 0; |
| 1667 | for (int i = 0; i <(int)x.size(); i++) { |
| 1668 | if (cluster[i] == ploidy){ |
| 1669 | xy(pos,degree) = y[i]; |
| 1670 | xy(pos,degree-1) = x[i];; |
| 1671 | for (int j = degree-2; j>=0; j--) { |
| 1672 | xy(pos,j)=xy(pos,j+1)*x[i]; |
| 1673 | } |
| 1674 | //xy(pos,0) = x[i]*x[i]*x[i]; |
| 1675 | //xy(pos,1) = x[i]*x[i]; |
| 1676 | //xy(pos,2) = x[i]; |
| 1677 | //xy(pos,3) = y[i]; |
| 1678 | pos++; |
| 1679 | } |
| 1680 | } |
| 1681 | |
| 1682 | linearmodel lm; |
| 1683 | int info; |
| 1684 | lrreport ar; |
| 1685 | |
| 1686 | if (intercept) |
| 1687 | lrbuild(xy,npoints,nvars,info,lm,ar); |
| 1688 | else |
| 1689 | lrbuildz(xy,npoints,nvars,info,lm,ar); |
| 1690 | if (info != 1) { |
| 1691 | cerr << "Error in linear regression, code: " << info <<"\n"; |
| 1692 | break; |
| 1693 | } |
no test coverage detected