using PyPlot
using NtToolBox

n = 500; #1000 points leads to very slow computation

x = rand(2,n);

v = 3*pi/2*(.1 + 2*x[1,:])
X  = zeros(3,n)
X[2,:] = 20*x[2,:]
X[1,:] = - cos(v).*v
X[3,:] = sin(v).*v;

ms = 200
elev = 20; azim = -110;

fig = figure(figsize=(15,11))
ax = gca(projection="3d")

#swiss roll
scatter3D(X[1,:], X[2,:], X[3,:], c=get_cmap("jet")((X[1,:].^2+X[3,:].^2)/100), s=ms, lw=0, alpha=1)

#params
xlim(minimum(X[1,:]),maximum(X[1,:]))
ylim(minimum(X[2,:]),maximum(X[2,:]))
zlim(minimum(X[3,:]),maximum(X[3,:]))
axis("off")
ax[:view_init](elev, azim)

D1 = repeat(sum(X.^2, 1), outer=(n,1))
D1 = D1 + D1' - 2*X'*X
D1[D1.<0] = NaN
D1=sqrt(D1);

k = 6;

DNN, NN = zeros(size(D1)), zeros(size(D1))
for i in 1:size(D1,2)
    DNN[i,:], NN[i,:] = sort(D1[i,:]), sortperm(D1[i,:])
end
NN = Int.(NN[:,2:k+1])
DNN = DNN[:,2:k+1];

B = repeat((1:n)', outer=(k,1))
A = sparse(vec(B), vec(NN'), ones(k*n));

W = sparse(vec(B),vec(NN'), vec(DNN'));

fig = figure(figsize=(15,11))
ax = gca(projection="3d")

#swiss roll
scatter3D(X[1,:], X[2,:], X[3,:], c=get_cmap("jet")((X[1,:].^2+X[3,:].^2)/100), s=ms, lw=0, alpha=1)

#graph
(I,J) = findn(A)
#(I,J) = (vec(B), vec(NN))
xx = hcat(X[1,I],X[1,J])
yy = hcat(X[2,I],X[2,J])
zz = hcat(X[3,I],X[3,J])

for i in 1:length(I)
    plot3D(xx[i,:], yy[i,:], zz[i,:], color="black")
end
    
#params
xlim(minimum(X[1,:]),maximum(X[1,:]))
ylim(minimum(X[2,:]),maximum(X[2,:]))
zlim(minimum(X[3,:]),maximum(X[3,:]))
axis("off")
ax[:view_init](elev, azim)

#W = readdlm("W.txt")
D = full(W)
D = (D + D')/2.;

D[D .== 0] = Inf;

D = D - diagm(diag(D))
D[isnan(D)] = Inf;

include("NtSolutions/shapes_7_isomap/exo1.jl")

## Insert your code here.

Iremove = find(D[:,1].==Inf);

D[D .== Inf] = 0;

include("NtSolutions/shapes_7_isomap/exo2.jl");

(L,U) = eig(Xstrain*Xstrain' / n);
Xstrain1 = U'*Xstrain;

Xstrain1[:,Iremove] = Inf;

#plot size
figure(figsize = (15,6))

#plot points
scatter(Xstrain1[2,:], Xstrain1[1,:], c=get_cmap("jet")((X[1,:].^2+X[3,:].^2)/100), s=ms, lw=0, alpha=1)

#plot vertices
(I,J) = findn(A)
xx = hcat(Xstrain1[2,I], Xstrain1[2,J])
yy = hcat(Xstrain1[1,I], Xstrain1[1,J])

for i in 1:length(I)
    plot(xx[i,:], yy[i,:], color="black")
end

#params
xlim(minimum(Xstrain1[2,:]-1),maximum(Xstrain1[2,:])+1)
ylim(minimum(Xstrain1[1,:]-1),maximum(Xstrain1[1,:])+1)
axis("off");

Y = [v'; X[2,:]']
Y[1,:] = rescale(Y[1,:], minimum(Xstrain[1,:]), maximum(Xstrain[1,:]))
Y[2,:] = rescale(Y[2,:], minimum(Xstrain[2,:]), maximum(Xstrain[2,:]));

#plot size
figure(figsize = (15,6))

#plot points
scatter(Y[1,:], Y[2,:], c=get_cmap("jet")((X[1,:].^2+X[3,:].^2)/100), s=ms, lw=0, alpha=1)

#plot vertices
(I,J) = findn(A)
xx = hcat(Y[1,I],Y[1,J])
yy = hcat(Y[2,I],Y[2,J])

for i in 1:length(I)
    plot(xx[i,:], yy[i,:], color="black")
end 

#params
axis("off")
xlim(minimum(Y[1,:]-1), maximum(Y[1,:])+1)
ylim(minimum(Y[2,:]-1), maximum(Y[2,:])+1);

include("NtSolutions/shapes_7_isomap/exo3.jl")

## Insert your code here.

figure(figsize=(10,7))
plot(stress, ".-");

(L,U)=eig(Xstress*Xstress' / n)
(L,I) = (sort(L), sortperm(L))
U = U[:,I[2:3]];

Xstress1 = U'*Xstress;

Xstress1[:,Iremove] = Inf;

#plot size
figure(figsize = (15,6))

#plot points
scatter(Xstress1[2,:], Xstress1[1,:], c=get_cmap("jet")((X[1,:].^2+X[3,:].^2)/100), s=ms, lw=0, alpha=1)

#plot vertices
(I,J) = findn(A)
xx = [Xstress1[2,I] Xstress1[2,J]]
yy = [Xstress1[1,I] Xstress1[1,J]]

for i in 1:length(I)
    plot(xx[i,:], yy[i,:], color="black")
end

#params
axis("off")
xlim(minimum(Xstress1[2,:]-1),maximum(Xstress1[2,:])+1)
ylim(minimum(Xstress1[1,:]-1),maximum(Xstress1[1,:])+1);
<script>
  $(document).ready(function(){
      $('div.prompt').hide();
  });
</script>