function [ result ] = toy_sample( mani_type, N, ExParam )
%TOY_SAMPLE Generate samples from toy manifold
%
% Input:
%  - mani_type: manifold type
%    0 : 2D Plane
%    1 : Swiss Role
%    2 : Swiss Role with hole
%    3 : Corner Planes
%    4 : Punctured Sphere
%    5 : Twin Peaks
%    6 : 3D Clusters
%    7 : Toroidal Helix
%    8 : Gaussian
%    9 : Plane with 3 Peaks
%   10 : Full sphere
%   11 : Half sphere
%   12 : Cylinder
%   13 : Torus
%  - N: number of samples to generate
%  - ExParam: optional, have different meaning for different manifold types
%
% Output:
%  - result: a struct containing the following fields
%    - X: N-by-3 sample points
%    - C: Color vector, usually indicate the intrinsic generating parameter
%    - N: Number of samples, note this might not equal to the input N for
%         some types of the manifold
%    - name: A human-readable name of the manifold
%  The following fields might not be available for some manifold types
%    - IsometricEmbedding: N-by-2 or N-by-1 isometric embedding
%    - TangentSpace: 3-by-2-by-N or 3-by-1-by-N tangent spaces
%
% Examples:
%   X = toy_sample(2,1000);

if ~exist('ExParam', 'var')
    ExParam = 1.0;
end

X = zeros(N, 3);
ColorVector = zeros(N, 1);
if mani_type == 0 % 2D Plane
    k = 1;
    xMax = floor(sqrt(N));
    yMax = ceil(N/xMax);
    for x = 1 : xMax
        for y = 1 : yMax
            X(k, :) = [x,y,0];
            ColorVector(k) = y;
            k = k+1;            
        end        
    end
    result = [];
    result.X = X;
    result.C = ColorVector;
    result.N = length(ColorVector);
    result.name = '2D Plane';
elseif mani_type == 1 % Swiss Role
    t = (3*pi/2) * (1+2*rand(N, 1));
    height = 21*rand(N, 1);
    X = [t.*cos(t), height, ExParam*t.*sin(t)];
    
    % compute embedding
    t0 = 3*pi/2;
    s0 = 0.5*t0*sqrt(1+t0*t0)+asinh(t0);
    s = 0.5*t.*sqrt(1+t.*t)+asinh(t) - s0;
    
    % ground truth of tangent space
    TM = zeros(3, 2, N);
    for i = 1:N
        TM(:, 1, i) = [cos(t(i))-t(i)*sin(t(i)) 0 sin(t(i))+t(i)*cos(t(i))];
        TM(:, 1, i) = TM(:, 1, i) / norm(TM(:, 1, i));
        TM(:, 2, i) = [0 1 0];
    end
    
    result = [];
    result.X = X;
    result.C = t;
    result.N = N;
    result.name = 'Swiss Roll';
    result.IsometricEmbedding = [s height];
    result.TangentSpace = TM;
elseif mani_type == 2 % swiss roll with hole
    res = toy_sample(1, 2*N, ExParam);
    I = (res.C > 9) & (res.C < 12) & (res.X(:,2) > 9) & (res.X(:,2) < 14);
    I2 = 1:2*N; I2(I) = []; I2(N+1:end) = [];
    
    result = [];
    result.X = res.X(I2,:);
    result.C = res.C(I2);
    result.N = length(I2);
    result.name = 'Swiss Roll with Hole';
    result.IsometricEmbedding = res.IsometricEmbedding(I2,:);
    result.TangentSpace = res.TangentSpace(:,:,I2);
elseif mani_type == 3
    % Corner Planes
    k = 1;
    xMax = floor(sqrt(N));
    yMax = ceil(N/xMax);
    cornerPoint = floor(yMax/2);
    for x = 0:xMax
        for y = 0:yMax
            if y <= cornerPoint
                X(k,:) = [x,y,0];
                ColorVector(k) = y;
            else
                X(k,:) = [x,cornerPoint+(y-cornerPoint)*cos(pi*ExParam/180),(y-cornerPoint)*sin(pi*ExParam/180)];
                ColorVector(k) = y;
            end;
            k = k+1;
        end;
    end;
    
    result = [];
    result.X = X;
    result.C = ColorVector;
    result.N = length(ColorVector);
    result.name = 'Corner Planes';
elseif mani_type == 4 % Punctured Sphere by Saul & Roweis
    inc = 9/sqrt(N);   %inc = 1/4;
    [xx,yy] = meshgrid(-5:inc:5);
    rr2 = xx(:).^2 + yy(:).^2;
    [tmp ii] = sort(rr2);
    Y = [xx(ii(1:N))'; yy(ii(1:N))'];
    a = 4./(4+sum(Y.^2));
    X = [a.*Y(1,:); a.*Y(2,:); ExParam*2*(1-a)]';
    
    result = [];
    result.X = X;
    result.C = X(:,3);
    result.N = N;
    result.name =  'Punctured Sphere';
