diff --git a/README.md b/README.md
index e83577c903e98b9762e9f07ec361ff2f887005d5..300942dad00d0f5af4e2d520806555f8d60c122a 100644
--- a/README.md
+++ b/README.md
@@ -233,12 +233,15 @@ All the examples are located in the `demos` folder, and the functions are locate
 | [demo_search01.m](./demos/demo_search01.m) | [[2]](#ref-2) | EM-based stochastic optimization |
 | [demo_spring01.m](./demos/demo_spring01.m) | [[10]](#ref-10) | Influence of the damping ratio in mass-spring-damper systems |
 | [demo_stdPGMM01.m](./demos/demo_stdPGMM01.m) | [[1]](#ref-1) | Parametric Gaussian mixture model (PGMM) used as a task-parameterized model, with DS-GMR employed to retrieve continuous movements |
+| [demo_TPHDDC01.m](./demos/demo_TPHDDC01.m) | [[1]](#ref-1) | Task-parameterized high dimensional data clustering (TP-HDDC) |
 | [demo_TPGMM01.m](./demos/demo_TPGMM01.m) | [[1]](#ref-1) | Task-parameterized Gaussian mixture model (TP-GMM) encoding |
 | [demo_TPGMM_bimanualReaching01.m](./demos/demo_TPGMM_bimanualReaching01.m) | [[1]](#ref-1) | Time-invariant task-parameterized GMM applied to a bimanual reaching task |
 | [demo_TPGMM_teleoperation01.m](./demos/demo_TPGMM_teleoperation01.m) | [[1]](#ref-1) | Time-invariant task-parameterized GMM applied to a teleoperation task (position only) |
 | [demo_TPGMR01.m](./demos/demo_TPGMR01.m) | [[1]](#ref-1) | Task-parameterized Gaussian mixture model (TP-GMM), with GMR used for reproduction (without dynamical system) |
 | [demo_TPGP01.m](./demos/demo_TPGP01.m) | [[1]](#ref-1) | Task-parameterized Gaussian process regression (TP-GPR) |
-| [demo_TPMPC01.m](./demos/demo_TPMPC01.m) | [[1]](#ref-1) | Task-parameterized probabilistic model encoding position data, with MPC (batch version of LQR) used to track the associated stepwise reference path |
+| [demo_TPLQT01.m](./demos/demo_TPLQT01.m) | [[1]](#ref-1) | Task-parameterized probabilistic model encoding position data, with LQT used to track the associated stepwise reference path |
+| [demo_TPMFA01.m](./demos/demo_TPMFA01.m) | [[1]](#ref-1) | Task-parameterized mixture of factor analyzers (TP-MFA) |
+| [demo_TPMPPCA01.m](./demos/demo_TPMPPCA01.m) | [[1]](#ref-1) | Task-parameterized mixture of probabilistic principal component analyzers (TP-MPPCA) |
 | [demo_TPproMP01.m](./demos/demo_TPproMP01.m) | [[1]](#ref-1) | Task-parameterized probabilistic movement primitives (TP-ProMP) |
 | [demo_TPtrajDistrib01.m](./demos/demo_TPtrajDistrib01.m) | [[1]](#ref-1) | Task-parameterized model with trajectory distribution and eigendecomposition |
 | [demo_TPtrajGMM01.m](./demos/demo_TPtrajGMM01.m) | [[1]](#ref-1) | Task-parameterized model with trajectory-GMM encoding |
diff --git a/demos/demo_OC_DDP_car01.m b/demos/demo_OC_DDP_car01.m
index d94bd55600a7babaed09fc19ab006ba49081af84..62ba642acb930ad21d1f35b385266c8b4ce462a2 100644
--- a/demos/demo_OC_DDP_car01.m
+++ b/demos/demo_OC_DDP_car01.m
@@ -27,7 +27,6 @@ function demo_OC_DDP_car01
 % You should have received a copy of the GNU General Public License
 % along with PbDlib. If not, see <https://www.gnu.org/licenses/>.
 
-
 addpath('./m_fcts/');
 
 
diff --git a/demos/demo_TPHDDC01.m b/demos/demo_TPHDDC01.m
new file mode 100644
index 0000000000000000000000000000000000000000..ee1a1a29475c2576d99ff24dd2b328d4bc2b496c
--- /dev/null
+++ b/demos/demo_TPHDDC01.m
@@ -0,0 +1,109 @@
+function demo_TPHDDC01
+% Task-parameterized high dimensional data clustering (TP-HDDC).
+%
+% If this code is useful for your research, please cite the related publication:
+% @article{Calinon16JIST,
+%   author="Calinon, S.",
+%   title="A Tutorial on Task-Parameterized Movement Learning and Retrieval",
+%   journal="Intelligent Service Robotics",
+%		publisher="Springer Berlin Heidelberg",
+%		doi="10.1007/s11370-015-0187-9",
+%		year="2016",
+%		volume="9",
+%		number="1",
+%		pages="1--29"
+% }
+% 
+% Copyright (c) 2019 Idiap Research Institute, http://idiap.ch/
+% Written by Sylvain Calinon, http://calinon.ch/
+% 
+% This file is part of PbDlib, http://www.idiap.ch/software/pbdlib/
+% 
+% PbDlib 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.
+% 
+% PbDlib 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 PbDlib. If not, see <http://www.gnu.org/licenses/>.
+
+addpath('./m_fcts/');
+
+
+%% Parameters
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+model.nbStates = 3; %Number of Gaussians in the GMM
+model.nbFrames = 2; %Number of candidate frames of reference
+model.nbVar = 2; %Dimension of the datapoints in the dataset (here: x1,x2)
+model.nbFA = 1; %Dimension of the subspace
+
+
+%% Load 3rd order tensor data
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+disp('Load 3rd order tensor data...');
+% The MAT file contains a structure 's' with the multiple demonstrations. 's(n).Data' is a matrix data for
+% sample n (with 's(n).nbData' datapoints). 's(n).p(m).b' and 's(n).p(m).A' contain the position and
+% orientation of the m-th candidate coordinate system for this demonstration. 'Data' contains the observations
+% in the different frames. It is a 3rd order tensor of dimension D x P x N, with D=2 the dimension of a
+% datapoint, P=2 the number of candidate frames, and N=TM the number of datapoints in a trajectory (T=200)
+% multiplied by the number of demonstrations (M=5).
+load('data/Data01.mat');
+
+
+%% TP-HDDC learning
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+fprintf('Parameters estimation of TP-HDDC with EM:');
+model = init_tensorGMM_kmeans(Data, model); 
+model = EM_tensorHDGMM(Data, model);
+
+%Reconstruct GMM for each demonstration
+for n=1:nbSamples
+	[s(n).Mu, s(n).Sigma] = productTPGMM0(model, s(n).p); %See Eq. (6)
+end
+
+
+%% Plots
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+figure('position',[20,50,1300,500]);
+xx = round(linspace(1,64,nbSamples));
+clrmap = colormap('jet');
+clrmap = min(clrmap(xx,:),.95);
+limAxes = [-1.2 0.8 -1.1 0.9];
+colPegs = [[.9,.5,.9];[.5,.9,.5]];
+
+%DEMOS
+subplot(1,model.nbFrames+1,1); hold on; box on; title('Demonstrations');
+for n=1:nbSamples
+	%Plot frames
+	for m=1:model.nbFrames
+		plot([s(n).p(m).b(1) s(n).p(m).b(1)+s(n).p(m).A(1,2)], [s(n).p(m).b(2) s(n).p(m).b(2)+s(n).p(m).A(2,2)], '-','linewidth',6,'color',colPegs(m,:));
+		plot(s(n).p(m).b(1), s(n).p(m).b(2),'.','markersize',30,'color',colPegs(m,:)-[.05,.05,.05]);
+	end
+	%Plot trajectories
+	plot(s(n).Data(2,1), s(n).Data(3,1),'.','markersize',15,'color',clrmap(n,:));
+	plot(s(n).Data(2,:), s(n).Data(3,:),'-','linewidth',1.5,'color',clrmap(n,:));
+	%Plot Gaussians
+	plotGMM(s(n).Mu, s(n).Sigma, [.5 .5 .5],.8);
+end
+axis(limAxes); axis square; set(gca,'xtick',[],'ytick',[]);
+
+%FRAMES
+for m=1:model.nbFrames
+	subplot(1,model.nbFrames+1,1+m); hold on; grid on; box on; title(['Frame ' num2str(m)]);
+	for n=1:nbSamples
+		plot(squeeze(Data(1,m,(n-1)*s(1).nbData+1)), ...
+			squeeze(Data(2,m,(n-1)*s(1).nbData+1)), '.','markersize',15,'color',clrmap(n,:));
+		plot(squeeze(Data(1,m,(n-1)*s(1).nbData+1:n*s(1).nbData)), ...
+			squeeze(Data(2,m,(n-1)*s(1).nbData+1:n*s(1).nbData)), '-','linewidth',1.5,'color',clrmap(n,:));
+	end
+	plotGMM(squeeze(model.Mu(:,m,:)), squeeze(model.Sigma(:,:,m,:)), [.5 .5 .5],.8);
+	axis square; set(gca,'xtick',[0],'ytick',[0]);
+end
+
+%print('-dpng','graphs/demo_TPHDDC01.png');
+pause;
+close all;
\ No newline at end of file
diff --git a/demos/demo_TPMPC01.m b/demos/demo_TPLQT01.m
similarity index 99%
rename from demos/demo_TPMPC01.m
rename to demos/demo_TPLQT01.m
index 0bb50a12752c45f3feb6d23d76e7989b4619b8ce..3cd83baf49489f3e7ce6ae4ec669a4cb1f485647 100644
--- a/demos/demo_TPMPC01.m
+++ b/demos/demo_TPLQT01.m
@@ -1,5 +1,5 @@
-function demo_TPMPC01
-% Linear quadratic control (unconstrained linear MPC) acting in multiple frames,
+function demo_TPLQT01
+% Linear quadratic control acting in multiple frames,
 % which is equivalent to a product of Gaussian controllers through a TP-GMM representation
 %
 % If this code is useful for your research, please cite the related publication:
diff --git a/demos/demo_TPMFA01.m b/demos/demo_TPMFA01.m
new file mode 100644
index 0000000000000000000000000000000000000000..c67b3b57c2812ded20061e732a1ad5d6ad43d2e3
--- /dev/null
+++ b/demos/demo_TPMFA01.m
@@ -0,0 +1,109 @@
+function demo_TPMFA01
+% Task-parameterized mixture of factor analyzers (TP-MFA).
+%
+% If this code is useful for your research, please cite the related publication:
+% @article{Calinon16JIST,
+%   author="Calinon, S.",
+%   title="A Tutorial on Task-Parameterized Movement Learning and Retrieval",
+%   journal="Intelligent Service Robotics",
+%		publisher="Springer Berlin Heidelberg",
+%		doi="10.1007/s11370-015-0187-9",
+%		year="2016",
+%		volume="9",
+%		number="1",
+%		pages="1--29"
+% }
+% 
+% Copyright (c) 2019 Idiap Research Institute, http://idiap.ch/
+% Written by Sylvain Calinon, http://calinon.ch/
+% 
+% This file is part of PbDlib, http://www.idiap.ch/software/pbdlib/
+% 
+% PbDlib 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.
+% 
+% PbDlib 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 PbDlib. If not, see <http://www.gnu.org/licenses/>.
+
+addpath('./m_fcts/');
+
+
+%% Parameters
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+model.nbStates = 3; %Number of Gaussians in the GMM
+model.nbFrames = 2; %Number of candidate frames of reference
+model.nbVar = 2; %Dimension of the datapoints in the dataset (here: x1,x2)
+model.nbFA = 1; %Dimension of the subspace (number of factor analyzers)
+
+
+%% Load 3rd order tensor data
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+disp('Load 3rd order tensor data...');
+% The MAT file contains a structure 's' with the multiple demonstrations. 's(n).Data' is a matrix data for
+% sample n (with 's(n).nbData' datapoints). 's(n).p(m).b' and 's(n).p(m).A' contain the position and
+% orientation of the m-th candidate coordinate system for this demonstration. 'Data' contains the observations
+% in the different frames. It is a 3rd order tensor of dimension D x P x N, with D=2 the dimension of a
+% datapoint, P=2 the number of candidate frames, and N=TM the number of datapoints in a trajectory (T=200)
+% multiplied by the number of demonstrations (M=5).
+load('data/Data01.mat');
+
+
+%% TP-MFA learning
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+fprintf('Parameters estimation of TP-MFA with EM:');
+model = init_tensorGMM_kmeans(Data, model); 
+model = EM_tensorMFA(Data, model);
+
+%Reconstruct GMM for each demonstration
+for n=1:nbSamples
+	[s(n).Mu, s(n).Sigma] = productTPGMM0(model, s(n).p); %See Eq. (6)
+end
+
+
+%% Plots
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+figure('position',[20,50,1300,500]);
+xx = round(linspace(1,64,nbSamples));
+clrmap = colormap('jet');
+clrmap = min(clrmap(xx,:),.95);
+limAxes = [-1.2 0.8 -1.1 0.9];
+colPegs = [[.9,.5,.9];[.5,.9,.5]];
+
+%DEMOS
+subplot(1,model.nbFrames+1,1); hold on; box on; title('Demonstrations');
+for n=1:nbSamples
+	%Plot frames
+	for m=1:model.nbFrames
+		plot([s(n).p(m).b(1) s(n).p(m).b(1)+s(n).p(m).A(1,2)], [s(n).p(m).b(2) s(n).p(m).b(2)+s(n).p(m).A(2,2)], '-','linewidth',6,'color',colPegs(m,:));
+		plot(s(n).p(m).b(1), s(n).p(m).b(2),'.','markersize',30,'color',colPegs(m,:)-[.05,.05,.05]);
+	end
+	%Plot trajectories
+	plot(s(n).Data(2,1), s(n).Data(3,1),'.','markersize',15,'color',clrmap(n,:));
+	plot(s(n).Data(2,:), s(n).Data(3,:),'-','linewidth',1.5,'color',clrmap(n,:));
+	%Plot Gaussians
+	plotGMM(s(n).Mu, s(n).Sigma, [.5 .5 .5],.8);
+end
+axis(limAxes); axis square; set(gca,'xtick',[],'ytick',[]);
+
+%FRAMES
+for m=1:model.nbFrames
+	subplot(1,model.nbFrames+1,1+m); hold on; grid on; box on; title(['Frame ' num2str(m)]);
+	for n=1:nbSamples
+		plot(squeeze(Data(1,m,(n-1)*s(1).nbData+1)), ...
+			squeeze(Data(2,m,(n-1)*s(1).nbData+1)), '.','markersize',15,'color',clrmap(n,:));
+		plot(squeeze(Data(1,m,(n-1)*s(1).nbData+1:n*s(1).nbData)), ...
+			squeeze(Data(2,m,(n-1)*s(1).nbData+1:n*s(1).nbData)), '-','linewidth',1.5,'color',clrmap(n,:));
+	end
+	plotGMM(squeeze(model.Mu(:,m,:)), squeeze(model.Sigma(:,:,m,:)), [.5 .5 .5],.8);
+	axis square; set(gca,'xtick',[0],'ytick',[0]);
+end
+
+%print('-dpng','graphs/demo_TPMFA01.png');
+pause;
+close all;
\ No newline at end of file
diff --git a/demos/demo_TPMPPCA01.m b/demos/demo_TPMPPCA01.m
new file mode 100644
index 0000000000000000000000000000000000000000..f14ec97ac5b0a92ae9699b0688b20d4c8685f5e9
--- /dev/null
+++ b/demos/demo_TPMPPCA01.m
@@ -0,0 +1,109 @@
+function demo_TPMPPCA01
+% Task-parameterized mixture of probabilistic principal component analyzers (TP-MPPCA). 
+%
+% If this code is useful for your research, please cite the related publication:
+% @article{Calinon16JIST,
+%   author="Calinon, S.",
+%   title="A Tutorial on Task-Parameterized Movement Learning and Retrieval",
+%   journal="Intelligent Service Robotics",
+%		publisher="Springer Berlin Heidelberg",
+%		doi="10.1007/s11370-015-0187-9",
+%		year="2016",
+%		volume="9",
+%		number="1",
+%		pages="1--29"
+% }
+% 
+% Copyright (c) 2019 Idiap Research Institute, http://idiap.ch/
+% Written by Sylvain Calinon, http://calinon.ch/
+% 
+% This file is part of PbDlib, http://www.idiap.ch/software/pbdlib/
+% 
+% PbDlib 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.
+% 
+% PbDlib 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 PbDlib. If not, see <http://www.gnu.org/licenses/>.
+
+addpath('./m_fcts/');
+
+
+%% Parameters
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+model.nbStates = 3; %Number of Gaussians in the GMM
+model.nbFrames = 2; %Number of candidate frames of reference
+model.nbVar = 2; %Dimension of the datapoints in the dataset (here: x1,x2)
+model.nbFA = 1; %Dimension of the subspace (number of principal components)
+
+
+%% Load 3rd order tensor data
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+disp('Load 3rd order tensor data...');
+% The MAT file contains a structure 's' with the multiple demonstrations. 's(n).Data' is a matrix data for
+% sample n (with 's(n).nbData' datapoints). 's(n).p(m).b' and 's(n).p(m).A' contain the position and
+% orientation of the m-th candidate coordinate system for this demonstration. 'Data' contains the observations
+% in the different frames. It is a 3rd order tensor of dimension D x P x N, with D=2 the dimension of a
+% datapoint, P=2 the number of candidate frames, and N=TM the number of datapoints in a trajectory (T=200)
+% multiplied by the number of demonstrations (M=5).
+load('data/Data01.mat');
+
+
+%% TP-MPPCA learning
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+fprintf('Parameters estimation of TP-MPPCA with EM:');
+model = init_tensorGMM_kmeans(Data, model); 
+model = EM_tensorMPPCA(Data, model);
+
+%Reconstruct GMM for each demonstration
+for n=1:nbSamples
+	[s(n).Mu, s(n).Sigma] = productTPGMM0(model, s(n).p); %See Eq. (6)
+end
+
+
+%% Plots
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+figure('position',[20,50,1300,500]);
+xx = round(linspace(1,64,nbSamples));
+clrmap = colormap('jet');
+clrmap = min(clrmap(xx,:),.95);
+limAxes = [-1.2 0.8 -1.1 0.9];
+colPegs = [[.9,.5,.9];[.5,.9,.5]];
+
+%DEMOS
+subplot(1,model.nbFrames+1,1); hold on; box on; title('Demonstrations');
+for n=1:nbSamples
+	%Plot frames
+	for m=1:model.nbFrames
+		plot([s(n).p(m).b(1) s(n).p(m).b(1)+s(n).p(m).A(1,2)], [s(n).p(m).b(2) s(n).p(m).b(2)+s(n).p(m).A(2,2)], '-','linewidth',6,'color',colPegs(m,:));
+		plot(s(n).p(m).b(1), s(n).p(m).b(2),'.','markersize',30,'color',colPegs(m,:)-[.05,.05,.05]);
+	end
+	%Plot trajectories
+	plot(s(n).Data(2,1), s(n).Data(3,1),'.','markersize',15,'color',clrmap(n,:));
+	plot(s(n).Data(2,:), s(n).Data(3,:),'-','linewidth',1.5,'color',clrmap(n,:));
+	%Plot Gaussians
+	plotGMM(s(n).Mu, s(n).Sigma, [.5 .5 .5],.8);
+end
+axis(limAxes); axis square; set(gca,'xtick',[],'ytick',[]);
+
+%FRAMES
+for m=1:model.nbFrames
+	subplot(1,model.nbFrames+1,1+m); hold on; grid on; box on; title(['Frame ' num2str(m)]);
+	for n=1:nbSamples
+		plot(squeeze(Data(1,m,(n-1)*s(1).nbData+1)), ...
+			squeeze(Data(2,m,(n-1)*s(1).nbData+1)), '.','markersize',15,'color',clrmap(n,:));
+		plot(squeeze(Data(1,m,(n-1)*s(1).nbData+1:n*s(1).nbData)), ...
+			squeeze(Data(2,m,(n-1)*s(1).nbData+1:n*s(1).nbData)), '-','linewidth',1.5,'color',clrmap(n,:));
+	end
+	plotGMM(squeeze(model.Mu(:,m,:)), squeeze(model.Sigma(:,:,m,:)), [.5 .5 .5],.8);
+	axis square; set(gca,'xtick',[0],'ytick',[0]);
+end
+
+%print('-dpng','graphs/demo_TPMPPCA01.png');
+pause;
+close all;
\ No newline at end of file