Skip to content
Snippets Groups Projects
Commit 86424e69 authored by Blaß, Michael's avatar Blaß, Michael :speech_balloon:
Browse files

Finished

parent 2ca013d5
No related branches found
No related tags found
No related merge requests found
#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);
}
...@@ -4,42 +4,32 @@ ...@@ -4,42 +4,32 @@
static PyObject * static PyObject *
hmm_poisson_fit(PyObject* self, PyObject* args) hmm_poisson_fit_em (PyObject *self, PyObject *args)
{ {
/* 1. alloc alpha, beta, pprob PyObject *py_Xtrain = NULL;
* 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_lambda = NULL; PyObject *py_lambda = NULL;
PyObject *py_gamma = NULL; PyObject *py_gamma = NULL;
PyObject *py_delta = NULL; PyObject *py_delta = NULL;
npy_intp m_states = 0;
npy_intp max_iter = 0; npy_intp max_iter = 0;
double tol = 0.0; double tol = 0.0;
PyArrayObject *X = NULL; PyArrayObject *X_t = NULL;
PyArrayObject *_lambda = NULL; PyArrayObject *_lambda = NULL;
PyArrayObject *_gamma = NULL; PyArrayObject *_gamma = NULL;
PyArrayObject *_delta = 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; &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); _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); _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); _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 (_lambda);
Py_XDECREF (_gamma); Py_XDECREF (_gamma);
Py_XDECREF (_delta); Py_XDECREF (_delta);
...@@ -47,16 +37,15 @@ hmm_poisson_EM(PyObject* self, PyObject* args) ...@@ -47,16 +37,15 @@ hmm_poisson_EM(PyObject* self, PyObject* args)
Py_RETURN_NONE; Py_RETURN_NONE;
} }
npy_intp m = PyArray_SIZE (_lambda); PoisHmm *ph = PoisHmm_FromData ((size_t) m_states,
PoisHmm *ph = PoisHmm_FromData ((size_t) m,
PyArray_DATA (_lambda), PyArray_DATA (_lambda),
PyArray_DATA (_gamma), PyArray_DATA (_gamma),
PyArray_DATA (_delta), PyArray_DATA (_delta),
(size_t) max_iter, (scalar) tol); (size_t) max_iter, (scalar) tol);
if (ph == NULL) if (ph == NULL)
{ {
Py_XDECREF (X); Py_XDECREF (X_t);
Py_XDECREF (_lambda); Py_XDECREF (_lambda);
Py_XDECREF (_gamma); Py_XDECREF (_gamma);
Py_XDECREF (_delta); Py_XDECREF (_delta);
...@@ -64,15 +53,15 @@ hmm_poisson_EM(PyObject* self, PyObject* args) ...@@ -64,15 +53,15 @@ hmm_poisson_EM(PyObject* self, PyObject* args)
Py_RETURN_NONE; 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_1d[] = { m_states };
npy_intp dims_2d[] = { m, m }; npy_intp dims_2d[] = { m_states, m_states };
size_t vector_s = (size_t) m * sizeof (scalar); size_t vector_s = (size_t) m_states * sizeof (scalar);
size_t matrix_s = (size_t) m * vector_s; size_t matrix_s = (size_t) m_states * vector_s;
PyArrayObject *lambda_ = Apollon_NewPyArray1d (dims_1d); PyArrayObject *lambda_ = Apollon_NewPyArray1d (dims_1d);
PyArrayObject *gamma_ = Apollon_NewPyArray2d (dims_2d); PyArrayObject *gamma_ = Apollon_NewPyArray2d (dims_2d);
...@@ -90,7 +79,7 @@ hmm_poisson_EM(PyObject* self, PyObject* args) ...@@ -90,7 +79,7 @@ hmm_poisson_EM(PyObject* self, PyObject* args)
(double) ph->aic, (double) ph->bic, (double) ph->aic, (double) ph->bic,
(double) ph->nll, ph->n_iter); (double) ph->nll, ph->n_iter);
Py_DECREF (X); Py_DECREF (X_t);
Py_DECREF (_lambda); Py_DECREF (_lambda);
Py_DECREF (_gamma); Py_DECREF (_gamma);
Py_DECREF (_delta); Py_DECREF (_delta);
...@@ -104,35 +93,15 @@ hmm_poisson_EM(PyObject* self, PyObject* args) ...@@ -104,35 +93,15 @@ hmm_poisson_EM(PyObject* self, PyObject* args)
static PyObject * static PyObject *
hmm_poisson_fwbw(PyObject *self, PyObject *args) hmm_poisson_fwbw(PyObject *self, PyObject *args)
{ {
PyArrayObject *x = NULL; PyErr_SetString (PyExc_NotImplementedError, "");
PyArrayObject *lambda_ = NULL; Py_RETURN_NONE;
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);
/* static PyObject *
if (success == 1) hmm_poisson_em (PyObject *self, PyObject *args)
{
PyObject *out = PyTuple_Pack(3, alpha, beta, lp_prob);
Py_INCREF(out);
return out;
}
else
{ {
PyErr_SetString(PyExc_RuntimeError, "Training failed."); PyErr_SetString (PyExc_NotImplementedError, "");
Py_RETURN_NONE;
}
*/
Py_RETURN_NONE; Py_RETURN_NONE;
} }
...@@ -140,18 +109,20 @@ hmm_poisson_fwbw(PyObject *self, PyObject *args) ...@@ -140,18 +109,20 @@ hmm_poisson_fwbw(PyObject *self, PyObject *args)
static PyObject * static PyObject *
hmm_poisson_viterbi (PyObject *self, PyObject *args) hmm_poisson_viterbi (PyObject *self, PyObject *args)
{ {
PyErr_SetString (PyExc_NotImplementedError, "");
Py_RETURN_NONE; Py_RETURN_NONE;
} }
static PyMethodDef static PyMethodDef
HMM_Methods[] = { HMM_Methods[] = {
{"hmm_poisson_fit", hmm_poisson_fit, METH_VARARGS, {"hmm_poisson_fit_em", hmm_poisson_fit_em, METH_VARARGS,
"hmm_poisson_fit(x, m, _lambda, _gamma, _delta, max_iter, tol)"}, "hmm_poisson_fit_em (x, m, _lambda, _gamma, _delta, max_iter, tol)"},
{"hmm_poisson_EM", hmm_poisson_EM, METH_VARARGS, {"hmm_poisson_fwbw", hmm_poisson_fwbw, METH_VARARGS,
"docstring"}, "docstring"},
{"hmm_poisson_fwbw", hmm_poisson_fwbw, METH_VARARGS, {"hmm_poisson_em", hmm_poisson_em, METH_VARARGS,
"docstring"}, "docstring"},
{"hmm_poisson_viterbi", hmm_poisson_viterbi, METH_VARARGS, {"hmm_poisson_viterbi", hmm_poisson_viterbi, METH_VARARGS,
......
...@@ -23,4 +23,22 @@ ...@@ -23,4 +23,22 @@
2, shape, NULL, NULL, 0, NULL)); 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 */ #endif /* hmm_module_h */
...@@ -2,17 +2,17 @@ import numpy as np ...@@ -2,17 +2,17 @@ import numpy as np
import hmm import hmm
x = np.random.randint(0, 500, 1000) x = np.random.randint(0, 500, 1000)
m = 3
l = np.array([10, 250, 450]) l = np.array([10, 250, 450])
g = np.array([[.8, .1, .1], [.1, .8, .1], [.1, .1, .8]]) g = np.array([[.8, .1, .1], [.1, .8, .1], [.1, .1, .8]])
d = np.array([.3, .3, .3]) 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) (success, lambda_, gamma_, delta_, aic, bic, nll, niter) = hmm.hmm_poisson_fit_em(x, m, l, g, d, 1000, 1e-6)
print(s, aic, bic, nll, niter) print(success, aic, bic, nll, niter)
print(l_) print(lambda_)
print("\n") print("\n")
print(g_) print(gamma_)
print("\n") print("\n")
print(d_) print(delta_)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment