[core] misc bug fixes

Author: terhorst

smcpp/_smcpp.pyx

@@ -170,6 +170,8 @@ cdef class _PyInferenceManager:
             self._im.setRho(rho)
 
     def E_step(self, forward_backward_only=False):
+        if None in (self.theta, self.rho):
+            raise RuntimeError("theta and rho must be set")
         cdef bool fbOnly = forward_backward_only
         with nogil:
             self._im.Estep(fbOnly)

smcpp/commands/posterior.py

@@ -96,6 +96,7 @@ def main(self, args):
             im = _smcpp.PyTwoPopInferenceManager(
                 *n, *a, all_obs, hidden_states, contig.key, args.polarization_error)
         im.theta = j['theta']
+        im.rho = j['rho']
         im.save_gamma = True
         im.model = m
         im.E_step()
@@ -111,6 +112,7 @@ def main(self, args):
         kwargs.update({path + "_sites": obs[:, 0] for path, obs in zip(args.data, all_obs)})
         np.savez_compressed(args.output, hidden_states=hidden_states, **kwargs)
         if args.heatmap:
+            obs = all_obs[0]
             if len(args.data) > 1:
                 logger.error("--heatmap is only supported for one data set")
                 sys.exit(1)

src/conditioned_sfs.cpp

@@ -104,6 +104,8 @@ std::vector<Matrix<T> > incorporate_theta(const std::vector<Matrix<T> > &csfs, d
     std::vector<Matrix<T> > ret(csfs.size());
     for (unsigned int i = 0; i < csfs.size(); ++i)
     {
+        assert(csfs[i](0, 0) == 0.);
+        assert(csfs[i](2, n) == 0.);
         T tauh = csfs[i].sum();
         if (toDouble(tauh) > 1.0 / theta)
             throw improper_sfs_exception();
@@ -112,7 +114,8 @@ std::vector<Matrix<T> > incorporate_theta(const std::vector<Matrix<T> > &csfs, d
             CHECK_NAN(tauh);
             ret[i] = csfs[i] * -expm1(-theta * tauh) / tauh;
             CHECK_NAN(ret[i]);
-        } catch (std::runtime_error)
+        }
+        catch (std::runtime_error)
         {
             std::cout << i << std::endl << csfs[i].template cast<double>() << std::endl;
             std::cout << tauh << std::endl;

src/inference_manager.cpp

@@ -422,7 +422,7 @@ void NPopInferenceManager<P>::recompute_emission_probs()
         }
         else
         {
-            e2(m, 1) = 2. * theta * avg_ct.at(m);
+            e2(m, 1) = -expm1(-2. * theta * avg_ct.at(m));
             e2(m, 0) = 1. - e2(m, 1);
         }
         // CHECK_NAN(e2(m, 1));