Refactored probability check; expokit test removed from active code

barparser.c

@@ -310,7 +310,7 @@ my_getline(FILE *fp)
 }
 
 // visualisation
-void VisulizeRates(char *filename, double *R, BarData *Data, int dim)
+void VisualizeRates(char *filename, double *R, BarData *Data, int dim)
 {
   FILE *dotf;
 

barparser.h

@@ -50,7 +50,7 @@ int  ParseSaddleFile(TypeDegSaddle **my_saddle);
 int  ParseRatesFile(FILE *rates_FP, double **Raten, int nstates, int max);
 
 // visualisation
-void VisulizeRates(char *filename, double *R, BarData *Data, int dim);
+void VisualizeRates(char *filename, double *R, BarData *Data, int dim);
 
 #endif
 

calc.c

@@ -395,7 +395,7 @@ MxIterate (double *p0, double *p8, double *S)
     tmpVec2 = exp(time * EV) *tmpVec
     p(t)    = S * tmpVec2
   */
-  int i,  count = 0, pdiff_counter = 0;
+  int i,  count = 0;
   double time, check = 0.;
   double *CL = NULL, *CR, *exptL, *tmpVec, *tmpVec2, *pt, *St, *pdiff;
   double *ptFULL = NULL;  /* prob dist 4 of the effective lmins of the tree at time t */
@@ -479,41 +479,11 @@ MxIterate (double *p0, double *p8, double *S)
     else                 check = PrintProb(pt, dim, time);
 
     //PrintProbNR(p8FULL, lmins, -1);
-
-    if ( ((check-1) < -0.05) || ((check-1) > 0.05) ) {
-      fprintf(stderr, "overall probability at time %e is %e != 1. ! exiting\n", time,check );
-      if (opt.num_err == 'H') exit(EXIT_FAILURE);
-    }
-    check = 0.;
-    /* now check if we have converged yet */
-    if(opt.method=='F') {
-      for(i = 1; i <= lmins; i++) {
-        pdiff[i] = p8FULL[i] - ptFULL[i];
-        if (fabs(pdiff[i]) >= 0.000001) {
-          pdiff_counter++;
-          break;
-        }
-      }
-    }
-    else {
-      for(i = 0; i < dim; i++) {
-        pdiff[i] = p8[i] - pt[i];
-
-        if (fabs(pdiff[i]) >= 0.000001) {
-          pdiff_counter++;
-          break;
-        }
-      }}
-              //PrintProb(p8, dim, -1);
-        //PrintProb(pt, dim, -2);
-    //}
-    if (pdiff_counter < 1) /* all mins' pdiff lies within threshold */
-      break;
-    pdiff_counter = 0;
-    memset(pdiff, 0, (lmins+1)*sizeof(double));
-    /* end check of convergence */
-
+    int reached;
+    if (opt.method=='F') reached = ConvergenceReached(p8FULL, ptFULL, lmins, 1);
+    else                 reached = ConvergenceReached(p8, pt, dim, 0);
     fflush(stdout);
+    if (reached) break;
   }
   if (time < opt.t8) {
     if (opt.method=='F') PrintProbFull(pt, dim, opt.t8, lmins);

calcpp.cpp

@@ -19,7 +19,8 @@ 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);
+extern "C" int ConvergenceReached(double *p8, double *pt, int dim, int full);
+extern "C" void TestExpokit(double *R, int dim, double *p0, double t_start, double t_end, double t_inc);
 
 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;
@@ -131,6 +132,13 @@ double PrintProbFull(double *line, int dim, double time, int lmins)
     check += fabs(ptFULL[i]);
   }
   printf("\n");
+
+  // check for overall propability
+  if ( ((check-1) < -0.05) || ((check-1) > 0.05) ) {
+    fprintf(stderr, "overall probability at time %e is %e != 1. ! exiting\n", time,check );
+    if (opt.num_err == 'H') exit(EXIT_FAILURE);
+  }
+
   return check;
 }
 
@@ -151,6 +159,13 @@ double PrintProbNR(double *line, int dim, double time)
   }
 
   printf("\n");
+
+  // check for overall propability
+  if ( ((check-1) < -0.05) || ((check-1) > 0.05) ) {
+    fprintf(stderr, "overall probability at time %e is %e != 1. ! exiting\n", time,check );
+    if (opt.num_err == 'H') exit(EXIT_FAILURE);
+  }
+
   return check;
 }
 
