Commit 90a72cc7 authored by Emmanuel PIGNAT's avatar Emmanuel PIGNAT
Browse files

changing name

(cherry picked from commit c9d3f966)
parent 43e9f4e3
......@@ -2,9 +2,9 @@ import numpy as np
from model import *
from functions import multi_variate_normal
from scipy.linalg import block_diag
import scipy.sparse as ss
from termcolor import colored
from mvn import MVN
from mvn import MVN, SparseMVN
class GMM(Model):
......@@ -136,7 +136,7 @@ class GMM(Model):
return gmm
def concatenate_gaussian(self, q, get_mvn=True):
def concatenate_gaussian(self, q, get_mvn=True, sparse=False):
"""
Get a concatenated-block-diagonal replication of the GMM with sequence of state
given by q.
......@@ -150,14 +150,20 @@ class GMM(Model):
if not get_mvn:
return np.concatenate([self.mu[i] for i in q]), block_diag(*[self.sigma[i] for i in q])
else:
mvn = MVN()
mvn.mu = np.concatenate([self.mu[i] for i in q])
mvn._sigma = block_diag(*[self.sigma[i] for i in q])
mvn._lmbda = block_diag(*[self.lmbda[i] for i in q])
if not sparse:
mvn = MVN()
mvn.mu = np.concatenate([self.mu[i] for i in q])
mvn._sigma = block_diag(*[self.sigma[i] for i in q])
mvn._lmbda = block_diag(*[self.lmbda[i] for i in q])
else:
mvn = SparseMVN()
mvn.mu = np.concatenate([self.mu[i] for i in q])
mvn._sigma = ss.block_diag([self.sigma[i] for i in q])
mvn._lmbda = ss.block_diag([self.lmbda[i] for i in q])
return mvn
def compute_resp(self, demo=None, dep=None, table=None, marginal=None, ):
def compute_resp(self, demo=None, dep=None, table=None, marginal=None, norm=True):
sample_size = demo.shape[0]
B = np.ones((self.nb_states, sample_size))
......@@ -180,15 +186,22 @@ class GMM(Model):
B[[i], :] *= multi_variate_normal(demo, mu[i, d],
sigma[dGrid][:, :, 0], log=False)
B *= self.priors[:, None]
return B/np.sum(B, axis=0)
if norm:
return B/np.sum(B, axis=0)
else:
return B
def init_params_scikit(self, data):
def init_params_scikit(self, data, cov_type='full'):
from sklearn.mixture import BayesianGaussianMixture, GaussianMixture
gmm_init = GaussianMixture(self.nb_states, 'full', n_init=10, init_params='random')
gmm_init = GaussianMixture(self.nb_states, cov_type, n_init=5, init_params='random')
gmm_init.fit(data)
self.mu = gmm_init.means_
self.sigma = gmm_init.covariances_
if cov_type == 'diag':
self.sigma = np.array([np.diag(gmm_init.covariances_[i]) for i in range(self.nb_states)])
else:
self.sigma = gmm_init.covariances_
self.priors = gmm_init.weights_
self.Trans = np.ones((self.nb_states, self.nb_states)) * 0.01
......@@ -221,7 +234,8 @@ class GMM(Model):
self.priors = np.ones(self.nb_states) / self.nb_states
def em(self, data, reg=1e-8, maxiter=100, minstepsize=1e-5, diag=False, reg_finish=False,
kmeans_init=False, random_init=True, dep_mask=None, verbose=False):
kmeans_init=False, random_init=True, dep_mask=None, verbose=False, only_scikit=False,
no_init=False):
"""
:param data: [np.array([nb_timesteps, nb_dim])]
......@@ -250,15 +264,18 @@ class GMM(Model):
nb_samples = data.shape[0]
if not no_init:
if random_init and not only_scikit:
self.init_params_random(data)
elif kmeans_init and not only_scikit:
self.init_params_kmeans(data)
else:
if diag:
self.init_params_scikit(data, 'diag')
else:
self.init_params_scikit(data, 'full')
if random_init:
self.init_params_random(data)
elif kmeans_init:
self.init_params_kmeans(data)
else:
self.init_params_scikit(data)
if only_scikit: return
data = data.T
LL = np.zeros(nb_max_steps)
......
......@@ -140,7 +140,22 @@ class MVN(object):
self.sigma = np.cov(data.T)
self.lmbda = np.linalg.inv(self.sigma)
def log_prob(self, x):
def log_prob(self, x, marginal=None, reg=None):
"""
:param x:
:param marginal:
:type marginal: slice
:return:
"""
if marginal is not None:
_mu = self.mu[marginal]
_sigma = self.sigma[marginal, marginal]
if reg is not None:
_sigma += np.eye(marginal.stop-marginal.start) * reg
return multi_variate_normal(x, _mu, _sigma)
return multi_variate_normal(x, self.mu, self.sigma)
def transform(self, A, b=None, dA=None, db=None):
......
......@@ -269,6 +269,36 @@ def plot_linear_system(K, b=None, name=None, nb_sub=10, ax0=None, xlim=[-1, 1],
return [strm]
def plot_function_map(f, nb_sub=10, ax0=None, xlim=[-1, 1], ylim=[-1, 1], opp=False):
"""
:param f: [function]
A function to plot that can take an array((N, nb_dim)) as input
:param nb_sub:
:param ax0:
:param xlim:
:param ylim:
:return:
"""
x = np.linspace(*xlim)
y = np.linspace(*ylim)
xx, yy = np.meshgrid(x, y)
# Y, X = np.mgrid[ylim[0]:ylim[1]:complex(nb_sub), xlim[0]:xlim[1]:complex(nb_sub)]
mesh_data = np.concatenate([np.atleast_2d(xx.ravel()), np.atleast_2d(yy.ravel())]).T
try:
zz = f(mesh_data)
except: # if function cannot take a vector as input
zz = np.array([f(_x) for _x in mesh_data])
z = zz.reshape(xx.shape)
CS = plt.contour(xx, yy, z, cmap='viridis')
plt.clabel(CS, inline=1, fontsize=10)
if opp: z = -z
plt.imshow(np.exp(z), interpolation='bilinear', origin='lower', extent=xlim + ylim,
alpha=0.5, cmap='viridis')
def plot_mixture_linear_system(model, mode='glob', nb_sub=20, gmm=True, min_alpha=0.,
cmap=plt.cm.jet, A=None,b=None, gmr=False, return_strm=False,
**kwargs):
......@@ -351,7 +381,7 @@ def plot_mixture_linear_system(model, mode='glob', nb_sub=20, gmm=True, min_alph
def plot_gmm(Mu, Sigma, dim=None, color=[1, 0, 0], alpha=0.5, linewidth=1, markersize=6,
ax=None, empty=False, edgecolor=None, edgealpha=None,
ax=None, empty=False, edgecolor=None, edgealpha=None, priors=None,
border=False, nb=1, swap=True, center=True):
''' This function displays the parameters of a Gaussian Mixture Model (GMM).
......@@ -365,7 +395,8 @@ def plot_gmm(Mu, Sigma, dim=None, color=[1, 0, 0], alpha=0.5, linewidth=1, marke
Note- Daniel Berio, switched matrix layout to be consistent with pbdlib matlab,
probably breaks with gmm now.
'''
Mu = np.array(Mu)
Sigma = np.array(Sigma)
if (Mu.ndim == 1):
if not swap:
Mu = Mu[:, np.newaxis]
......@@ -393,6 +424,9 @@ def plot_gmm(Mu, Sigma, dim=None, color=[1, 0, 0], alpha=0.5, linewidth=1, marke
sl = np.ix_(dim, dim)
Sigma = Sigma[sl]
if priors is not None:
priors /= np.max(priors)
priors = np.clip(priors, 0.1, 1.)
nbDrawingSeg = 35;
t = np.linspace(-np.pi, np.pi, nbDrawingSeg);
......@@ -423,6 +457,8 @@ def plot_gmm(Mu, Sigma, dim=None, color=[1, 0, 0], alpha=0.5, linewidth=1, marke
if edgecolor is None:
edgecolor = c
if priors is not None: a *= priors[i]
polygon = plt.Polygon(points.transpose().tolist(), facecolor=c, alpha=a,
linewidth=linewidth, zorder=20, edgecolor=edgecolor)
......@@ -487,6 +523,7 @@ def plot_gmm(Mu, Sigma, dim=None, color=[1, 0, 0], alpha=0.5, linewidth=1, marke
def plot_gaussian(mu, sigma, dim=None, color='r', alpha=0.5, lw=1, markersize=6,
ax=None, plots=None, nb_segm=24, **kwargs):
mu, sigma = np.array(mu), np.array(sigma)
t = np.linspace(-np.pi, np.pi, nb_segm)
R = np.real(sp.linalg.sqrtm(1.0 * sigma))
......@@ -507,6 +544,78 @@ def plot_gaussian(mu, sigma, dim=None, color='r', alpha=0.5, lw=1, markersize=6,
return center, line
def plot_y_gaussian(x, mu, sigma, dim=0, alpha=1., alpha_fill=None, color='r', lw=1.,
ax=None):
"""
:param mu: [n_states]
:param mu: [n_states, n_dim]
:param sigma: [n_states, n_dim, n_dim]
:param dim:
:return:
"""
if x.ndim == 2:
x = x[:, 0]
if alpha_fill is None:
alpha_fill = 0.4 * alpha
if ax is None:
ax = plt
ax.plot(x, mu[:, dim], alpha=alpha, color=color)
ax.fill_between(x,
mu[:, dim] - sigma[:, dim, dim] ** 0.5,
mu[:, dim] + sigma[:, dim, dim] ** 0.5,
alpha=alpha_fill, color=color)
def plot_dynamic_system(f, nb_sub=10, ax=None, xlim=[-1, 1], ylim=[-1, 1], scale=0.01,
name=None, equal=False, **kwargs):
"""
Plot a dynamical system dx = f(x)
:param f: a function that takes as input x as [N,2] and return dx [N, 2]
:param nb_sub:
:param ax0:
:param xlim:
:param ylim:
:param scale:
:param kwargs:
:return:
"""
Y, X = np.mgrid[ylim[0]:ylim[1]:complex(nb_sub), xlim[0]:xlim[1]:complex(nb_sub)]
mesh_data = np.vstack([X.ravel(), Y.ravel()])
field = f(mesh_data.T).T
U = field[0]
V = field[1]
U = U.reshape(nb_sub, nb_sub)
V = V.reshape(nb_sub, nb_sub)
speed = np.sqrt(U * U + V * V)
if name is not None:
plt.suptitle(name)
if ax is not None:
strm = ax.streamplot(X, Y, U, V, linewidth=scale* speed,**kwargs)
ax.set_xlim(xlim)
ax.set_ylim(ylim)
if equal:
ax.set_aspect('equal')
else:
strm = plt.streamplot(X, Y, U, V, linewidth=scale* speed, **kwargs)
plt.xlim(xlim)
plt.ylim(ylim)
if equal:
plt.axes().set_aspect('equal')
return [strm]
def plot_trajdist(td, ix=0, iy=1, covScale=1, color=[1, 0, 0], alpha=0.1, linewidth=0.1):
'''Plot 2D representation of a trajectory distribution'''
......
from setuptools import setup, find_packages
import os
# Pybdlib can be installed with demonstration data. The demonstration data should be
# located in the following folder:
data_path = '%s/pbdlib/data/'%(os.curdir)
# The script below will search this folder recursively for files that have the
# following extensions:
data_ext = ['mat']
# The folder structure is copied when installing the data with the pybdlib module.
# Gather data files:
def gatherfiles(path,ext_list,fname):
dlist = []
# gather files located direclty in path:
print('path: ', path)
file_names = ['%s%s'%(path,fn) for fn in os.listdir((path))
if any(fn.endswith(ext) for ext in ext_list)]
if len(file_names)>0:
dlist.append((fname,file_names))
# collect files located in subdirectories of path:
for fn in os.listdir(path):
if os.path.isdir('%s%s'%(path,fn)):
# select first element in the list
tmp = gatherfiles('%s%s/'%(path,fn),ext_list,'%s/%s'%(fname,fn))
if len(tmp) >0:
dlist.append(tmp[0])
return dlist
dlist = gatherfiles(data_path,data_ext,'pybdlib/data')
# Setup:
setup(name='pbdlib',
......@@ -9,5 +40,6 @@ setup(name='pbdlib',
author_email='emmanuel.pignat@idiap.ch',
license='MIT',
packages=find_packages(),
install_requires = ['ipython==5','numpy','scipy','matplotlib','sklearn', 'dtw', 'jupyter', 'enum', 'termcolor'],
data_files = dlist,
install_requires = ['numpy','scipy','matplotlib','sklearn', 'dtw', 'jupyter', 'enum', 'termcolor'],
zip_safe=False)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment