% [tr, anom] = create_ts2(ntraj, nclass, nanom, nanomclass, trlen, randangle)
% create data sets
%
% tr and anom will contain the normal and anomalous trajectories, one
% trajectory for each line, in the form [x(1), ..., x(trlen), y(1), ..., y(trlen) ]
%
% ntraj: number of trajectories in each cluster of normal data
% nclass: number of clusters of normal data
% nanom: number of trajectories in each anomalous 'cluster' (generally 1 or 2,
% depending if the data set is split in a train and a test set or not)
% nanomclass: number of 'clusters' of anomalies
% trlen: trajectory length
% randangle: trajectory randomness
%
% typical usage:
%
%        [TS,anom] = create_ts2(100 * 2, 5, 2, 10, 16, 0.7);
%        trainset = TS(1:2:end, :);  
%        testset = TS(2:2:end, :); 
%        trainset = [ trainset; anom(1:2:end,:) ];
%        testset  = [ testset; anom(2:2:end,:) ]; 
%
% Author: 
% Claudio Piciarelli <claudio.piciarelli at dimi.uniud.it>
%


function [tr, anom] = create_ts2(ntraj, nclass, nanom, nanomclass, trlen, randangle)

    % discretization intervals
    gridsize = 4;   % gridsize x gridsize cells for starting position
    nangles = 5;    % divide the -90:+90 interval in nangles slices
    
    %%%
    %create a grid of all the possible trajecotry parameters
    featspace = zeros( (gridsize-1) * 4 * nangles^3, 1 );

    % choose nclass+nanomclass trajectories
    for i=1:nclass+nanomclass
        idx = find( ~featspace );  
        chosen = floor( rand() * length(idx) + 1 );
        featspace( idx( chosen ) ) = 1;  % the choose a free trajectory prototype in the tr. feature space

        tmp = ceil ( idx(chosen) / nangles^3 );   % convert the global index into startside/startcell/angle1/2/3 parameters
        startside(i) = ceil(tmp / (gridsize-1));
        startcell(i) = mod( tmp, (gridsize-1) ) + 1;
        tmp = mod( idx(chosen), nangles^3 ) + 1;
        angle1(i) = ceil(tmp / nangles^2);
        tmp = mod( tmp, nangles^2 ) + 1;
        angle2(i) = ceil( tmp / nangles);
        angle3(i) = mod(tmp, nangles) + 1;
    end
    
    % create normal trajectories
    tr = zeros( ntraj*nclass, trlen*2);
    for i=1:nclass
        tmp = create_traj( ntraj, startside(i), startcell(i), angle1(i), angle2(i), angle3(i), trlen, randangle, gridsize, nangles );
        tr(ntraj*(i-1)+1:ntraj*i, :) = tmp;
    end
    
    % create anomalies
    anom = zeros( nanom * nanomclass, trlen*2);
    for j=1:nanomclass
        tmp = create_traj( nanom, startside(j+nclass), startcell(j+nclass), angle1(j+nclass), angle2(j+nclass), angle3(j+nclass), trlen, randangle, gridsize, nangles );
        anom(nanom*(j-1)+1:nanom*j, :) = tmp;
    end

    % shift from [0-1] range to [-1,1] 
    tr = tr .* 2;
    tr = tr - 1;
    anom = anom .* 2;
    anom = anom - 1;

end


% create a single trajectory with prototype number idx (startside(idx),
% startcell(idx), angle1(idx), angle2(idx), angle3(idx)) and put it in
% tr(count,:). The trajectory will have trlen length, the perturbation
% angle is randangle. Gridsize and nangles are the discretization steps of
% the trajectory prototype feature space.
function tr = create_traj( trnum, startside, startcell, angle1, angle2, angle3, trlen, randangle, gridsize, nangles )    
    
    % set the boudaries of the starting cell
    if( startside == 1) % bot
        miny = 0;
        minx = (startcell-1) / gridsize;
        dir = pi / 2;
    elseif( startside == 2) % right
        minx = (gridsize-1)/gridsize;
        miny = (startcell-1) / gridsize;
        dir = pi;
    elseif( startside == 3) % top
        miny = (gridsize-1) / gridsize;
        minx = 1/gridsize + ((startcell-1) / gridsize);
        dir = pi * 3 / 2;
    else % left
        minx = 0;
        miny = 1/gridsize + ((startcell-1) / gridsize);
        dir = 0;
    end

    % set the starting angles 
    startdir  = dir + (((angle1-1) * (pi/nangles)) + rand() * (pi/nangles) - pi/2) /2; % /2 --> smaller variance in the starting angle
    seconddir = (angle2-1) * (pi/nangles) + rand() * (pi/nangles) - pi/2;
    thirddir  = (angle3-1) * (pi/nangles) + rand() * (pi/nangles) - pi/2;

    % rotation matrices for the three angles
    mat1 = [ cos(startdir), -sin(startdir); sin(startdir), cos(startdir) ];
    mat2 = [ cos(seconddir), -sin(seconddir); sin(seconddir), cos(seconddir) ];
    mat3 = [ cos(thirddir), -sin(thirddir); sin(thirddir), cos(thirddir) ];
    
    tr = zeros(trnum, trlen*2);
    for j = 1:trnum
        % starting pos
        posx = minx + 1/(gridsize*4) + rand() / (gridsize*2);
        posy = miny + 1/(gridsize*4) + rand() / (gridsize*2);
        % starting dir
        d = mat1 * [1/trlen, 0]';
        % create the trajectory!
        for i=1:trlen
            if( i == floor(trlen/3) )
                d = mat2 * d;
            elseif ( i == floor(trlen*2/3) )
                d = mat3 * d;
            end
            theta1 = (rand()-0.5)*(randangle/10);
            theta2 = (rand()-0.5)*randangle;
            randmat1 = [ cos(theta1), -sin(theta1); sin(theta1), cos(theta1) ];
            randmat2 = [ cos(theta2), -sin(theta2); sin(theta2), cos(theta2) ];
            d = randmat1 * d;  % cumulative randomness
            d2 = randmat2 * d; % non-cumulative randomness
            % add the point
            posx = posx + d2(1) ;
            posy = posy + d2(2);
            tr(j,i) = posx;
            tr(j,i+trlen) = posy;
        end
    end
end
        
    