"""
VStatePic.py data_dir y_name

Creates varg_stateN (N in 0..11) in the directory named by data

y[t][0] is a numpy.array containing the observation
y[t][1] is a numpy.array containing the context
context = (y[t-1][0],y[t-2][0],...,y[t-taumax][0],1.0)

Copyright (c) 2005, 2007 Andrew Fraser
This file is part of HMM_DS_Code.

HMM_DS_Code is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.

HMM_DS_Code 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
this program; if not, write to the Free Software Foundation, Inc.,
59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
"""

import sys, os.path, math, random, scipy
data_dir,y_name = sys.argv[1:3]
import VARG

def randomP(A):
    """ Fill allocated array A with random normalized probability
    """
    sum = 0
    for i in xrange(len(A)):
        x = random.random()
        sum += x
        A[i] = x
    A /= sum
    return A

def MakeVARG_HMM(Nstates,Odim,Cdim):
    # Returns a normalized random initial model
    P_S0 = randomP(scipy.zeros(Nstates))
    P_S0_ergodic = randomP(scipy.zeros(Nstates))
    P_ScS = scipy.zeros((Nstates,Nstates))
    for row in P_ScS:
        randomP(row)
    As = Nstates*[None]
    ICovs = Nstates*[None]
    for i in xrange(Nstates):
        ICovs[i] = scipy.eye(Odim)/100
        As[i] = scipy.zeros((Odim,Cdim))
        for row in As[i]:
            randomP(row)
    return VARG.VARG_HMM(P_S0,P_S0_ergodic,P_ScS,As,ICovs,a=30.,b=.1)
# Recall:   Cov = (b * scipy.eye(dim_Y) + ZZT)/(a + sum_w[i])

# Read data and construct Y.
Odim,Cdim,N_states = (3,4,12)
Y = []
f = file(os.path.join(data_dir, y_name), 'r')
tail_one = scipy.ones(1,scipy.float32) # Because Odim = 3 and Cdim = 4
data = scipy.zeros(Odim,scipy.float32)
for line in f.xreadlines():
    context = scipy.concatenate((data,tail_one))
    data = scipy.array(map(float,line.split()),scipy.float32)
    Y.append([data,context])
f.close()
Y.pop(0) # first y has no context, delete it

# Now for each t, y[t] = [observation,context] where observation is
# the 3-d state at time t and context is 4-d; the 3-d state at time
# t-1 with 1.0 appended.  1.0 lets the state use a column of A to
# provide a fixed offset
random.seed(6) # 96 is interesting too
model = MakeVARG_HMM(N_states,Odim,Cdim) # Make a random initial model
model.train(Y, 25)
states = model.decode(Y)                 # Do Viterbi decoding
vs = N_states*[None]
for s in xrange(N_states):
    vs[s] = []  # vs[s] is a list of vectors that belong in state s.
for t in range(0,len(states)):
    vs[states[t]].append(Y[t][0])
for s in xrange(N_states):
    f = file(os.path.join(data_dir, 'varg_state'+str(s)), 'w')
    for v in vs[s]:
        print >>f, v[0], v[1], v[2]
                                 
#Local Variables:
#mode:python
#End: