Merge branch 'master' of github.com:marcelTBI/treekin

Conflicts: calcpp.cpp
ViennaRNA · Dec 5, 2013 · 91be572 · 91be572
2 parents b61d60d + 40f0adc
commit 91be572
Show file tree

Hide file tree

Showing 6 changed files with 74 additions and 37 deletions.
diff --git a/calc.c b/calc.c
@@ -311,6 +311,7 @@ MxDiagonalize ( double *U, double **_S, double *P8)
     if (opt.want_verbose) MxPrint (U, "force symmetrized U", 'm');
   }
 
+  // get eigenv*
   if(opt.absrb) {
     MxEVLapackNonSym(U);
   } else {
@@ -436,8 +437,12 @@ MxIterate (double *p0, double *p8, double *S)
     for (i = 0; i < dim; i++)   p8FULL[E[i].ag] += p8[i];
     for (i = 0; i < lmins; i++) check += fabs(p8FULL[i]);
     if ( ((check-1) < -0.1) || ((check-1) > 0.1) ) {
-      fprintf(stderr, "overall equilibrium probability is %e != 1. ! exiting\n", check);
-      exit(EXIT_FAILURE);
+      fprintf(stderr, "overall equilibrium probability is %e != 1. !\n", check);
+      if (opt.num_err == 'H') exit(EXIT_FAILURE);
+      else if (opt.num_err == 'R') {
+        for (i=0; i<dim; i++) p8[i] /= check;
+        check = 1.0;
+      }
     }
   }
   check = 0.;
@@ -477,7 +482,7 @@ MxIterate (double *p0, double *p8, double *S)
 
     if ( ((check-1) < -0.05) || ((check-1) > 0.05) ) {
       fprintf(stderr, "overall probability at time %e is %e != 1. ! exiting\n", time,check );
-      exit(EXIT_FAILURE);
+      if (opt.num_err == 'H') exit(EXIT_FAILURE);
     }
     check = 0.;
     /* now check if we have converged yet */
@@ -1262,7 +1267,7 @@ MxExponent(double *p0, double *p8, double *U)
 
     if ( ((check-1) < -0.01) || ((check-1) > 0.01) ) {
       fprintf(stderr, "overall probability at time %e is %e != 1. ! exiting\n", time,check );
-      exit(EXIT_FAILURE);
+      if (opt.num_err == 'H') exit(EXIT_FAILURE);
     }
     memset(pt,   0, dim*sizeof(double));
     memset(pdiff, 0, dim*sizeof(double));
@@ -1583,9 +1588,7 @@ MxEVLapackNonSym(double *origU)
   int one, ilo, ihi, lwork, *iwork, nfo;
   //int dimx2;
   /* for sorting */
-  int i, j;
-  //float p;
-  double tmp;
+  int i;
 
   dim = dim+500;
   one = 1;
@@ -1607,12 +1610,13 @@ MxEVLapackNonSym(double *origU)
   }
 
   /* instead of more fiddling, we transpose the input */
