function [A,b]=mrm_solver(Y,W,L,m,numOuterIters)
%[A,b]=mrm_solver(Y,W,L,m,numOuterIters)
%Solves sum (max(0, sign(L(i,j))*[W(:,i)'AY(:,j)-b(i)] )*abs(L(i,j) ) 
%Solution via gradient descend on convex. Gradients are updated easily
%each calculation of objective we rebalnce biases. 
dimY = size(Y,1);
dimX = size(W,1);
disp('begining hinge loss solver, using linear kernel'); 
N = size(Y,2);
K = size(W,2);
EPS = 1.0e-10; 
WtW = (W'*W);
YtY = (Y'*Y); 
M = repmat(m,1,N); 
b = zeros(size(m));
if(nargout>1)
    b = optim_otzu_locally( M , L , zeros(size(m)), m*0.02 ); %one shot tought solution
end;
[kVecPlus,jVecPlus] = find(L>=EPS);
[kVecMinu,jVecMinu] = find(L<=-EPS);

T = zeros(K,N); %outer product weights matrix

kVecPlusFull = kVecPlus;
jVecPlusFull = jVecPlus;
kVecMinuFull = kVecMinu;
jVecMinuFull = jVecMinu;

for outerIdx=1:numOuterIters

kVecPlus = kVecPlusFull;
jVecPlus = jVecPlusFull;
kVecMinu = kVecMinuFull;
jVecMinu = jVecMinuFull;
%preallocation for xsi,psi will be updated later
xsi = zeros(length(kVecPlus),1);
psi = zeros(length(kVecMinu),1);
dXsi = zeros(length(kVecPlus),1);
dPsi = zeros(length(kVecMinu),1);
weightXsi = dXsi;
weightPsi = dPsi;
%if no bias then at zero all hinges are in the money. 
dT = zeros(K,N); %outer product weights matrix

TYtY = T*YtY;
WtWTYtY = WtW*TYtY; 
if(exist('backupT')) %not first iteration
    if(nargout>1)
        b = optim_otzu_locally ( M-WtWTYtY , L , b , m*0.01) ; %change in margin one percent up or down towards optimum
    end;
end;
xsi = m(kVecPlus)-WtWTYtY(sub2ind(size(WtWTYtY),kVecPlus,jVecPlus))-b(kVecPlus);
psi = m(kVecMinu)+WtWTYtY(sub2ind(size(WtWTYtY),kVecMinu,jVecMinu))+b(kVecMinu); 
weightXsi = abs(L(sub2ind(size(L),kVecPlus,jVecPlus)));
weightPsi = abs(L(sub2ind(size(L),kVecMinu,jVecMinu)));
posXsiIndices = find(xsi>EPS);
posXSIdTind = sub2ind(size(dT),kVecPlus(posXsiIndices),jVecPlus(posXsiIndices));
dT(posXSIdTind) = dT(posXSIdTind) - abs(L(posXSIdTind)); 
posPsiIndices = find(psi>EPS);
posPSIdTind = sub2ind(size(dT),kVecMinu(posPsiIndices),jVecMinu(posPsiIndices));
dT(posPSIdTind) = dT(posPSIdTind) + abs(L(posPSIdTind)); 

dTYtY = dT*YtY;
WtWdTYtY = WtW*dTYtY; 
dXsi = -WtWdTYtY(sub2ind(size(WtWdTYtY),kVecPlus,jVecPlus));
dPsi =  WtWdTYtY(sub2ind(size(WtWdTYtY),kVecMinu,jVecMinu));

currentObj = sum(weightPsi.*max(psi,0))+sum(weightXsi.*max(xsi,0)); 

if(exist('backupT'))
    if(currentObj>=(backupObj-EPS))
        T = backupT;
        dT = backupDT;
        b = backupB; 
        psi = backupPsi ; xsi = backupXsi ; 
        dPsi = dBackupPsi ; dXsi = dBackupXsi ;     
    end;
end;

backupT = T;
backupDT = dT; 
backupB = b;
backupPsi = psi ; 
backupXsi = xsi; 
dBackupPsi = dPsi ; 
dBackupXsi = dXsi; 

backupObj = sum(weightPsi.*max(psi,0))+sum(weightXsi.*max(xsi,0)); 
if(outerIdx>=numOuterIters)
    break;
end;

dXsiFull = dXsi;
dPsiFull = dPsi;
xsiFull = xsi;
psiFull = psi;
    
indices4xsi = randperm(length(xsiFull));
indices4xsi = indices4xsi(1:min(100,length(indices4xsi)));
indices4psi = randperm(length(psiFull));
indices4psi = indices4psi(1:min(100,length(indices4psi))); 

xsi = xsiFull(indices4xsi);
psi = psiFull(indices4psi);
dXsi = dXsiFull(indices4xsi);
dPsi = dPsiFull(indices4psi);

kVecPlus = kVecPlusFull(indices4xsi);
jVecPlus = jVecPlusFull(indices4xsi);
kVecMinu = kVecMinuFull(indices4psi);
jVecMinu = jVecMinuFull(indices4psi);

numInnerIter = ceil(exp(rand(1)*5));
for iter=1:numInnerIter
    %gradient descend, so same sign means go down (we do xsi-=dXsi) 
    xsi4quit = find((xsi>0)&(dXsi>0));
    psi4quit = find((psi>0)&(dPsi>0));
    xsi4enter = find((xsi<0)&(dXsi<0));
    psi4enter = find((psi<0)&(dPsi<0)); 
    if((length(xsi4quit)+length(psi4quit)+length(xsi4enter)+length(psi4enter))<1)
        break;
    end;
    %all those who collided
    [minXsiQuit,idxXsiQuit] = min([xsi(xsi4quit)./dXsi(xsi4quit) ; Inf]);
    [minPsiQuit,idxPsiQuit] = min([psi(psi4quit)./dPsi(psi4quit) ; Inf]);
    ttcQuit = min(minXsiQuit,minPsiQuit); %get into first change, hinge hit. 
    [minXsiEnter,idxXsiEnter] = min([xsi(xsi4enter)./dXsi(xsi4enter) ; Inf]);
    [minPsiEnter,idxPsiEnter] = min([psi(psi4enter)./dPsi(psi4enter) ; Inf]);
    ttcEnter = min(minXsiEnter,minPsiEnter); %get into first change, hinge hit. 
    ttc = 1.01*min(ttcQuit,ttcEnter); %overcome hinge struddle. 
    %update - gradient descend itself
    T = T-(ttc*dT) ; xsi = xsi-(ttc*dXsi) ; psi = psi-(ttc*dPsi);  %reducing gradient (minimization prob.) 
    %Detect what have changed sign. a xsi or a psi, in or out
    %update dA accordingly by subtracting relevant outer-product
    s0 = 0; %sign undefined at first
    j0 = -1 ; k0 = -1; %undefined
    if(ttcQuit<ttcEnter)
        if(minXsiQuit<minPsiQuit) % a xsi quit
            s0 = -1;j0 = jVecPlus(idxXsiQuit) ; k0 = kVecPlus(idxXsiQuit) ;
        else %psi quit
            s0 = 1;j0 = jVecMinu(idxPsiQuit) ; k0 = kVecMinu(idxPsiQuit) ;
        end;
    else
        if(minXsiQuit<minPsiQuit) % a xsi eneterd
            s0 = 1;j0 = jVecPlus(idxXsiQuit) ; k0 = kVecPlus(idxXsiQuit) ;
        else %psi eneterd
            if((idxPsiQuit>length(kVecMinu))|(idxPsiQuit>length(jVecMinu)))
                ttcQuit
                ttcEnter
                disp('found bugubug');
                return;
            end;
            s0 = -1; j0 = jVecMinu(idxPsiQuit) ; k0 = kVecMinu(idxPsiQuit) ;
        end;
    end; %hinge change catch
    s0 = s0*abs(L(k0,j0));
    dT(k0,j0) = dT(k0,j0) + s0;
    wtw = WtW(:,k0) ; yty = YtY(:,j0) ; 
    dXsi = dXsi + (  s0*(wtw(kVecPlus)).*(yty(jVecPlus)) ) ; 
    dPsi = dPsi -  (  s0*(wtw(kVecMinu)).*(yty(jVecMinu)) ) ; 
end;%main inner loop
%check if consider the inner iteration or not.
    
end;%main outer loop
WT = W*T;
A = WT*(Y');

disp('hinge-loss optimized'); 


function b = optim_otzu_locally( M , S ,b0,steps)
%b = optim_otzu_locally(m,s,b0)
%optimize per k sum_j max(( M(k,j)-sign(S(k,j))*b(k) ),0)abs(S(k,j))
%This optimized by testing b0+m9k)*0.01 and - to check direction.
%Then run in steps of 1% of margin to finish
b = b0;
for k=1:size(M,1)
    Sline = S(k,:);
    sg = sign(Sline);
    weight = (abs(Sline))';
    mu = M(k,:);
    bk = b(k);
    dk = steps(k);
    bk = reshape(bk+(-10:1:10)*dk,21,1);
    
    [val,idx] = min ( (max((repmat(mu,21,1)-bk*sg),0))*weight ) ;
    b(k) = bk(idx);
end;