@@ -188,32 +203,88 @@ double PrintProb(double *line, int dim, double time)
     }
   }
   printf("\n");
+
+  // check for overall propability
+  if ( ((check-1) < -0.05) || ((check-1) > 0.05) ) {
+    fprintf(stderr, "overall probability at time %e is %e != 1. ! exiting\n", time,check );
+    if (opt.num_err == 'H') exit(EXIT_FAILURE);
+  }
+
   return check;
 }
 
-double TestExpokit(double *R, int dim, double *V)
+int ConvergenceReached(double *p8, double *pt, int dim, int full) {
+
+  int pdiff_counter = 0;
+  full = (full?1:0);
+  /* now check if we have converged yet */
+  for(int i = 0; i < dim; i++) {
+    if (fabs(p8[i] - pt[i]) >= 0.000001) {
+      pdiff_counter++;
+      break;
+    }
+  }
+
+  if (pdiff_counter < 1) /* all mins' pdiff lies within threshold */
+    return 1;
+  pdiff_counter = 0;
+  /* end check of convergence */
+  return false;
+}
+
+void TestExpokit(double *R, int dim, double *p0, double t_start, double t_end, double t_inc)
 {
 
-  int n = 5;
+  int n = dim;
   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]; //??
+  double tol = 0.01; // change to something better maybe?
+  double anorm[n*n]; //??
   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);
+  // change the R into ia/ja matrix:
+
+    // now get the non-zero ones
+  vector<double> non_zero;
+  vector<int> ia_vec;
+  vector<int> ja_vec;
+  for (int i=0; i<dim; i++) {
+    for (int j=0; j<dim; j++) {
+      if (fabs(R[i*dim+j]-(i==j?1.0:0.0)) > 0.0) {
+        non_zero.push_back(R[i*dim+j]-(i==j?1.0:0.0));
+        ia_vec.push_back(i);
+        ja_vec.push_back(j);
+      }
+    }
+  }
+
+    // convert them into C code
+  int ia[non_zero.size()];
+  int ja[non_zero.size()];
+  double a[non_zero.size()];
+  int nz = non_zero.size();
+  for (unsigned int i=0; i<non_zero.size(); i++) {
+    ia[i] = ia_vec[i];
+    ja[i] = ja_vec[i];
+    a[i] = non_zero[i];
+  }
+
+/*
+  const int nzc = 9;
+  int ia[nzc] = {1,1,1,2,2,2,3,3,3};
+  int ja[nzc] = {1,2,3,1,2,3,1,2,3};
+  double a[nzc] = {-0.9, 0.5, 0.5, 0.5, -1.0, 0.5, 0.4, 0.5, -1.0};
+  int nz = nzc;
+*/
+  // main loop
+  for (double time=t_start; time<=t_end; time *= t_inc) {
+    wrapsingledgexpv_(&n, &m, &time, p0, w, &tol, anorm, wsp, &lwsp, iwsp, &liwsp, &itrace, &iflag, ia, ja, a, &nz, res);
+    PrintProb(w, dim, time);
+  }
 
-  return 0.0;
 }

calcpp.h

@@ -11,8 +11,10 @@ double PrintProb(double *line, int dim, double time);
 double PrintProbFull(double *line, int dim, double time, int lmins);
 double PrintProbNR(double *line, int dim, double time);
 
+int ConvergenceReached(double *p8, double *pt, int dim, int full);
 
-double TestExpokit(double *R, int dim, double *V);
+
+void TestExpokit(double *R, int dim, double *p0, double t_start, double t_end, double t_inc);
 
 
 #endif

expokit_wrappers.h

