More tests with ml2 to find optimal results

1 files changed, 93 insertions, 110 deletions

diff --git a/experiments/ml.c b/experiments/ml.c
index b350aaf..64d4e3a 100644
--- a/experiments/ml.c
+++ b/experiments/ml.c
@@ -983,8 +983,6 @@ static double __Pyx__PyObject_AsDouble(PyObject* obj);
  PyFloat_AS_DOUBLE(obj) : __Pyx__PyObject_AsDouble(obj))
 #endif
 
-static CYTHON_INLINE void __Pyx_RaiseUnboundLocalError(const char *varname);
-
 static CYTHON_INLINE void __Pyx_ErrRestore(PyObject *type, PyObject *value, PyObject *tb);
 static CYTHON_INLINE void __Pyx_ErrFetch(PyObject **type, PyObject **value, PyObject **tb);
 
@@ -1475,32 +1473,32 @@ static __pyx_t_2ml_DTYPE_t __pyx_f_2ml_weight_success(int __pyx_v_dist, int __py
   /* "ml.pyx":16
  *     """weight for successful infection, exponential time model"""
  *     cdef DTYPE_t structural, temporal, result
- *     structural = delta ** dist             # <<<<<<<<<<<<<<
+ *     structural = dist * log(delta)             # <<<<<<<<<<<<<<
  *     # structural = plogis(w1,delta) * plogis(w2,delta) * plogis(w3,delta)
- *     temporal = exp(-alpha*dt) * (exp(alpha)-1.)
+ *     temporal = log(exp(alpha)-1.) - alpha*dt
  */
-  __pyx_v_structural = pow(__pyx_v_delta, ((__pyx_t_2ml_DTYPE_t)__pyx_v_dist));
+  __pyx_v_structural = (__pyx_v_dist * log(__pyx_v_delta));
 
   /* "ml.pyx":18
- *     structural = delta ** dist
+ *     structural = dist * log(delta)
  *     # structural = plogis(w1,delta) * plogis(w2,delta) * plogis(w3,delta)
- *     temporal = exp(-alpha*dt) * (exp(alpha)-1.)             # <<<<<<<<<<<<<<
+ *     temporal = log(exp(alpha)-1.) - alpha*dt             # <<<<<<<<<<<<<<
  *     # temporal = 1 - exp(-alpha*dt)
  *     # if exp(-alpha*dt)==0.: print 'UNDERFLOW ERROR'
  */
-  __pyx_v_temporal = (exp(((-__pyx_v_alpha) * __pyx_v_dt)) * (exp(__pyx_v_alpha) - 1.));
+  __pyx_v_temporal = (log((exp(__pyx_v_alpha) - 1.)) - (__pyx_v_alpha * __pyx_v_dt));
 
   /* "ml.pyx":22
  *     # if exp(-alpha*dt)==0.: print 'UNDERFLOW ERROR'
  *     # temporal = 1. / (1. + (dt - 1.)/alpha)**0.01 - 1. / (1. + dt/alpha)**0.01
- *     result = log(structural * temporal)             # <<<<<<<<<<<<<<
+ *     result = structural + temporal             # <<<<<<<<<<<<<<
  *     # print 'st', structural, temporal
  *     return result
  */
-  __pyx_v_result = log((__pyx_v_structural * __pyx_v_temporal));
+  __pyx_v_result = (__pyx_v_structural + __pyx_v_temporal);
 
   /* "ml.pyx":24
- *     result = log(structural * temporal)
+ *     result = structural + temporal
  *     # print 'st', structural, temporal
  *     return result             # <<<<<<<<<<<<<<
  * 
@@ -1851,14 +1849,14 @@ static PyObject *__pyx_pf_2ml_ml(CYTHON_UNUSED PyObject *__pyx_self, PyObject *_
   int __pyx_v_dt;
   __pyx_t_2ml_DTYPE_t __pyx_v_beta;
   __pyx_t_2ml_DTYPE_t __pyx_v_ll;
+  __pyx_t_2ml_DTYPE_t __pyx_v_beta_add;
+  __pyx_t_2ml_DTYPE_t __pyx_v_max_beta;
+  __pyx_t_2ml_DTYPE_t __pyx_v_max_beta_add;
   PyObject *__pyx_v_parents = 0;
   PyArrayObject *__pyx_v_probs = 0;
   PyArrayObject *__pyx_v_probs_fail = 0;
   PyArrayObject *__pyx_v_probs_nv = 0;
-  PyObject *__pyx_v_max_beta_add = NULL;
   double __pyx_v_thresh;
-  PyObject *__pyx_v_beta_add = NULL;
-  PyObject *__pyx_v_max_beta = NULL;
   int __pyx_v_max_roots;
   PyObject *__pyx_v_w1 = NULL;
   PyObject *__pyx_v_w2 = NULL;
@@ -1924,7 +1922,7 @@ static PyObject *__pyx_pf_2ml_ml(CYTHON_UNUSED PyObject *__pyx_self, PyObject *_
   __pyx_pybuffernd_probs_nv.rcbuffer = &__pyx_pybuffer_probs_nv;
 
   /* "ml.pyx":44
- *         DTYPE_t beta, ll
+ *         DTYPE_t beta, ll, beta_add, max_beta, max_beta_add
  *         list parents, failures, successes
  *     n_roots, n_victims = len(root_victims), len(victims)             # <<<<<<<<<<<<<<
  *     n_nodes = 148152
@@ -2703,16 +2701,13 @@ static PyObject *__pyx_pf_2ml_ml(CYTHON_UNUSED PyObject *__pyx_self, PyObject *_
  *     for beta in np.arange(0.001, .2, .002):
  */
   __pyx_t_27 = __Pyx_PyObject_AsDouble(__pyx_kp_s_inf); if (unlikely(__pyx_t_27 == ((double)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 84; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-  __pyx_t_5 = PyFloat_FromDouble(__pyx_t_27); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 84; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-  __Pyx_GOTREF(__pyx_t_5);
-  __pyx_v_max_beta_add = __pyx_t_5;
-  __pyx_t_5 = 0;
+  __pyx_v_max_beta_add = __pyx_t_27;
 
   /* "ml.pyx":86
  *     max_beta_add = float('-inf')
  *     # iterate over all victim nodes to find the optimal threshold
  *     for beta in np.arange(0.001, .2, .002):             # <<<<<<<<<<<<<<
- *         thresh = log(beta/(3012*(1.-beta)))
+ *         thresh = log(beta/(3012.*(1.-beta)))
  *         # print 'beta:', beta, 'thresh:', thresh, 'infected:', len(probs[probs>=thresh])
  */
   __pyx_t_5 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 86; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
@@ -2768,18 +2763,18 @@ static PyObject *__pyx_pf_2ml_ml(CYTHON_UNUSED PyObject *__pyx_self, PyObject *_
     /* "ml.pyx":87
  *     # iterate over all victim nodes to find the optimal threshold
  *     for beta in np.arange(0.001, .2, .002):
- *         thresh = log(beta/(3012*(1.-beta)))             # <<<<<<<<<<<<<<
+ *         thresh = log(beta/(3012.*(1.-beta)))             # <<<<<<<<<<<<<<
  *         # print 'beta:', beta, 'thresh:', thresh, 'infected:', len(probs[probs>=thresh])
  *         roots = n_roots + len(probs[probs<thresh])
  */
-    __pyx_v_thresh = log((__pyx_v_beta / (3012.0 * (1. - __pyx_v_beta))));
+    __pyx_v_thresh = log((__pyx_v_beta / (3012. * (1. - __pyx_v_beta))));
 
     /* "ml.pyx":89
- *         thresh = log(beta/(3012*(1.-beta)))
+ *         thresh = log(beta/(3012.*(1.-beta)))
  *         # print 'beta:', beta, 'thresh:', thresh, 'infected:', len(probs[probs>=thresh])
  *         roots = n_roots + len(probs[probs<thresh])             # <<<<<<<<<<<<<<
  * 
- *         # add probability for realized edges and subtract probability these edges fail
+ *         beta_add = 0.
  */
     __pyx_t_5 = PyFloat_FromDouble(__pyx_v_thresh); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 89; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
     __Pyx_GOTREF(__pyx_t_5);
@@ -2792,84 +2787,89 @@ static PyObject *__pyx_pf_2ml_ml(CYTHON_UNUSED PyObject *__pyx_self, PyObject *_
     __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
     __pyx_v_roots = (__pyx_v_n_roots + __pyx_t_1);
 
-    /* "ml.pyx":92
+    /* "ml.pyx":91
+ *         roots = n_roots + len(probs[probs<thresh])
  * 
+ *         beta_add = 0.             # <<<<<<<<<<<<<<
  *         # add probability for realized edges and subtract probability these edges fail
- *         beta_add = (probs[probs>=thresh]).sum()             # <<<<<<<<<<<<<<
+ *         beta_add += (probs[probs>=thresh]).sum()
+ */
+    __pyx_v_beta_add = 0.;
+
+    /* "ml.pyx":93
+ *         beta_add = 0.
+ *         # add probability for realized edges and subtract probability these edges fail
+ *         beta_add += (probs[probs>=thresh]).sum()             # <<<<<<<<<<<<<<
  *         # add probability for the seeds and non-seeds
  *         beta_add += roots * log(beta) + (n_nodes-roots) * log(1. - beta)
  */
-    __pyx_t_3 = PyFloat_FromDouble(__pyx_v_thresh); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 92; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-    __Pyx_GOTREF(__pyx_t_3);
-    __pyx_t_6 = PyObject_RichCompare(((PyObject *)__pyx_v_probs), __pyx_t_3, Py_GE); __Pyx_XGOTREF(__pyx_t_6); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 92; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
-    __pyx_t_3 = PyObject_GetItem(((PyObject *)__pyx_v_probs), __pyx_t_6); if (unlikely(__pyx_t_3 == NULL)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 92; __pyx_clineno = __LINE__; goto __pyx_L1_error;};
-    __Pyx_GOTREF(__pyx_t_3);
+    __pyx_t_5 = PyFloat_FromDouble(__pyx_v_beta_add); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_5);
+    __pyx_t_6 = PyFloat_FromDouble(__pyx_v_thresh); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_6);
+    __pyx_t_18 = PyObject_RichCompare(((PyObject *)__pyx_v_probs), __pyx_t_6, Py_GE); __Pyx_XGOTREF(__pyx_t_18); if (unlikely(!__pyx_t_18)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
     __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
-    __pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_sum); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 92; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __pyx_t_6 = PyObject_GetItem(((PyObject *)__pyx_v_probs), __pyx_t_18); if (unlikely(__pyx_t_6 == NULL)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __pyx_clineno = __LINE__; goto __pyx_L1_error;};
     __Pyx_GOTREF(__pyx_t_6);
-    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
-    __pyx_t_3 = NULL;
-    if (CYTHON_COMPILING_IN_CPYTHON && likely(PyMethod_Check(__pyx_t_6))) {
-      __pyx_t_3 = PyMethod_GET_SELF(__pyx_t_6);
-      if (likely(__pyx_t_3)) {
-        PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_6);
-        __Pyx_INCREF(__pyx_t_3);
+    __Pyx_DECREF(__pyx_t_18); __pyx_t_18 = 0;
+    __pyx_t_18 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_n_s_sum); if (unlikely(!__pyx_t_18)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_18);
+    __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
+    __pyx_t_6 = NULL;
+    if (CYTHON_COMPILING_IN_CPYTHON && likely(PyMethod_Check(__pyx_t_18))) {
+      __pyx_t_6 = PyMethod_GET_SELF(__pyx_t_18);
+      if (likely(__pyx_t_6)) {
+        PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_18);
+        __Pyx_INCREF(__pyx_t_6);
         __Pyx_INCREF(function);
-        __Pyx_DECREF_SET(__pyx_t_6, function);
+        __Pyx_DECREF_SET(__pyx_t_18, function);
       }
     }
-    if (__pyx_t_3) {
-      __pyx_t_5 = __Pyx_PyObject_CallOneArg(__pyx_t_6, __pyx_t_3); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 92; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-      __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+    if (__pyx_t_6) {
+      __pyx_t_3 = __Pyx_PyObject_CallOneArg(__pyx_t_18, __pyx_t_6); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
     } else {
-      __pyx_t_5 = __Pyx_PyObject_CallNoArg(__pyx_t_6); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 92; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __pyx_t_3 = __Pyx_PyObject_CallNoArg(__pyx_t_18); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
     }
-    __Pyx_GOTREF(__pyx_t_5);
-    __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
-    __Pyx_XDECREF_SET(__pyx_v_beta_add, __pyx_t_5);
-    __pyx_t_5 = 0;
+    __Pyx_GOTREF(__pyx_t_3);
+    __Pyx_DECREF(__pyx_t_18); __pyx_t_18 = 0;
+    __pyx_t_18 = PyNumber_InPlaceAdd(__pyx_t_5, __pyx_t_3); if (unlikely(!__pyx_t_18)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_18);
+    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
+    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+    __pyx_t_24 = __pyx_PyFloat_AsDouble(__pyx_t_18); if (unlikely((__pyx_t_24 == (npy_double)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_DECREF(__pyx_t_18); __pyx_t_18 = 0;
+    __pyx_v_beta_add = __pyx_t_24;
 
-    /* "ml.pyx":94
- *         beta_add = (probs[probs>=thresh]).sum()
+    /* "ml.pyx":95
+ *         beta_add += (probs[probs>=thresh]).sum()
  *         # add probability for the seeds and non-seeds
  *         beta_add += roots * log(beta) + (n_nodes-roots) * log(1. - beta)             # <<<<<<<<<<<<<<
  * 
  *         if beta_add > max_beta_add:
  */
-    __pyx_t_5 = PyFloat_FromDouble(((__pyx_v_roots * log(__pyx_v_beta)) + ((__pyx_v_n_nodes - __pyx_v_roots) * log((1. - __pyx_v_beta))))); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 94; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-    __Pyx_GOTREF(__pyx_t_5);
-    __pyx_t_6 = PyNumber_InPlaceAdd(__pyx_v_beta_add, __pyx_t_5); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 94; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-    __Pyx_GOTREF(__pyx_t_6);
-    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
-    __Pyx_DECREF_SET(__pyx_v_beta_add, __pyx_t_6);
-    __pyx_t_6 = 0;
+    __pyx_v_beta_add = (__pyx_v_beta_add + ((__pyx_v_roots * log(__pyx_v_beta)) + ((__pyx_v_n_nodes - __pyx_v_roots) * log((1. - __pyx_v_beta)))));
 
-    /* "ml.pyx":96
+    /* "ml.pyx":97
  *         beta_add += roots * log(beta) + (n_nodes-roots) * log(1. - beta)
  * 
  *         if beta_add > max_beta_add:             # <<<<<<<<<<<<<<
  *             max_beta = beta
  *             max_roots = roots
  */
-    __pyx_t_6 = PyObject_RichCompare(__pyx_v_beta_add, __pyx_v_max_beta_add, Py_GT); __Pyx_XGOTREF(__pyx_t_6); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 96; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-    __pyx_t_29 = __Pyx_PyObject_IsTrue(__pyx_t_6); if (unlikely(__pyx_t_29 < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 96; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-    __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
+    __pyx_t_29 = ((__pyx_v_beta_add > __pyx_v_max_beta_add) != 0);
     if (__pyx_t_29) {
 
-      /* "ml.pyx":97
+      /* "ml.pyx":98
  * 
  *         if beta_add > max_beta_add:
  *             max_beta = beta             # <<<<<<<<<<<<<<
  *             max_roots = roots
  *             max_beta_add = beta_add
  */
-      __pyx_t_6 = PyFloat_FromDouble(__pyx_v_beta); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 97; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-      __Pyx_GOTREF(__pyx_t_6);
-      __Pyx_XDECREF_SET(__pyx_v_max_beta, __pyx_t_6);
-      __pyx_t_6 = 0;
+      __pyx_v_max_beta = __pyx_v_beta;
 
-      /* "ml.pyx":98
+      /* "ml.pyx":99
  *         if beta_add > max_beta_add:
  *             max_beta = beta
  *             max_roots = roots             # <<<<<<<<<<<<<<
@@ -2878,15 +2878,14 @@ static PyObject *__pyx_pf_2ml_ml(CYTHON_UNUSED PyObject *__pyx_self, PyObject *_
  */
       __pyx_v_max_roots = __pyx_v_roots;
 
-      /* "ml.pyx":99
+      /* "ml.pyx":100
  *             max_beta = beta
  *             max_roots = roots
  *             max_beta_add = beta_add             # <<<<<<<<<<<<<<
  *             # print 'beta:', max_beta, 'add:', max_beta_add, 'roots:', max_roots
  * 
  */
-      __Pyx_INCREF(__pyx_v_beta_add);
-      __Pyx_DECREF_SET(__pyx_v_max_beta_add, __pyx_v_beta_add);
+      __pyx_v_max_beta_add = __pyx_v_beta_add;
       goto __pyx_L21;
     }
     __pyx_L21:;
@@ -2895,29 +2894,22 @@ static PyObject *__pyx_pf_2ml_ml(CYTHON_UNUSED PyObject *__pyx_self, PyObject *_
  *     max_beta_add = float('-inf')
  *     # iterate over all victim nodes to find the optimal threshold
  *     for beta in np.arange(0.001, .2, .002):             # <<<<<<<<<<<<<<
- *         thresh = log(beta/(3012*(1.-beta)))
+ *         thresh = log(beta/(3012.*(1.-beta)))
  *         # print 'beta:', beta, 'thresh:', thresh, 'infected:', len(probs[probs>=thresh])
  */
   }
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
 
-  /* "ml.pyx":102
+  /* "ml.pyx":103
  *             # print 'beta:', max_beta, 'add:', max_beta_add, 'roots:', max_roots
  * 
  *     ll += max_beta_add             # <<<<<<<<<<<<<<
  *     roots = max_roots
  *     beta = max_beta
  */
-  __pyx_t_4 = PyFloat_FromDouble(__pyx_v_ll); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 102; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-  __Pyx_GOTREF(__pyx_t_4);
-  __pyx_t_6 = PyNumber_InPlaceAdd(__pyx_t_4, __pyx_v_max_beta_add); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 102; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-  __Pyx_GOTREF(__pyx_t_6);
-  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
-  __pyx_t_24 = __pyx_PyFloat_AsDouble(__pyx_t_6); if (unlikely((__pyx_t_24 == (npy_double)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 102; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
-  __pyx_v_ll = __pyx_t_24;
+  __pyx_v_ll = (__pyx_v_ll + __pyx_v_max_beta_add);
 
-  /* "ml.pyx":103
+  /* "ml.pyx":104
  * 
  *     ll += max_beta_add
  *     roots = max_roots             # <<<<<<<<<<<<<<
@@ -2926,42 +2918,40 @@ static PyObject *__pyx_pf_2ml_ml(CYTHON_UNUSED PyObject *__pyx_self, PyObject *_
  */
   __pyx_v_roots = __pyx_v_max_roots;
 
-  /* "ml.pyx":104
+  /* "ml.pyx":105
  *     ll += max_beta_add
  *     roots = max_roots
  *     beta = max_beta             # <<<<<<<<<<<<<<
  *     # print n_nodes, n_roots, n_victims, max_i, roots
  *     return (beta, roots, ll)
  */
-  if (unlikely(!__pyx_v_max_beta)) { __Pyx_RaiseUnboundLocalError("max_beta"); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 104; __pyx_clineno = __LINE__; goto __pyx_L1_error;} }
-  __pyx_t_24 = __pyx_PyFloat_AsDouble(__pyx_v_max_beta); if (unlikely((__pyx_t_24 == (npy_double)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 104; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-  __pyx_v_beta = __pyx_t_24;
+  __pyx_v_beta = __pyx_v_max_beta;
 
-  /* "ml.pyx":106
+  /* "ml.pyx":107
  *     beta = max_beta
  *     # print n_nodes, n_roots, n_victims, max_i, roots
  *     return (beta, roots, ll)             # <<<<<<<<<<<<<<
  */
   __Pyx_XDECREF(__pyx_r);
-  __pyx_t_6 = PyFloat_FromDouble(__pyx_v_beta); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 106; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-  __Pyx_GOTREF(__pyx_t_6);
-  __pyx_t_4 = __Pyx_PyInt_From_int(__pyx_v_roots); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 106; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __pyx_t_4 = PyFloat_FromDouble(__pyx_v_beta); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 107; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
   __Pyx_GOTREF(__pyx_t_4);
-  __pyx_t_5 = PyFloat_FromDouble(__pyx_v_ll); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 106; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
-  __Pyx_GOTREF(__pyx_t_5);
-  __pyx_t_3 = PyTuple_New(3); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 106; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __pyx_t_18 = __Pyx_PyInt_From_int(__pyx_v_roots); if (unlikely(!__pyx_t_18)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 107; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_18);
+  __pyx_t_3 = PyFloat_FromDouble(__pyx_v_ll); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 107; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
   __Pyx_GOTREF(__pyx_t_3);
-  PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_6);
-  __Pyx_GIVEREF(__pyx_t_6);
-  PyTuple_SET_ITEM(__pyx_t_3, 1, __pyx_t_4);
+  __pyx_t_5 = PyTuple_New(3); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 107; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_5);
+  PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_4);
   __Pyx_GIVEREF(__pyx_t_4);
-  PyTuple_SET_ITEM(__pyx_t_3, 2, __pyx_t_5);
-  __Pyx_GIVEREF(__pyx_t_5);
-  __pyx_t_6 = 0;
+  PyTuple_SET_ITEM(__pyx_t_5, 1, __pyx_t_18);
+  __Pyx_GIVEREF(__pyx_t_18);
+  PyTuple_SET_ITEM(__pyx_t_5, 2, __pyx_t_3);
+  __Pyx_GIVEREF(__pyx_t_3);
   __pyx_t_4 = 0;
-  __pyx_t_5 = 0;
-  __pyx_r = __pyx_t_3;
+  __pyx_t_18 = 0;
   __pyx_t_3 = 0;
+  __pyx_r = __pyx_t_5;
+  __pyx_t_5 = 0;
   goto __pyx_L0;
 
   /* "ml.pyx":38
@@ -3002,9 +2992,6 @@ static PyObject *__pyx_pf_2ml_ml(CYTHON_UNUSED PyObject *__pyx_self, PyObject *_
   __Pyx_XDECREF((PyObject *)__pyx_v_probs);
   __Pyx_XDECREF((PyObject *)__pyx_v_probs_fail);
   __Pyx_XDECREF((PyObject *)__pyx_v_probs_nv);
-  __Pyx_XDECREF(__pyx_v_max_beta_add);
-  __Pyx_XDECREF(__pyx_v_beta_add);
-  __Pyx_XDECREF(__pyx_v_max_beta);
   __Pyx_XDECREF(__pyx_v_w1);
   __Pyx_XDECREF(__pyx_v_w2);
   __Pyx_XDECREF(__pyx_v_w3);
@@ -5351,7 +5338,7 @@ static int __Pyx_InitCachedConstants(void) {
  *     max_beta_add = float('-inf')
  *     # iterate over all victim nodes to find the optimal threshold
  *     for beta in np.arange(0.001, .2, .002):             # <<<<<<<<<<<<<<
- *         thresh = log(beta/(3012*(1.-beta)))
+ *         thresh = log(beta/(3012.*(1.-beta)))
  *         # print 'beta:', beta, 'thresh:', thresh, 'infected:', len(probs[probs>=thresh])
  */
   __pyx_tuple_ = PyTuple_Pack(3, __pyx_float_0_001, __pyx_float__2, __pyx_float__002); if (unlikely(!__pyx_tuple_)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 86; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
@@ -5431,7 +5418,7 @@ static int __Pyx_InitCachedConstants(void) {
  *        DTYPE_t alpha, DTYPE_t delta):
  *     cdef:
  */
-  __pyx_tuple__8 = PyTuple_Pack(33, __pyx_n_s_root_victims, __pyx_n_s_victims, __pyx_n_s_non_victims, __pyx_n_s_age, __pyx_n_s_alpha, __pyx_n_s_delta, __pyx_n_s_n_roots, __pyx_n_s_n_victims, __pyx_n_s_n_nodes, __pyx_n_s_roots, __pyx_n_s_i, __pyx_n_s_dist, __pyx_n_s_dt, __pyx_n_s_t, __pyx_n_s_l, __pyx_n_s_beta, __pyx_n_s_ll, __pyx_n_s_parents, __pyx_n_s_failures, __pyx_n_s_successes, __pyx_n_s_probs, __pyx_n_s_probs_fail, __pyx_n_s_probs_nv, __pyx_n_s_max_beta_add, __pyx_n_s_thresh, __pyx_n_s_beta_add, __pyx_n_s_max_beta, __pyx_n_s_max_roots, __pyx_n_s_w1, __pyx_n_s_w2, __pyx_n_s_w3, __pyx_n_s_genexpr, __pyx_n_s_genexpr); if (unlikely(!__pyx_tuple__8)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 38; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __pyx_tuple__8 = PyTuple_Pack(33, __pyx_n_s_root_victims, __pyx_n_s_victims, __pyx_n_s_non_victims, __pyx_n_s_age, __pyx_n_s_alpha, __pyx_n_s_delta, __pyx_n_s_n_roots, __pyx_n_s_n_victims, __pyx_n_s_n_nodes, __pyx_n_s_roots, __pyx_n_s_i, __pyx_n_s_dist, __pyx_n_s_dt, __pyx_n_s_t, __pyx_n_s_l, __pyx_n_s_beta, __pyx_n_s_ll, __pyx_n_s_beta_add, __pyx_n_s_max_beta, __pyx_n_s_max_beta_add, __pyx_n_s_parents, __pyx_n_s_failures, __pyx_n_s_successes, __pyx_n_s_probs, __pyx_n_s_probs_fail, __pyx_n_s_probs_nv, __pyx_n_s_thresh, __pyx_n_s_max_roots, __pyx_n_s_w1, __pyx_n_s_w2, __pyx_n_s_w3, __pyx_n_s_genexpr, __pyx_n_s_genexpr); if (unlikely(!__pyx_tuple__8)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 38; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
   __Pyx_GOTREF(__pyx_tuple__8);
   __Pyx_GIVEREF(__pyx_tuple__8);
   __pyx_codeobj__9 = (PyObject*)__Pyx_PyCode_New(6, 0, 33, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__8, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_ben_Documents_Cascade_Pro, __pyx_n_s_ml, 38, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__9)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 38; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
@@ -6843,10 +6830,6 @@ bad:
     return (double)-1;
 }
 
-static CYTHON_INLINE void __Pyx_RaiseUnboundLocalError(const char *varname) {
-    PyErr_Format(PyExc_UnboundLocalError, "local variable '%s' referenced before assignment", varname);
-}
-
 static CYTHON_INLINE void __Pyx_ErrRestore(PyObject *type, PyObject *value, PyObject *tb) {
 #if CYTHON_COMPILING_IN_CPYTHON
     PyObject *tmp_type, *tmp_value, *tmp_tb;