How to use the hmmlearn.hmm function in hmmlearn

__author__ = 'conghuai'

import numpy as np
from hmmlearn import hmm

np.random.seed(42)
model = hmm.GaussianHMM(n_components=3, covariance_type='full')
model.startprob_ = np.array([0.6, 0.3, 0.1])
model.transmat_ = np.array([[0.7, 0.2, 0.1],
                           [0.3, 0.5, 0.2],
                           [0.3, 0.3, 0.4]])
model.means_ = np.array([[0.0, 0.0], [3.0, -3.0], [5.0, 10.0]])
model.covars_ = np.tile(np.identity(2), (3, 1, 1))
X, Z = model.sample(100)

remodel = hmm.GaussianHMM(n_components=3, covariance_type='full', n_iter=100)
remodel.fit(X)
Z2 = remodel.predict(X)

if method == 'hmm-tumor':
        state_means = np.column_stack([[-4, -1, 0, .585, 1]])
    else:
        state_means = np.column_stack([[-1, 0, .585]])
    kwargs = {
        'n_components': len(state_means),
        'means_prior': state_means,
        'means_weight': .5,
    }
    if method == 'hmm-germline':
        # Don't alter the state means when training the model
        kwargs.update({
            'init_params': 'stc',
            'params': 'stc',
        })
    model = hmm.GaussianHMM(covariance_type='diag', n_iter=100,
                            random_state=0xA5EED, **kwargs)
    if method == 'germline':
        model.means_ = state_means

    # Create observation matrix
    # TODO incorporate weights -- currently handled by smoothing
    # TODO incorporate inter-bin distances
    cnarr = cnarr.autosomes()
    freqs = as_observation_matrix(cnarr)

    if cnarr.chromosome.nunique() == 1:
        model.fit(freqs)
    else:
        model.fit(freqs, lengths=chrom_arm_lengths(cnarr))
    return model

# -*- coding:utf-8 -*-  
__author__ = 'conghuai'

import numpy as np
from hmmlearn import hmm

np.random.seed(42)
model = hmm.GaussianHMM(n_components=3, covariance_type='full')
model.startprob_ = np.array([0.6, 0.3, 0.1])
model.transmat_ = np.array([[0.7, 0.2, 0.1],
                           [0.3, 0.5, 0.2],
                           [0.3, 0.3, 0.4]])
model.means_ = np.array([[0.0, 0.0], [3.0, -3.0], [5.0, 10.0]])
model.covars_ = np.tile(np.identity(2), (3, 1, 1))
X, Z = model.sample(100)

remodel = hmm.GaussianHMM(n_components=3, covariance_type='full', n_iter=100)
remodel.fit(X)
Z2 = remodel.predict(X)

def hmmpred(obs_seq, na3, forw_rerv, hmmoptions, commonOptions):
	obs_seq = obs_seq.replace('-', '')
	#obs_seq = obs_seq.replace('N', ''); obs_seq = obs_seq.replace('n', '');
	
	bp = getBasePair()
	len_repPat = printHMMmatrix.get_len_repPat(na3, commonOptions)

	trainsmat, startprob, emisionmat, obs_symbols, states, numStates, numSymbols, state3class, tol_info = hmmoptions
	hmmmodel = hmm.MultinomialHMM(numStates)
	hmmmodel.transmat_ = trainsmat;
	hmmmodel.startprob_ = startprob;
	hmmmodel.emissionprob_ = emisionmat;
	hmmmodel.n_features = numSymbols

	myobs = []
	for osi in range(len(obs_seq)): 
		myobs.append((np.where(obs_symbols==obs_seq[osi]))[0][0]) 

	logprob, predstats = hmmmodel.decode(np.array([myobs]).T, algorithm="viterbi")

	newstr, ststar, pre0 = getPred(predstats, obs_seq, state3class, len_repPat)
	if myheader.cur_M_STAT <= myheader.M_DEBUG: #int(len(newstr)/float(len_repPat)+0.5)<14: #False: #True: #False: #int(len(newstr)/float(len_repPat)) in [8,13]:
		print ('hmmB:', obs_seq, int(len(newstr)/float(len_repPat)+0.5))
		psstr = []
		for ps in predstats: psstr.append(str(ps))

#    
    # Find the frames that match to each clustered 3DMM state and set a uniform PDF for these frames as the emission probabilities
