Skip to content
Snippets Groups Projects
Commit e9adb51d authored by Sylvain CALINON's avatar Sylvain CALINON
Browse files

Gradient added

parent 9561c765
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
demos/demo_GPR_closedShape01.m
+ 2
− 2
View file @ e9adb51d
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
......
demos/demo_GPR_closedShape02.m
+ 81
− 39
View file @ e9adb51d
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment