From 86424e696cf19d4d44c9a01ab68ea4c028daaf42 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Michael=20Bla=C3=9F?= <michael.blass@uni-hamburg.de>
Date: Sat, 2 Mar 2019 10:26:11 +0100
Subject: [PATCH] Finished
---
hmm/hmm_module.bckp.c | 161 ------------------------------------------
hmm/hmm_module.c | 105 ++++++++++-----------------
include/hmm_module.h | 18 +++++
test.py | 12 ++--
4 files changed, 62 insertions(+), 234 deletions(-)
delete mode 100644 hmm/hmm_module.bckp.c
diff --git a/hmm/hmm_module.bckp.c b/hmm/hmm_module.bckp.c
deleted file mode 100644
index bf53d68..0000000
--- a/hmm/hmm_module.bckp.c
+++ /dev/null
@@ -1,161 +0,0 @@
-#define NPY_NO_DEPRECATED_API NPY_1_8_API_VERSION
-
-#include "hmm_module.h"
-
-
-static PyObject *
-hmm_poisson_fit(PyObject* self, PyObject* args)
-{
- /* 1. alloc alpha, beta, pprob
- * 2. call fwbw
- * 3. wrap alpha, beta, pprob in arrays
- * 4. return the arrays */
- Py_RETURN_NONE;
-}
-
-
-static PyObject *
-hmm_poisson_EM(PyObject* self, PyObject* args)
-{
- PyArrayObject *x = NULL;
- PyArrayObject *_lambda= NULL;
- PyArrayObject *_gamma = NULL;
- PyArrayObject *_delta = NULL;
- npy_intp max_iter= 0;
- double tol = 0;
-
- if (!PyArg_ParseTuple(args, "OOOOld", &x, &_lambda, &_gamma,
- &_delta, &max_iter, &tol)) return NULL;
-
- npy_intp n = PyArray_SIZE(x);
- npy_intp m = PyArray_SIZE(_lambda);
-
- PoissonHMM *hmm = NewPoissonHMM( m,
- PyArray_DATA(_lambda),
- PyArray_DATA(_gamma),
- PyArray_DATA(_delta),
- (size_t) max_iter, (scalar) tol);
-
- int success = poisson_expectation_maximization(PyArray_DATA(x),
- (size_t) n, hmm);
-
- if (success == 1)
- {
- npy_intp dims_1d[] = { m };
- npy_intp dims_2d[] = { m, m };
- size_t v_size = (size_t) m * sizeof(scalar);
- size_t m_size = (size_t) m * v_size;
-
- PyArrayObject *lambda_ = Apollon_NewPyArray1d(dims_1d);
- PyArrayObject *gamma_ = Apollon_NewPyArray2d(dims_2d);
- PyArrayObject *delta_ = Apollon_NewPyArray1d(dims_1d);
-
- scalar *data = PyArray_DATA(lambda_);
- memcpy (PyArray_DATA(lambda_), hmm->lambda_, v_size);
- memcpy (PyArray_DATA(gamma_), hmm->gamma_, m_size);
- memcpy (PyArray_DATA(delta_), hmm->delta_, v_size);
-
- for (npy_intp i = 0; i < m; i++)
- fprintf(stdout, "%Lf\t", data[i]);
-
- hmm->aic = compute_aic(hmm->nll, m, n);
- hmm->bic = compute_bic(hmm->nll, m, n);
-
- PyObject *out = PyTuple_Pack(7, lambda_, gamma_, delta_,
- PyFloat_FromDouble(hmm->aic), PyFloat_FromDouble(hmm->bic),
- PyFloat_FromDouble(hmm->nll), PyLong_FromSize_t(hmm->n_iter));
-
- Py_INCREF(out);
- DeletePoissonHMM(hmm);
- return out;
- }
- else
- {
- DeletePoissonHMM(hmm);
- PyErr_SetString(PyExc_ValueError, "Training failed.");
- Py_RETURN_NONE;
- }
-}
-
-
-static PyObject *
-hmm_poisson_fwbw(PyObject* self, PyObject* args)
-{
- PyArrayObject *x = NULL;
- PyArrayObject *lambda_= NULL;
- PyArrayObject *gamma_ = NULL;
- PyArrayObject *delta_ = NULL;
-
- if (!PyArg_ParseTuple(args, "OOOO", &x, &lambda_, &gamma_, &delta_))
- return NULL;
-
- npy_intp n = PyArray_SIZE(x);
- npy_intp m = PyArray_SIZE(lambda_);
- npy_intp dims[] = { n, m };
-
- PyArrayObject *alpha = Apollon_NewPyArray2d(dims);
- PyArrayObject *beta = Apollon_NewPyArray2d(dims);
- PyArrayObject *lp_prob = Apollon_NewPyArray2d(dims);
-
- int success = log_poisson_forward_backward (PyArray_DATA(x),
- (size_t) n, (size_t) m,
- PyArray_DATA(lambda_),
- PyArray_DATA(gamma_),
- PyArray_DATA(delta_),
- PyArray_DATA(alpha),
- PyArray_DATA(beta),
- PyArray_DATA(lp_prob));
-
- if (success == 1)
- {
- PyObject *out = PyTuple_Pack(3, alpha, beta, lp_prob);
- Py_INCREF(out);
- return out;
- }
- else
- {
- PyErr_SetString(PyExc_RuntimeError, "Training failed.");
- Py_RETURN_NONE;
- }
-}
-
-
-static PyObject *
-hmm_poisson_viterbi(PyObject* self, PyObject* args)
-{
- Py_RETURN_NONE;
-}
-
-static PyMethodDef
-HMM_Methods[] = {
- {"hmm_poisson_fit", hmm_poisson_fit, METH_VARARGS,
- "hmm_poisson_fit(x, m, _lambda, _gamma, _delta, max_iter, tol)"},
-
- {"hmm_poisson_EM", hmm_poisson_EM, METH_VARARGS,
- "docstring"},
-
- {"hmm_poisson_fwbw", hmm_poisson_fwbw, METH_VARARGS,
- "docstring"},
-
- {"hmm_poisson_viterbi", hmm_poisson_viterbi, METH_VARARGS,
- "docstring"},
-
- {NULL, NULL, 0, NULL}
-};
-
-static struct PyModuleDef
-hmm_module = {
- PyModuleDef_HEAD_INIT,
- "hmm",
- NULL,
- -1,
- HMM_Methods
-};
-
-PyMODINIT_FUNC
-PyInit_hmm(void)
-{
- import_array();
- return PyModule_Create (&hmm_module);
-}
-
diff --git a/hmm/hmm_module.c b/hmm/hmm_module.c
index 045177f..96a63f3 100644
--- a/hmm/hmm_module.c
+++ b/hmm/hmm_module.c
@@ -4,42 +4,32 @@
static PyObject *
-hmm_poisson_fit(PyObject* self, PyObject* args)
+hmm_poisson_fit_em (PyObject *self, PyObject *args)
{
- /* 1. alloc alpha, beta, pprob
- * 2. call fwbw
- * 3. wrap alpha, beta, pprob in arrays
- * 4. return the arrays */
- Py_RETURN_NONE;
-}
-
-
-static PyObject *
-hmm_poisson_EM(PyObject* self, PyObject* args)
-{
- PyObject *py_X = NULL;
+ PyObject *py_Xtrain = NULL;
PyObject *py_lambda = NULL;
PyObject *py_gamma = NULL;
PyObject *py_delta = NULL;
+ npy_intp m_states = 0;
npy_intp max_iter = 0;
double tol = 0.0;
- PyArrayObject *X = NULL;
+ PyArrayObject *X_t = NULL;
PyArrayObject *_lambda = NULL;
PyArrayObject *_gamma = NULL;
PyArrayObject *_delta = NULL;
- if (!PyArg_ParseTuple(args, "OOOOld", &py_X, &py_lambda, &py_gamma,
+ if (!PyArg_ParseTuple(args, "OlOOOld", &py_Xtrain, &m_states, &py_lambda, &py_gamma,
&py_delta, &max_iter, &tol)) return NULL;
- X = (PyArrayObject *) PyArray_FROM_OTF(py_X, NPY_LONG, NPY_ARRAY_IN_ARRAY);
+ X_t = (PyArrayObject *) PyArray_FROM_OTF(py_Xtrain, NPY_LONG, NPY_ARRAY_IN_ARRAY);
_lambda = (PyArrayObject *) PyArray_FROM_OTF(py_lambda, NPY_LONGDOUBLE, NPY_ARRAY_IN_ARRAY);
- _gamma = (PyArrayObject *) PyArray_FROM_OTF(py_gamma, NPY_LONGDOUBLE, NPY_ARRAY_IN_ARRAY);
- _delta = (PyArrayObject *) PyArray_FROM_OTF(py_delta, NPY_LONGDOUBLE, NPY_ARRAY_IN_ARRAY);
+ _gamma = (PyArrayObject *) PyArray_FROM_OTF(py_gamma, NPY_LONGDOUBLE, NPY_ARRAY_IN_ARRAY);
+ _delta = (PyArrayObject *) PyArray_FROM_OTF(py_delta, NPY_LONGDOUBLE, NPY_ARRAY_IN_ARRAY);
- if (X == NULL || _lambda == NULL || _gamma == NULL || _delta == NULL)
+ if (X_t == NULL || _lambda == NULL || _gamma == NULL || _delta == NULL)
{
- Py_XDECREF (X);
+ Py_XDECREF (X_t);
Py_XDECREF (_lambda);
Py_XDECREF (_gamma);
Py_XDECREF (_delta);
@@ -47,16 +37,15 @@ hmm_poisson_EM(PyObject* self, PyObject* args)
Py_RETURN_NONE;
}
- npy_intp m = PyArray_SIZE (_lambda);
-
- PoisHmm *ph = PoisHmm_FromData ((size_t) m,
+ PoisHmm *ph = PoisHmm_FromData ((size_t) m_states,
PyArray_DATA (_lambda),
PyArray_DATA (_gamma),
PyArray_DATA (_delta),
(size_t) max_iter, (scalar) tol);
+
if (ph == NULL)
{
- Py_XDECREF (X);
+ Py_XDECREF (X_t);
Py_XDECREF (_lambda);
Py_XDECREF (_gamma);
Py_XDECREF (_delta);
@@ -64,15 +53,15 @@ hmm_poisson_EM(PyObject* self, PyObject* args)
Py_RETURN_NONE;
}
- DataSet X_train = {PyArray_DATA (X), (size_t) PyArray_SIZE (X)};
+ DataSet X_train = {PyArray_DATA (X_t), (size_t) PyArray_SIZE (X_t)};
- int success = PoisHmm_EM (&X_train, ph);
+ int success = PoisHmm_EM (&X_train, ph);
{
- npy_intp dims_1d[] = { m };
- npy_intp dims_2d[] = { m, m };
- size_t vector_s = (size_t) m * sizeof (scalar);
- size_t matrix_s = (size_t) m * vector_s;
+ npy_intp dims_1d[] = { m_states };
+ npy_intp dims_2d[] = { m_states, m_states };
+ size_t vector_s = (size_t) m_states * sizeof (scalar);
+ size_t matrix_s = (size_t) m_states * vector_s;
PyArrayObject *lambda_ = Apollon_NewPyArray1d (dims_1d);
PyArrayObject *gamma_ = Apollon_NewPyArray2d (dims_2d);
@@ -86,11 +75,11 @@ hmm_poisson_EM(PyObject* self, PyObject* args)
ph->bic = compute_bic(ph->nll, ph->m, X_train.size);
PyObject *out = NULL;
- out = Py_BuildValue("iNNNdddk", success, lambda_, gamma_, delta_,
+ out = Py_BuildValue("iNNNdddk", success, lambda_, gamma_, delta_,
(double) ph->aic, (double) ph->bic,
(double) ph->nll, ph->n_iter);
- Py_DECREF (X);
+ Py_DECREF (X_t);
Py_DECREF (_lambda);
Py_DECREF (_gamma);
Py_DECREF (_delta);
@@ -104,56 +93,38 @@ hmm_poisson_EM(PyObject* self, PyObject* args)
static PyObject *
hmm_poisson_fwbw(PyObject *self, PyObject *args)
{
- PyArrayObject *x = NULL;
- PyArrayObject *lambda_ = NULL;
- PyArrayObject *gamma_ = NULL;
- PyArrayObject *delta_ = NULL;
-
- if (!PyArg_ParseTuple(args, "OOOO", &x, &lambda_, &gamma_, &delta_))
- return NULL;
-
- npy_intp n = PyArray_SIZE (x);
- npy_intp m = PyArray_SIZE (lambda_);
- npy_intp dims[] = { n, m };
-
- PyArrayObject *alpha = Apollon_NewPyArray2d (dims);
- PyArrayObject *beta = Apollon_NewPyArray2d (dims);
- PyArrayObject *pois_pr = Apollon_NewPyArray2d (dims);
-
- /*
- if (success == 1)
- {
- PyObject *out = PyTuple_Pack(3, alpha, beta, lp_prob);
- Py_INCREF(out);
- return out;
- }
- else
- {
- PyErr_SetString(PyExc_RuntimeError, "Training failed.");
- Py_RETURN_NONE;
- }
- */
+ PyErr_SetString (PyExc_NotImplementedError, "");
+ Py_RETURN_NONE;
+}
+
+
+static PyObject *
+hmm_poisson_em (PyObject *self, PyObject *args)
+{
+ PyErr_SetString (PyExc_NotImplementedError, "");
Py_RETURN_NONE;
}
static PyObject *
-hmm_poisson_viterbi(PyObject* self, PyObject* args)
+hmm_poisson_viterbi (PyObject *self, PyObject *args)
{
+ PyErr_SetString (PyExc_NotImplementedError, "");
Py_RETURN_NONE;
}
+
static PyMethodDef
HMM_Methods[] = {
- {"hmm_poisson_fit", hmm_poisson_fit, METH_VARARGS,
- "hmm_poisson_fit(x, m, _lambda, _gamma, _delta, max_iter, tol)"},
-
- {"hmm_poisson_EM", hmm_poisson_EM, METH_VARARGS,
- "docstring"},
+ {"hmm_poisson_fit_em", hmm_poisson_fit_em, METH_VARARGS,
+ "hmm_poisson_fit_em (x, m, _lambda, _gamma, _delta, max_iter, tol)"},
{"hmm_poisson_fwbw", hmm_poisson_fwbw, METH_VARARGS,
"docstring"},
+ {"hmm_poisson_em", hmm_poisson_em, METH_VARARGS,
+ "docstring"},
+
{"hmm_poisson_viterbi", hmm_poisson_viterbi, METH_VARARGS,
"docstring"},
diff --git a/include/hmm_module.h b/include/hmm_module.h
index de3ff10..d07097a 100644
--- a/include/hmm_module.h
+++ b/include/hmm_module.h
@@ -23,4 +23,22 @@
2, shape, NULL, NULL, 0, NULL));
+/** Perform EM for given parameters. Return only fitted params. */
+static PyObject *hmm_poisson_em (PyObject* self, PyObject* args);
+
+
+/** Fit HMM for given parameters using EM. Return all params
+ * and quality measures.
+ */
+static PyObject *hmm_poisson_fit_em (PyObject* self, PyObject* args);
+
+
+/** Return forward backward, and state-dependent probabilites. */
+static PyObject *hmm_poisson_fwbw (PyObject *self, PyObject *args);
+
+
+/** Calculate viterbi path given HMM. */
+static PyObject *hmm_poisson_viterbi (PyObject *self, PyObject *args);
+
+
#endif /* hmm_module_h */
diff --git a/test.py b/test.py
index bfdb69f..e338101 100644
--- a/test.py
+++ b/test.py
@@ -2,17 +2,17 @@ import numpy as np
import hmm
x = np.random.randint(0, 500, 1000)
-
+m = 3
l = np.array([10, 250, 450])
g = np.array([[.8, .1, .1], [.1, .8, .1], [.1, .1, .8]])
d = np.array([.3, .3, .3])
-(s, l_, g_, d_, aic, bic, nll, niter) = hmm.hmm_poisson_EM(x, l, g, d, 1000, 1e-6)
-print(s, aic, bic, nll, niter)
+(success, lambda_, gamma_, delta_, aic, bic, nll, niter) = hmm.hmm_poisson_fit_em(x, m, l, g, d, 1000, 1e-6)
+print(success, aic, bic, nll, niter)
-print(l_)
+print(lambda_)
print("\n")
-print(g_)
+print(gamma_)
print("\n")
-print(d_)
+print(delta_)
--
GitLab