added -e option to start estimation with custom initial parameters

panicmage.c

@@ -50,7 +50,7 @@ int g_includecoreflag = 0;
 //parse the commandline for options
 
 
-  int get_args(int argc, char** argv, int* notest_flag, int* samplingbias_flag,  float* theta_value, float* rho_value, float* core_value, int* estimate_flag, int* details_flag, int* skipall_flag, int* runs_input, int* printtree_flag, int* pansize_flag, float* millgenerationstoMRCA, int* scale_flag)
+  int get_args(int argc, char** argv, int* notest_flag, int* samplingbias_flag,  float* theta_value, float* rho_value, float* core_value, int* estimate_flag, int* details_flag, int* skipall_flag, int* runs_input, int* printtree_flag, int* pansize_flag, float* millgenerationstoMRCA, int* scale_flag, float* customstarttheta, float* customstartrho, int* customstart_flag)
 {
     int i;
 
@@ -110,6 +110,19 @@ int g_includecoreflag = 0;
 				*millgenerationstoMRCA = atof(argv[++i]);
 				*pansize_flag = 1;
 				break;
+        case 'e':   if(i+2 >= argc) return -1;
+                *customstarttheta = atof(argv[++i]);
+                if (argv[i][0] == '-'){ 
+                    fprintf(stderr, "Missing parameter\n");
+                    return -1;
+                }
+                *customstartrho = atof(argv[++i]);
+                if (argv[i][0] == '-'){ 
+                    fprintf(stderr, "Missing second parameter\n");
+                    return -1;
+                }
+                *customstart_flag = 1;
+                break;
 
 		default:	fprintf(stderr,
 				"Unknown option %s\n", argv[i]);
@@ -222,15 +235,16 @@ printf("Filename Gene Frequency Data:\t %s\n", argv[2]);
 //   -a skip all but the results for a typical population
 //   -p print tree
 //   -b scaling of the tree is skipped
+//   -e custom start values for the estimation are used (need two )
 
 
     // Set defaults for all parameters
-    int notest_flag = 0, samplingbias_flag = 0, estimate_flag = 1, skipall_flag = 0, details_flag = 1, printtree_flag = 0, pansize_flag = 0, oldalgo_flag = 0, scale_flag = 1;
-    float theta_input = 0.0, rho_input = 0.0, core_input = 0.0, millgenstoMRCA_input = 0.0;
+    int notest_flag = 0, samplingbias_flag = 0, estimate_flag = 1, skipall_flag = 0, details_flag = 1, printtree_flag = 0, pansize_flag = 0, oldalgo_flag = 0, scale_flag = 1, customstart_flag = 0;
+    float theta_input = 0.0, rho_input = 0.0, core_input = 0.0, millgenstoMRCA_input = 0.0, customstarttheta = 0.0, customstartrho = 0.0;
     int runs_input = 0;
 
 
-    if ( get_args(argc, argv, &notest_flag, &samplingbias_flag, &theta_input, &rho_input, &core_input, &estimate_flag, &details_flag, &skipall_flag, &runs_input, &printtree_flag, &pansize_flag, &millgenstoMRCA_input, &scale_flag ) == -1 ){
+    if ( get_args(argc, argv, &notest_flag, &samplingbias_flag, &theta_input, &rho_input, &core_input, &estimate_flag, &details_flag, &skipall_flag, &runs_input, &printtree_flag, &pansize_flag, &millgenstoMRCA_input, &scale_flag, &customstarttheta, &customstartrho, &customstart_flag ) == -1 ){
       printf("something was wrong with your options, maybe you wrote '-t200' instead of '-t 200' ?\n");
       return -1;
     }
@@ -538,6 +552,9 @@ if (skipall_flag == 0) {
   if (estimate_flag == 1){
     printf("Estimating theta and rho...this may take some time\n");
     if (oldalgo_flag == 1){
+        if(customstart_flag == 1){
+            printf("WARNING: custom start values are not available in the non-numeric estimation procedure!\n");
+        }
       estimate(&theta_hat,&rho_hat,para);
     }
     else{
@@ -587,7 +604,7 @@ if (skipall_flag == 0) {
   // // // // //    
 
 //  cout << "We successfully called treegfs_symbolic_fast\n";        
-      estimate_numeric(&theta_hat,&rho_hat,para);
+      estimate_numeric(&theta_hat,&rho_hat,para,customstart_flag,customstarttheta,customstartrho);
     }
   }
   //////////////////* END estimate theta and rho*//////////////////////////////////////////

source/treenumeric.h

@@ -206,7 +206,7 @@ double my_f_numeric (const gsl_vector *v, void *params)
 // para->tree = tree;
 // }
 
-void estimate_numeric(float *theta_hat, float *rho_hat, Params *para){   
+void estimate_numeric(float *theta_hat, float *rho_hat, Params *para, int customstart_flag,float customstarttheta, float customstartrho){   
 
        const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex2;
        gsl_multimin_fminimizer *s = NULL;
@@ -236,6 +236,13 @@ void estimate_numeric(float *theta_hat, float *rho_hat, Params *para){
        if (firsttheta < 50.) firsttheta = 50;
 
 
+       if (customstart_flag == 1){
+        firsttheta = customstarttheta;
+        firstrho = customstartrho;
+        printf("Using custom start values for numeric parameter estimation:\n");
+        printf("theta start value = %.0f\t rho start value = %.2f\n",firsttheta, firstrho);
+        }
+
        /* Starting point */
        x = gsl_vector_alloc (2);
 
@@ -245,7 +252,7 @@ void estimate_numeric(float *theta_hat, float *rho_hat, Params *para){
        gsl_vector_set (x, 0, firstrho);  // estimated starting point rho
        gsl_vector_set (x, 1, firsttheta); // estimated starting point theta
 
-       //        printf("firsttheta = %.0f\tfirstrho = %.2f\n",firsttheta, firstrho);
+//        printf("firsttheta = %.0f\tfirstrho = %.2f\n",firsttheta, firstrho);
 
 
        /* Set initial step sizes to 1 */