Smooth probability function in MC optimization

lafita · lafita · commit e95ffa3fcae5 · 2015-06-26T15:11:18.000+02:00
diff --git a/biojava-structure-gui/src/main/java/demo/DemoMultipleMC.java b/biojava-structure-gui/src/main/java/demo/DemoMultipleMC.java
@@ -32,7 +32,7 @@ public static void main(String[] args) throws IOException, StructureException, I
 		//Protein Kinases (CEMC paper)
 		//List<String> names = Arrays.asList("1cdk.A", "1cja.A", "1csn", "1b6c.B", "1ir3.A", "1fgk.A", "1byg.A", "1hck", "1blx.A", "3erk", "1bmk.A", "1kob.A", "1tki.A", "1phk", "1a06");
 		//DHFR (Gerstein 1998 paper)
-		List<String> names = Arrays.asList("d1dhfa_", "8dfr", "d4dfra_", "3dfr");
+		//List<String> names = Arrays.asList("d1dhfa_", "8dfr", "d4dfra_", "3dfr");
 		//TIM barrels (MUSTA paper)
 		//List<String> names = Arrays.asList("1tim.A", "1vzw", "1nsj", "3tha.A", "4enl", "2mnr", "7tim.A", "1tml", "1btc", "a1piia1", "6xia", "5rub.A", "2taa.B");
 		//Calcium Binding (MUSTA paper)
@@ -44,7 +44,7 @@ public static void main(String[] args) throws IOException, StructureException, I
 		//GPCRs
 		//List<String> names = Arrays.asList("2z73.A", "1u19.A", "4ug2.A", "4xt3", "4or2.A", "3odu.A");
 		//Immunoglobulins (MAMMOTH paper)
-		//List<String> names = Arrays.asList("2hla.B", "3hla.B", "1cd8", "2rhe", "1tlk", "1ten", "1ttf");
+		List<String> names = Arrays.asList("2hla.B", "3hla.B", "1cd8", "2rhe", "1tlk", "1ten", "1ttf");
 		//Globins (MAMMOTH, POSA, Gerstein&Levitt and MUSTA papers)
 		//List<String> names = Arrays.asList("1mbc", "1hlb", "1thb.A", "1ith.A", "1idr.A", "1dlw", "1kr7.A", "1ew6.A", "1it2.A", "1eco", "3sdh.A", "1cg5.B", "1fhj.B", "1ird.A", "1mba", "2gdm", "1b0b", "1h97.A", "1ash.A", "1jl7.A");
 		//Rossman-Fold (POSA paper)
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/mc/MultipleAlignmentOptimizerMC.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/mc/MultipleAlignmentOptimizerMC.java
@@ -38,15 +38,15 @@
  */
 public class MultipleAlignmentOptimizerMC implements Callable<MultipleAlignment> {
 	
-	private static final boolean debug = false;  //Prints the optimization moves and saves a file of the history
+	private static final boolean debug = true;  //Prints the optimization moves and saves a file of the history
 	private Random rnd;
 	private MultipleSuperimposer imposer;
 	
 	//Optimization parameters
 	private int Rmin;   	//Minimum number of aligned structures without a gap (33% of initial)
 	private int Lmin;   	//Minimum alignment length of a Block
 	private int convergenceSteps; //Steps without a score change before stopping the optimization
-	private static final double C = 20; //Probability function constant (probability of acceptance for bad moves)
+	private double C; //Probability function constant (probability of acceptance for bad moves)
 	
 	//Score function parameters - they are defined by the user
 	private double Gopen; //Penalty for opening gap
@@ -90,6 +90,7 @@ public MultipleAlignmentOptimizerMC(MultipleAlignment seedAln, MultipleMcParamet
 		if (params.getMinAlignedStructures() == 0) Rmin = Math.max(size/3, 2);
 		else Rmin = Math.min(Math.max(params.getMinAlignedStructures(),2),size);  //has to be in the range [2-size]
 		Lmin = params.getMinBlockLen();
+		C = 10*size;
 	}
 
 	@Override
@@ -192,15 +193,15 @@ private void optimizeMC(int maxIter) throws StructureException {
 				//Randomly select one of the steps to modify the alignment. 
 				//Because two moves are biased, the probabilities are not the same
 				double move = rnd.nextDouble();
-				if (move < 0.4){
+				if (move < 0.5){
 						moved = shiftRow();
 						if (debug) System.out.println("did shift");
 				}
-				else if (move < 0.7){
+				else if (move < 0.8){
 						moved = expandBlock();
 						if (debug) System.out.println("did expand");
 				}
-				else if (move < 0.85){
+				else if (move < 0.9){
 						moved = shrinkBlock();
 						if (debug) System.out.println("did shrink");
 				}
@@ -640,12 +641,12 @@ private boolean shrinkBlock(){
 	/**
 	 *  Calculates the probability of accepting a bad move given the iteration step and the score change.
 	 *  
-	 *  Function: p=(C-AS)/m   *slightly different from the CEMC algorithm.
+	 *  Function: p=(C-AS)/(m*C)   *slightly different from the CEMC algorithm.
 	 *  Added a normalization factor so that the probability approaches 0 when the maxIter is reached.
 	 */
 	private double probabilityFunction(double AS, int m, int maxIter) {
 		
-		double prob = (C+AS)/(m);
+		double prob = (C+AS)/(m*C);
 		double norm = (1-(m*1.0)/maxIter);  //Normalization factor
 		return Math.min(Math.max(prob*norm,0.0),1.0);
 	}
