From 22d6020c609ea5e75af0f1720bceb01c1fddb095 Mon Sep 17 00:00:00 2001 From: Jonathan Jesuraj Date: Wed, 26 Oct 2016 20:49:52 -0400 Subject: [PATCH 01/15] updating with latest changes to solver and new branch --- mflib/CMakeLists.txt | 2 + mflib/MFCpgSolverGradient.cpp | 164 ++++++++++++++++++++++++++++++++++ mflib/MFCpgSolverGradient.hpp | 35 ++++++++ mflib/MFSolver.hpp | 4 +- 4 files changed, 203 insertions(+), 2 deletions(-) create mode 100644 mflib/MFCpgSolverGradient.cpp create mode 100644 mflib/MFCpgSolverGradient.hpp diff --git a/mflib/CMakeLists.txt b/mflib/CMakeLists.txt index c9f536a..f7d5e8b 100644 --- a/mflib/CMakeLists.txt +++ b/mflib/CMakeLists.txt @@ -14,6 +14,8 @@ ADD_LIBRARY(mflib MFRegionSolver.hpp MFCpgSolver.cpp MFCpgSolver.hpp + MFCpgSolverGradient.hpp + MFCpgSolverGradient.cpp MethylRead.cpp MethylRead.hpp MFRegionPrinter.cpp diff --git a/mflib/MFCpgSolverGradient.cpp b/mflib/MFCpgSolverGradient.cpp new file mode 100644 index 0000000..e9214b7 --- /dev/null +++ b/mflib/MFCpgSolverGradient.cpp @@ -0,0 +1,164 @@ +#include "MFCpgSolverGradient.hpp" +#include "MFCpgEstimator.hpp" +#include "MFSolver.hpp" + +namespace methylFlow { + + MFCpgSolverGradient::MFCpgSolverGradient(MFGraph* mfobj, const float length_mult) : MFSolver(mfobj), + estimator(mfobj, &std::cout, length_mult), + alpha_lambda(mfobj->get_graph()), + beta_lambda(mfobj->get_graph()) + { + estimator.computeNormalized(); + } + + MFCpgSolverGradient::~MFCpgSolverGradient() + { + } + + float MFCpgSolverGradient::score(const float lambda) + { + float obj = lp->primal(); + for (ListDigraph::InArcIt arc(mf->mfGraph, mf->sink); arc != INVALID; ++arc) { + obj -= lambda * lp->dual(rows[arc]); + } + return obj; + } + + int MFCpgSolverGradient::add_cols() + { + for (MFCpgEstimator::CpgMap::iterator it = estimator.normalized_map.begin(); + it != estimator.normalized_map.end(); ++it) { + int pos = it->first; + alpha_y[pos] = lp->addCol(); + lp->colLowerBound(alpha_y[pos], 0.); + lp->colUpperBound(alpha_y[pos], 1.); + + beta_y[pos] = lp->addCol(); + lp->colLowerBound(beta_y[pos], 0.); + lp->colUpperBound(beta_y[pos], 1.); + + alpha_m[pos] = lp->addCol(); + lp->colLowerBound(alpha_m[pos], 0.); + lp->colUpperBound(alpha_m[pos], 1.); + + beta_m[pos] = lp->addCol(); + lp->colLowerBound(beta_m[pos], 0.); + lp->colUpperBound(beta_m[pos], 1.); + + } + + // now add alpha and beta for lambda nodes + for (ListDigraph::InArcIt arc(mf->mfGraph, mf->sink); arc != INVALID; ++arc) { + ListDigraph::Node v = mf->mfGraph.source(arc); + + alpha_lambda[v] = lp->addCol(); + lp->colLowerBound(alpha_lambda[v], 0.); + lp->colUpperBound(alpha_lambda[v], 1.); + + beta_lambda[v] = lp->addCol(); + lp->colLowerBound(beta_lambda[v], 0.); + lp->colUpperBound(beta_lambda[v], 1.); + } + return 0; + } + + int MFCpgSolverGradient::make_deviance_objective(Lp::Expr &obj) + { + for (MFCpgEstimator::CpgMap::iterator it = estimator.normalized_map.begin(); + it != estimator.normalized_map.end(); ++it) { + int pos = it->first; + MFCpgEstimator::CpgEntry entry = it->second; + + obj += entry.Cov * ( beta_y[pos] - alpha_y[pos]); + obj += entry.Meth * ( beta_m[pos] - alpha_m[pos]); + } + return 0; + } + + int MFCpgSolverGradient::make_lambda_objective(const float lambda, Lp::Expr &obj) + { + const ListDigraph &mfGraph = mf->get_graph(); + + // add terms to objective for lambda nodes + for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { + ListDigraph::Node v = mfGraph.source(arc); + obj += lambda * (beta_lambda[v] - alpha_lambda[v]); + } + return 0; + } + + int MFCpgSolverGradient::add_constraints() + { + const ListDigraph &mfGraph = mf->get_graph(); + + //add sink constraints + for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { + ListDigraph::Node v = mfGraph.source(arc); + + rows[arc] = lp->addRow(scaled_length[arc] * beta_lambda[v] - + scaled_length[arc] * alpha_lambda[v] - nu[v] <= -CONSISTENCY_FACTOR); + } + + //add remaining constraints (if not childless) + for (IterableBoolMap::FalseIt v(mf->fake); + v != INVALID; ++v) { + if (mf->childless[v]) continue; + + for (ListDigraph::OutArcIt arc(mfGraph, v); arc != INVALID; ++arc) { + ListDigraph::Node u = mfGraph.target(arc); + if (u == INVALID) { + std::cerr << "error getting target from arc" << std::endl; + return -1; + } + + MethylRead *m = mf->read(v); + if (!m) continue; + + int rPos = m->start(); + Lp::Expr expr; + + for (std::vector::iterator it = m->cpgs.begin(); + it != m->cpgs.end(); ++it) { + MethylRead::CpgEntry entry = *it; + int loc = rPos + entry.offset - 1; + expr += beta_y[loc] - alpha_y[loc]; + if (entry.methyl) { + expr += beta_m[loc] - alpha_m[loc]; + } + } + expr += nu[u] - nu[v]; + rows[arc] = lp->addRow(expr <= -CONSISTENCY_FACTOR); + } + } + return 0; + } + + int MFCpgSolverGradient::modify_lambda_constraints(const float lambda) + { + const ListDigraph &mfGraph = mf->get_graph(); + + for (ListDigraph::InArcIt arc(mf->mfGraph, mf->sink); arc != INVALID; ++arc) { + ListDigraph::Node v = mfGraph.source(arc); + Lp::Row row = rows[arc]; + lp->row(row, -lambda * beta_lambda[v] - (-lambda * alpha_lambda[v]) - nu[v] <= -CONSISTENCY_FACTOR); + } + return 0; + } + + void MFCpgSolverGradient::print_primal() + { + for (MFCpgEstimator::CpgMap::iterator it = estimator.normalized_map.begin(); + it != estimator.normalized_map.end(); ++it) { + int pos = it->first; + std::cout << pos << ": ay=" << lp->primal(alpha_y[pos]); + std::cout << " by=" << lp->primal(beta_y[pos]); + std::cout << " am=" << lp->primal(alpha_m[pos]); + std::cout << " bm=" << lp->primal(beta_m[pos]); + + MFCpgEstimator::CpgEntry entry = it->second; + std::cout << " y=" << entry.Cov << " my=" << entry.Meth << std::endl; + } + print_nus(); + } +} // namespace methylFlow diff --git a/mflib/MFCpgSolverGradient.hpp b/mflib/MFCpgSolverGradient.hpp new file mode 100644 index 0000000..5d29003 --- /dev/null +++ b/mflib/MFCpgSolverGradient.hpp @@ -0,0 +1,35 @@ +#include "MFSolver.hpp" +#include "MFCpgEstimator.hpp" + +#ifndef MFCPGSOLVERGRADIENT_H +#define MFCPGSOLVERGRADIENT_H + +namespace methylFlow { + + class MFCpgSolverGradient : public MFSolver { + friend class MFCpgEstimator; + + public: + MFCpgSolverGradient(MFGraph *mfobj, const float length_mult); + ~MFCpgSolverGradient(); + + private: + MFCpgEstimator estimator; + std::map alpha_y; + std::map beta_y; + std::map alpha_m; + std::map beta_m; + ListDigraph::NodeMap alpha_lambda; + ListDigraph::NodeMap beta_lambda; + + protected: + float score(const float lambda); + int add_cols(); + int make_deviance_objective(Lp::Expr &obj); + int make_lambda_objective(const float lambda, Lp::Expr &obj); + int add_constraints(); + int modify_lambda_constraints(const float lambda); + void print_primal(); + }; +} // namespace methylFlow +#endif // MFCPGSOLVER_H diff --git a/mflib/MFSolver.hpp b/mflib/MFSolver.hpp index 7b06499..23a98bf 100644 --- a/mflib/MFSolver.hpp +++ b/mflib/MFSolver.hpp @@ -34,7 +34,7 @@ namespace methylFlow { // make the LP object // virtual LP creator - int make_lp(const float length_mult); + virtual int make_lp(const float length_mult)=0; virtual int add_cols() =0; virtual int make_deviance_objective(Lp::Expr &obj) =0; @@ -46,7 +46,7 @@ namespace methylFlow { // solve the optimization problem // lambda: penalty parameter - int solve_for_lambda(const float lambda); + virtual int solve_for_lambda(const float lambda)=0; // find the best lambda int search_lambda(const float epislon, float &best_lambda, const float scale_mult, const bool verbose, const bool verboseTime); From a7d62205a2ebaddfa439c90d1032be9c7ae0fa77 Mon Sep 17 00:00:00 2001 From: Jonathan Jesuraj Date: Mon, 31 Oct 2016 19:41:47 -0400 Subject: [PATCH 02/15] forgot to commit these --- mflib/MFCpgSolverGradient.cpp | 103 ++++++++++++++++++++++++++++++++++ mflib/MFCpgSolverGradient.hpp | 1 + 2 files changed, 104 insertions(+) diff --git a/mflib/MFCpgSolverGradient.cpp b/mflib/MFCpgSolverGradient.cpp index e9214b7..c1fb8b3 100644 --- a/mflib/MFCpgSolverGradient.cpp +++ b/mflib/MFCpgSolverGradient.cpp @@ -16,6 +16,109 @@ namespace methylFlow { { } + + //LF Calculation + private float LofF(int estimate_f){ + double sum = 0; + for(int i = 0; i <= l; i++){ + sum += abs(takeLog(ml, ul) - log(numerator(estimate)/denominator(estimate))) + } + return sum; + } + + private void numerator(int estimate_f){ + double sum = 0; + for(int i = 0; i <= v; i++){ + if (p(v) == m){ + for(int j = 0; j <= u; j++){ + sum += estimate_f/lvu; + } + } + } + } + + private void denominator(int estimate_f){ + double sum = 0; + for(int i = 0; i <= v; i++){ + if (p(v) == u){ + for(int j = 0; j <= u; j++){ + sum += estimate_f/lvu; + } + } + } + } + + private float takeLog(double ml, double ul){ + return log(ml/ul); + } + + //L Prime F Calculation + private float LPrimeOfF(int estimate_f){ + double sum = 0; + for(int i = 0; i <= l; i++){ + sum += signFunction(estimate_f)*(firstNumerator/firstDenominator - secondNumerator/secondDenominator) + } + } + + private double firstNumerator(int estimate_f){ + if(p(v) == u){ + return 1/lvu; + } + } + private double firstDenominator(int estimate_f){ + double sum = 0; + for(int i =0; i <= v; i++){ + if(pv==u){ + int sum2 = 0; + for(int j =0; j <=u; j++){ + sum2+= estimate_f * vu/lvu; + } + } + } + return log(sum); + } + + + private double secondNumerator(int estimate_f){ + if(p(v) == m){ + return 1/lvu; + } + } + + private double secondDenominator(int estimate_f){ + double sum = 0; + for(int i =0; i <= v; i++){ + if(pv==m){ + int sum2 = 0; + for(int j =0; j <=u; j++){ + sum2+= estimate_f/lvu; + } + } + } + return log(sum); + } + + + private int signFunction(int estimate_f){ + if(estimate_f < 0){ + return -1; + }else if (estimate_f > 0){ + return 1; + }else + return 0; + } + + virtual int solve_for_lambda(const float lambda){ + + double f0 = 0; + double fn = 0; + for(int i = 0; i < 5; i ++{ + fn = f0 - LPrimeOfF(f0); + f0 = fn; + } + return LofF(fn); + } + float MFCpgSolverGradient::score(const float lambda) { float obj = lp->primal(); diff --git a/mflib/MFCpgSolverGradient.hpp b/mflib/MFCpgSolverGradient.hpp index 5d29003..65bee6a 100644 --- a/mflib/MFCpgSolverGradient.hpp +++ b/mflib/MFCpgSolverGradient.hpp @@ -30,6 +30,7 @@ namespace methylFlow { int add_constraints(); int modify_lambda_constraints(const float lambda); void print_primal(); + virtual int solve_for_lambda(const float lambda); }; } // namespace methylFlow #endif // MFCPGSOLVER_H From b37d248cf5acef38012f2c8070e2e5f1eed34c71 Mon Sep 17 00:00:00 2001 From: Jonathan Jesuraj Date: Mon, 31 Oct 2016 20:41:09 -0400 Subject: [PATCH 03/15] using cpg entry --- mflib/MFCpgSolverGradient.cpp | 98 +++++++++++++++++++++++++---------- 1 file changed, 72 insertions(+), 26 deletions(-) diff --git a/mflib/MFCpgSolverGradient.cpp b/mflib/MFCpgSolverGradient.cpp index c1fb8b3..4515e2f 100644 --- a/mflib/MFCpgSolverGradient.cpp +++ b/mflib/MFCpgSolverGradient.cpp @@ -18,34 +18,63 @@ namespace methylFlow { //LF Calculation + //Loop through the CPGs (number of Ls?) private float LofF(int estimate_f){ double sum = 0; - for(int i = 0; i <= l; i++){ + /*for(int i = 0; i <= l; i++){ sum += abs(takeLog(ml, ul) - log(numerator(estimate)/denominator(estimate))) } - return sum; + return sum;*/ + + + for (std::vector::iterator it = m->cpgs.begin(); + it != m->cpgs.end(); ++it) { + int pos = it->first; + MFCpgEstimator::CpgEntry entry = it->second; + + float u = entry.Cov; + float m = entry.Meth; + float logmu = takeLog(m,u); + + sum += abs(takeLog(m, u) - log(numerator(estimate,m,u)/denominator(estimate,u))) + } + } - private void numerator(int estimate_f){ + + //What is v? Is it list digraph (like i have it here?) + //what is P(v) -- how do we calculate it? + // what is L(Vu), where do we calculate it? + private double numerator(int estimate_f, int m, int u) { double sum = 0; - for(int i = 0; i <= v; i++){ - if (p(v) == m){ + + const ListDigraph &mfGraph = mf->get_graph(); + for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { + ListDigraph::Node v = mfGraph.source(arc); + if (p(v) == m){ for(int j = 0; j <= u; j++){ sum += estimate_f/lvu; } } - } + } + + return sum; } - private void denominator(int estimate_f){ + private double denominator(int estimate_f, int u){ double sum = 0; - for(int i = 0; i <= v; i++){ - if (p(v) == u){ + + const ListDigraph &mfGraph = mf->get_graph(); + for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { + ListDigraph::Node v = mfGraph.source(arc); + if (p(v) == u){ for(int j = 0; j <= u; j++){ sum += estimate_f/lvu; } } - } + } + + return sum; } private float takeLog(double ml, double ul){ @@ -58,43 +87,60 @@ namespace methylFlow { for(int i = 0; i <= l; i++){ sum += signFunction(estimate_f)*(firstNumerator/firstDenominator - secondNumerator/secondDenominator) } + + for (std::vector::iterator it = m->cpgs.begin(); + it != m->cpgs.end(); ++it) { + int pos = it->first; + MFCpgEstimator::CpgEntry entry = it->second; + + float u = entry.Cov; + float m = entry.Meth; + + sum += signFunction(estimate_f) * (firstNumerator(u)/firstDenominator(estimate_f, u) - secondNumerator(m)/secondDenominator(estimate_f,m,u)); + + } + } - private double firstNumerator(int estimate_f){ + private double firstNumerator(int u){ if(p(v) == u){ return 1/lvu; } } - private double firstDenominator(int estimate_f){ + + //LVu, vu, pv + private double firstDenominator(int estimate_f, int u){ double sum = 0; - for(int i =0; i <= v; i++){ - if(pv==u){ - int sum2 = 0; - for(int j =0; j <=u; j++){ - sum2+= estimate_f * vu/lvu; + for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { + ListDigraph::Node v = mfGraph.source(arc); + if (p(v) == u){ + for(int j = 0; j <= u; j++){ + sum += estimate_f * vu /lvu; } } - } + } + return log(sum); } - private double secondNumerator(int estimate_f){ + private double secondNumerator(int m){ if(p(v) == m){ return 1/lvu; } } - private double secondDenominator(int estimate_f){ + private double secondDenominator(int estimate_f, int m, int u){ double sum = 0; - for(int i =0; i <= v; i++){ - if(pv==m){ - int sum2 = 0; - for(int j =0; j <=u; j++){ - sum2+= estimate_f/lvu; + for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { + ListDigraph::Node v = mfGraph.source(arc); + if (p(v) == m){ + for(int j = 0; j <= u; j++){ + sum += estimate_f * vu /lvu; } } - } + } + return log(sum); } From bc18172e6506bf862da789ce59b7065db0d992fb Mon Sep 17 00:00:00 2001 From: Jonathan Jesuraj Date: Tue, 1 Nov 2016 20:58:34 -0400 Subject: [PATCH 04/15] Cleaning up solver --- mflib/MFCpgSolverGradient.cpp | 111 +++++++--------------------------- 1 file changed, 21 insertions(+), 90 deletions(-) diff --git a/mflib/MFCpgSolverGradient.cpp b/mflib/MFCpgSolverGradient.cpp index 4515e2f..34fcf0d 100644 --- a/mflib/MFCpgSolverGradient.cpp +++ b/mflib/MFCpgSolverGradient.cpp @@ -17,76 +17,9 @@ namespace methylFlow { } - //LF Calculation - //Loop through the CPGs (number of Ls?) - private float LofF(int estimate_f){ - double sum = 0; - /*for(int i = 0; i <= l; i++){ - sum += abs(takeLog(ml, ul) - log(numerator(estimate)/denominator(estimate))) - } - return sum;*/ - - - for (std::vector::iterator it = m->cpgs.begin(); - it != m->cpgs.end(); ++it) { - int pos = it->first; - MFCpgEstimator::CpgEntry entry = it->second; - - float u = entry.Cov; - float m = entry.Meth; - float logmu = takeLog(m,u); - - sum += abs(takeLog(m, u) - log(numerator(estimate,m,u)/denominator(estimate,u))) - } - - } - - - //What is v? Is it list digraph (like i have it here?) - //what is P(v) -- how do we calculate it? - // what is L(Vu), where do we calculate it? - private double numerator(int estimate_f, int m, int u) { - double sum = 0; - - const ListDigraph &mfGraph = mf->get_graph(); - for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { - ListDigraph::Node v = mfGraph.source(arc); - if (p(v) == m){ - for(int j = 0; j <= u; j++){ - sum += estimate_f/lvu; - } - } - } - - return sum; - } - - private double denominator(int estimate_f, int u){ - double sum = 0; - - const ListDigraph &mfGraph = mf->get_graph(); - for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { - ListDigraph::Node v = mfGraph.source(arc); - if (p(v) == u){ - for(int j = 0; j <= u; j++){ - sum += estimate_f/lvu; - } - } - } - - return sum; - } - - private float takeLog(double ml, double ul){ - return log(ml/ul); - } - //L Prime F Calculation private float LPrimeOfF(int estimate_f){ double sum = 0; - for(int i = 0; i <= l; i++){ - sum += signFunction(estimate_f)*(firstNumerator/firstDenominator - secondNumerator/secondDenominator) - } for (std::vector::iterator it = m->cpgs.begin(); it != m->cpgs.end(); ++it) { @@ -94,54 +27,51 @@ namespace methylFlow { MFCpgEstimator::CpgEntry entry = it->second; float u = entry.Cov; - float m = entry.Meth; + float m = entry.Meth - sum += signFunction(estimate_f) * (firstNumerator(u)/firstDenominator(estimate_f, u) - secondNumerator(m)/secondDenominator(estimate_f,m,u)); + sum += signFunction(estimate_f) * (-(firstNumerator(m)/MLofF(estimate_f) + secondNumerator(m)/ULofF(f))); } } - private double firstNumerator(int u){ - if(p(v) == u){ + private double firstNumerator(int m){ + if(p(v) == m){ return 1/lvu; } } - //LVu, vu, pv - private double firstDenominator(int estimate_f, int u){ + private int MLofF(double estimate_f){ double sum = 0; for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { ListDigraph::Node v = mfGraph.source(arc); - if (p(v) == u){ + if (p(v) == m){ for(int j = 0; j <= u; j++){ - sum += estimate_f * vu /lvu; + sum += estimate_f /lvu; } } } - - return log(sum); - } - - - private double secondNumerator(int m){ - if(p(v) == m){ - return 1/lvu; - } + return sum; } - private double secondDenominator(int estimate_f, int m, int u){ + private int ULofF(double estimate_f){ double sum = 0; for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { ListDigraph::Node v = mfGraph.source(arc); - if (p(v) == m){ + if (p(v) == u){ for(int j = 0; j <= u; j++){ - sum += estimate_f * vu /lvu; + sum += estimate_f /lvu; } } } + } + + - return log(sum); + private double secondNumerator(int u){ + if(p(v) == u){ + return 1/lvu; + } } @@ -158,11 +88,12 @@ namespace methylFlow { double f0 = 0; double fn = 0; + double gamma = 0.01; for(int i = 0; i < 5; i ++{ - fn = f0 - LPrimeOfF(f0); f0 = fn; + fn += -gamma*LPrimeOfF(f0); } - return LofF(fn); + return fn; } float MFCpgSolverGradient::score(const float lambda) From fe5890efc823a9498cc57978cedb627844a6181e Mon Sep 17 00:00:00 2001 From: jonjesuraj Date: Sun, 6 Nov 2016 13:39:03 -0500 Subject: [PATCH 05/15] Trying to get it to compile in linux --- mflib/MFSolver.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mflib/MFSolver.hpp b/mflib/MFSolver.hpp index 23a98bf..fc95a19 100644 --- a/mflib/MFSolver.hpp +++ b/mflib/MFSolver.hpp @@ -34,7 +34,7 @@ namespace methylFlow { // make the LP object // virtual LP creator - virtual int make_lp(const float length_mult)=0; + int make_lp(const float length_mult); virtual int add_cols() =0; virtual int make_deviance_objective(Lp::Expr &obj) =0; @@ -46,7 +46,7 @@ namespace methylFlow { // solve the optimization problem // lambda: penalty parameter - virtual int solve_for_lambda(const float lambda)=0; + int solve_for_lambda(const float lambda); // find the best lambda int search_lambda(const float epislon, float &best_lambda, const float scale_mult, const bool verbose, const bool verboseTime); From 0705f99f71d66346c7aeba263e180723be0c26c1 Mon Sep 17 00:00:00 2001 From: jonjesuraj Date: Sun, 6 Nov 2016 17:30:38 -0500 Subject: [PATCH 06/15] Fixing the build in linux --- mflib/CMakeLists.txt | 4 +- mflib/MFCpgEstimator.cpp | 1 + mflib/MFCpgEstimator.hpp | 8 +- mflib/MFCpgSolverGradient.cpp | 39 +++++--- mflib/MFCpgSolverGradient.hpp | 1 - mflib/MFCpgSolverSuccessive.cpp | 164 ++++++++++++++++++++++++++++++++ mflib/MFCpgSolverSuccessive.hpp | 35 +++++++ mflib/MFGraph.hpp | 2 +- mflib/MFGraph_solve.cpp | 1 + mflib/MethylRead.hpp | 2 + 10 files changed, 237 insertions(+), 20 deletions(-) create mode 100644 mflib/MFCpgSolverSuccessive.cpp create mode 100644 mflib/MFCpgSolverSuccessive.hpp diff --git a/mflib/CMakeLists.txt b/mflib/CMakeLists.txt index f7d5e8b..36a08d6 100644 --- a/mflib/CMakeLists.txt +++ b/mflib/CMakeLists.txt @@ -14,8 +14,8 @@ ADD_LIBRARY(mflib MFRegionSolver.hpp MFCpgSolver.cpp MFCpgSolver.hpp - MFCpgSolverGradient.hpp - MFCpgSolverGradient.cpp + MFCpgSolverSuccessive.hpp + MFCpgSolverSuccessive.cpp MethylRead.cpp MethylRead.hpp MFRegionPrinter.cpp diff --git a/mflib/MFCpgEstimator.cpp b/mflib/MFCpgEstimator.cpp index c88a3ba..13212a6 100644 --- a/mflib/MFCpgEstimator.cpp +++ b/mflib/MFCpgEstimator.cpp @@ -1,5 +1,6 @@ #include "MFCpgEstimator.hpp" #include "MFCpgSolver.hpp" +#include "MFCpgSolverSuccessive.hpp" #include diff --git a/mflib/MFCpgEstimator.hpp b/mflib/MFCpgEstimator.hpp index 97db57e..ac3b511 100644 --- a/mflib/MFCpgEstimator.hpp +++ b/mflib/MFCpgEstimator.hpp @@ -7,11 +7,12 @@ using namespace lemon; namespace methylFlow { class MFCpgSolver; - + class MFCpgSolverSuccessive; + class MFCpgEstimator { friend class MFCpgSolver; friend class MFGraph; - + friend class MFCpgSolverSuccessive; public: MFCpgEstimator(MFGraph *graph, std::ostream * ostream, const float scale_mult); ~MFCpgEstimator(); @@ -33,7 +34,8 @@ namespace methylFlow { class CpgEntry { friend class MFCpgEstimator; friend class MFCpgSolver; - public: + friend class MFCpgSolverSuccessive; + public: CpgEntry(T cov, T meth) : Cov(cov), Meth(meth), Beta(0) {}; CpgEntry() : Cov(0), Meth(0), Beta(0) {}; protected: diff --git a/mflib/MFCpgSolverGradient.cpp b/mflib/MFCpgSolverGradient.cpp index 34fcf0d..9e2bacd 100644 --- a/mflib/MFCpgSolverGradient.cpp +++ b/mflib/MFCpgSolverGradient.cpp @@ -18,8 +18,10 @@ namespace methylFlow { //L Prime F Calculation - private float LPrimeOfF(int estimate_f){ - double sum = 0; +/* + float LPrimeOfF(int estimate_f){ + + double sum = 0; for (std::vector::iterator it = m->cpgs.begin(); it != m->cpgs.end(); ++it) { @@ -31,18 +33,23 @@ namespace methylFlow { sum += signFunction(estimate_f) * (-(firstNumerator(m)/MLofF(estimate_f) + secondNumerator(m)/ULofF(f))); - } + }/ + +return 0; } - private double firstNumerator(int m){ + double firstNumerator(int m){ + /* if(p(v) == m){ return 1/lvu; } +return 0; } - private int MLofF(double estimate_f){ - double sum = 0; + int MLofF(double estimate_f){ +/* + double sum = 0; for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { ListDigraph::Node v = mfGraph.source(arc); if (p(v) == m){ @@ -52,10 +59,13 @@ namespace methylFlow { } } return sum; + +return 0; } - private int ULofF(double estimate_f){ - double sum = 0; + int ULofF(double estimate_f){ +/* + double sum = 0; for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { ListDigraph::Node v = mfGraph.source(arc); if (p(v) == u){ @@ -64,18 +74,21 @@ namespace methylFlow { } } } +return 0; } - private double secondNumerator(int u){ + double secondNumerator(int u){ + /* if(p(v) == u){ return 1/lvu; } + return 0; } - private int signFunction(int estimate_f){ + int signFunction(int estimate_f){ if(estimate_f < 0){ return -1; }else if (estimate_f > 0){ @@ -83,13 +96,13 @@ namespace methylFlow { }else return 0; } - - virtual int solve_for_lambda(const float lambda){ +*/ + int solve_for_lambda(const float lambda){ double f0 = 0; double fn = 0; double gamma = 0.01; - for(int i = 0; i < 5; i ++{ + for(int i = 0; i < 5; i ++){ f0 = fn; fn += -gamma*LPrimeOfF(f0); } diff --git a/mflib/MFCpgSolverGradient.hpp b/mflib/MFCpgSolverGradient.hpp index 65bee6a..5d29003 100644 --- a/mflib/MFCpgSolverGradient.hpp +++ b/mflib/MFCpgSolverGradient.hpp @@ -30,7 +30,6 @@ namespace methylFlow { int add_constraints(); int modify_lambda_constraints(const float lambda); void print_primal(); - virtual int solve_for_lambda(const float lambda); }; } // namespace methylFlow #endif // MFCPGSOLVER_H diff --git a/mflib/MFCpgSolverSuccessive.cpp b/mflib/MFCpgSolverSuccessive.cpp new file mode 100644 index 0000000..8aa54ec --- /dev/null +++ b/mflib/MFCpgSolverSuccessive.cpp @@ -0,0 +1,164 @@ +#include "MFCpgSolverSuccessive.hpp" +#include "MFCpgEstimator.hpp" +#include "MFSolver.hpp" + +namespace methylFlow { + + MFCpgSolverSuccessive::MFCpgSolverSuccessive(MFGraph* mfobj, const float length_mult) : MFSolver(mfobj), + estimator(mfobj, &std::cout, length_mult), + alpha_lambda(mfobj->get_graph()), + beta_lambda(mfobj->get_graph()) + { + estimator.computeNormalized(); + } + + MFCpgSolverSuccessive::~MFCpgSolverSuccessive() + { + } + + float MFCpgSolverSuccessive::score(const float lambda) + { + float obj = lp->primal(); + for (ListDigraph::InArcIt arc(mf->mfGraph, mf->sink); arc != INVALID; ++arc) { + obj -= lambda * lp->dual(rows[arc]); + } + return obj; + } + + int MFCpgSolverSuccessive::add_cols() + { + for (MFCpgEstimator::CpgMap::iterator it = estimator.normalized_map.begin(); + it != estimator.normalized_map.end(); ++it) { + int pos = it->first; + alpha_y[pos] = lp->addCol(); + lp->colLowerBound(alpha_y[pos], 0.); + lp->colUpperBound(alpha_y[pos], 1.); + + beta_y[pos] = lp->addCol(); + lp->colLowerBound(beta_y[pos], 0.); + lp->colUpperBound(beta_y[pos], 1.); + + alpha_m[pos] = lp->addCol(); + lp->colLowerBound(alpha_m[pos], 0.); + lp->colUpperBound(alpha_m[pos], 1.); + + beta_m[pos] = lp->addCol(); + lp->colLowerBound(beta_m[pos], 0.); + lp->colUpperBound(beta_m[pos], 1.); + + } + + // now add alpha and beta for lambda nodes + for (ListDigraph::InArcIt arc(mf->mfGraph, mf->sink); arc != INVALID; ++arc) { + ListDigraph::Node v = mf->mfGraph.source(arc); + + alpha_lambda[v] = lp->addCol(); + lp->colLowerBound(alpha_lambda[v], 0.); + lp->colUpperBound(alpha_lambda[v], 1.); + + beta_lambda[v] = lp->addCol(); + lp->colLowerBound(beta_lambda[v], 0.); + lp->colUpperBound(beta_lambda[v], 1.); + } + return 0; + } + + int MFCpgSolverSuccessive::make_deviance_objective(Lp::Expr &obj) + { + for (MFCpgEstimator::CpgMap::iterator it = estimator.normalized_map.begin(); + it != estimator.normalized_map.end(); ++it) { + int pos = it->first; + MFCpgEstimator::CpgEntry entry = it->second; + + obj += entry.Cov * ( beta_y[pos] - alpha_y[pos]); + obj += entry.Meth * ( beta_m[pos] - alpha_m[pos]); + } + return 0; + } + + int MFCpgSolverSuccessive::make_lambda_objective(const float lambda, Lp::Expr &obj) + { + const ListDigraph &mfGraph = mf->get_graph(); + + // add terms to objective for lambda nodes + for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { + ListDigraph::Node v = mfGraph.source(arc); + obj += lambda * (beta_lambda[v] - alpha_lambda[v]); + } + return 0; + } + + int MFCpgSolverSuccessive::add_constraints() + { + const ListDigraph &mfGraph = mf->get_graph(); + + //add sink constraints + for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { + ListDigraph::Node v = mfGraph.source(arc); + + rows[arc] = lp->addRow(scaled_length[arc] * beta_lambda[v] - + scaled_length[arc] * alpha_lambda[v] - nu[v] <= -CONSISTENCY_FACTOR); + } + + //add remaining constraints (if not childless) + for (IterableBoolMap::FalseIt v(mf->fake); + v != INVALID; ++v) { + if (mf->childless[v]) continue; + + for (ListDigraph::OutArcIt arc(mfGraph, v); arc != INVALID; ++arc) { + ListDigraph::Node u = mfGraph.target(arc); + if (u == INVALID) { + std::cerr << "error getting target from arc" << std::endl; + return -1; + } + + MethylRead *m = mf->read(v); + if (!m) continue; + + int rPos = m->start(); + Lp::Expr expr; + + for (std::vector::iterator it = m->cpgs.begin(); + it != m->cpgs.end(); ++it) { + MethylRead::CpgEntry entry = *it; + int loc = rPos + entry.offset - 1; + expr += beta_y[loc] - alpha_y[loc]; + if (entry.methyl) { + expr += beta_m[loc] - alpha_m[loc]; + } + } + expr += nu[u] - nu[v]; + rows[arc] = lp->addRow(expr <= -CONSISTENCY_FACTOR); + } + } + return 0; + } + + int MFCpgSolverSuccessive::modify_lambda_constraints(const float lambda) + { + const ListDigraph &mfGraph = mf->get_graph(); + + for (ListDigraph::InArcIt arc(mf->mfGraph, mf->sink); arc != INVALID; ++arc) { + ListDigraph::Node v = mfGraph.source(arc); + Lp::Row row = rows[arc]; + lp->row(row, -lambda * beta_lambda[v] - (-lambda * alpha_lambda[v]) - nu[v] <= -CONSISTENCY_FACTOR); + } + return 0; + } + + void MFCpgSolverSuccessive::print_primal() + { + for (MFCpgEstimator::CpgMap::iterator it = estimator.normalized_map.begin(); + it != estimator.normalized_map.end(); ++it) { + int pos = it->first; + std::cout << pos << ": ay=" << lp->primal(alpha_y[pos]); + std::cout << " by=" << lp->primal(beta_y[pos]); + std::cout << " am=" << lp->primal(alpha_m[pos]); + std::cout << " bm=" << lp->primal(beta_m[pos]); + + MFCpgEstimator::CpgEntry entry = it->second; + std::cout << " y=" << entry.Cov << " my=" << entry.Meth << std::endl; + } + print_nus(); + } +} // namespace methylFlow diff --git a/mflib/MFCpgSolverSuccessive.hpp b/mflib/MFCpgSolverSuccessive.hpp new file mode 100644 index 0000000..4900179 --- /dev/null +++ b/mflib/MFCpgSolverSuccessive.hpp @@ -0,0 +1,35 @@ +#include "MFSolver.hpp" +#include "MFCpgEstimator.hpp" + +#ifndef MFCPGSOLVERSUCCESSIVE_H +#define MFCPGSOLVERSUCCESSIVE_H + +namespace methylFlow { + + class MFCpgSolverSuccessive : public MFSolver { + friend class MFCpgEstimator; + + public: + MFCpgSolverSuccessive(MFGraph *mfobj, const float length_mult); + ~MFCpgSolverSuccessive(); + + private: + MFCpgEstimator estimator; + std::map alpha_y; + std::map beta_y; + std::map alpha_m; + std::map beta_m; + ListDigraph::NodeMap alpha_lambda; + ListDigraph::NodeMap beta_lambda; + + protected: + float score(const float lambda); + int add_cols(); + int make_deviance_objective(Lp::Expr &obj); + int make_lambda_objective(const float lambda, Lp::Expr &obj); + int add_constraints(); + int modify_lambda_constraints(const float lambda); + void print_primal(); + }; +} // namespace methylFlow +#endif // MFCPGSOLVER_H diff --git a/mflib/MFGraph.hpp b/mflib/MFGraph.hpp index 101f6da..59e3d64 100644 --- a/mflib/MFGraph.hpp +++ b/mflib/MFGraph.hpp @@ -19,7 +19,7 @@ namespace methylFlow { friend class MFSolver; friend class MFRegionSolver; friend class MFCpgSolver; - + friend class MFCpgSolverSuccessive; public: MFGraph(); ~MFGraph(); diff --git a/mflib/MFGraph_solve.cpp b/mflib/MFGraph_solve.cpp index ae44f21..469f860 100644 --- a/mflib/MFGraph_solve.cpp +++ b/mflib/MFGraph_solve.cpp @@ -10,6 +10,7 @@ #include "MFGraph.hpp" #include "MFRegionSolver.hpp" #include "MFCpgSolver.hpp" +#include "MFCpgSolver.hpp" namespace methylFlow { void MFGraph::add_terminals() diff --git a/mflib/MethylRead.hpp b/mflib/MethylRead.hpp index 8582290..9f42aa5 100644 --- a/mflib/MethylRead.hpp +++ b/mflib/MethylRead.hpp @@ -14,6 +14,8 @@ namespace methylFlow { private: friend class MFCpgEstimator; friend class MFCpgSolver; + friend class MFCpgSolverSuccessive; + typedef struct {int offset; bool methyl; } CpgEntry; typedef struct {int start; int length; char indel; } CigarEntry; From a150b644bbdddd6669c4afafde80cbd0d7437f4a Mon Sep 17 00:00:00 2001 From: Jonathan Jesuraj Date: Sun, 6 Nov 2016 19:29:45 -0500 Subject: [PATCH 07/15] Adding back what was in gradient --- mflib/MFCpgSolverSuccessive.cpp | 74 ++++++++++++++++++++++++++++++++- 1 file changed, 73 insertions(+), 1 deletion(-) diff --git a/mflib/MFCpgSolverSuccessive.cpp b/mflib/MFCpgSolverSuccessive.cpp index 8aa54ec..a00a18d 100644 --- a/mflib/MFCpgSolverSuccessive.cpp +++ b/mflib/MFCpgSolverSuccessive.cpp @@ -14,7 +14,79 @@ namespace methylFlow { MFCpgSolverSuccessive::~MFCpgSolverSuccessive() { - } + + //L Prime F Calculation + float LPrimeOfF(int estimate_f){ + float sum = 0; + for (std::vector::iterator it = m->cpgs.begin(); + it != m->cpgs.end(); ++it) { + int pos = it->first; + MFCpgEstimator::CpgEntry entry = it->second; + float u = entry.Cov; + float m = entry.Meth + sum += signFunction(estimate_f) * (-(firstNumerator(m)/MLofF(estimate_f) + secondNumerator(m)/ULofF(f))); + } + return sum; + } + + double firstNumerator(int m){ + /* + if(p(v) == m){ + return 1/lvu; + }*/ + return 0; + } + + int MLofF(double estimate_f){ + /* + double sum = 0; + for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { + ListDigraph::Node v = mfGraph.source(arc); + if (p(v) == m){ + for(int j = 0; j <= u; j++){ + sum += estimate_f /lvu; + } + } + } + return sum; + */ + return 0; + } + + int ULofF(double estimate_f){ + /* + double sum = 0; + for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { + ListDigraph::Node v = mfGraph.source(arc); + if (p(v) == u){ + for(int j = 0; j <= u; j++){ + sum += estimate_f /lvu; + } + } + } + */ + return 0; + } + + + + double secondNumerator(int u){ + /* + if(p(v) == u){ + return 1/lvu; + }*/ + return 0; + } + + + int signFunction(int estimate_f){ + if(estimate_f < 0){ + return -1; + }else if (estimate_f > 0){ + return 1; + }else + return 0; + } float MFCpgSolverSuccessive::score(const float lambda) { From a52761e0d08caba6a0b8d0f20b6f43f992fbd857 Mon Sep 17 00:00:00 2001 From: jonjesuraj Date: Sun, 6 Nov 2016 19:35:08 -0500 Subject: [PATCH 08/15] Fixing --- mflib/MFCpgSolverSuccessive.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mflib/MFCpgSolverSuccessive.cpp b/mflib/MFCpgSolverSuccessive.cpp index a00a18d..14739ac 100644 --- a/mflib/MFCpgSolverSuccessive.cpp +++ b/mflib/MFCpgSolverSuccessive.cpp @@ -14,7 +14,7 @@ namespace methylFlow { MFCpgSolverSuccessive::~MFCpgSolverSuccessive() { - + } //L Prime F Calculation float LPrimeOfF(int estimate_f){ float sum = 0; @@ -23,8 +23,8 @@ namespace methylFlow { int pos = it->first; MFCpgEstimator::CpgEntry entry = it->second; float u = entry.Cov; - float m = entry.Meth - sum += signFunction(estimate_f) * (-(firstNumerator(m)/MLofF(estimate_f) + secondNumerator(m)/ULofF(f))); + float m = entry.Meth; + sum += signFunction(estimate_f) * (-(firstNumerator(m)/MLofF(estimate_f) + secondNumerator(m)/ULofF(estimate_f))); } return sum; } From 1ed7f9f311eb314bd7e9e78f7db810810679ab89 Mon Sep 17 00:00:00 2001 From: Jonathan Jesuraj Date: Sun, 6 Nov 2016 19:38:07 -0500 Subject: [PATCH 09/15] fixing some things to help with compilation --- mflib/MFCpgSolverSuccessive.cpp | 12 ++++++------ mflib/MFCpgSolverSuccessive.hpp | 6 ++++++ 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/mflib/MFCpgSolverSuccessive.cpp b/mflib/MFCpgSolverSuccessive.cpp index 14739ac..e6a6bae 100644 --- a/mflib/MFCpgSolverSuccessive.cpp +++ b/mflib/MFCpgSolverSuccessive.cpp @@ -16,7 +16,7 @@ namespace methylFlow { { } //L Prime F Calculation - float LPrimeOfF(int estimate_f){ + float MFCpgSolverSuccessive::LPrimeOfF(int estimate_f){ float sum = 0; for (std::vector::iterator it = m->cpgs.begin(); it != m->cpgs.end(); ++it) { @@ -29,7 +29,7 @@ namespace methylFlow { return sum; } - double firstNumerator(int m){ + double MFCpgSolverSuccessive::firstNumerator(int m){ /* if(p(v) == m){ return 1/lvu; @@ -37,7 +37,7 @@ namespace methylFlow { return 0; } - int MLofF(double estimate_f){ + int MFCpgSolverSuccessive::MLofF(double estimate_f){ /* double sum = 0; for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { @@ -53,7 +53,7 @@ namespace methylFlow { return 0; } - int ULofF(double estimate_f){ + int MFCpgSolverSuccessive::ULofF(double estimate_f){ /* double sum = 0; for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { @@ -70,7 +70,7 @@ namespace methylFlow { - double secondNumerator(int u){ + double MFCpgSolverSuccessive::secondNumerator(int u){ /* if(p(v) == u){ return 1/lvu; @@ -79,7 +79,7 @@ namespace methylFlow { } - int signFunction(int estimate_f){ + int MFCpgSolverSuccessive::signFunction(int estimate_f){ if(estimate_f < 0){ return -1; }else if (estimate_f > 0){ diff --git a/mflib/MFCpgSolverSuccessive.hpp b/mflib/MFCpgSolverSuccessive.hpp index 4900179..d106ec1 100644 --- a/mflib/MFCpgSolverSuccessive.hpp +++ b/mflib/MFCpgSolverSuccessive.hpp @@ -21,6 +21,12 @@ namespace methylFlow { std::map beta_m; ListDigraph::NodeMap alpha_lambda; ListDigraph::NodeMap beta_lambda; + float LPrimeOfF(int estimate_f); + double firstNumerator(int m); + int MLofF(double estimate_f); + int ULofF(double estimate_f); + double secondNumerator(int u); + int signFunction(int estimate_f); protected: float score(const float lambda); From f2317c6b42278b1880a2f633beb1c721ba829fba Mon Sep 17 00:00:00 2001 From: Jonathan Jesuraj Date: Sun, 6 Nov 2016 19:46:19 -0500 Subject: [PATCH 10/15] fixing another error build --- mflib/MFCpgSolverSuccessive.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/mflib/MFCpgSolverSuccessive.cpp b/mflib/MFCpgSolverSuccessive.cpp index e6a6bae..30fdc4c 100644 --- a/mflib/MFCpgSolverSuccessive.cpp +++ b/mflib/MFCpgSolverSuccessive.cpp @@ -18,9 +18,10 @@ namespace methylFlow { //L Prime F Calculation float MFCpgSolverSuccessive::LPrimeOfF(int estimate_f){ float sum = 0; - for (std::vector::iterator it = m->cpgs.begin(); - it != m->cpgs.end(); ++it) { - int pos = it->first; + + for (MFCpgEstimator::CpgMap::iterator it = estimator.normalized_map.begin(); + it != estimator.normalized_map.end(); ++it){ + int pos = it->first; MFCpgEstimator::CpgEntry entry = it->second; float u = entry.Cov; float m = entry.Meth; From 2a877ca54874ed3291aef9a1008a38eb0c0bb33d Mon Sep 17 00:00:00 2001 From: Jonathan Jesuraj Date: Sun, 6 Nov 2016 20:53:20 -0500 Subject: [PATCH 11/15] fixing a few things --- mflib/MFCpgSolverSuccessive.cpp | 43 ++++++++++++++++----------------- mflib/MFCpgSolverSuccessive.hpp | 10 ++++---- 2 files changed, 26 insertions(+), 27 deletions(-) diff --git a/mflib/MFCpgSolverSuccessive.cpp b/mflib/MFCpgSolverSuccessive.cpp index 30fdc4c..f32dea2 100644 --- a/mflib/MFCpgSolverSuccessive.cpp +++ b/mflib/MFCpgSolverSuccessive.cpp @@ -25,57 +25,56 @@ namespace methylFlow { MFCpgEstimator::CpgEntry entry = it->second; float u = entry.Cov; float m = entry.Meth; - sum += signFunction(estimate_f) * (-(firstNumerator(m)/MLofF(estimate_f) + secondNumerator(m)/ULofF(estimate_f))); + sum += signFunction(estimate_f) * (-(firstNumerator(pos, m)/MLofF(pos,estimate_f) + secondNumerator(pos,m)/ULofF(pos,estimate_f))); } return sum; } - double MFCpgSolverSuccessive::firstNumerator(int m){ - /* - if(p(v) == m){ - return 1/lvu; - }*/ + double MFCpgSolverSuccessive::firstNumerator(int ell, int m){ + + if(ell == m){ + return 1/ell; + } return 0; } - int MFCpgSolverSuccessive::MLofF(double estimate_f){ - /* + int MFCpgSolverSuccessive::MLofF(int ell, double estimate_f){ + double sum = 0; for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { ListDigraph::Node v = mfGraph.source(arc); - if (p(v) == m){ + if (ell == m){ for(int j = 0; j <= u; j++){ sum += estimate_f /lvu; } } } return sum; - */ - return 0; + } - int MFCpgSolverSuccessive::ULofF(double estimate_f){ - /* + int MFCpgSolverSuccessive::ULofF(int ell, double estimate_f){ + double sum = 0; for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { ListDigraph::Node v = mfGraph.source(arc); - if (p(v) == u){ + if (ell == u){ for(int j = 0; j <= u; j++){ - sum += estimate_f /lvu; + sum += estimate_f /ell; } } } - */ - return 0; + + return sum; } - double MFCpgSolverSuccessive::secondNumerator(int u){ - /* - if(p(v) == u){ - return 1/lvu; - }*/ + double MFCpgSolverSuccessive::secondNumerator(int ell, int u){ + + if(ell == u){ + return 1/ell; + } return 0; } diff --git a/mflib/MFCpgSolverSuccessive.hpp b/mflib/MFCpgSolverSuccessive.hpp index d106ec1..fdc99f0 100644 --- a/mflib/MFCpgSolverSuccessive.hpp +++ b/mflib/MFCpgSolverSuccessive.hpp @@ -22,11 +22,11 @@ namespace methylFlow { ListDigraph::NodeMap alpha_lambda; ListDigraph::NodeMap beta_lambda; float LPrimeOfF(int estimate_f); - double firstNumerator(int m); - int MLofF(double estimate_f); - int ULofF(double estimate_f); - double secondNumerator(int u); - int signFunction(int estimate_f); + double firstNumerator(int ell, int m); + int MLofF(int ell, double estimate_f); + int ULofF(int ell, double estimate_f); + double secondNumerator(int ell, int u); + int signFunction(int ell, int estimate_f); protected: float score(const float lambda); From 72cffa8359b59417f1adbd11c3558510f4922517 Mon Sep 17 00:00:00 2001 From: Jonathan Jesuraj Date: Sun, 6 Nov 2016 20:56:26 -0500 Subject: [PATCH 12/15] another fix --- mflib/MFCpgSolverSuccessive.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mflib/MFCpgSolverSuccessive.hpp b/mflib/MFCpgSolverSuccessive.hpp index fdc99f0..44ed2e7 100644 --- a/mflib/MFCpgSolverSuccessive.hpp +++ b/mflib/MFCpgSolverSuccessive.hpp @@ -26,7 +26,7 @@ namespace methylFlow { int MLofF(int ell, double estimate_f); int ULofF(int ell, double estimate_f); double secondNumerator(int ell, int u); - int signFunction(int ell, int estimate_f); + int signFunction(int estimate_f); protected: float score(const float lambda); From 873af7ff37a9a7b3c1cbffab5160a184e6c4de13 Mon Sep 17 00:00:00 2001 From: Jonathan Jesuraj Date: Sun, 6 Nov 2016 21:03:48 -0500 Subject: [PATCH 13/15] fixing more build issues --- mflib/MFCpgSolverSuccessive.cpp | 10 +++++----- mflib/MFCpgSolverSuccessive.hpp | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/mflib/MFCpgSolverSuccessive.cpp b/mflib/MFCpgSolverSuccessive.cpp index f32dea2..ae73a30 100644 --- a/mflib/MFCpgSolverSuccessive.cpp +++ b/mflib/MFCpgSolverSuccessive.cpp @@ -25,7 +25,7 @@ namespace methylFlow { MFCpgEstimator::CpgEntry entry = it->second; float u = entry.Cov; float m = entry.Meth; - sum += signFunction(estimate_f) * (-(firstNumerator(pos, m)/MLofF(pos,estimate_f) + secondNumerator(pos,m)/ULofF(pos,estimate_f))); + sum += signFunction(estimate_f) * (-(firstNumerator(pos, m)/MLofF(pos,estimate_f,u, m) + secondNumerator(pos,m)/ULofF(pos,estimate_f, u, m))); } return sum; } @@ -38,10 +38,10 @@ namespace methylFlow { return 0; } - int MFCpgSolverSuccessive::MLofF(int ell, double estimate_f){ + int MFCpgSolverSuccessive::MLofF(int ell, double estimate_f, int u, int m){ double sum = 0; - for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { + for (ListDigraph::InArcIt arc(mf->mfGraph, mf->sink); arc != INVALID; ++arc) { ListDigraph::Node v = mfGraph.source(arc); if (ell == m){ for(int j = 0; j <= u; j++){ @@ -53,10 +53,10 @@ namespace methylFlow { } - int MFCpgSolverSuccessive::ULofF(int ell, double estimate_f){ + int MFCpgSolverSuccessive::ULofF(int ell, double estimate_f, int u, int m){ double sum = 0; - for (ListDigraph::InArcIt arc(mfGraph, mf->get_sink()); arc != INVALID; ++arc) { + for (ListDigraph::InArcIt arc(mf->mfGraph, mf->sink); arc != INVALID; ++arc) { ListDigraph::Node v = mfGraph.source(arc); if (ell == u){ for(int j = 0; j <= u; j++){ diff --git a/mflib/MFCpgSolverSuccessive.hpp b/mflib/MFCpgSolverSuccessive.hpp index 44ed2e7..34fa4fd 100644 --- a/mflib/MFCpgSolverSuccessive.hpp +++ b/mflib/MFCpgSolverSuccessive.hpp @@ -23,8 +23,8 @@ namespace methylFlow { ListDigraph::NodeMap beta_lambda; float LPrimeOfF(int estimate_f); double firstNumerator(int ell, int m); - int MLofF(int ell, double estimate_f); - int ULofF(int ell, double estimate_f); + int MLofF(int ell, double estimate_f, int u, int m); + int ULofF(int ell, double estimate_f, int u, int m); double secondNumerator(int ell, int u); int signFunction(int estimate_f); From 5693cd957a60118e4d97b8fe3dfe1052f1cc762f Mon Sep 17 00:00:00 2001 From: Jonathan Jesuraj Date: Sun, 6 Nov 2016 21:10:53 -0500 Subject: [PATCH 14/15] cleaning up again --- mflib/MFCpgSolverSuccessive.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mflib/MFCpgSolverSuccessive.cpp b/mflib/MFCpgSolverSuccessive.cpp index ae73a30..49b849e 100644 --- a/mflib/MFCpgSolverSuccessive.cpp +++ b/mflib/MFCpgSolverSuccessive.cpp @@ -45,7 +45,7 @@ namespace methylFlow { ListDigraph::Node v = mfGraph.source(arc); if (ell == m){ for(int j = 0; j <= u; j++){ - sum += estimate_f /lvu; + sum += estimate_f /ell; } } } From fb7c3dffe6196bf9a0f53a893102bdef9d0fcd98 Mon Sep 17 00:00:00 2001 From: Jonathan Jesuraj Date: Sun, 6 Nov 2016 21:14:45 -0500 Subject: [PATCH 15/15] another fix (sorry for all the commits) --- mflib/MFCpgSolverSuccessive.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mflib/MFCpgSolverSuccessive.cpp b/mflib/MFCpgSolverSuccessive.cpp index 49b849e..a73eb61 100644 --- a/mflib/MFCpgSolverSuccessive.cpp +++ b/mflib/MFCpgSolverSuccessive.cpp @@ -42,7 +42,7 @@ namespace methylFlow { double sum = 0; for (ListDigraph::InArcIt arc(mf->mfGraph, mf->sink); arc != INVALID; ++arc) { - ListDigraph::Node v = mfGraph.source(arc); + ListDigraph::Node v = mf->mfGraph.source(arc); if (ell == m){ for(int j = 0; j <= u; j++){ sum += estimate_f /ell; @@ -57,7 +57,7 @@ namespace methylFlow { double sum = 0; for (ListDigraph::InArcIt arc(mf->mfGraph, mf->sink); arc != INVALID; ++arc) { - ListDigraph::Node v = mfGraph.source(arc); + ListDigraph::Node v = mf->mfGraph.source(arc); if (ell == u){ for(int j = 0; j <= u; j++){ sum += estimate_f /ell;