elseif mani_type == 5 % Twin Peaks by Saul & Roweis
    xy = 1-2*rand(2,N);
    X = [xy; sin(pi*xy(1,:)).*tanh(3*xy(2,:))]';
    X(:,3) = ExParam * X(:,3);
    
    result = [];
    result.X = X;
    result.C = X(:,3);
    result.N = N;
    result.name = 'Twin Peaks';
elseif mani_type == 6 % 3D Clusters
    numClusters = ExParam;
    numClusters = max(1,numClusters);
    Centers = 10*rand(numClusters,3);
    D = squareform(pdist(Centers'));
    minDistance = min(D(D>0));
    k = 1;
    N2 = N - (numClusters-1)*9;
    for i = 1:numClusters
        for j = 1:ceil(N2/numClusters)
            X(k,1:3) = Centers(i,1:3)+(rand(1,3)-0.5)*minDistance/sqrt(12);
            ColorVector(k) = i;
            k = k + 1;
        end;
        % Connect clusters with straight line.
        if i < numClusters
            for t = 0.1:0.1:0.9
                X(k,1:3) = Centers(i,1:3) + (Centers(i+1,1:3)-Centers(i,1:3))*t;
                ColorVector(k) = 0;
                k = k+1;
            end;
        end;
    end;
    
    result = [];
    result.X = X;
    result.C = ColorVector;
    result.N = length(ColorVector);
    result.name = '3D Clusters';
elseif mani_type == 7 % Toroidal Helix by Coifman & Lafon
     noiseSigma=0.05;   %noise parameter
     t = (1:N)'/N;
     t = t.^(ExParam)*2*pi;
     X = [(2+cos(8*t)).*cos(t) (2+cos(8*t)).*sin(t) sin(8*t)]+noiseSigma*randn(N,3);
     
     result = [];
     result.X = X;
     result.C = t;
     result.N = N;
     result.name = 'Toroidal Helix';
elseif mani_type == 8 % Gaussian Randomly Sampled
    X = ExParam * randn(N,3);
    X(:,3) = 1 / (ExParam^2 * 2 * pi) * exp ( (-X(:,1).^2 - X(:,2).^2) / (2*ExParam^2) );
    
    result = [];
    result.X = X;
    result.C = X(:,3);
    result.N = N;
    result.name = 'Gaussian';
elseif mani_type == 9 % Plane with 3 Peaks (See paper MLLE), Approximately isometric
    Nt = ceil(sqrt(N));
    [t s] = meshgrid(linspace(-1.5,1.5,Nt), linspace(-1.5,1.5,Nt));
    t = reshape(t, Nt*Nt, 1);
    s = reshape(s, Nt*Nt, 1);
    h = exp(-10*((t-0.5).^2 + (s-0.5).^2)) - exp(-10*(t.^2+(s+1).^2)) - exp(-10*((1+t).^2+s.^2));
    X = [t s h];
    X = X(1:N, :); % make sure we return N points
    
    result = [];
    result.X = X;
    result.C = X(:,3);
    result.N = N;
    result.name = 'Plane with 3 Peaks';
elseif mani_type == 10 % full sphere
    theta = pi*rand(N, 1);
    phi = 2*pi*rand(N, 1);
    X = [sin(theta).*cos(phi) sin(theta).*sin(phi) cos(theta)];
    
    result = [];
    result.X = X;
    result.C = cos(theta);
    result.N = N;
    result.name = 'Full Sphere';
elseif mani_type == 11 % half sphere
    theta = pi*rand(N, 1);
    phi = 2*pi*rand(N, 1);
    X = [cos(theta).*cos(phi) cos(theta).*sin(phi) sin(theta)];
    
    result = [];
    result.X = X;
    result.C = sin(theta);
    result.N = N;
    result.name = 'Half Sphere';
elseif mani_type == 12 % cylender
    theta = 2*pi*rand(N, 1);
    height = 0.2*pi*rand(N, 1);
    X = [cos(theta) height sin(theta)];
    
    result = [];
    result.X = X;
    result.C = theta;
    result.N = N;
    result.name = 'Cylinder';
elseif mani_type == 13 % Torus
    t1 = 2*pi*rand(N,1);
    t2 = 2*pi*rand(N,1);
    X = [(2+cos(t1)).*cos(t2) (2+cos(t1)).*sin(t2) sin(t1)];
    
    result = [];
    result.X = X;
    result.C = t2;
    result.N = N;
    result.name = 'Torus';
else
    error('Unknown manifold type: %d', mani_type);
end

end

