######################################################################
# Copyright (C) 2011,2012 Jaakko Luttinen
#
# This file is licensed under Version 3.0 of the GNU General Public
# License. See LICENSE for a text of the license.
######################################################################

######################################################################
# This file is part of BayesPy.
#
# BayesPy is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License version 3 as
# published by the Free Software Foundation.
#
# BayesPy is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with BayesPy.  If not, see <http://www.gnu.org/licenses/>.
######################################################################


import numpy as np
import matplotlib.pyplot as plt
import time

from bayespy.utils import utils
from bayespy.plot import plotting as myplt
import bayespy.inference.vmp.nodes as EF

import imp
imp.reload(utils)
imp.reload(myplt)
imp.reload(EF)

def categorical_model(M, D):

    p = EF.Dirichlet(1*np.ones(D), name='p')
    z = EF.Categorical(p, plates=(M,), name='z')
    return (z, p)


def run(M=30, D=5):

    # Generate data
    y = np.random.randint(D, size=(M,))

    # Construct model
    (z, p) = categorical_model(M, D)

    # Initialize nodes
    p.update()
    z.update()

    # Observe the data with randomly missing values
    mask = np.random.rand(M) < 0.5 # randomly missing
    z.observe(y, mask)

    # Inference loop.
    L_last = -np.inf
    for i in range(10):
        t = time.clock()

        # Update nodes
        p.update()
        z.update()

        # Compute lower bound
        L_p = p.lower_bound_contribution()
        L_z = z.lower_bound_contribution()
        L = L_p + L_z

        # Check convergence
        print("Iteration %d: loglike=%e (%.3f seconds)" % (i+1, L, time.clock()-t))
        if L_last > L:
            L_diff = (L_last - L)

        if L - L_last < 1e-12:
            print("Converged.")

        L_last = L


    z.show()
    p.show()

if __name__ == '__main__':
    run()


@@ -1030,17 +1030,16 @@
 ######%0A%0A
-%0A
 import n
@@ -1052,52 +1052,8 @@
 s np
-%0Aimport matplotlib.pyplot as plt%0Aimport time
 %0A%0Afr
@@ -1076,21 +1076,22 @@
  import 
-utils
+random
 %0Afrom ba
@@ -1099,21 +1099,16 @@
 espy
-.plot
  import 
 plot
@@ -1107,273 +1107,55 @@
 ort 
-plotting as myplt%0Aimport bayespy.inference.vmp.nodes as EF%0A%0Aimport imp%0Aimp.reload(utils)%0Aimp.reload(myplt)%0Aimp.reload(EF)%0A%0Adef categorical_model(M, D):%0A%0A    p = EF.Dirichlet(1*np.ones(D), name='p')%0A    z = EF.Categorical(p, plates=(M,), name='z')%0A    return (z, p)%0A
+nodes%0A%0Afrom bayespy.inference.vmp.vmp import VB
 %0A%0Ade
@@ -1264,774 +1264,364 @@
     
-(z, p) = categorical_model(M, D)%0A
+p = nodes.Dirichlet(1*np.ones(D),
 %0A    
-# Initialize nodes%0A    p.update()%0A    z.update()%0A%0A    # Observe the data with randomly missing values%0A    mask = np.random.rand(M) %3C 0.5 # randomly missing%0A    z.observe(y, mask)%0A%0A    # Inference loop.%0A    L_last = -np.inf%0A    for i in range(10):%0A        t = time.clock()%0A%0A        # Update nodes%0A        p.update()%0A        z.update()%0A%0A        # Compute lower bound%0A        L_p = p.lower_bound_contribution()%0A        L_z = z.lower_bound_contribution()%0A        L = L_p + L_z%0A%0A        # Check convergence%0A        print(%22Iteration %25d: loglike=%25e (%25.3f seconds)%22 %25 (i+1, L, time.clock()-t))%0A        if L_last %3E L:%0A            L_diff = (L_last - L
+                    name='p')%0A    z = nodes.Categorical(p, %0A                          plates=(M,), %0A                          name='z')%0A%0A    # Observe the data with randomly missing values%0A    mask = random.mask(M, p=0.5)%0A    z.observe(y, mask=mask
 )%0A%0A    
-    if L - L_last %3C 1e-12:%0A            print(%22Converged.%22)%0A%0A        L_last = L%0A%0A
+# Run VB-EM%0A    Q = VB(p, z)%0A    Q.update()%0A%0A    # Show results
 %0A   

Update demo for categorical distribution