-  for (i=0; i<dim; i++)
+  trnm(origU, dim);
+  /*for (i=0; i<dim; i++)
     for (j=i+1; j<dim; j++) {
       tmp = origU[dim*i+j];
       origU[dim*i+j]=origU[dim*j+i];
       origU[dim*j+i]=tmp;
-    }
+    }*/
 
   dgeevx_("B", "N", "V", "V", &dim, origU, &dim, evals_re, evals_im, NULL, &one,\
         evecs ,&dim, &ilo, &ihi, scale, &abnrm, rconde, rcondv, work,\
@@ -1707,6 +1711,17 @@ int MxShorten(double **shorten, int nstates, int my_dim, int max) {
   if (my_dim>nstates) {
     if (!opt.quiet) fprintf(stderr, "decreasing %d to %d\n", my_dim, nstates);
 
+    // first we need to fix the diagonal entries tmp_rates[i][i] = sum_j tmp_rates[i][j]
+    double *tmp_rates = *shorten;
+    int i, j;
+    for (i = 0; i < my_dim; i++) tmp_rates[my_dim*i+i] = 0.0;
+    for (i = 0; i < my_dim; i++) {
+      double tmp = 0.00;
+      // calculate column sum
+      for(j = 0; j < my_dim; j++)  tmp += tmp_rates[my_dim*j+i];
+      tmp_rates[my_dim*i+i] = -tmp;
+    }
+
     // shorten by some value max
     if (max!=1) {
       while (my_dim-max > nstates) {
@@ -1720,6 +1735,7 @@ int MxShorten(double **shorten, int nstates, int my_dim, int max) {
       while (my_dim!=nstates) {
         MxOneShorten(shorten, my_dim);
         my_dim--;
+        if (!opt.quiet && my_dim%100==0 && my_dim>0) fprintf(stderr, "%d done...\n", my_dim);
       }
     }
   }
@@ -1737,23 +1753,13 @@ void MxRShorten(double **shorten, int fulldim, int gdim)
 
   int bdim = fulldim - gdim;
   int i,j;
+  double *tmp_rates = *shorten;
 
   double *gg = (double *)calloc(gdim*gdim,sizeof(double));
   double *bg = (double *)calloc(bdim*gdim,sizeof(double));
   double *bb = (double *)calloc(bdim*bdim,sizeof(double));
   double *gb = (double *)calloc(gdim*bdim,sizeof(double));
 
-  // first we need to fix the diagonal entries tmp_rates[i][i] = sum_j tmp_rates[i][j]
-  double *tmp_rates = *shorten;
-  for (i = 0; i < fulldim; i++) tmp_rates[fulldim*i+i] = 0.0;
-  for (i = 0; i < fulldim; i++) {
-    double tmp = 0.00;
-    // calculate column sum
-    for(j = 0; j < fulldim; j++)  tmp += tmp_rates[fulldim*j+i];
-    tmp_rates[fulldim*i+i] = -tmp;
-  }
-
-
   // fill the matrices: (row = i; column = j)
   for (i=0; i<gdim; i++) {
     for (j=0; j<gdim; j++) {
@@ -1832,12 +1838,12 @@ void MxOneShorten(double **shorten, int fulldim)
   for (i=0; i<gdim; i++) {
     for (j=0; j<gdim; j++) {
       // just x - a*c^-1*b
-      result[gdim*i+j] = tmp_rates[fulldim*i+j] - c*tmp_rates[fulldim*gdim + i]*tmp_rates[fulldim*j + gdim];
+      result[gdim*i+j] = tmp_rates[fulldim*i+j] - c*tmp_rates[fulldim*gdim + j]*tmp_rates[fulldim*i + gdim];
     }
   }
 
-  MxFPrintD(tmp_rates, "Q", fulldim, fulldim, stderr);
-  MxFPrintD(result, "Q-1", gdim, gdim, stderr);
+  //MxFPrintD(tmp_rates, "Q", fulldim, fulldim, stderr);
+  //MxFPrintD(result, "Q-1", gdim, gdim, stderr);
 
   free(*shorten);
   *shorten = result;

diff --git a/calcpp.cpp b/calcpp.cpp
@@ -8,6 +8,7 @@
 
 extern "C" {
   #include "calc.h"
+  #include "expokit_wrappers.h"
 }
 
 #include "expokit_wrappers.h"
@@ -18,6 +19,7 @@ extern "C" void MxEgro(double **U, double **p0, int dim);
 extern "C" double PrintProb(double *line, int dim, double time);
 extern "C" double PrintProbNR(double *line, int dim, double time);
 extern "C" double PrintProbFull(double *line, int dim, double time, int lmins);
+extern "C" double TestExpokit(double *R, int dim, double *V);
 
 vector<int> reorganize; // reorganize array (so if LM 0 1 3 were reachable and 2 not, reorganize will contain r[0]=0 r[1]=1 r[2]=3), so r[x] = old position of x
 int last_dim;
@@ -121,10 +123,11 @@ double PrintProbFull(double *line, int dim, double time, int lmins)
   for (int i=0; i<lmins; i++) {
     if(ptFULL[i] < -0.01) {
       fprintf(stderr, "prob of lmin %i at time %e has become negative: %e \n", i+1, time, ptFULL[i]);
-      exit(EXIT_FAILURE);
+      if (opt.num_err == 'H') exit(EXIT_FAILURE);
+      else if (opt.num_err == 'R') ptFULL[i] = 0.0;
     }
     /* map individual structure -> gradient basins */
-    else   printf("%e ", fabs(ptFULL[i]));
+    printf("%e ", fabs(ptFULL[i]));
     check += fabs(ptFULL[i]);
   }
   printf("\n");
@@ -139,10 +142,11 @@ double PrintProbNR(double *line, int dim, double time)
   for (int i=0; i<dim; i++) {
     if(line[i] < -0.01) {
       fprintf(stderr, "prob of lmin %i at time %e has become negative: %e \n", i+1, time, line[i]);
-      exit(EXIT_FAILURE);
+      if (opt.num_err == 'H') exit(EXIT_FAILURE);
+      else if (opt.num_err == 'R') line[i] = 0.0;
     }
     /* map individual structure -> gradient basins */
-    else printf("%e ", fabs(line[i]));
+    printf("%e ", fabs(line[i]));
     check += fabs(line[i]);
   }
 
@@ -159,10 +163,11 @@ double PrintProb(double *line, int dim, double time)
     for (int i=0; i<dim; i++) {
       if(line[i] < -0.01) {
         fprintf(stderr, "prob of lmin %i at time %e has become negative: %e \n", i+1, time, line[i]);
-        exit(EXIT_FAILURE);
+        if (opt.num_err == 'H') exit(EXIT_FAILURE);
+        else if (opt.num_err == 'R') line[i] = 0.0;
       }
       /* map individual structure -> gradient basins */
-      else printf("%e ", fabs(line[i]));
+      printf("%e ", fabs(line[i]));
       check += fabs(line[i]);
     }
   } else {
@@ -173,7 +178,8 @@ double PrintProb(double *line, int dim, double time)
       } else {
         if(line[j] < -0.01) {
           fprintf(stderr, "prob of lmin %i at time %e has become negative: %e \n", i+1, time, line[i]);
-          exit(EXIT_FAILURE);
+          if (opt.num_err == 'H') exit(EXIT_FAILURE);
+          else if (opt.num_err == 'R') line[j] = 0.0;
         }
         printf("%e ", fabs(line[j]));
         check += fabs(line[j]);
@@ -187,12 +193,27 @@ double PrintProb(double *line, int dim, double time)
 
 double TestExpokit(double *R, int dim, double *V)
 {
-  void wrapsingledgexpv_(int * n,int* m,double * t,double* v,double * w,double* tol,
-	double* anorm,double* wsp,int * lwsp,int* iwsp,int *liwsp, int * itrace,int *iflag,
-			int *ia, int *ja, double *a, int *nz, double * res);
 
-
- return 0.0;
+  int n = 5;
+  int m = n-1;
+  double t = 2.1;
+  double v[n];
+  v[0] = 1.0;
+  double w[n];
+  double tol = 0.01;
+  double anorm[25]; //??
+  int lwsp = n*(m+2)+5*(m+2)*(m+2)+7;
+  int liwsp = max(m+2, 7);
+  int iwsp[liwsp];
+  int itrace = 1, iflag = 1;
+  double wsp[lwsp];
+  int ia[10] = {1,2,3,4,5,1,2,3,4,5};
+  int ja[10] = {1,2,3,4,5,2,3,4,5,2};
+  double a[10] = {0.5, 0.5, 0.5, 0.5, 0.5, 0.1, 0.2, 0.3, 0.4, 0.1};
+  int nz = 10;
+  double res[n*n];
+
+  wrapsingledmexpv_(&n, &m, &t, v, w, &tol, anorm, wsp, &lwsp, iwsp, &liwsp, &itrace, &iflag, ia, ja, a, &nz, res);
+
+  return 0.0;
 }
-
-
diff --git a/globals.c b/globals.c
@@ -102,6 +102,11 @@ set_parameters(void)
   else if (strncmp(args_info.method_arg, "A", 1)==0)
     opt.method = 'A';
 
+  if(strncmp(args_info.num_err_arg, "I", 1)==0)
+    opt.num_err = 'I';
+  else if (strncmp(args_info.num_err_arg, "R", 1)==0)
+    opt.num_err = 'R';
+
   if (args_info.fptfile_given) {
     opt.fpt_file = strdup(args_info.fptfile_arg);
   } else opt.fpt_file = NULL;
@@ -252,6 +257,7 @@ ini_globs(void)
   opt.basename        =          NULL;
   opt.vis_file        =          NULL;
   opt.just_sh         =          0;
+  opt.num_err         =          'H';
 }
 
 /*==============================*/

diff --git a/globals.h b/globals.h
@@ -47,6 +47,7 @@ typedef struct {         /* command-line options */
   int quiet;             // be quiet?
   int just_sh;           // just shorten
   int max_decrease;      // how many states to decrease at once
+  char num_err;          // method for numerical error handling
 } treekin_options;
 
 treekin_options opt;

diff --git a/main.c b/main.c
@@ -30,6 +30,8 @@ main (int argc, char **argv)
 {
   clock_t clck1 = clock();
 
+  TestExpokit(NULL, 5, NULL);
+
   BarData *Data=NULL;
   /*  U - matrix (Q+I)^T, where Q is infetisimal generator (^T - transposed)
       S - eigenvectors of U

diff --git a/treekin.ggo b/treekin.ggo
@@ -11,6 +11,7 @@ args "--file-name=treekin_cmdline --include-getopt --default-optional --unamed-o
 # Options
 option "absorb"             a "Make a state absorbing" typestr="state" int no
 option "method"             m "Select method to build transition matrix:\nA ==> Arrhenius-like kinetics\nF ==> Full process kinetics (whole subopt)\nI ==> use rates from barriers" values="A","F","I" default="I" no
+option "num-err"            - "Method how to deal wih numerical errors in probability:\nI ==> Ignore\nH ==> Halt the program\nR ==> Rescale the probability" values="I","H","R" default="H" no
 option "t0"                 - "Start time" typestr="time" default="0.1" double no
 option "t8"                 - "Stop time" typestr="time" default="1E9" double no
 option "Temp"               T "Temperature in Celsius" default="37.0" double no