#    state2frame = [None] * N
    B = 0.01*np.ones((numFramesSiro, N))
#    for i in range(N):
#        state2frame[i] = np.nonzero(stateLabels == i)[0].tolist()
#        B[state2frame[i], i] = 1. / len(state2frame[i])
#    
    B[np.arange(numFramesSiro), stateLabels] = 1
    B /= B.sum(1)[:, np.newaxis]
#    
    # Use a uniform PDF over all frames as the initial distribution
    pi = np.ones(numFramesSiro) / numFramesSiro
    
    # Set up the HMM
    model = hmm.MultinomialHMM(n_components = numFramesSiro)
    model.startprob_ = pi
    model.transmat_ = A.toarray()
    model.emissionprob_ = B
    
    frames = model.predict(stateSeq_kuro.reshape(-1, 1))
#    np.save('kuroSelectedFrames2.npy', frames)
    
    # Animations
    s = stateShapes.reshape((N, 3, mouthIdx.size), order = 'F')
    v = mouthVertices.reshape((numFramesSiro, 3, mouthIdx.size), order = 'F')
    animate(v[:240], mouthFace, 'color/1_origSiroFrames', m.texMean[:, mouthIdx])
    animate(s[stateLabels[:240]], mouthFace, 'color/2_siroKmeansStateSeq', m.texMean[:, mouthIdx])
    animate(s[stateSeq_siro[:240]], mouthFace, 'color/3_siroHMMStateSeq', m.texMean[:, mouthIdx])
    animate(s[stateSeq_kuro[:240]], mouthFace, 'color/4_kuroHMMStateSeq', m.texMean[:, mouthIdx])
    animate(v[selectedFrames[:240]], mouthFace, 'color/5_closestFramesToKuroFromSiroViaRTS', m.texMean[:, mouthIdx])
    animate(v[frames[:240]], mouthFace, 'color/6_closestFramesToKuroFromSiroViaHMM', m.texMean[:, mouthIdx])

def test_anomaly_detection_through_hidden_state(model, threshold_dict):
    normial_data_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'tag_7')
    normial_data_group_by_tag = util.get_anomaly_data_for_labelled_case(training_config, normial_data_path)
    temp = normial_data_group_by_tag.values()
    lengths = []
    for i in range(len(temp)):
        lengths.append(temp[i][1].shape[0])
        if i == 0:
            x = temp[i][1]
        else:
            x = np.concatenate((x, temp[i][1]), axis = 0)
            
    if issubclass(type(model), hmmlearn.hmm._BaseHMM):
        zhats = model.predict(x)
        logliks = model.score(x)
        for i in range(len(x)):
            zhat = model.predict(x[i].reshape(1,-1))
            curr_loglik = model.score(x[i].reshape(1,-1))
            ipdb.set_trace()
            if curr_loglik < threshold_dict[zhat[0]]:
                print 'anomaly'
            else:
                print 'success'
                
    else:
        prev_sample = None
        for isampe in range(len(x)):
            if prev_sample is None:
                prev_sample = x[isampe]

min_sm = min(lenF, lenL)
        F = F[:, 0:min_sm]
        flags = flags[0:min_sm]

        flagsNew = []
        for j, fl in enumerate(flags):      # append features and labels
            flagsNew.append(classes_all.index(class_names[flags[j]]))

        flags_all = numpy.append(flags_all, numpy.array(flagsNew))

        if i == 0:
            f_all = F
        else:
            f_all = numpy.concatenate((f_all, F), axis=1)
    start_prob, transmat, means, cov = trainHMM_computeStatistics(f_all, flags_all)        # compute HMM statistics
    hmm = hmmlearn.hmm.GaussianHMM(start_prob.shape[0], "diag")      # train HMM
    hmm.startprob_ = start_prob
    hmm.transmat_ = transmat        
    hmm.means_ = means
    hmm.covars_ = cov

    fo = open(hmm_model_name, "wb")   # save HMM model
    cPickle.dump(hmm, fo, protocol=cPickle.HIGHEST_PROTOCOL)
    cPickle.dump(classes_all, fo, protocol=cPickle.HIGHEST_PROTOCOL)
    cPickle.dump(mt_win, fo, protocol=cPickle.HIGHEST_PROTOCOL)
    cPickle.dump(mt_step, fo, protocol=cPickle.HIGHEST_PROTOCOL)
    fo.close()

    return hmm, classes_all