List of usage examples for org.apache.commons.math3.optimization GoalType MAXIMIZE
GoalType MAXIMIZE
To view the source code for org.apache.commons.math3.optimization GoalType MAXIMIZE.
Click Source Link
From source file:com.itemanalysis.psychometrics.kernel.LikelihoodCrossValidation.java
public double value() { BrentOptimizer brent = new BrentOptimizer(1e-10, 1e-14); UnivariatePointValuePair result;/*from w w w. j av a 2s . c o m*/ result = brent.optimize(400, this, GoalType.MAXIMIZE, 0.001, max); return result.getPoint(); }
From source file:edu.rice.cs.bioinfo.programs.phylonet.algos.network.NetworkLikelihoodFromGTT.java
protected double findOptimalBranchLength(final Network<Object> speciesNetwork, final Map<String, List<String>> species2alleles, final List distinctTrees, final List gtCorrespondence, final Set<String> singleAlleleSpecies) { boolean continueRounds = true; // keep trying to improve network for (NetNode<Object> node : speciesNetwork.dfs()) { for (NetNode<Object> parent : node.getParents()) { node.setParentDistance(parent, 1.0); if (node.isNetworkNode()) { node.setParentProbability(parent, 0.5); }/* w w w . j a v a2 s . c o m*/ } } Set<NetNode> node2ignoreForBL = findEdgeHavingNoBL(speciesNetwork, singleAlleleSpecies); double initalProb = computeProbabilityForCached(speciesNetwork, distinctTrees, species2alleles, gtCorrespondence); if (_printDetails) System.out.println(speciesNetwork.toString() + " : " + initalProb); final Container<Double> lnGtProbOfSpeciesNetwork = new Container<Double>(initalProb); // records the GTProb of the network at all times int roundIndex = 0; for (; roundIndex < _maxRounds && continueRounds; roundIndex++) { /* * Prepare a random ordering of network edge examinations each of which attempts to change a branch length or hybrid prob to improve the GTProb score. */ double lnGtProbLastRound = lnGtProbOfSpeciesNetwork.getContents(); List<Proc> assigmentActions = new ArrayList<Proc>(); // store adjustment commands here. Will execute them one by one later. for (final NetNode<Object> parent : edu.rice.cs.bioinfo.programs.phylonet.structs.network.util.Networks .postTraversal(speciesNetwork)) { for (final NetNode<Object> child : parent.getChildren()) { if (node2ignoreForBL.contains(child)) { continue; } assigmentActions.add(new Proc() { public void execute() { UnivariateFunction functionToOptimize = new UnivariateFunction() { public double value(double suggestedBranchLength) { double incumbentBranchLength = child.getParentDistance(parent); child.setParentDistance(parent, suggestedBranchLength); double lnProb = updateProbabilityForCached(speciesNetwork, distinctTrees, gtCorrespondence, child, parent); //System.out.println(speciesNetwork + ": " + lnProb); if (lnProb > lnGtProbOfSpeciesNetwork.getContents()) // did improve, keep change { lnGtProbOfSpeciesNetwork.setContents(lnProb); } else // didn't improve, roll back change { child.setParentDistance(parent, incumbentBranchLength); } return lnProb; } }; BrentOptimizer optimizer = new BrentOptimizer(_Brent1, _Brent2); // very small numbers so we control when brent stops, not brent. try { optimizer.optimize(_maxTryPerBranch, functionToOptimize, GoalType.MAXIMIZE, Double.MIN_VALUE, _maxBranchLength); } catch (TooManyEvaluationsException e) // _maxAssigmentAttemptsPerBranchParam exceeded { } updateProbabilityForCached(speciesNetwork, distinctTrees, gtCorrespondence, child, parent); if (_printDetails) System.out.println( speciesNetwork.toString() + " : " + lnGtProbOfSpeciesNetwork.getContents()); } }); } } for (final NetNode<Object> child : speciesNetwork.getNetworkNodes()) // find every hybrid node { Iterator<NetNode<Object>> hybridParents = child.getParents().iterator(); final NetNode hybridParent1 = hybridParents.next(); final NetNode hybridParent2 = hybridParents.next(); assigmentActions.add(new Proc() { public void execute() { UnivariateFunction functionToOptimize = new UnivariateFunction() { public double value(double suggestedProb) { double incumbentHybridProbParent1 = child.getParentProbability(hybridParent1); child.setParentProbability(hybridParent1, suggestedProb); child.setParentProbability(hybridParent2, 1.0 - suggestedProb); double lnProb = updateProbabilityForCached(speciesNetwork, distinctTrees, gtCorrespondence, child, null); //System.out.println(speciesNetwork + ": " + lnProb); if (lnProb > lnGtProbOfSpeciesNetwork.getContents()) // change improved GTProb, keep it { lnGtProbOfSpeciesNetwork.setContents(lnProb); } else // change did not improve, roll back { child.setParentProbability(hybridParent1, incumbentHybridProbParent1); child.setParentProbability(hybridParent2, 1.0 - incumbentHybridProbParent1); } return lnProb; } }; BrentOptimizer optimizer = new BrentOptimizer(_Brent1, _Brent2); // very small numbers so we control when brent stops, not brent. try { if (child.getName().equals("Y")) optimizer.optimize(_maxTryPerBranch, functionToOptimize, GoalType.MAXIMIZE, 0.6, 0.8); else optimizer.optimize(_maxTryPerBranch, functionToOptimize, GoalType.MAXIMIZE, 0, 1.0); } catch (TooManyEvaluationsException e) // _maxAssigmentAttemptsPerBranchParam exceeded { } updateProbabilityForCached(speciesNetwork, distinctTrees, gtCorrespondence, child, null); if (_printDetails) System.out.println( speciesNetwork.toString() + " : " + lnGtProbOfSpeciesNetwork.getContents()); } }); } // add hybrid probs to hybrid edges Collections.shuffle(assigmentActions); for (Proc assigment : assigmentActions) // for each change attempt, perform attempt { assigment.execute(); } if (_printDetails) { System.out.println("Round end ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"); System.out .println(speciesNetwork.toString() + "\n" + lnGtProbOfSpeciesNetwork.getContents() + "\n"); } if (((double) lnGtProbOfSpeciesNetwork.getContents()) == lnGtProbLastRound) // if no improvement was made wrt to last around, stop trying to find a better assignment { continueRounds = false; } else if (lnGtProbOfSpeciesNetwork.getContents() > lnGtProbLastRound) // improvement was made, ensure it is large enough wrt to improvement threshold to continue searching { double improvementPercentage = Math.pow(Math.E, (lnGtProbOfSpeciesNetwork.getContents() - lnGtProbLastRound)) - 1.0; // how much did we improve over last round if (improvementPercentage < _improvementThreshold) // improved, but not enough to keep searching { continueRounds = false; } } else { throw new IllegalStateException("Should never have decreased prob."); } } //System.out.println("\n" + lnGtProbOfSpeciesNetwork.getContents() + ": " + speciesNetwork); return lnGtProbOfSpeciesNetwork.getContents(); }
From source file:edu.rice.cs.bioinfo.programs.phylonet.commands.SearchBranchLengthsMaxGTProb.java
@Override protected String produceResult() { StringBuffer result = new StringBuffer(); final List<Tree> geneTrees = new ArrayList<Tree>(); final List<Integer> counter = new ArrayList<Integer>(); for (NetworkNonEmpty geneTree : _geneTrees) { String phylonetGeneTree = NetworkTransformer.toENewickTree(geneTree); NewickReader nr = new NewickReader(new StringReader(phylonetGeneTree)); STITree<Double> newtr = new STITree<Double>(true); try {// w ww .ja v a2s .c o m nr.readTree(newtr); } catch (Exception e) { errorDetected.execute(e.getMessage(), this._motivatingCommand.getLine(), this._motivatingCommand.getColumn()); } boolean found = false; int index = 0; for (Tree tr : geneTrees) { if (Trees.haveSameRootedTopology(tr, newtr)) { found = true; break; } index++; } if (found) { counter.set(index, counter.get(index) + 1); } else { geneTrees.add(newtr); counter.add(1); } } NetworkFactoryFromRNNetwork transformer = new NetworkFactoryFromRNNetwork(); final Network<Double> speciesNetwork = transformer.makeNetwork(_speciesNetwork); /* * Make the branch length of every edge initially 1 if no initial user value is specified. * Make the hybrid prob of every hybrid edge initially .5 if no initial user value is specified. */ for (NetNode<Double> parent : speciesNetwork.bfs()) { for (NetNode<Double> child : parent.getChildren()) { double initialBL = child.getParentDistance(parent); if (initialBL == NetNode.NO_DISTANCE || Double.isNaN(initialBL)) // no specification from user on branch length initialBL = 1.0; child.setParentDistance(parent, initialBL); if (child.getParentCount() == 2) { for (NetNode<Double> hybridParent : child.getParents()) { if (child.getParentProbability(hybridParent) == 1.0) // no specification from user on hybrid prob { child.setParentProbability(parent, 0.5); } } } } } /* * Try to assign branch lengths and hybrid probs to increase GTProb from the initial network. * Except branch lengths of leaf edges. They don't impact GTProb. */ // def: a round is an attempt to tweak each branch length and each hybrid prob. boolean continueRounds = true; // keep trying to improve network final Container<Double> lnGtProbOfSpeciesNetwork = new Container<Double>( _computeGTProbStrategy.execute(speciesNetwork, geneTrees, counter)); // records the GTProb of the network at all times int assigmentRound = 0; for (; assigmentRound < _assigmentRounds && continueRounds; assigmentRound++) { /* * Prepare a random ordering of network edge examinations each of which attempts to change a branch length or hybrid prob to improve the GTProb score. */ double lnGtProbLastRound = lnGtProbOfSpeciesNetwork.getContents(); List<Proc> assigmentActions = new ArrayList<Proc>(); // store adjustment commands here. Will execute them one by one later. // add branch length adjustments to the list for (final NetNode<Double> parent : speciesNetwork.bfs()) { for (final NetNode<Double> child : parent.getChildren()) { if (!parent.isRoot()) // jd tmp continue; if (child.isLeaf()) // leaf edge, skip continue; assigmentActions.add(new Proc() { public void execute() { UnivariateFunction functionToOptimize = new UnivariateFunction() { public double value(double suggestedBranchLength) { // brent suggests a new branch length double incumbentBranchLength = child.getParentDistance(parent); // mutate and see if it yields an improved network child.setParentDistance(parent, suggestedBranchLength); double lnProb = _computeGTProbStrategy.execute(speciesNetwork, geneTrees, counter); RnNewickPrinter<Double> rnNewickPrinter = new RnNewickPrinter<Double>(); StringWriter sw = new StringWriter(); rnNewickPrinter.print(speciesNetwork, sw); // String inferredNetwork = sw.toString(); // System.out.println(inferredNetwork + "\t" + lnProb); if (lnProb > lnGtProbOfSpeciesNetwork.getContents()) // did improve, keep change { lnGtProbOfSpeciesNetwork.setContents(lnProb); // System.out.println("(improved)"); } else // didn't improve, roll back change { child.setParentDistance(parent, incumbentBranchLength); } return lnProb; } }; BrentOptimizer optimizer = new BrentOptimizer(.000000000001, .0000000000000001); // very small numbers so we control when brent stops, not brent. try { optimizer.optimize(_maxAssigmentAttemptsPerBranchParam, functionToOptimize, GoalType.MAXIMIZE, Double.MIN_VALUE, _maxBranchLength); } catch (TooManyEvaluationsException e) // _maxAssigmentAttemptsPerBranchParam exceeded { } // System.out.println("-----------------------------------------------------------------------"); } }); } } // add hybrid probs to hybrid edges for (final NetNode<Double> child : speciesNetwork.bfs()) // find every hybrid node { if (child.isRoot()) // calling getParentNumber on root causes NPE. Bug workaround. continue; if (child.getParentCount() == 2) // hybrid node { Iterator<NetNode<Double>> hybridParents = child.getParents().iterator(); final NetNode<Double> hybridParent1 = hybridParents.next(); final NetNode<Double> hybridParent2 = hybridParents.next(); assigmentActions.add(new Proc() { public void execute() { UnivariateFunction functionToOptimize = new UnivariateFunction() { public double value(double suggestedProb) { double incumbentHybridProbParent1 = child.getParentProbability(hybridParent1); // try new pair of hybrid probs child.setParentProbability(hybridParent1, suggestedProb); child.setParentProbability(hybridParent2, 1.0 - suggestedProb); double lnProb = _computeGTProbStrategy.execute(speciesNetwork, geneTrees, counter); if (lnProb > lnGtProbOfSpeciesNetwork.getContents()) // change improved GTProb, keep it { lnGtProbOfSpeciesNetwork.setContents(lnProb); } else // change did not improve, roll back { child.setParentProbability(hybridParent1, incumbentHybridProbParent1); child.setParentProbability(hybridParent2, 1.0 - incumbentHybridProbParent1); } return lnProb; } }; BrentOptimizer optimizer = new BrentOptimizer(.000000000001, .0000000000000001); // very small numbers so we control when brent stops, not brent. try { optimizer.optimize(_maxAssigmentAttemptsPerBranchParam, functionToOptimize, GoalType.MAXIMIZE, 0, 1.0); } catch (TooManyEvaluationsException e) // _maxAssigmentAttemptsPerBranchParam exceeded { } } }); } } // Collections.shuffle(assigmentActions); // randomize the order we will try to adjust network edge properties for (Proc assigment : assigmentActions) // for each change attempt, perform attempt { assigment.execute(); } if (((double) lnGtProbOfSpeciesNetwork.getContents()) == lnGtProbLastRound) // if no improvement was made wrt to last around, stop trying to find a better assignment { continueRounds = false; } else if (lnGtProbOfSpeciesNetwork.getContents() > lnGtProbLastRound) // improvement was made, ensure it is large enough wrt to improvement threshold to continue searching { double improvementPercentage = Math.pow(Math.E, (lnGtProbOfSpeciesNetwork.getContents() - lnGtProbLastRound)) - 1.0; // how much did we improve over last round if (improvementPercentage < _improvementThreshold) // improved, but not enough to keep searching { continueRounds = false; } } else { throw new IllegalStateException("Should never have decreased prob."); } } RnNewickPrinter<Double> rnNewickPrinter = new RnNewickPrinter<Double>(); StringWriter sw = new StringWriter(); rnNewickPrinter.print(speciesNetwork, sw); String inferredNetwork = sw.toString(); this.richNewickGenerated(inferredNetwork); result.append( "\nTotal log probability: " + lnGtProbOfSpeciesNetwork.getContents() + ": " + inferredNetwork); return result.toString(); }
From source file:edu.rice.cs.bioinfo.programs.phylonet.algos.network.NetworkLikelihoodFromGTTBL.java
protected double findOptimalBranchLength(final Network<Object> speciesNetwork, final Map<String, List<String>> species2alleles, final List gts, final List gtCorrespondence, final Set<String> singleAlleleSpecies) { boolean continueRounds = true; // keep trying to improve network if (_pair2time == null) { computePairwiseCoalesceTime(gts, species2alleles); }/* w w w . jav a 2 s . co m*/ //System.out.println("\n"+speciesNetwork); final Map<NetNode, Double> node2constraints = new Hashtable<NetNode, Double>(); computeNodeHeightUpperbound(speciesNetwork, node2constraints); final Map<NetNode<Object>, Double> node2height = new Hashtable<NetNode<Object>, Double>(); initializeNetwork(speciesNetwork, node2constraints, node2height); double initialProb = computeProbability(speciesNetwork, gts, species2alleles, gtCorrespondence); final Container<Double> lnGtProbOfSpeciesNetwork = new Container<Double>(initialProb); // records the GTProb of the network at all times final Container<Map<NetNode<Object>, Double>> node2heightContainer = new Container<Map<NetNode<Object>, Double>>( node2height); int roundIndex = 0; for (; roundIndex < _maxRounds && continueRounds; roundIndex++) { double lnGtProbLastRound = lnGtProbOfSpeciesNetwork.getContents(); List<Proc> assigmentActions = new ArrayList<Proc>(); // store adjustment commands here. Will execute them one by one later. for (final NetNode<Object> child : speciesNetwork.getNetworkNodes()) // find every hybrid node { Iterator<NetNode<Object>> hybridParents = child.getParents().iterator(); final NetNode hybridParent1 = hybridParents.next(); final NetNode hybridParent2 = hybridParents.next(); assigmentActions.add(new Proc() { public void execute() { UnivariateFunction functionToOptimize = new UnivariateFunction() { public double value(double suggestedProb) { double incumbentHybridProbParent1 = child.getParentProbability(hybridParent1); child.setParentProbability(hybridParent1, suggestedProb); child.setParentProbability(hybridParent2, 1.0 - suggestedProb); double lnProb = computeProbability(speciesNetwork, gts, species2alleles, gtCorrespondence); if (lnProb > lnGtProbOfSpeciesNetwork.getContents()) // change improved GTProb, keep it { lnGtProbOfSpeciesNetwork.setContents(lnProb); } else // change did not improve, roll back { child.setParentProbability(hybridParent1, incumbentHybridProbParent1); child.setParentProbability(hybridParent2, 1.0 - incumbentHybridProbParent1); } return lnProb; } }; BrentOptimizer optimizer = new BrentOptimizer(_Brent1, _Brent2); // very small numbers so we control when brent stops, not brent. try { optimizer.optimize(_maxTryPerBranch, functionToOptimize, GoalType.MAXIMIZE, 0, 1.0); } catch (TooManyEvaluationsException e) // _maxAssigmentAttemptsPerBranchParam exceeded { } } }); } for (final NetNode<Object> node : Networks.postTraversal(speciesNetwork)) { if (node.isLeaf()) { continue; } assigmentActions.add(new Proc() { public void execute() { final Container<Double> minHeight = new Container<Double>(0.0); final Container<Double> maxHeight = new Container<Double>(Double.MAX_VALUE); for (NetNode<Object> child : node.getChildren()) { double childHeight = node2heightContainer.getContents().get(child); minHeight.setContents(Math.max(minHeight.getContents(), childHeight)); } double minParent = Double.MAX_VALUE; if (!node.isRoot()) { for (NetNode<Object> parent : node.getParents()) { double parentHeight = node2heightContainer.getContents().get(parent); minParent = Math.min(minParent, parentHeight); } } else { minParent = minHeight.getContents() + _maxBranchLength; } maxHeight.setContents(Math.min(minParent, node2constraints.get(node))); //System.out.println("\nChanging node " + node.getName() + " from " + minHeight.getContents() + " to " + maxHeight.getContents()); UnivariateFunction functionToOptimize = new UnivariateFunction() { public double value(double suggestedHeight) { // brent suggests a new branch length double incumbentHeight = node2heightContainer.getContents().get(node); for (NetNode<Object> child : node.getChildren()) { child.setParentDistance(node, suggestedHeight - node2heightContainer.getContents().get(child)); } if (!node.isRoot()) { for (NetNode<Object> parent : node.getParents()) { node.setParentDistance(parent, node2heightContainer.getContents().get(parent) - suggestedHeight); } } double lnProb = computeProbability(speciesNetwork, gts, species2alleles, gtCorrespondence); //System.out.print("suggest: "+ suggestedHeight + " " + lnProb + " vs. " + lnGtProbOfSpeciesNetwork.getContents() + ": "); if (lnProb > lnGtProbOfSpeciesNetwork.getContents()) // did improve, keep change { lnGtProbOfSpeciesNetwork.setContents(lnProb); node2heightContainer.getContents().put(node, suggestedHeight); //System.out.println( " better "); } else // didn't improve, roll back change { for (NetNode<Object> child : node.getChildren()) { child.setParentDistance(node, incumbentHeight - node2heightContainer.getContents().get(child)); } if (!node.isRoot()) { for (NetNode<Object> parent : node.getParents()) { node.setParentDistance(parent, node2heightContainer.getContents().get(parent) - incumbentHeight); } } //System.out.println( " worse "); } return lnProb; } }; BrentOptimizer optimizer = new BrentOptimizer(_Brent1, _Brent2); // very small numbers so we control when brent stops, not brent. try { optimizer.optimize(_maxTryPerBranch, functionToOptimize, GoalType.MAXIMIZE, minHeight.getContents(), maxHeight.getContents()); } catch (TooManyEvaluationsException e) // _maxAssigmentAttemptsPerBranchParam exceeded { } //System.out.println(network2String(speciesNetwork) + " " + lnGtProbOfSpeciesNetwork.getContents()); } }); } Collections.shuffle(assigmentActions); for (Proc assigment : assigmentActions) // for each change attempt, perform attempt { assigment.execute(); } if (((double) lnGtProbOfSpeciesNetwork.getContents()) == lnGtProbLastRound) // if no improvement was made wrt to last around, stop trying to find a better assignment { continueRounds = false; } else if (lnGtProbOfSpeciesNetwork.getContents() > lnGtProbLastRound) // improvement was made, ensure it is large enough wrt to improvement threshold to continue searching { double improvementPercentage = Math.pow(Math.E, (lnGtProbOfSpeciesNetwork.getContents() - lnGtProbLastRound)) - 1.0; // how much did we improve over last round //System.out.println(improvementPercentage + " vs. " + _improvementThreshold); if (improvementPercentage < _improvementThreshold) // improved, but not enough to keep searching { continueRounds = false; } } else { throw new IllegalStateException("Should never have decreased prob."); } } //System.out.print("\n" + lnGtProbOfSpeciesNetwork.getContents() + ": " + speciesNetwork); return lnGtProbOfSpeciesNetwork.getContents(); }
From source file:edu.rice.cs.bioinfo.programs.phylonet.algos.network.NetworkPseudoLikelihoodFromGTT.java
protected double findOptimalBranchLength(final Network<Object> speciesNetwork, final Map<String, List<String>> species2alleles, final List tripleFrequencies, final List gtCorrespondence, final Set<String> singleAlleleSpecies) { boolean continueRounds = true; // keep trying to improve network for (NetNode<Object> node : speciesNetwork.dfs()) { for (NetNode<Object> parent : node.getParents()) { node.setParentDistance(parent, 1.0); if (node.isNetworkNode()) { node.setParentProbability(parent, 0.5); }// www . j a v a 2 s . c o m } } Set<NetNode> node2ignoreForBL = findEdgeHavingNoBL(speciesNetwork); double initalProb = computeProbability(speciesNetwork, tripleFrequencies, species2alleles, gtCorrespondence); if (_printDetails) System.out.println(speciesNetwork.toString() + " : " + initalProb); final Container<Double> lnGtProbOfSpeciesNetwork = new Container<Double>(initalProb); // records the GTProb of the network at all times int roundIndex = 0; for (; roundIndex < _maxRounds && continueRounds; roundIndex++) { /* * Prepare a random ordering of network edge examinations each of which attempts to change a branch length or hybrid prob to improve the GTProb score. */ double lnGtProbLastRound = lnGtProbOfSpeciesNetwork.getContents(); List<Proc> assigmentActions = new ArrayList<Proc>(); // store adjustment commands here. Will execute them one by one later. for (final NetNode<Object> parent : edu.rice.cs.bioinfo.programs.phylonet.structs.network.util.Networks .postTraversal(speciesNetwork)) { for (final NetNode<Object> child : parent.getChildren()) { if (node2ignoreForBL.contains(child)) { continue; } assigmentActions.add(new Proc() { public void execute() { UnivariateFunction functionToOptimize = new UnivariateFunction() { public double value(double suggestedBranchLength) { double incumbentBranchLength = child.getParentDistance(parent); child.setParentDistance(parent, suggestedBranchLength); double lnProb = computeProbability(speciesNetwork, tripleFrequencies, species2alleles, gtCorrespondence); //System.out.println(speciesNetwork + ": " + lnProb); if (lnProb > lnGtProbOfSpeciesNetwork.getContents()) // did improve, keep change { lnGtProbOfSpeciesNetwork.setContents(lnProb); } else // didn't improve, roll back change { child.setParentDistance(parent, incumbentBranchLength); } return lnProb; } }; BrentOptimizer optimizer = new BrentOptimizer(_Brent1, _Brent2); // very small numbers so we control when brent stops, not brent. try { optimizer.optimize(_maxTryPerBranch, functionToOptimize, GoalType.MAXIMIZE, Double.MIN_VALUE, _maxBranchLength); } catch (TooManyEvaluationsException e) // _maxAssigmentAttemptsPerBranchParam exceeded { } if (_printDetails) System.out.println( speciesNetwork.toString() + " : " + lnGtProbOfSpeciesNetwork.getContents()); } }); } } for (final NetNode<Object> child : speciesNetwork.getNetworkNodes()) // find every hybrid node { Iterator<NetNode<Object>> hybridParents = child.getParents().iterator(); final NetNode hybridParent1 = hybridParents.next(); final NetNode hybridParent2 = hybridParents.next(); assigmentActions.add(new Proc() { public void execute() { UnivariateFunction functionToOptimize = new UnivariateFunction() { public double value(double suggestedProb) { double incumbentHybridProbParent1 = child.getParentProbability(hybridParent1); child.setParentProbability(hybridParent1, suggestedProb); child.setParentProbability(hybridParent2, 1.0 - suggestedProb); double lnProb = computeProbability(speciesNetwork, tripleFrequencies, species2alleles, gtCorrespondence); //System.out.println(speciesNetwork + ": " + lnProb); if (lnProb > lnGtProbOfSpeciesNetwork.getContents()) // change improved GTProb, keep it { lnGtProbOfSpeciesNetwork.setContents(lnProb); } else // change did not improve, roll back { child.setParentProbability(hybridParent1, incumbentHybridProbParent1); child.setParentProbability(hybridParent2, 1.0 - incumbentHybridProbParent1); } return lnProb; } }; BrentOptimizer optimizer = new BrentOptimizer(_Brent1, _Brent2); // very small numbers so we control when brent stops, not brent. try { optimizer.optimize(_maxTryPerBranch, functionToOptimize, GoalType.MAXIMIZE, 0, 1.0); } catch (TooManyEvaluationsException e) // _maxAssigmentAttemptsPerBranchParam exceeded { } if (_printDetails) System.out.println( speciesNetwork.toString() + " : " + lnGtProbOfSpeciesNetwork.getContents()); } }); } // add hybrid probs to hybrid edges Collections.shuffle(assigmentActions); for (Proc assigment : assigmentActions) // for each change attempt, perform attempt { assigment.execute(); } if (_printDetails) { System.out.println("Round end ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"); System.out .println(speciesNetwork.toString() + "\n" + lnGtProbOfSpeciesNetwork.getContents() + "\n"); } if (((double) lnGtProbOfSpeciesNetwork.getContents()) == lnGtProbLastRound) // if no improvement was made wrt to last around, stop trying to find a better assignment { continueRounds = false; } else if (lnGtProbOfSpeciesNetwork.getContents() > lnGtProbLastRound) // improvement was made, ensure it is large enough wrt to improvement threshold to continue searching { double improvementPercentage = Math.pow(Math.E, (lnGtProbOfSpeciesNetwork.getContents() - lnGtProbLastRound)) - 1.0; // how much did we improve over last round if (improvementPercentage < _improvementThreshold) // improved, but not enough to keep searching { continueRounds = false; } } else { throw new IllegalStateException("Should never have decreased prob."); } } //System.out.println(speciesNetwork + " " + lnGtProbOfSpeciesNetwork.getContents()); return lnGtProbOfSpeciesNetwork.getContents(); }
From source file:edu.rice.cs.bioinfo.programs.phylonet.algos.network.InferMLNetworkFromSequences.java
protected double findNonUltrametricOptimalBranchLength(final String[] gtTaxa, final Network<Object> speciesNetwork, final Map<String, String> allele2species, final List<Tuple<char[], Integer>> sequences, final RateModel rModel, final double theta) { boolean continueRounds = true; // keep trying to improve network for (NetNode<Object> node : speciesNetwork.dfs()) { for (NetNode<Object> parent : node.getParents()) { node.setParentDistance(parent, theta); if (node.isNetworkNode()) { //node.setParentDistance(parent,0.0); node.setParentProbability(parent, 0.5); }/*from w w w .jav a2 s.co m*/ } } //long start = System.currentTimeMillis(); double initalProb = computeProbability(gtTaxa, speciesNetwork, allele2species, sequences, rModel, theta); //System.out.println(initalProb); //System.out.print("\n"+(System.currentTimeMillis()-start)); //computeProbability(speciesNetwork, distinctTrees, species2alleles, nbTreeAndCountAndBinaryIDList); final Container<Double> lnGtProbOfSpeciesNetwork = new Container<Double>(initalProb); // records the GTProb of the network at all times //final Container<Integer> callCount = new Container<Integer>(0); int roundIndex = 0; for (; roundIndex < _maxRounds && continueRounds; roundIndex++) { /* * Prepare a random ordering of network edge examinations each of which attempts to change a branch length or hybrid prob to improve the GTProb score. */ double lnGtProbLastRound = lnGtProbOfSpeciesNetwork.getContents(); List<Proc> assigmentActions = new ArrayList<Proc>(); // store adjustment commands here. Will execute them one by one later. for (final NetNode<Object> parent : edu.rice.cs.bioinfo.programs.phylonet.structs.network.util.Networks .postTraversal(speciesNetwork)) { for (final NetNode<Object> child : parent.getChildren()) { assigmentActions.add(new Proc() { public void execute() { UnivariateFunction functionToOptimize = new UnivariateFunction() { public double value(double suggestedBranchLength) { //System.out.print(" l/"+child.getName()+" "); //callCount.setContents(callCount.getContents()+1); double incumbentBranchLength = child.getParentDistance(parent); child.setParentDistance(parent, suggestedBranchLength); double lnProb = computeProbability(gtTaxa, speciesNetwork, allele2species, sequences, rModel, theta); //System.out.println("Changing branch ("+parent.getName()+","+child.getName()+") to " + suggestedBranchLength); //System.out.println(network2String(speciesNetwork)+": " + lnProb); if (lnProb > lnGtProbOfSpeciesNetwork.getContents()) // did improve, keep change { lnGtProbOfSpeciesNetwork.setContents(lnProb); } else // didn't improve, roll back change { child.setParentDistance(parent, incumbentBranchLength); } return lnProb; } }; BrentOptimizer optimizer = new BrentOptimizer(_Brent1, _Brent2); // very small numbers so we control when brent stops, not brent. try { optimizer.optimize(_maxTryPerBranch, functionToOptimize, GoalType.MAXIMIZE, Double.MIN_VALUE, _maxBranchLength); } catch (TooManyEvaluationsException e) // _maxAssigmentAttemptsPerBranchParam exceeded { } } }); } } for (final NetNode<Object> child : speciesNetwork.getNetworkNodes()) // find every hybrid node { Iterator<NetNode<Object>> hybridParents = child.getParents().iterator(); final NetNode hybridParent1 = hybridParents.next(); final NetNode hybridParent2 = hybridParents.next(); assigmentActions.add(new Proc() { public void execute() { UnivariateFunction functionToOptimize = new UnivariateFunction() { public double value(double suggestedProb) { //callCount.setContents(callCount.getContents()+1); //System.out.print(" p/"+child.getName()+" "); double incumbentHybridProbParent1 = child.getParentProbability(hybridParent1); child.setParentProbability(hybridParent1, suggestedProb); child.setParentProbability(hybridParent2, 1.0 - suggestedProb); double lnProb = computeProbability(gtTaxa, speciesNetwork, allele2species, sequences, rModel, theta); //System.out.println("Changing node probability to "+ suggestedProb); //System.out.println(network2String(speciesNetwork)+": " + lnProb); //System.out.println(Math.abs(computeProbability(speciesNetwork, distinctTrees, species2alleles, nbTreeAndCountAndBinaryIDList) - lnProb)); if (lnProb > lnGtProbOfSpeciesNetwork.getContents()) // change improved GTProb, keep it { lnGtProbOfSpeciesNetwork.setContents(lnProb); } else // change did not improve, roll back { child.setParentProbability(hybridParent1, incumbentHybridProbParent1); child.setParentProbability(hybridParent2, 1.0 - incumbentHybridProbParent1); } return lnProb; } }; BrentOptimizer optimizer = new BrentOptimizer(_Brent1, _Brent2); // very small numbers so we control when brent stops, not brent. try { optimizer.optimize(_maxTryPerBranch, functionToOptimize, GoalType.MAXIMIZE, 0, 1.0); } catch (TooManyEvaluationsException e) // _maxAssigmentAttemptsPerBranchParam exceeded { } //System.out.println(network2String(speciesNetwork) + " : " + lnGtProbOfSpeciesNetwork.getContents()); } }); } // add hybrid probs to hybrid edges if (_seed == null) { //Collections.shuffle(assigmentActions); } for (Proc assigment : assigmentActions) // for each change attempt, perform attempt { assigment.execute(); } if (((double) lnGtProbOfSpeciesNetwork.getContents()) == lnGtProbLastRound) // if no improvement was made wrt to last around, stop trying to find a better assignment { continueRounds = false; } else if (lnGtProbOfSpeciesNetwork.getContents() > lnGtProbLastRound) // improvement was made, ensure it is large enough wrt to improvement threshold to continue searching { double improvementPercentage = Math.pow(Math.E, (lnGtProbOfSpeciesNetwork.getContents() - lnGtProbLastRound)) - 1.0; // how much did we improve over last round if (improvementPercentage < _improvementThreshold) // improved, but not enough to keep searching { continueRounds = false; } } else { throw new IllegalStateException("Should never have decreased prob."); } } //System.out.println(callCount.getContents()); //System.out.println(computeProbability(speciesNetwork, distinctTrees, species2alleles, nbTreeAndCountAndBinaryIDList) + " vs. " + lnGtProbOfSpeciesNetwork.getContents()); return lnGtProbOfSpeciesNetwork.getContents(); }