1 files changed, 46 insertions, 0 deletions
diff --git a/gi/pf/quasi_model2.h b/gi/pf/quasi_model2.h
new file mode 100644
index 00000000..0095289f
--- /dev/null
+++ b/gi/pf/quasi_model2.h
@@ -0,0 +1,46 @@
+#ifndef _QUASI_MODEL2_H_
+#define _QUASI_MODEL2_H_
+
+#include <vector>
+#include <cmath>
+#include "prob.h"
+#include "array2d.h"
+
+struct QuasiModel2 {
+  explicit QuasiModel2(double alpha, double pnull = 0.1) :
+      alpha_(alpha),
+      pnull_(pnull),
+      pnotnull_(1 - pnull),
+      z_(1000,1000) {}
+  // a_j = 0 => NULL; src_len does *not* include null
+  prob_t Pa_j(unsigned a_j, unsigned j, unsigned src_len, unsigned trg_len) const {
+    if (!a_j) return pnull_;
+    std::vector<prob_t>& zv = z_(src_len, trg_len);
+    if (zv.size() == 0)
+      zv.resize(trg_len);
+    
+    prob_t& z = zv[j];
+    if (z.is_0()) z = ComputeZ(j, src_len, trg_len);
+
+    prob_t p;
+    p.logeq(-fabs(double(a_j - 1) / src_len - double(j) / trg_len) * alpha_);
+    p *= pnotnull_;
+    p /= z;
+    return p;
+  }
+ private:
+  prob_t ComputeZ(unsigned j, unsigned src_len, unsigned trg_len) const {
+    prob_t p, z = prob_t::Zero();
+    for (int a_j = 1; a_j <= src_len; ++a_j) {
+      p.logeq(-fabs(double(a_j - 1) / src_len - double(j) / trg_len) * alpha_);
+      z += p;
+    }
+    return z;
+  }
+  double alpha_;
+  const prob_t pnull_;
+  const prob_t pnotnull_;
+  mutable Array2D<std::vector<prob_t> > z_;
+};
+
+#endif