diff --git a/demos/demo_GPR_closedShape01.m b/demos/demo_GPR_closedShape01.m
index 60f5c294dd09712b6218cbcd21750cd2ca1d9705..e7e786d244f1a5800fc941098a912293a5993340 100644
--- a/demos/demo_GPR_closedShape01.m
+++ b/demos/demo_GPR_closedShape01.m
@@ -1,5 +1,5 @@
function demo_GPR_closedShape01
-% Closed shape modeling with sequential points and RBFs in Gaussian process regression (GPR)
+% Closed shape modeling with sequential points and periodic RBF kernel in Gaussian process regression (GPR)
%
% If this code is useful for your research, please cite the related publication:
% @incollection{Calinon19chapter,
@@ -170,7 +170,7 @@ close all;
end
% %User-defined distance function
-% function d = distfun(x1,x2)
+% function d = distfun(x1, x2)
% % d = (x1*x2').^4 + (x1*x2').^3 + (x1*x2').^2 + x1*x2' + 1E-1;
% d = exp(-1E1 .* (x1-x2).^2);
% end
diff --git a/demos/demo_GPR_closedShape02.m b/demos/demo_GPR_closedShape02.m
index 9f6a5e2554a3cf76b8ff14c14bb9532d4bd17406..ac6642463ae0390305070d9d7f738914cfd2ef58 100644
--- a/demos/demo_GPR_closedShape02.m
+++ b/demos/demo_GPR_closedShape02.m
@@ -39,24 +39,31 @@ function demo_GPR_closedShape02
% addpath('./m_fcts/');
+
%% Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-nbVarX = 2; %Dimension of x
-nbVarY = 1; %Dimension of y
nbDataRepro = 50^2; %Number of datapoints in a reproduction
+nbVarX = 2; %Dimension of x
+% nbVarY = 1; %Dimension of y
% nbRepros = 6; %Number of randomly sampled reproductions
+% SigmaY = eye(nbVarY);
p(1)=1E0; p(2)=1E-5; %Thin-plate covariance function parameters
-SigmaY = eye(nbVarY);
%% Generate data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-x = [-3 -2 -1 0 1 2 3 3 3 2 1 1 1 0 -1 -1 -1 -2 -3 -3;
- 2 2 2 2 2 2 2 1 0 0 0 -1 -2 -2 -2 -1 0 0 0 1];
-x = [x, mean(x,2), [-4, -4, 4, 4; 4 -4 -4 4]] * 1E-1;
-x = x(:,6:end); %Simulate missing datapoints
-nbData = size(x,2);
-y = [zeros(1,nbData-5), 1, -1, -1, -1, -1]; %0, 1 and -1 represent border, interior, and exterior points
+% x = [-3 -2 -1 0 1 2 3 3 3 2 1 1 1 0 -1 -1 -1 -2 -3 -3;
+% 2 2 2 2 2 2 2 1 0 0 0 -1 -2 -2 -2 -1 0 0 0 1];
+% x = [x, mean(x,2), [-4, -4, 4, 4; 4 -4 -4 4]] * 1E-1;
+% x = x(:,6:end); %Simulate missing datapoints
+% nbData = size(x,2);
+% y = [zeros(1,nbData-5), 1, -1, -1, -1, -1]; %0, 1 and -1 represent border, interior, and exterior points
+%mean(x,2)+[-1;0], mean(x,2)+[1;0], mean(x,2)+[0;-1], mean(x,2)+[0;1]
+
+x0 = [3 3 2 2 0 -1 -1 -3 -3;
+ 2 1 0 -1 -2 -2 -1 2 1];
+x = [x0, mean(x0,2), [-3;3], [-2;-3], [4;-2], [-3;3], [4;2], [2;-2] [4;4]] * 1E-1;
+y = [zeros(1,size(x0,2)), 1, -1, -1, -1, -1, -1, -1, -1]; %0, 1 and -1 represent border, interior, and exterior points
Data = [x; y];
@@ -75,10 +82,12 @@ Mu = .5 * rc * diag(1 - x' * S * x)';
K = covFct(x, x, p, 1); %Inclusion of noise on the inputs for the computation of K
[Ks, dKs] = covFct(xs, x, p);
% r(1).Data = [xs; (Ks / K * y')']; %GPR with Mu=0
-r(1).Data = [xs; MuS + (Ks / K * (y - Mu)')'];
-r(1).dData = [xs; (dKs / K * (y - Mu)')'];
-% size(r(1).dData)
-% return
+r(1).Data = [xs; MuS + (Ks / K * (y - Mu)')'];
+r(1).dData = [xs; (dKs(:,:,1) / K * (y - Mu)')'; (dKs(:,:,2) / K * (y - Mu)')']; %Gradients
+a = min(r(1).Data(3,:),0); %Amplitude
+for t=1:size(r(1).dData,2)
+ r(1).dData(3:4,t) = - exp(a(t)) * r(1).dData(3:4,t) / norm(r(1).dData(3:4,t)); %Vector moving away from contour, with stronger influence when close to the border
+end
% %Uncertainty evaluation
% Kss = covFct(xs, xs, p);
@@ -103,50 +112,57 @@ r(1).dData = [xs; (dKs / K * (y - Mu)')'];
% end
-%% Spatial plots
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-figure('PaperPosition',[0 0 14 5],'position',[10 10 2500 900]);
-colormap([.7 .7 .7]);
+% %% Spatial plots
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+idi = (Data(3,:)==1);
+idb = (Data(3,:)==0);
+ide = (Data(3,:)==-1);
+
+figure('position',[10 10 2300 1300],'color',[1,1,1]); hold on; axis off; %'PaperPosition',[0 0 14 5]
+set(0,'DefaultAxesLooseInset',[0,0,0,0]);
+set(gca,'LooseInset',[0,0,0,0]);
+colormap(repmat(linspace(1,0,64),3,1)');
%Plot center of GP (acting as prior if points on the contours are missing)
-subplot(1,3,1); hold on; axis off; rotate3d on; title('Prior');
+subplot(2,3,1); hold on; axis off; rotate3d on; title('Prior (circular contour)','fontsize',20);
% for n=2:nbRepros/2
% coltmp = [.5 .5 .5] + [.5 .5 .5].*rand(1);
% mesh(Xm, Ym, reshape(r(n).Data(3,:), nbDataRepro^.5, nbDataRepro^.5), 'facealpha',.4,'edgealpha',.4,'facecolor',coltmp,'edgecolor',coltmp); %Prior samples
% end
-mesh(Xm, Ym, reshape(MuS, nbDataRepro^.5, nbDataRepro^.5), 'facealpha',.8,'edgealpha',.8);
-mesh(Xm, Ym, zeros(nbDataRepro^.5, nbDataRepro^.5), 'facealpha',.3,'edgealpha',.3,'facecolor',[0 0 0],'edgecolor',[0 0 0]);
+mesh(Xm, Ym, reshape(MuS, nbDataRepro^.5, nbDataRepro^.5), 'facealpha',.8,'edgealpha',.8,'edgecolor',[.7 .7 .7]);
+mesh(Xm, Ym, zeros(nbDataRepro^.5, nbDataRepro^.5), 'facealpha',.3,'edgealpha',.6,'facecolor',[0 0 0],'edgecolor','none');
tl = linspace(0,2*pi,100);
plot(cos(tl)*rc, sin(tl)*rc, 'k-','linewidth',2);
view(3); axis vis3d;
-%Plot posterior distribution
-subplot(1,3,2); hold on; axis off; rotate3d on; title('Posterior');
+%Plot posterior distribution (3D)
+subplot(2,3,4); hold on; axis off; rotate3d on; title('Posterior (3d plot)','fontsize',20);
% for n=nbRepros/2+1:nbRepros
% coltmp = [.5 .5 .5] + [.5 .5 .5].*rand(1);
% mesh(Xm, Ym, reshape(r(n).Data(3,:), nbDataRepro.^.5, nbDataRepro.^.5), 'facealpha',.4,'edgealpha',.4,'facecolor',coltmp,'edgecolor',coltmp); %posterior samples
% end
-mesh(Xm, Ym, reshape(r(1).Data(3,:), nbDataRepro^.5, nbDataRepro^.5), 'facealpha',.8,'edgealpha',.8);
-mesh(Xm, Ym, zeros(nbDataRepro^.5, nbDataRepro^.5), 'facealpha',.3,'edgealpha',.3,'facecolor',[0 0 0],'edgecolor',[0 0 0]);
-plot3(Data(1,1:end-5), Data(2,1:end-5), Data(3,1:end-5), '.','markersize',18,'color',[.8 0 0]);
-plot3(Data(1,end-4), Data(2,end-4), Data(3,end-4), '.','markersize',18,'color',[0 0 .8]);
-plot3(Data(1,end-3:end), Data(2,end-3:end), Data(3,end-3:end), '.','markersize',18,'color',[0 .7 0]);
+mesh(Xm, Ym, reshape(r(1).Data(3,:), nbDataRepro^.5, nbDataRepro^.5), 'facealpha',.8,'edgealpha',.8,'edgecolor',[.7 .7 .7]);
+mesh(Xm, Ym, zeros(nbDataRepro^.5, nbDataRepro^.5), 'facealpha',.3,'edgealpha',.6,'facecolor',[0 0 0],'edgecolor','none');
contour(Xm, Ym, reshape(r(1).Data(3,:), nbDataRepro^.5, nbDataRepro^.5), [0,0], 'linewidth',2,'color',[0 0 0]);
+plot3(Data(1,idi), Data(2,idi), Data(3,idi)+.04, '.','markersize',18,'color',[.8 0 0]); %Interior points
+plot3(Data(1,idb), Data(2,idb), Data(3,idb)+.04, '.','markersize',18,'color',[.8 .4 0]); %Border points
+plot3(Data(1,ide), Data(2,ide), Data(3,ide)+.04, '.','markersize',18,'color',[0 .6 0]); %Exterior points
view(3); axis vis3d;
-subplot(1,3,3); hold on; axis off; title('Posterior');
-colormap(repmat(linspace(1,0,64),3,1)');
+%Plot posterior distribution (2D)
+subplot(2,3,[2,3,5,6]); hold on; axis off; title('Posterior (2d plot)','fontsize',20);
surface(Xm, Ym, reshape(r(1).Data(3,:), nbDataRepro^.5, nbDataRepro^.5)-2, 'FaceColor','interp','EdgeColor','interp');
-% imagesc(Xm, Ym, reshape(r(1).Data(3,:), nbDataRepro^.5, nbDataRepro^.5));
-h(1) = plot(Data(1,1:end-5), Data(2,1:end-5), '.','markersize',18,'color',[1 0 0]);
-h(2) = plot(Data(1,end-4), Data(2,end-4), '.','markersize',18,'color',[0 0 .8]);
-h(3) = plot(Data(1,end-3:end), Data(2,end-3:end), '.','markersize',18,'color',[0 .7 0]);
+quiver(r(1).dData(1,:), r(1).dData(2,:), r(1).dData(3,:), r(1).dData(4,:),'color',[0 0 .8]);
+h(1) = plot(Data(1,idi), Data(2,idi), '.','markersize',28,'color',[.8 0 0]); %Interior points
+h(2) = plot(Data(1,idb), Data(2,idb), '.','markersize',28,'color',[.8 .4 0]); %Border points
+h(3) = plot(Data(1,ide), Data(2,ide), '.','markersize',28,'color',[0 .6 0]); %Exterior points
+contour(Xm, Ym, reshape(r(1).Data(3,:), nbDataRepro.^.5, nbDataRepro.^.5), [0,0], 'linewidth',2,'color',[0 0 0]);
[c, h(4)] = contour(Xm, Ym, reshape(r(1).Data(3,:), nbDataRepro.^.5, nbDataRepro.^.5), [0,0], 'linewidth',2,'color',[0 0 0]);
% h(4) = plot3(c(1,2:end), c(2,2:end), ones(1,size(c,2)-1), '-','linewidth',2,'color',[0 0 0]);
% contour(Xm, Ym, reshape(r(1).Data(3,:)+2E0.*diag(S).^.5', nbDataRepro.^.5, nbDataRepro.^.5), [0,0], 'linewidth',1,'color',[.6 .6 .6]);
% contour(Xm, Ym, reshape(r(1).Data(3,:)-2E0.*diag(S).^.5', nbDataRepro.^.5, nbDataRepro.^.5), [0,0], 'linewidth',1,'color',[.6 .6 .6]);
-axis tight; axis equal; axis ij;
-legend(h,{'Demonstrated contour points (y=0)','Demonstrated interior points (y=1)','Demonstrated exterior points (y=-1)','Estimated contour'});
+axis tight; axis equal; axis([-.5 .5 -.5 .5]); axis ij;
+legend(h,{'Demonstrated interior points (y=1)','Demonstrated contour points (y=0)','Demonstrated exterior points (y=-1)','Estimated contour'},'fontsize',20);
% legend('boxoff');
% print('-dpng','graphs/GPR_GPIS01.png');
@@ -154,6 +170,9 @@ pause;
close all;
end
+function d = substr(x1, x2)
+ d = x1 - x2;
+end
function [K, dK] = covFct(x1, x2, p, flag_noiseObs)
if nargin<4
@@ -161,14 +180,37 @@ function [K, dK] = covFct(x1, x2, p, flag_noiseObs)
end
%Thin plate covariance function (for 3D implicit shape)
- K = 12^-1 .* (2 * pdist2(x1',x2').^3 - 3 * p(1) * pdist2(x1',x2').^2 + p(1)^3);
+ K = 12^-1 * (2 * pdist2(x1',x2').^3 - 3 * p(1) * pdist2(x1',x2').^2 + p(1)^3);
+ dK(:,:,1) = 12^-1 * (6 * pdist2(x1',x2') .* substr(x1(1,:)',x2(1,:)) - 6 * p(1) * substr(x1(1,:)',x2(1,:))); %Derivatives along x1
+ dK(:,:,2) = 12^-1 * (6 * pdist2(x1',x2') .* substr(x1(2,:)',x2(2,:)) - 6 * p(1) * substr(x1(2,:)',x2(2,:))); %Derivatives along x2
+% for i=1:size(x1,2)
+% for j=1:size(x2,2)
+% e = x1(:,i) - x2(:,j);
+% % K(i,j) = 12^-1 * (2 * pdist2(x1(:,i)',x2(:,j)')^3 - 3 * p(1) * pdist2(x1(:,i)',x2(:,j)')^2 + p(1)^3);
+% % K(i,j) = 12^-1 * (2 * norm(e)^3 - 3 * p(1) * e'*e + p(1)^3);
+% K(i,j) = 12^-1 * (2 * (e'*e)^1.5 - 3 * p(1) * e'*e + p(1)^3); %Kernel (slow one by one computation)
+% dK(i,j,:) = 12^-1 * (6 * (e'*e)^.5 * e - 6 * p(1) * e); %Derivatives (slow one by one computation)
+% end
+% end
+
% %Thin plate covariance function (for 2D implicit shape -> does not seem to work)
% K = 2 * pdist2(x1',x2').^2 .* log(pdist2(x1',x2')) - (1 + 2*log(p(1))) .* pdist2(x1',x2').^2 + p(1)^2;
-
- %RBF covariance function
+
+
+% %RBF covariance function
% K = 1E-1 * exp(-1E-1^-1 * pdist2(x1',x2').^2);
- dK = 1E-1 * exp(-1E-1^-1 * pdist2(x1',x2').^2) .* pdist2(x1',x2'); %Derivative
+% % size(x1)
+% % size(x2)
+% % dK(:,:,1) = 1E-1 * exp(-1E-1^-1 * pdist2(x1',x2').^2) .* repmat(x1(1,:)-x2(1,:), size(x1,2), 1); %Derivative
+% % dK(:,:,2) = 1E-1 * exp(-1E-1^-1 * pdist2(x1',x2').^2) .* repmat(x1(2,:)-x2(2,:), size(x1,2), 1); %Derivative
+% for i=1:size(x1,2)
+% for j=1:size(x2,2)
+% dK(i,j,:) = 1E-1 * exp(-1E-1^-1 * pdist2(x1(:,i)',x2(:,j)').^2) * (x1(:,i) - x2(:,j));
+% % e = x1(:,i) - x2(:,j);
+% % dK(i,j,:) = 1E-1 * exp(-1E-1^-1 * e'*e) * e;
+% end
+% end
if flag_noiseObs==1
K = K + p(2) * eye(size(x1,2),size(x2,2)); %Consideration of noisy observation y