Commit 894d5307 authored by Sylvain CALINON's avatar Sylvain CALINON

PbDlib v1.0

parent dae71d94
......@@ -155,4 +155,4 @@ view(180,-90);
%print('-dpng','-r600','graphs/demo_DMP01.png');
pause;
close all;
close all;
\ No newline at end of file
......@@ -35,7 +35,7 @@ addpath('./m_fcts/');
%% Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model.nbStates = 5; %Number of activation functions (i.e., number of states in the GMM)
model.nbStates = 8; %Number of activation functions (i.e., number of states in the GMM)
model.nbVar = 1; %Number of variables for the radial basis functions [s] (decay term)
model.nbVarPos = 2; %Number of motion variables [x1,x2]
model.kP = 50; %Stiffness gain
......@@ -80,7 +80,7 @@ model = init_GMM_timeBased(sIn, model);
%Set Sigma
for i=1:model.nbStates
model.Sigma(:,:,i) = 1E-2; %Heuristic setting of covariance
model.Sigma(:,:,i) = 1E-3; %Heuristic setting of covariance
end
%Compute activation
......@@ -243,4 +243,4 @@ view(180,-90);
%print('-dpng','-r100','graphs/demo_DMP02.png');
pause;
close all;
close all;
\ No newline at end of file
......@@ -141,4 +141,4 @@ view(180,-90);
%print('-dpng','graphs/demo_DMP_GMR01.png');
pause;
close all;
close all;
\ No newline at end of file
function demo_DMP_GMR02
% Enhanced dynamic movement primitive (DMP) model trained with EM by using a Gaussian mixture
% model (GMM) representation, with full covariance matrices coordinating the different variables
% in the feature space. After learning (i.e., autonomous organization of the basis functions (position
% and spread), Gaussian mixture regression (GMR) is used to regenerate the nonlinear force profile.
% in the feature space. After learning (i.e., autonomous organization of the basis functions
% (position and spread), Gaussian mixture regression (GMR) is used to regenerate the nonlinear force profile.
%
% If this code is useful for your research, please cite the related publication:
% @incollection{Calinon19chapter,
......@@ -139,4 +139,4 @@ view(180,-90);
%print('-dpng','graphs/demo_DMP_GMR02.png');
pause;
close all;
close all;
\ No newline at end of file
......@@ -139,4 +139,4 @@ view(180,-90);
%print('-dpng','graphs/demo_DMP_GMR03.png');
pause;
close all;
close all;
\ No newline at end of file
......@@ -154,4 +154,4 @@ view(180,-90);
%print('-dpng','graphs/demo_DMP_GMR04.png');
pause;
close all;
close all;
\ No newline at end of file
......@@ -54,8 +54,8 @@ nbData = 200; %Number of datapoints in a trajectory
nbSamples = 4; %Number of demonstrations
%Canonical system parameters
A = kron([0 1; 0 0], eye(model.nbVarPos)); %See Eq. (5.1.1) in doc/TechnicalReport.pdf
B = kron([0; 1], eye(model.nbVarPos)); %See Eq. (5.1.1) in doc/TechnicalReport.pdf
A = kron([0 1; 0 0], eye(model.nbVarPos));
B = kron([0; 1], eye(model.nbVarPos));
C = kron([1,0],eye(model.nbVarPos));
%Discretize system (Euler method)
Ad = A*model.dt + eye(size(A));
......@@ -64,7 +64,7 @@ Bd = B*model.dt;
R = eye(model.nbVar) * model.rFactor;
R = kron(eye(nbData),R);
%Create transformation matrix to compute xhat = x + dx*kV/kP + ddx/kP, see Eq. (4.0.2) in doc/TechnicalReport.pdf
%Create transformation matrix to compute xhat = x + dx*kV/kP + ddx/kP
K1d = [1, model.kV/model.kP, 1/model.kP];
K = kron(K1d,eye(model.nbVarPos));
......@@ -103,10 +103,10 @@ model = EM_tensorGMM(DataDMP, model);
%Task-adaptive spring-damper attractor path retrieval
r(n).p(1).A = eye(model.nbVar);
r(n).p(1).b = s(n).p(1).b + [0; 0; 5]; %Offset added to test generalization capability
[r(n).Mu, r(n).Sigma] = productTPGMM0(model, r(n).p); %See Eq. (6.0.5), (6.0.6) and (6.0.7) in doc/TechnicalReport.pdf
[r(n).Mu, r(n).Sigma] = productTPGMM0(model, r(n).p);
r(n).Priors = model.Priors;
r(n).nbStates = model.nbStates;
[r(1).currTar, r(1).currSigma] = GMR(r(n), sIn, 1, 2:model.nbVar); %See Eq. (3.0.2) to (3.0.5) in doc/TechnicalReport.pdf
[r(1).currTar, r(1).currSigma] = GMR(r(n), sIn, 1, 2:model.nbVar);
%LQR tracking (discrete version)
Q = zeros(model.nbVarPos*2, model.nbVarPos*2);
......@@ -175,4 +175,4 @@ view(180,-90);
%print('-dpng','graphs/demo_DMP_GMR_LQR01.png');
pause;
close all;
close all;
\ No newline at end of file
......@@ -93,7 +93,7 @@ end
%% Plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure('PaperPosition',[0 0 10 4.0],'position',[10,10,1300,450],'color',[1 1 1]);
figure('PaperPosition',[0 0 10 4.0],'position',[10,10,2300,900],'color',[1 1 1]);
clrmap = lines(model.nbStates);
%Activation of the basis functions
......@@ -149,7 +149,7 @@ set(gca,'xtick',[],'ytick',[sIn(end),sIn(1)],'yticklabel',{'0','1'});
xlabel('$t$','fontsize',16,'interpreter','latex');
ylabel('$s$','fontsize',16,'interpreter','latex');
%print('-dpng','-r300','graphs/DMP_GMR01.png');
%print('-dpng','graphs/DMP_GMR01.png');
pause;
close all;
end
......@@ -184,4 +184,4 @@ function hg = displaySpring(x,colTmp)
%h = [h draw2DArrow(msh(:,25)+[0.5;-0.5], [1.2;1.1], [.7 0 0])];
hg = [hg plot(x(1,1),x(2,1),'.','markersize',15,'color',colTmp)];
hg = [hg plot(x(1,2),x(2,2),'.','markersize',15,'color',colTmp)];
end
end
\ No newline at end of file
......@@ -70,4 +70,4 @@ axis equal; set(gca,'Xtick',[]); set(gca,'Ytick',[]);
%print('-dpng','graphs/demo_GMM01.png');
%pause;
%close all;
%close all;
\ No newline at end of file
......@@ -80,6 +80,6 @@ for i=1:2
xlabel(['x_' num2str((i-1)*2+1)]); ylabel(['x_' num2str(i*2)]);
end
%print('-dpng','graphs/demo_MFA01.png');
%pause;
%close all;
%print('-dpng','graphs/demo_GMM_MFA01.png');
pause;
close all;
\ No newline at end of file
......@@ -78,6 +78,6 @@ for i=1:2
xlabel(['x_' num2str((i-1)*2+1)]); ylabel(['x_' num2str(i*2)]);
end
%print('-dpng','graphs/demo_MPPCA01.png');
%pause;
%close all;
%print('-dpng','graphs/demo_GMM_MPPCA01.png');
pause;
close all;
\ No newline at end of file
function demo_GMM_logGMM01
% Multivariate lognormal mixture model parameters estimation with EM algorithm
% Multivariate lognormal mixture model parameters estimation with EM algorithm.
%
% If this code is useful for your research, please cite the related publication:
% @article{Calinon16JIST,
......@@ -31,7 +31,6 @@ function demo_GMM_logGMM01
% You should have received a copy of the GNU General Public License
% along with PbDlib. If not, see <http://www.gnu.org/licenses/>.
addpath('./m_fcts/');
......@@ -42,9 +41,6 @@ model.nbVar = 2; %Number of variables [x1,x2]
nbData = 500; %Number of datapoints
figure('PaperPosition',[0 0 24 4],'position',[10,10,1000,600]); hold on;
% for nb=1:20
%% Generate data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model.Mu(:,1) = [1; 0.5];
......@@ -69,12 +65,8 @@ end
%% Parameters estimation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model.Mu
% model.Sigma
model = init_GMM_kmeans(Data, model);
model = EM_logGMM(Data, model);
% model.Mu
% model.Sigma
nbX = 100;
[X0,Y0] = meshgrid(linspace(min(Data(1,:)),max(Data(1,:)),nbX), linspace(min(Data(2,:)),max(Data(2,:)),nbX));
......@@ -91,17 +83,11 @@ hmix = sum(h,1);
%% Plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf; hold on; axis off;
for i=1:model.nbStates
%plot(DataIn, h(i,:), '-','color',[.7 .7 .7]);
end
figure('PaperPosition',[0 0 24 4],'position',[10,10,1600,900]); hold on; hold on; axis off;
pcolor(X0,Y0,reshape(hmix,nbX,nbX)); shading flat;
plot(Data(1,:), Data(2,:), '.','markersize',6,'color',[0 0 0]);
axis([min(Data(1,:)), max(Data(1,:)), min(Data(2,:)), max(Data(2,:))]);
% pause
% end %nb
%print('-dpng','graphs/demo_logGMM01.png');
pause;
close all;
close all;
\ No newline at end of file
function demo_GMM_logNormal01
%Conditional probability with multivariate log-normal distribution
% Conditional probability with multivariate log-normal distribution
%
% If this code is useful for your research, please cite the related publication:
% @article{Calinon16JIST,
......@@ -42,21 +42,15 @@ model0.nbStates = 1;
model0.Priors = 1;
model = model0;
figure('position',[10 10 1300 500],'color',[1 1 1]);
% disp('Press enter to regenerate data...');
% for nb=1:5
%% Multivariate normal distribution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model0.Mu = zeros(model0.nbVar,1);
%model0.Mu = rand(model0.nbVar,1)*10;
%if nb<2
V0 = rand(model0.nbVar,1);
model0.Sigma = V0*V0' + diag(rand(model0.nbVar,1))*1E-1;
[V,D] = eigs(model0.Sigma);
%end
Data0 = V*D^.5 * randn(model0.nbVar,nbData) + repmat(model0.Mu,1,nbData);
......@@ -86,7 +80,8 @@ DataOut = exp(DataOut0);
%% Plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf;
figure('position',[10 10 2300 1200],'color',[1 1 1]);
%Normal distribution
subplot(1,2,1); hold on; grid on; title('Normal distribution');
plot(Data0(1,:), Data0(2,:), '.','markersize',8,'color',[.7 .7 .7]);
......@@ -102,10 +97,8 @@ set(gca,'xtick',model0.Mu(1),'ytick',model0.Mu(2));
% X = R * [cos(t); sin(t)] + repmat(model0.Mu, 1, nbDrawingSeg);
% gaussPDF(X, model0.Mu, model0.Sigma)
%Log-normal distribution
subplot(1,2,2); hold on; grid on; title('Log-normal distribution');
nbX = 50;
[X0,Y0] = meshgrid(linspace(min(Data(1,:)),max(Data(1,:)),nbX), linspace(min(Data(2,:)),max(Data(2,:)),nbX));
X = reshape(X0,1,nbX^2);
......@@ -129,9 +122,7 @@ plot(model.Mu(1,:), model.Mu(2,:), '.', 'markersize', 6, 'color', [.4 0 0]);
axis equal; axis([min(Data(1,:)), max(Data(1,:)), min(Data(2,:)), max(Data(2,:))]);
xlabel('exp(x_1)'); ylabel('exp(x_2)');
set(gca,'xtick',exp(model0.Mu(1)),'ytick',exp(model0.Mu(2)));
%print('-dpng','graphs/demo_GMM_logNormal01.png');
% pause;
% end %nb
%print('-dpng','graphs/demo_logNormal01.png');
close all;
pause;
close all;
\ No newline at end of file
function demo_GMM_profileGMM01
% Example of univariate velocity profile fitting with a Gaussian mixture model (GMM) and a weighted EM algorithm.
% Univariate velocity profile fitting with a Gaussian mixture model (GMM) and a weighted EM algorithm.
% The approach shares links with radial basis function networks (RBFN), but adapts the Gaussian basis functions with EM
% based on the distribution profiles instead of considering a simple clustering of the input space.
%
......@@ -41,35 +41,19 @@ addpath('./m_fcts/');
model.nbStates = 40; %Number of states in the GMM
model.nbVar = 1; %Number of variables [x1]
nbData = 1000; %Length of each trajectory
nbSamples = 1; %Number of samples
%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% demos=[];
% load('data/2Dletters/W.mat');
% Data=[]; w=[];
% for n=1:nbSamples
% s(n).w = spline(1:size(demos{n}.vel(1:model.nbVar,:),2), demos{n}.vel(1:model.nbVar,:), linspace(1,size(demos{n}.pos,2),nbData)); %Resampling
% Data = [Data, 1:nbData];
% w = [w, s(n).w];
% end
% w = w-min(w);
% w = w/max(w);
% Data = 1:nbData;
% load('data/lognormal03.mat');
load('data/velprofile01.mat');
Data = spline(1:size(Data,2),Data,linspace(1,size(Data,2),nbData)); %Resample data
tlist = Data(1,:) + 1E-8; %The first value should be non-zero
w = Data(2,:);
Data = tlist;
%Rescale data and make it positive to represent a probability density function
w = w-min(w);
w = w/max(w);
w = w - min(w);
w = w / max(w);
%% Parameters estimation
......@@ -81,8 +65,7 @@ model = EM_weighted_univariateGMM(Data, w, model);
for i=1:model.nbStates
Phi(:,i) = gaussPDF(Data, model.Mu(:,i), model.Sigma(:,:,i));
end
model.Priors = (Phi'*Phi)\Phi'*w';
%sum(model.Priors)
model.Priors = (Phi' * Phi) \ Phi' * w';
%Probability density function of normal distributions
for i=1:model.nbStates
......@@ -90,28 +73,23 @@ for i=1:model.nbStates
end
hmix = sum(h,1);
disp(['Error: ' num2str(norm(hmix-w))]);
%% Plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure('PaperPosition',[0 0 24 4],'position',[10,10,1200,200]); hold on; box on;
figure('position',[10,10,2300,500]); hold on; box on;
%[~,hmix] = plotGMM1D(model, [.8 0 0], [tlist(1) 0 tlist(end) 1], .8, nbData);
hf(1) = plot(tlist, w, '-','linewidth',2,'color',[.7 .7 .7]);
for i=1:model.nbStates
plot(tlist, h(i,:),'-','linewidth',4,'color',[1 .7 .7]);
plot(tlist, h(i,:),'-','linewidth',2,'color',[1 .7 .7]);
end
for n=1:nbSamples
plot(tlist, w((n-1)*nbData+1:n*nbData), '-','linewidth',6,'color',[.7 .7 .7]);
end
plot(tlist, hmix,'-','linewidth',2,'color',[.8 0 0]);
hf(2) = plot(tlist, hmix,'-','linewidth',2,'color',[.8 0 0]);
axis([tlist(1) tlist(end) 0 1.05]); set(gca,'Xtick',[]); set(gca,'Ytick',[]);
xlabel('$t$','fontsize',18,'interpreter','latex');
ylabel('$|\dot{x}|$','fontsize',18,'interpreter','latex');
legend(hf, {'Reference','Reconstructed'});
% print('-dpng','graphs/demo_GMM_profileGMM01.png');
% w=hmix;
% save('data/lognormal03.mat','w');
% return
disp(['Error: ' num2str(norm(hmix-w))]);
% print('-dpng','graphs/demo_profileGMM01.png');
pause;
close all;
close all;
\ No newline at end of file
function demo_GMM_profileGMM_multivariate01
% Example of multivariate velocity profile fitting with a Gaussian mixture model (GMM) and a weighted EM algorithm
% Multivariate velocity profile fitting with a Gaussian mixture model (GMM) and a weighted EM algorithm
%
% If this code is useful for your research, please cite the related publication:
% @article{Calinon16JIST,
......@@ -36,25 +36,13 @@ addpath('./m_fcts/');
%% Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model.nbStates = 4; %Number of states in the GMM
model.nbStates = 8; %Number of states in the GMM
model.nbVar = 2; %Number of variables [x1]
nbData = 200; %Length of each trajectory
nbSamples = 1; %Number of samples
%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% demos=[];
% load('data/2Dletters/W.mat');
% Data=[]; w=[];
% for n=1:nbSamples
% s(n).w = spline(1:size(demos{n}.vel(1:model.nbVar,:),2), demos{n}.vel(1:model.nbVar,:), linspace(1,size(demos{n}.pos,2),nbData)); %Resampling
% Data = [Data, repmat(1:nbData,model.nbVar,1)];
% w = [w, s(n).w];
% end
% w = w - repmat(min(w,[],2),1,nbData*nbSamples);
% w = w ./ repmat(max(w,[],2),1,nbData*nbSamples);
Data = repmat(1:nbData,model.nbVar,1);
load('data/lognormal05.mat'); %load 'w'
......@@ -68,24 +56,24 @@ model = EM_weighted_multivariateGMM(Data, w, model);
%% Plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure('position',[10,10,1000,500]);
figure('position',[10,10,2000,900]);
for k=1:model.nbVar
subplot(model.nbVar,1,k); hold on; box on;
mtmp.nbStates = model.nbStates;
mtmp.Priors = model.Priors(k,:);
mtmp.Mu(1,:) = model.Mu(k,:);
mtmp.Sigma(1,1,:) = model.Sigma(k,k,:);
hf(1) = plot(1:nbData, w(k,:), '-','linewidth',2,'color',[.7 .7 .7]);
[~,hmix(k,:)] = plotGMM1D(mtmp, [1 nbData 0 1], [.8 0 0], .2);
for n=1:nbSamples
plot(1:nbData, w(k,(n-1)*nbData+1:n*nbData), '-','linewidth',2,'color',[.7 .7 .7]);
end
hf(2) = plot(1:nbData, hmix(k,:),'-','linewidth',2,'color',[.8 0 0]);
axis([1 nbData 0 1.05]); set(gca,'Xtick',[]); set(gca,'Ytick',[]);
xlabel('$t$','fontsize',18,'interpreter','latex');
ylabel(['$\dot{x}_' num2str(k) '$'],'fontsize',18,'interpreter','latex');
end
legend(hf, {'Reference','Reconstructed'});
%print('-dpng','graphs/demo_GMM_profileGMM_multivariate01.png');
norm(hmix-w)
disp(['Error: ' num2str(norm(hmix-w))]);
%print('-dpng','graphs/demo_profileGMM_multivariate01.png');
pause;
close all;
close all;
\ No newline at end of file
function demo_GMM_profileLogGMM01
% Example of univariate velocity profile fitting with a lognormal mixture model and a weighted EM algorithm
% Univariate velocity profile fitting with a lognormal mixture model and a weighted EM algorithm
%
% If this code is useful for your research, please cite the related publication:
% @article{Calinon16JIST,
......@@ -39,40 +39,19 @@ addpath('./m_fcts/');
model.nbStates = 40; %Number of states in the GMM
model.nbVar = 1; %Number of variables [x1]
nbData = 1000; %Length of each trajectory
nbSamples = 1; %Number of samples
%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% demos=[];
% load('data/2Dletters/W.mat');
% Data=[]; w=[];
% for n=1:nbSamples
% s(n).w = spline(1:size(demos{n}.vel(1:model.nbVar,:),2), demos{n}.vel(1:model.nbVar,:), linspace(1,size(demos{n}.pos,2),nbData)); %Resampling
% Data = [Data, log(1:nbData)]; %log data are collected
% w = [w, s(n).w];
% end
% w = w-min(w);
% w = w/max(w);
% Data = log(1:nbData);
% load('data/lognormal03.mat');
% s = importdata('data/2cold_vel_orig.csv');
% Data = s.data';
% save('data/velprofile01.mat','Data');
% return
load('data/velprofile01.mat');
Data = spline(1:size(Data,2),Data,linspace(1,size(Data,2),nbData)); %Resample data
tlist = Data(1,:) + 1E-8; %The first value should be non-zero
w = Data(2,:);
%Data = tlist;
Data = log(tlist);
%Rescale data and make it positive to represent a probability density function
w = w-min(w);
w = w/max(w);
w = w - min(w);
w = w / max(w);
%% Parameters estimation
......@@ -96,27 +75,21 @@ for i=1:model.nbStates
end
hmix = sum(h,1);
disp(['Error: ' num2str(norm(hmix-w))]);
%% Plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure('PaperPosition',[0 0 24 4],'position',[10,10,1200,200]); hold on; box on;
figure('position',[10,10,2300,500]); hold on; box on;
hf(1) = plot(tlist, w, '-','linewidth',2,'color',[.7 .7 .7]);
for i=1:model.nbStates
plot(tlist, h(i,:),'-','linewidth',4,'color',[1 .7 .7]);
plot(tlist, h(i,:),'-','linewidth',2,'color',[1 .7 .7]);
end
for n=1:nbSamples
plot(tlist,w((n-1)*nbData+1:n*nbData), '-','linewidth',6,'color',[.7 .7 .7]);
end
plot(tlist, hmix,'-','linewidth',2,'color',[.8 0 0]);
hf(2) = plot(tlist, hmix,'-','linewidth',2,'color',[.8 0 0]);
axis([tlist(1) tlist(end) 0 1.05]); set(gca,'Xtick',[]); set(gca,'Ytick',[]);
xlabel('$t$','fontsize',18,'interpreter','latex');
ylabel('$|\dot{x}|$','fontsize',18,'interpreter','latex');
legend(hf, {'Reference','Reconstructed'});
%print('-dpng','graphs/demo_GMM_profileLogGMM01.png');
% w=hmix;
% save('data/lognormal04.mat','w');
% return
disp(['Error: ' num2str(norm(hmix-w))]);
%print('-dpng','graphs/demo_profileLogGMM01.png');
pause;
close all;
close all;
\ No newline at end of file
function demo_GMM_profileLogGMM_multivariate01
%Example of multivariate velocity profile fitting with a Gaussian mixture model (GMM) and a weighted EM algorithm
% Multivariate velocity profile fitting with a Gaussian mixture model (GMM) and a weighted EM algorithm
%
% If this code is useful for your research, please cite the related publication:
% @article{Calinon16JIST,
......@@ -36,10 +36,9 @@ addpath('./m_fcts/');
%% Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model.nbStates = 4; %Number of states in the GMM
model.nbStates = 8; %Number of states in the GMM
model.nbVar = 2; %Number of variables [x1]
nbData = 200; %Length of each trajectory
nbSamples = 1; %Number of samples
%% Load data
......@@ -78,28 +77,26 @@ hmix = hmix ./ repmat(max(hmix,[],2), [1,nbData]);
%% Plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure('position',[10,10,1000,500]);
figure('position',[10,10,2000,900]);
for k=1:model.nbVar
subplot(model.nbVar,1,k); hold on; box on;
mtmp.nbStates = model.nbStates;
mtmp.Priors = model.Priors(k,:);
mtmp.Mu(1,:) = model.Mu(k,:);
mtmp.Sigma(1,1,:) = model.Sigma(k,k,:);
for n=1:nbSamples
plot(1:nbData, w(k,(n-1)*nbData+1:n*nbData), '-','linewidth',2,'color',[.7 .7 .7]);
end
hf(1) = plot(1:nbData, w(k,:), '-','linewidth',2,'color',[.7 .7 .7]);
for i=1:model.nbStates
plot(1:nbData, squeeze(h(k,i,:)), '-','linewidth',4,'color',[1 .7 .7]);
plot(1:nbData, squeeze(h(k,i,:)), '-','linewidth',2,'color',[1 .7 .7]);
end
plot(1:nbData, hmix(k,:), '-','linewidth',2,'color',[.8 0 0]);
hf(2) = plot(1:nbData, hmix(k,:),'-','linewidth',2,'color',[.8 0 0]);
axis([1 nbData 0 1.05]); set(gca,'Xtick',[]); set(gca,'Ytick',[]);
xlabel('$t$','fontsize',18,'interpreter','latex');
ylabel(['$\dot{x}_' num2str(k) '$'],'fontsize',18,'interpreter','latex');
end
legend(hf, {'Reference','Reconstructed'});
%print('-dpng','graphs/demo_GMM_profileLogGMM_multivariate01.png');
norm(hmix-w)
disp(['Error: ' num2str(norm(hmix-w))]);
%print('-dpng','graphs/demo_profileLogGMM_multivariate01.png');
pause;
close all;
\ No newline at end of file
......@@ -72,6 +72,6 @@ for i=1:model.nbStates
end
view(-40,6); axis equal;
%print('-dpng','graphs/demo_semitiedGMM01.png');
%print('-dpng','graphs/demo_GMM_semitiedGMM01.png');
pause;
close all;
\ No newline at end of file
......@@ -102,7 +102,13 @@ end
set(gca,'xtick',[],'ytick',[]); axis equal; axis square;
xlabel(['$x_1$'],'interpreter','latex','fontsize',18);
ylabel(['$x_2$'],'interpreter','latex','fontsize',18);
%print('-dpng','graphs/demo_GPR01.png');
% figure; hold on; axis off;
% colormap(flipud(gray));
% imagesc(abs(Kdd));
% axis tight; axis square; axis ij;
% print('-dpng','graphs/GPR_kernel_RBF01.png');
pause;
close all;
close all;
\ No newline at end of file
......@@ -161,7 +161,4 @@ function K = covFct(x1, x2, p, flag_noiseObs)
if flag_noiseObs==1
K = K + eye(size(K)) .* p(3); %Consideration of noisy observation y
end
end
end
\ No newline at end of file
......@@ -83,9 +83,9 @@ end
%% Plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure('PaperPosition',[0 0 12 4],'position',[10 10 1300 600]);
figure('PaperPosition',[0 0 22 4],'position',[10 10 2300 600]);
%Prior samples
subplot(1,3,1); hold on; title('Samples from prior','fontsize',14);
subplot(1,4,1); hold on; title('Samples from prior','fontsize',14);
for n=2:nbRepros/2
plot(r(n).Data(1,:), r(n).Data(2,:), '-','lineWidth',3.5,'color',[.9 .9 .9]*rand(1));
end
......@@ -96,7 +96,7 @@ xlabel('$x^{\scriptscriptstyle\mathcal{I}}_1$','interpreter','latex','fontsize',
ylabel('$x^{\scriptscriptstyle\mathcal{O}}_1$','interpreter','latex','fontsize',18);
%Posterior samples
subplot(1,3,2); hold on; title('Samples from posterior','fontsize',14);
subplot(1,4,2); hold on; title('Samples from posterior','fontsize',14);
for n=nbRepros/2+1:nbRepros
plot(r(n).Data(1,:), r(n).Data(2,:), '-','lineWidth',3.5,'color',[.9 .9 .9]*rand(1));
end
......@@ -108,7 +108,7 @@ xlabel('$x^{\scriptscriptstyle\mathcal{I}}_1$','interpreter','latex','fontsize',
ylabel('$x^{\scriptscriptstyle\mathcal{O}}_1$','interpreter','latex','fontsize',18);
%Trajectory distribution
subplot(1,3,3); hold on; title('Trajectory distribution','fontsize',14);
subplot(1,4,3); hold on; title('Trajectory distribution','fontsize',14);
patch([r(1).Data(1,:), r(1).Data(1,end:-1:1)], ...
[r(1).Data(2,:)+squeeze(r(1).SigmaOut.^.5)', r(1).Data(2,end:-1:1)-squeeze(r(1).SigmaOut(:,:,end:-1:1).^.5)'], ...
[.8 .8 .8],'edgecolor','none');
......@@ -120,6 +120,20 @@ set(gca,'xtick',[],'ytick',[]); axis([0, 1, -0.7 0.7]);
xlabel('$x^{\scriptscriptstyle\mathcal{I}}_1$','interpreter','latex','fontsize',18);
ylabel('$x^{\scriptscriptstyle\mathcal{O}}_1$','interpreter','latex','fontsize',18);
%Plot covariance
subplot(1,4,4); hold on; axis off; title('Covariance','fontsize',14);
colormap(flipud(gray));
imagesc(abs(Kdd));
axis tight; axis square; axis ij;
%print('-dpng','-r300','graphs/GPR03.png');
% figure; hold on; axis off;
% colormap(flipud(gray));
% imagesc(abs(Kdd));
% axis tight; axis square; axis ij;
% print('-dpng','graphs/GPR_kernel_periodic01.png');
pause;
close all;
......@@ -82,10 +82,10 @@ end
%% Plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure('PaperPosition',[0 0 12 4],'position',[10 10 1300 600]);
figure('PaperPosition',[0 0 22 4],'position',[10 10 2300 600]);
limAxes = [0, 1, -.5 .5];
%Prior samples
subplot(1,3,1); hold on; title('Samples from prior','fontsize',14);
subplot(1,4,1); hold on; title('Samples from prior','fontsize',14);
for n=2:nbRepros/2
plot(r(n).Data(1,:), r(n).Data(2,:), '-','lineWidth',3.5,'color',[.9 .9 .9]*rand(1));
end
......@@ -96,7 +96,7 @@ xlabel('$x^{\scriptscriptstyle\mathcal{I}}_1$','interpreter','latex','fontsize',
ylabel('$x^{\scriptscriptstyle\mathcal{O}}_1$','interpreter','latex','fontsize',18);
%Posterior samples
subplot(1,3,2); hold on; title('Samples from posterior','fontsize',14);
subplot(1,4,2); hold on; title('Samples from posterior','fontsize',14);
for n=nbRepros/2+1:nbRepros
plot(r(n).Data(1,:), r(n).Data(2,:), '-','lineWidth',3.5,'color',[.9 .9 .9]*rand(1));
end
......@@ -108,7 +108,7 @@ xlabel('$x^{\scriptscriptstyle\mathcal{I}}_1$','interpreter','latex','fontsize',
ylabel('$x^{\scriptscriptstyle\mathcal{O}}_1$','interpreter','latex','fontsize',18);
%Trajectory distribution
subplot(1,3,3); hold on; title('Trajectory distribution','fontsize',14);
subplot(1,4,3); hold on; title('Trajectory distribution','fontsize',14);
patch([r(1).Data(1,:), r(1).Data(1,end:-1:1)], ... <