mmm.c

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 */


#include <jwsc/config.h>

#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <inttypes.h>

#include "jwsc/base/limits.h"
#include "jwsc/vector/vector.h"
#include "jwsc/vector/vector_math.h"
#include "jwsc/matrix/matrix.h"
#include "jwsc/matrix/matrix_math.h"
#include "jwsc/matblock/matblock.h"
#include "jwsc/prob/mv_pmf.h"
#include "jwsc/stat/stat.h"
#include "jwsc/stat/kmeans.h"
#include "jwsc/stat/mmm.h"


static void mmm_e_step_d
(
    const Matblock_d*   mu,
    const Vector_d*     pi,
    Matrix_d*           gamma,
    const Matblock_u32* x
)
{
    double   p;
    double   C_n;
    uint32_t k, K;
    uint32_t n, N;
    uint32_t d, D;
    uint32_t m, M;

    Matrix_u32* x_n  = 0;
    Matrix_d*   mu_k = 0;

    K = pi->num_elts;
    N = x->num_mats;
    D = x->num_rows;
    M = x->num_cols;

    create_zero_matrix_d(&gamma, N, K);

    for (n = 0; n < N; n++)
    {
        copy_matrix_from_matblock_u32(&x_n, x, MATBLOCK_MAT_MATRIX, n);
        C_n = 0;

        for (k = 0; k < K; k++)
        {
            copy_matrix_from_matblock_d(&mu_k, mu, MATBLOCK_MAT_MATRIX, k);

            p = 1.0;
            for (d = 0; d < D; d++)
            {
                for (m = 0; m < M; m++)
                {
                    if (x_n->elts[ d ][ m ])
                    {
                        p *= mu_k->elts[ d ][ m ];
                    }
                }
            }

            gamma->elts[ n ][ k ] = pi->elts[ k ] * p;

            C_n += gamma->elts[ n ][ k ];
        }

        assert(C_n > 0);

        for (k = 0; k < K; k++)
        {
            gamma->elts[ n ][ k ] /= C_n;
        }
    }

    free_matrix_u32(x_n);
    free_matrix_d(mu_k);
}


static void mmm_m_step_d
(
    Matblock_d*         mu,
    Vector_d*           pi,
    const Matrix_d*     gamma,
    const Matblock_u32* x
)
{
    uint32_t k, K;
    uint32_t n, N;
    uint32_t d, D;
    uint32_t m, M;

    Matrix_d* mu_k = NULL;

    K = pi->num_elts;
    N = x->num_mats;
    D = x->num_rows;
    M = x->num_cols;

    create_zero_vector_d(&pi, K);

    // mu
    create_matblock_d(&mu, K, D, M);
    for (k = 0; k < K; k++)
    {
        create_zero_matrix_d(&mu_k, D, M);

        for (n = 0; n < N; n++)
        {
            for (d = 0; d < D; d++)
            {
                for (m = 0; m < M; m++)
                {
                    mu_k->elts[ d ][ m ] += 
                        gamma->elts[ n ][ k ]*x->elts[ n ][ d ][ m ];
                }
            }

            pi->elts[ k ] += gamma->elts[ n ][ k ];
        }

        assert(pi->elts[ k ] > 0);

        multiply_matrix_by_scalar_d(&mu_k, mu_k, 1.0 / pi->elts[ k ]);
        copy_matrix_into_matblock_d(&mu, mu, mu_k, MATBLOCK_MAT_MATRIX, k);
    }

    // pi
    multiply_vector_by_scalar_d(&pi, pi, 1.0 / N);

    free_matrix_d(mu_k);
}


double mmm_log_likelihood_d
(
    const Matblock_d*   mu,
    const Vector_d*     pi,
    const Matblock_u32* x
)
{
    uint32_t k, K;
    uint32_t n, N;
    uint32_t d, D;
    uint32_t m, M;
    double   p, q;
    double   ll;

    Matrix_u32* x_n  = NULL;
    Matrix_d*   mu_k = NULL;

    K = pi->num_elts;
    N = x->num_mats;
    D = x->num_rows;
    M = x->num_cols;

    ll = 0;

    for (n = 0; n < N; n++)
    {
        p = 0;
        copy_matrix_from_matblock_u32(&x_n, x, MATBLOCK_MAT_MATRIX, n);

        for (k = 0; k < K; k++)
        {
            copy_matrix_from_matblock_d(&mu_k, mu, MATBLOCK_MAT_MATRIX, k);

            q = 1.0;
            for (d = 0; d < D; d++)
            {
                for (m = 0; m < M; m++)
                {
                    if (x_n->elts[ d ][ m ])
                    {
                        q *= mu_k->elts[ d ][ m ];
                    }
                }
            }

            p += pi->elts[ k ] * q;

#if defined(JWSC_HAVE_ISNAN) && defined(JWSC_HAVE_ISINF)
            assert(!isnan(p) && !isinf(p));
#endif
        }

        ll += log((p < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p);
    }

    free_matrix_u32(x_n);
    free_matrix_d(mu_k);

    return ll;
}


void mmm_kmeans_init_d
(
    Matblock_d**        mu_out,
    Vector_d**          pi_out,
    const Matblock_u32* x,
    uint32_t            K
)
{
    uint32_t d, D;
    uint32_t m, M;
    uint32_t n, N;
    uint32_t k;
    double   sum;

    Vector_i32* members = NULL;
    Vector_u32* counts  = NULL;
    Matrix_d*   x_d     = NULL;
    Matrix_d*   mu_d    = NULL;

    N = x->num_mats;
    D = x->num_rows;
    M = x->num_cols;

    create_matrix_d(&x_d, N, D*M);
    for (n = 0; n < N; n++)
    {
        for (d = 0; d < D; d++)
        {
            for (m = 0; m < M; m++)
            {
                x_d->elts[ n ][ d*M + m ] = x->elts[ n ][ d ][ m ];
            }
        }
    }

    init_kmeans_randomly_d(&mu_d, x_d, K);
    kmeans_d(&mu_d, &members, x_d, K, 2);

    create_zero_vector_u32(&counts, K);
    for (n = 0; n < N; n++)
    {
        k = members->elts[ n ];
        counts->elts[ k ]++;
    }

    create_vector_d(pi_out, K);
    create_matblock_d(mu_out, K, D, M);

    for (k = 0; k < K; k++)
    {
        for (d = 0; d < D; d++)
        {
            sum = 0;

            for (m = 0; m < M; m++)
            {
                sum += mu_d->elts[ k ][ d*M + m ];
                (*mu_out)->elts[ k ][ d ][ m ] = mu_d->elts[ k ][ d*M + m ];
            }

            assert(fabs(sum - 1.0) < 1.0e-8);
        }

        (*pi_out)->elts[ k ] = counts->elts[ k ] / (double)N;
    }

    free_vector_i32(members);
    free_vector_u32(counts);
    free_matrix_d(x_d);
    free_matrix_d(mu_d);
}


void train_mmm_d
(
    Matblock_d**        mu_out,
    Vector_d**          pi_out,
    Matrix_d**          gamma_out,
    const Matblock_u32* x,
    uint32_t            K
)
{
    double ll, ll_prev;
    double ll_delta;

    Matblock_d* mu;
    Vector_d*   pi;
    Matrix_d*   gamma = 0;

    mmm_kmeans_init_d(mu_out, pi_out, x, K);
    mu    = *mu_out;
    pi    = *pi_out;

    if (gamma_out)
    {
        create_matrix_d(gamma_out, x->num_rows, pi->num_elts);
        gamma = *gamma_out;
    }
    else
    {
        create_matrix_d(&gamma, x->num_rows, pi->num_elts);
    }

    ll       =  log(0);
    ll_delta = -log(0);

    while (ll_delta > 1.0e-10)
    {
        mmm_e_step_d(mu, pi, gamma, x);
        mmm_m_step_d(mu, pi, gamma, x);

        ll_prev = ll;
        ll = mmm_log_likelihood_d(mu, pi, x);

        ll_delta = ll - ll_prev;
        fprintf(stderr, "%e\n", ll_delta);
    }

    if (!gamma_out)
    {
        free_matrix_d(gamma);
    }
}