% Example 7.8 Polar Smoothing

load L1bpm
% Pick the topics we will look at.
t1 = 8;
t2 = 11;
% Reduce the dimensionality using Isomap.
options.dims = 1:10;    % These are for ISOMAP.
options.display = 0;
[Yiso, Riso, Eiso] = isomap(L1bpm, 'k', 7, options);
% Get the data out.
X = Yiso.coords{2}';
% Get the data for each one and plot
ind1 = find(classlab == t1);
ind2 = find(classlab == t2);
plot(X(ind1,1),X(ind1,2),'.',X(ind2,1),X(ind2,2),'o')
title('Scatterplot of Topic 8 and Topic 11')
legend({'Topic 8';'Topic 11'})

% First lets do the polar smoothing just for
% the first topic. Get the x and y values.
x = X(ind1,1);
y = X(ind1,2);
% Step 1.
% Normalize using the median absolute deviation.
% We will use the Matlab 'inline' functionality.
md = inline('median(abs(x - median(x)))');
xstar = (x - median(x))/md(x);
ystar = (y - median(y))/md(y);
% Step 2.
s = ystar + xstar;
d = ystar - xstar;
% Step 3. Normalize these values.
sstar = s/md(s);
dstar = d/md(d);
% Step 4. Convert to polar coordinates.
[th,m] = cart2pol(sstar,dstar);
% Step 5. Transform radius m.
z = m.^(2/3);
% Step 6. Smooth z given theta.
n = length(x);
J = ceil(n/2);
% Get the temporary data for loess.
tx = -2*pi + th((n-J+1):n);
% So we can get the values back, find this.
ntx = length(tx);  
tx = [tx; th];
tx = [tx; th(1:J)];
ty = z((n-J+1):n);
ty = [ty; z];
ty = [ty; z(1:J)];
tyhat = loess(tx,ty,tx,0.5,1);
% Step 7. Transform the values back.
% Note that we only need the middle values.
tyhat(1:ntx) = [];
mhat = tyhat(1:n).^(3/2);
% Step 8. Convert back to Cartesian.
[shatstar,dhatstar] = pol2cart(th,mhat);
% Step 9. Transform to original scales.
shat = shatstar*md(s);
dhat = dhatstar*md(d);
xhat = ((shat-dhat)/2)*md(x) + median(x);
yhat = ((shat+dhat)/2)*md(y) + median(y);
% Step 10. Plot the smooth.
% We use the convex hull to make it easier
% for plotting.
K = convhull(xhat,yhat);
plot(X(ind1,1),X(ind1,2),'.',X(ind2,1),X(ind2,2),'o')
hold on
plot(xhat(K),yhat(K))

% Now use the polarloess function to get the
% other one.
[xhat2,yhat2] = ...
    polarloess(X(ind2,1),X(ind2,2),0.5,1);
plot(xhat2,yhat2)
hold off
title('Scatterplot of Topic 8 and Topic 11')
legend({'Topic 8';'Topic 11'})