@@ -0,0 +1,160 @@
+
+// not needed apparently
+#ifndef EXPOKIT_WRAPPERS_H
+#define EXPOKIT_WRAPPERS_H
+
+/*
+#######################################################
+# expokit_wrappers.h
+# Wrapper functions for expokit functions.
+#######################################################
+# Assembled, not really authored, by
+# Nick Matzke, matzke@berkeley.edu
+# June 2012
+#######################################################
+*/
+
+/*
+Acknowledgements/sources:
+
+1. Copied in part from a file in Lagrange:
+http://code.google.com/p/lagrange/
+https://github.com/blackrim/lagrange
+
+Specifically:
+ * RateMatrix.cpp
+ *
+ *  Created on: Aug 14, 2009
+ *      Author: smitty
+ *
+
+...and the my_*.f wrappers for the EXPOKIT *.f code files.
+
+*/
+
+/*
+2. Also copied in part (to get the .h file) from:
+
+Python package "Pyprop":
+http://code.google.com/p/pyprop/
+http://pyprop.googlecode.com/svn/trunk/core/krylov/expokit/expokitpropagator.cpp
+http://www.koders.com/python/fidCA95B5A4B2FB77455A72B8A361CF684FFE48F4DC.aspx?s=fourier+transform
+
+Specifically:
+pyprop/core/krylov/expokit/f2c/expokit.h
+*/
+
+/*
+3. EXPOKIT package is available at:
+http://www.maths.uq.edu.au/expokit/
+
+Copyright:
+http://www.maths.uq.edu.au/expokit/copyright
+...or...
+expokit_copyright.txt in this install
+*/
+
+/*
+Note: Because these functions are defined with
+extern"C", they can be sourced from R with e.g.
+these commands:
+
+################################
+# Begin R code
+################################
+# Specify the *.so file (a 'shared object' file, in this case
+# rexpokit.so), which is produced by the compiler linking
+# the individual *.o files
+wrapper_so_fn = paste(wd, "rexpokit.so", sep="")
+
+# Unload the previous *.so file (if needed), then load the new one
+dyn.unload(wrapper_so_fn)
+dyn.load(wrapper_so_fn)
+
+# Alledgely, this will list the loaded routines, but
+# it doesn't seem to work
+(dlls = getLoadedDLLs())
+getDLLRegisteredRoutines(dll="rexpokit", addNames = TRUE)
+
+# However, they are available, as you can see:
+is.loaded("wrapalldmexpv_")
+is.loaded("wrapsingledmexpv_")
+is.loaded("wrapdgpadm_")
+###############################
+*/
+
+
+// not needed apparently
+//#include <core/common.h>
+
+// not needed apparently
+//namespace expokit_wrappers {
+
+extern"C" {
+// These functions are defined in my_matexp.f, which calls the following functions/
+
+// DMEXPV contains an additional check ensuring sums to 1, for Markov-chain applications
+
+// This exponentiates a large, sparse matrix (sparse = lots of zeros)
+// Before input, the matrix should be transposed and
+// put into coordinate list (COO) format:
+// http://en.wikipedia.org/wiki/Sparse_matrix#Coordinate_list_.28COO.29
+//
+// i.e.:
+// ia = row number
+// ja = col number
+// a  = corresponding value of that cell (zeros are excluded)
+void wrapalldmexpv_(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);
+
+// This returns just one row (?) of the transition matrix, useful in Lagrange;
+// same inputs as wrapalldmexpv_.
+void wrapsingledmexpv_(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);
+
+// DGEXPV contains an additional check ensuring sums to 1, for Markov-chain applications
+void wrapalldgexpv_(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);
+
+// This returns just one row (?) of the transition matrix, useful in Lagrange;
+// same inputs as wrapalldmexpv_.
+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);
+
+
+// The myDMEXPV etc. functions provide direct access to EXPOKIT functions;
+// This should be faster, especially for sparse matrices
+// Here, you input v (starting probabilities) and it fills in w, which are the
+// output probabilities (in output list item #5)
+void myDMEXPV_(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 );
+
+void myDGEXPV_(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 );
+
+
+// This exponentiates a small, dense (no zeros) matrix with the padm
+// approximation.  Here, the input matrix is just the
+// matrix, tranposed and then put into a list of numbers, e.g.:
+//
+/*
+# R code
+# transpose the matrix in R
+tmat = t(mat)
+# convert into list of numbers (as numeric or double)
+tmat_nums = as.double(tmat)
+*/
+void wrapdgpadm_(int * ideg,int * m,double * t,double * H,int * ldh,
+	double * wsp,int * lwsp,int * ipiv,int * iexph,int *ns,int *iflag );
+
+} // end extern "C"
+
+// not needed apparently
+//};
+#endif