// You are free to use, redistribute and modifiy this code if you include
// as a comment the author of the original code 
// (c) Aniko Ekart & Sandor Zoltan Nemeth, 2009.

//base1 will have the columns of ebase for which index(i)=1 and base2 those of ubase that have index(i)=0
function [base1,base2]=split_indices(ebase,ubase,index)
base1=[]
base2=[]
[n1,n2]=size(ebase)
[n3,n4]=size(index)
if (n1==n2 & n1==n4 & n3==1) then
   n=n1
   if (det(ebase)==0) then 
     mprintf('zero determinant\n')
     abort
   end
   else
   //not a square matrix
   mprintf('not square')
   abort
end

for i=1:n
  if (index(i)==1) then
       base1=[base1 ebase(:,i)]
  else
      base2=[base2 ubase(:,i)]
  end
end
endfunction

//function to solve the ecuation x transposed times e  
//equals alfa times e transposed times e 
//alfa has size m
//x has size n
//e has size nxm
function alfa=solve_system(x,e)
alfa=[]
[n,m]=size(e)
//if there is a system to solve
if ~(m==0) then
[nx,mx]=size(x)
if ~(mx==n) then 
  mprintf('size wrong')
  abort
end

// det(ee) must be non zero!
ee=e*e'
inverse=inv(e'*e)
alfa=x*e*inverse
else 
alfa=[]
end
endfunction


//function to solve both the ecuation in alfa and bbeta
// x has size n
// base has size n
//index has size n and contains zero for each component that is not used and one for each component that is used

function [alfa,bbeta]=alfa_bbeta(x,base,index)
alfa=[]
bbeta=[]
//calculate the polar base
ubase=-inv(base')

//produce the reduced base that only has the components given in index
[base1,base2]=split_indices(base,ubase,index)

//solve the ecuation for the cone
alfa=solve_system(x,base1)

//solve the ecuation for the polar cone
bbeta=solve_system(x,base2)

endfunction

//moves the indices corresponding to the negative alfas into base2 (swaps the zero into one in index)
//and the indices corresponding to the negative bbetas into base1 (swaps the one into zero in index)  
//returns in changed the value 1 if at least one swap was done
//returns the number of negative alfas and negative bbetas
function [new_index,changed,neg_alfa,neg_bbeta]=move_index(alfa,bbeta,index)
new_index=index
changed=0
neg_alfa=0
neg_bbeta=0
[ma,na]=size(alfa)
[mb,nb]=size(bbeta)
[m,n]=size(index)

//if alfa is not empty
if(na>0) then
j=1;
for i=1:n
 if index(i)==1 then
    if alfa(j)<0 then
         new_index(i)=0
         changed=1
         neg_alfa=neg_alfa+1
    end
    j=j+1
 end 
end
end

//if bbeta is not empty
if (nb>0) then
j=1
for i=1:n
 if index(i)==0 then
    if bbeta(j)<0 then
       new_index(i)=1
       changed=1
       neg_bbeta=neg_bbeta+1
    end
    j=j+1
 end
end 
end

endfunction


//calculates the projection based on the split between the alfas and bbetas
function [y]=calculate_projection(x,base,index,alfa)
[m,n]=size(index)
y=zeros(n,1)
j=1
for i=1:n
  if (index(i)==1)
    y=y+alfa(j)*base(:,i)
    j=j+1
  end
end
y=y'
endfunction  

//iterate until no indices change, ie a positive solution is found, returns alfa,bbeta, the index split, 
//the nr of iterations and all indices for all iterations and the actual projection
//no memory of previously visited index combinations!
//results written to testfile and loopfile
function [alfa,bbeta,new_index,iteration_count,all_indices,projected]=iterate(x,base,index,testfile,loopfile)
changed=1
new_index=index
all_indices(1,:)=index;
projected=[]
iteration_count=0
fd=mopen(testfile,'a')
//mfprintf(fd,'%d ',index')
//mfprintf(fd,'\n')

while (changed==1)
[alfa,bbeta]=alfa_bbeta(x,base,new_index)
[new_index,changed,neg_alfa, neg_bbeta]=move_index(alfa,bbeta,new_index)
iteration_count=iteration_count+1
if (changed==1) then 
  //mfprintf(fd,'%d ',new_index');
  //mfprintf(fd,'  Neg alfa: %d  Neg bbeta: %d  Changes: %d \n',neg_alfa,neg_bbeta,neg_alfa+neg_bbeta)
  for i=1:iteration_count
     if new_index==all_indices(i,:) then
      //a loop was found
        //mfprintf(fd,'loop back to: ');
        //mfprintf(fd,'%d ',all_indices(i,:)');
        //mfprintf(fd,'\n\n');
        //mclose(fd);
        fl=mopen(loopfile,'a');
        mfprintf(fl,'x: \n')
        mfprintf(fl,'%f ',x')
        mfprintf(fl,'\n');
        mfprintf(fl,'base: \n')
        [nn,mm]=size(base)
        for k=1:nn
           mfprintf(fl,'%f ',base(k,:)')
           mfprintf(fl,'\n\n')
        end
       
        for j=1:iteration_count
           mfprintf(fl,'%d ',all_indices(j,:)');
           mfprintf(fl,'\n');
        end
        mfprintf(fl,'%d ',new_index');
        mfprintf(fl,'\n alfas: ');
        if (alfa~=[]) then mfprintf(fl,'%f ',alfa');
        end
        mfprintf(fl,'\n bbetas: ');
        if (bbeta~=[]) then mfprintf(fl,'%f ',bbeta');
        end
        mfprintf(fl,'\n\n');
        mclose(fl);
        mclose(fd);
        return;
     end
  end
  all_indices(iteration_count+1,:)=new_index
end

end
projected=calculate_projection(x,base,new_index,alfa)

//mfprintf(fd,'x: \n')
//mfprintf(fd,'%f ',x')
//mfprintf(fd,'\n')
//mfprintf(fd,'base: \n')
//[nn,mm]=size(base)
//for k=1:nn
//     mfprintf(fd,'%f ',base(k,:)')
//     mfprintf(fd,'\n\n')
//end




//for j=1:iteration_count
//  mfprintf(fd,'%d ',all_indices(j,:)');
//  mfprintf(fd,'\n');
//end

//mfprintf(fd,'%d ',new_index');
//mfprintf(fd,'\n alfas: ');
//mfprintf(fd,'%f ',alfa');
//mfprintf(fd,'\n bbetas: ');
//mfprintf(fd,'%f ',bbeta');
//mfprintf(fd,'\n\n');

//calculate statistics on all_indices: 
[mm,nn]=size(x)
n=max(mm,nn)
[pp,qq]=size(alfa)
sizealfas=max(pp,qq)
changes=zeros(1,n) //the changes for each dimension
it_changes=zeros(1,iteration_count-1)//the changes for each iteration
for i=1:n //for each dimension
  for j=2:iteration_count //for each iteration
     if ~(all_indices(j,i)==all_indices(j-1,i)) then //a change was made
        changes(i)=changes(i)+1
        it_changes(j-1)=it_changes(j-1)+1
     end
  end
end

//how many dimensions indices changed 0,1,2,3,...,iteration_count-1 times
changebins=zeros(1,iteration_count);

for ii=1:n
    changebins(changes(ii)+1)=changebins(changes(ii)+1)+1
end

totalincreases=0
for k=1:iteration_count-2
     if (it_changes(k+1)>it_changes(k)) then
             totalincreases=totalincreases+1
     end
end

propalfas=sizealfas/n

totalchanges=sum(changes)
totaliterations=iteration_count-1 //the last one does not change anything, that is the loop stopping condition, hence the -1

//proportion of alfas in solution, total nr. of changes, total nr. of iterations, nr of iterations where nr. of changes increased
//nr. of indices that changed once, nr. of indices that changed twice, etc.
mfprintf(fd,'%f ',propalfas)
mfprintf(fd,' %d ',totalchanges)
mfprintf(fd,' %d ',totaliterations)
mfprintf(fd,' %d : ',totalincreases)
if (totaliterations>0) then mfprintf(fd, '%d ',changebins')
end  
mfprintf(fd,'\n');
mclose(fd)
endfunction

function [good]=check_solution(x,base,index,alfa,bbeta)
good=0
ubase=-inv(base')
[base1,base2]=split_indices(base,ubase,index)
x1=x*base1
x2=alfa*base1'*base1
x3=x*base2
x4=bbeta*base2'*base2
if (max(abs(x1-x2))<0.00001) & (max(abs(x3-x4))<0.00001) then
   good=1
end
endfunction

//generates non-negative monotone cone of size n
function [monotone]=gen_monotone(n)
for i=1:n
   for j=1:i
     monotone(i,j)=1
   end
   for j=i+1:n
     monotone(i,j)=0
   end
end
endfunction

//generates internal point of non-negative monotone cone of size n
function [x]=gen_x(n)
for i=1:n
 x(1,i)=i
end
endfunction

//generates vector x on a face of cone with given base)
function [x]=gen_face_x(base)
  [nn,n]=size(base)
  if ~(nn==n) then error('base must be a square matrix'); abort
  end 
  ubase=-inv(base') 
  n1=round(rand()*(n-1)+1)
  n2=n-n1
  alfas=rand(1,n1)*20-10
  alfas(1)=0
  bbetas=rand(1,n2)*20-10
  
  x=zeros(n,1)
  // alfas are used for indices 1:n1
  for i=1:n1
    x=x+alfas(i)*base(:,i)
  end  

  //bbetas are used for indices n1+1:n
  for i=1:n2
    x=x+bbetas(i)*ubase(:,n1+i)
  end
  x=x'
endfunction

function [base]=gen_random_base(n)
base=rand(n,n)*10-5
endfunction

function [x]=gen_random_x(n)
x=rand(1,n)*10-5
endfunction

//returns the iteration nr, percentage of alfa for m random examples of size n 
//results (indices for all iterations) written to testfile
function [ics,alfas]=generate_test(m,n,testfile,loopfile)
ics=[]
alfas=[]
for i=1:m
//   x=gen_random_x(n)
   base=gen_random_base(n)
   x=gen_face_x(base)
   
   index=ones(1,n)
   [alfa bbeta n_index ic all_indices projection]=iterate(x,base,index,testfile,loopfile)
   ics=[ics ic]
   alfa_count=0
   for j=1:n
     if n_index(j)==1 then alfa_count=alfa_count+1
     end
   end
   alfas=[alfas alfa_count/n]
   //mprintf('x: ')
   //mprintf('%f ',x')
   //mprintf('  proj: ')
   //[mp,np]=size(projection)
   //if (mp>0)& (np>0) mprintf('%f ',projection')
   //end
   
end
endfunction

//the same as generate_test but the base is fixed and given as parameter
function [ics,alfas]=gen_test_fix_base(m,base,testfile,loopfile)
[n1,n]=size(base)
if ~(n1==n) then error('base must be a square matrix'); abort
end 
ics=[]
alfas=[]

for i=1:m
   x=gen_random_x(n)
   index=ones(1,n)
   [alfa bbeta n_index ic all_indices projection]=iterate(x,base,index,testfile,loopfile)
   ics=[ics ic]
   alfa_count=0
   for j=1:n
     if n_index(j)==1 then alfa_count=alfa_count+1
     end
   end
   alfas=[alfas alfa_count/n]
   mprintf('x: ')
   mprintf('%f ',x')
   mprintf('  proj: ')
   [mp,np]=size(projection)
   if (mp>0) & (np>0) mprintf('%f ',projection')
   end
end
endfunction


//creates a random base of size n and random values 
//for alfas and bbetas and creates the x as the sum of alfa(i)*base(i)+bbeta(i)&ubase(i)
function [x,n1,proj1,projection]=test_proj_fix_base(base)
[nn,n]=size(base)
if ~(nn==n) then error('base must be a square matrix'); abort
end 
  ubase=-inv(base')
  n1=round(rand()*n)
  n2=n-n1
  alfas=rand(1,n1)*10
  bbetas=rand(1,n2)*10
  
  x=zeros(n,1)
  // alfas are used for indices 1:n1
  for i=1:n1
    x=x+alfas(i)*base(:,i)
  end  
  proj1=x'
  //bbetas are used for indices n1+1:n
  for i=1:n2
    x=x+bbetas(i)*ubase(:,n1+i)
  end
  x=x'
  
  index=ones(1,n)
  [alfa bbeta n_index ic all_indices projection]=iterate(x,base,index,'test','testloop')  
endfunction

function [x,n1,proj1,projection]=test_projection(n)
  base=gen_random_base(n)
  ubase=-inv(base')
  n1=round(rand()*n)
  n2=n-n1
  alfas=rand(1,n1)*10
  bbetas=rand(1,n2)*10
  
  x=zeros(n,1)
  // alfas are used for indices 1:n1
  for i=1:n1
    x=x+alfas(i)*base(:,i)
  end  
  proj1=x'
  //bbetas are used for indices n1+1:n
  for i=1:n2
    x=x+bbetas(i)*ubase(:,n1+i)
  end
  x=x'
  
  index=ones(1,n)
  [alfa bbeta n_index ic all_indices projection]=iterate(x,base,index,'test','testloop')  
endfunction


//runs the generate_test 10 times for 10000 random examples each for given size n
function run_10()
 //   generate_test(10000,10,'10000_10-1.test','10000_10-1.loop')
    generate_test(10000,10,'10000_10-2.test','10000_10-2.loop')
    generate_test(10000,10,'10000_10-3.test','10000_10-3.loop')
    generate_test(10000,10,'10000_10-4.test','10000_10-4.loop')
    generate_test(10000,10,'10000_10-5.test','10000_10-5.loop')
    generate_test(10000,10,'10000_10-6.test','10000_10-6.loop')
    generate_test(10000,10,'10000_10-7.test','10000_10-7.loop')
    generate_test(10000,10,'10000_10-8.test','10000_10-8.loop')
    generate_test(10000,10,'10000_10-9.test','10000_10-9.loop')
    generate_test(10000,10,'10000_10-10.test','10000_10-10.loop')
endfunction
