function [model,parms] = loreta_similarity(data,class_labels,init_model,parms)
% [model,parms] = loreta_similarity(data,class_labels,init_model,parms)
%
% Input:   
% -- data             d X N matrix (d = instance dimension, N = number of instances)
% -- class_labels     N X 1 vector
% -- init_model       structure with init_model.A, init_model.B being d X k matrices where k is the model rank. The model is W = A*B'.
%   parameters:
%   -- parms.step_size   (default 1)
%   -- parms.rseed       random seed (default 1)
%   -- parms.num_steps   user defined
%   -- parms.inner_loop  size of inner loop for counting losses and reporting progress (default 1000)
%   -- parms.verbose     report once every parms.inner_loop steps (default 0)
%   -- parms.use_psd     use PSD version of Loreta instead of general
%   low-rank. Uses only init_model.A. The model is W=A*A'
%
% Output:
% -- model            structure with model.A, model.B being d X k matrices
% with k the model rank. The model is W = A*B', unless parms.use_psd = 1 in which case the model is W=A*A'.
%                     model.loss_steps is an indicator vector containing all loss events = all model update events
% -- parms


    [d,N] = size(data);

    if ~isfield(parms,'step_size')
        parms.step_size = 1;
        fprintf('step size value set to 1 \n')
    end

    [dA,kA] = size(init_model.A);
    [dB,kB] = size(init_model.B);

    if dA~=d 
        error('data dim is %d while model.A dim is %d',d,dA);
    elseif dB~=d
        error('data dim is %d while model.B dim is %d',d,dB);
    end

    if kA~=kB
        error('model.A and model.B are of incompatible sizes %d and %d',kA,kB);
    end



    if ~isfield(parms,'rseed')
        parms.rseed = 1;
    end

    rand('twister',parms.rseed);
    if ~isfield(parms,'inner_loop')
        parms.inner_loop = 1000;
    end

    if ~isfield(parms,'verbose')
        parms.verbose = 0;
    end

    if ~isfield(parms,'use_psd')
        parms.use_psd = 0;
        fprintf('Using Loreta-1 low-rank model \n');
    else
        if parms.use_psd == 0
            fprintf('Using Loreta-1 low-rank model \n');
        else
            fprintf('Using Loreta PSD model \n');
        end
    end

    

    if ~isfield(parms,'num_steps')       
        fprintf('no number of steps provided, please set parms.num_steps \n');
    end

   

   % rearrange data
   %-----------------------
   [unused_values, inds] = sort(class_labels); %#ok
   class_labels = class_labels(inds);
   data = data(:,inds);
   classes = sort(unique(class_labels));
   num_classes = length(classes);
  
   % Translate class labels to serial numbers 1,2,...
   new_class_labels = zeros(size(class_labels));
   for i=1:length(classes)
     new_class_labels(class_labels==classes(i)) = i;
   end
   class_labels = new_class_labels;
   class_sizes = zeros(num_classes,1);
   class_start = zeros(num_classes,1);
   for k=1:num_classes
     class_sizes(k) = sum(class_labels==k);
     class_start(k) = find(class_labels==k, 1, 'first');
   end    
   %------------------------


    A = init_model.A;
    if parms.use_psd == 0
        B = init_model.B;
    end    
    clear init_model;
    
    % initialize pseudo inverses - assumes A and B are of rank k
    Apinv = (A'*A)\A';
    if parms.use_psd == 0
        Bpinv = (B'*B)\B';
    end
    %---------

    num_steps = parms.num_steps;
    step_size = parms.step_size;
    inner_loop_size = parms.inner_loop;
    outer_loss = zeros(1,floor(num_steps/inner_loop_size));
    inner_loss = zeros(1,inner_loop_size);
    inner_l = 0;
    ctime = clock;



    for i=1:num_steps

        %sample a query/positive/negative trio
        q_ind =floor(rand*N)+1;
        q_label = class_labels(q_ind);
        p_ind = class_start(q_label)+floor(rand*class_sizes(q_label));
        n_ind = floor(rand*N)+1;        
        n_label = class_labels(n_ind);
        while n_label==q_label
            n_ind = floor(rand*N)+1;
            n_label = class_labels(n_ind);
        end
        %---------------

        x1 = data(:,q_ind); 
        x2 = data(:,p_ind)-data(:,n_ind); 
        
        if parms.use_psd == 0
            pred = (x1'*A)*(B'*x2);    
            if pred<1
                inner_loss(i-inner_l*inner_loop_size) = 1;
                [A,B,Apinv,Bpinv] = online_step_loreta1(A,B,Apinv,Bpinv,x1,x2,1,step_size);            
            end
        else            
            pred = (x1'*A)*(A'*x2);    
            if pred<1
                inner_loss(i-inner_l*inner_loop_size) = 1;
                [A,Apinv] = online_step_loreta_psd(A,Apinv,x1,x2,1,step_size);            
            end
        end           

        if mod(i,inner_loop_size) == 0              
            inner_l = inner_l+1;
            outer_loss(inner_l) = mean(inner_loss);
            inner_loss = zeros(1,inner_loop_size);
            if parms.verbose == 1       
                c = clock;
                fprintf('iteration %d out of %d, time elapsed %g sec \n',i,num_steps,etime(c,ctime));
            end
        end
    end
   
           

    %set output
    model.A = A;
    if parms.use_psd == 0
        model.B = B;
    end
    model.loss_steps = outer_loss;