//----------------------- Random permutation of integers from 1 to N
function perm=permut(N,randOption)
  //perm=grand(1,'prm',1:N);  Standard SciLab instruction, but
  // NB.  For certain operating systems (Linux Ubuntu 14.10)
  // the "setsd" initialization does not work properly, preventing
  // exactly the same series of tests to be repeated.
  // In this case, it is better to use the code given below

  temp=1:N; tempSize=N;
  for n=1:N
    // Random rank in temp
    //	r=1+floor(tempSize*rand(1,"uniform"));
    r=1+floor(tempSize*alea(0,1,randOption));
    // Allocates a value to temp(r)
    perm(1,n)=temp(r);
    // Compact temp
    if r<tempSize then
      for d=r:tempSize-1 temp(d)=temp(d+1); end
    end
    tempSize=tempSize-1;
  end
endfunction
//----------------------------------------------------------
function pos=randAroundCombin(centre,position,randOption,distanceType)
  //distanceType=1; // 0 => Kendall-Tau
  // 1 => Cayley
  [dummy, D]=size(centre);
  // Distance 
  radius=permutDist(centre,position,distanceType);
  // A number of transpositions is chosen at random
  if radius>1 then
    nbTrans=1+floor(alea(0,radius,randOption));
  else
    nbTrans=1; radius=1;
  end

  // These transpositions are randomly applied to the "center"
  // With Kendall-Tau, the method is not guaranteed to work, 
  //so we continue to loop as long as the distance is greater than the radius
  dist=%inf;
  while dist>radius 
    pos=centre; // Initial position
    for n=1:nbTrans
      p1=0; p2=0;
      while p1==p2  // Selection of two different elements
        p1=1+floor(alea(0,D,randOption));
        p2=1+floor(alea(0,D,randOption));
      end
      // Transposition
      temp=pos(p1); pos(p1)=pos(p2); pos(p2)=temp;
    end

    if distanceType==0 then
      // Checks whether the distance is acceptable
      dist=permutDist(centre,pos,distanceType);
      //if dist>radius then printf("\n too far away"); end
    else
      dist=radius-1; // Random value, just to end the while loop
    end
  end 
endfunction
//----------------------------------------------------------
function pos=randLineCombin(centre,position,randOption)
  // Defines a random permutation “between” two others
  transpo=permutDecompCayley(centre,position,[]);
  [dummy, nTrans]=size(transpo);
  // Number of transpositions selected at random
  if nTrans>1 then
    n=floor(alea(1,nTrans,randOption));
  else
    n=1;
  end

  // These transpositions are applied to the "center", 
  pos=transpoApply(centre,transpo(:,n));
endfunction

//----------------------------------------------------------
function dist=permutDist(perm1,perm2,distanceType)
  // Distance between two permutations
  select distanceType
  case 0  // Kendall-Tau distance
    dist=permutDistKT(perm1,perm2);
  else  // Cayley distance
    dist=permutDistCayley(perm1,perm2);
  end
endfunction
//----------------------------------------------------------
function dist=permutDistKT(perm1,perm2)
  // Kendall-Tau distance
  // The number of pairs which are not in the same order
  // in the two permutations
  transpo=permutDecompKT(perm1,perm2);
  [dummy,dist]=size(transpo);
endfunction
//----------------------------------------------------------
function transpo=permutDecompKT(perm1,perm2)
  // Minimum decomposition of a permutation into a transposition sequence
  // for the Kendall-Tau distance
  [dummy,D]=size(perm1);
  transpo=[];
  dist=0;
  for d1=1:D-1  // For all pairs in perm1 ...
    for d2=d1+1:D 
      m1=perm1(d1);
      m2=perm1(d2);
      // ... we look for the rank of their first and second elements
      //      in perm2
      for d3=1:D
        if perm2(d3)==m1 then rank1=d3; 
        else 
          if perm2(d3)==m2 then rank2=d3; end
        end
      end

      if rank1>rank2 then
        dist=dist+1;
        transpo(1,dist)=d1;
        transpo(2,dist)=d2;
      end
    end
  end
endfunction 
//----------------------------------------------------------
function dist=permutDistCayley(permutInit,permutEnd)
  // Cayley distance
  // Minimum number of transpositions needed
  // to pass from one permutation to another
  // Takes longer to calculate, but is more precise than Kendall-Tau
  transpo=permutDecompCayley(permutInit,permutEnd,[]);
  [dummy,dist]=size(transpo);
endfunction
//----------------------------------------------------------
function transpo=permutDecompCayley(permutInit,permutEnd,transpo0)
  // Constructs a minimum sequence of transpositions 
  // to go from a permutation of 1:D to another permutation of 1:D
  // NB, this function is said to be recursive with transpo=[] 
  [dummy,D]=size(permutEnd);
  [dummy,p]=size(transpo0);
  permut=permutInit;
  transpo=transpo0;
  // If the permutations are identical, we do nothing
  ident=1;
  for k=1:D
    if permutInit(k) ~=permutEnd(k) then 
      ident=0;
      break;
    end
  end
  // Otherwise, k is the first non-fixed point, 
  // We wish to find a transposition which replaces permutInit(k) in the right position
  if ident==0 then
    p=p+1;
    d=k+1;
    while permutEnd(d)~=permutInit(k)
      d=d+1;
    end
    // Memorizes the transpositiona
    transpo(1,p)=k; 
    transpo(2,p)=d;
    // Applies the transposition
    temp=permut(d);
    permut(d)=permut(k);
    permut(k)=temp;
    //printf("\npermut "); for d=1:D printf(" %i",permut(d)); end
    // Recursive callup
    transpo=permutDecompCayley(permut,permutEnd,transpo);
  end
endfunction
//----------------------------------------------------------
function permut=transpoApply(permut,transpo)
  // Applies a sequence of transpositions to a permutation
  // Note: no checks of sizes and values, 
  // presumed to be OK
  [dummy2,nbTranspo]=size(transpo);
  if nbTranspo>0 then
    for n=1:nbTranspo
      r1=transpo(1,n);
      r2=transpo(2,n);
      temp=permut(r1);
      permut(r1)=permut(r2);
      permut(r2)=temp;
    end
  end
endfunction
//----------------------------------------------------------
function pr=probaTsp(D,noRepeat,t)
  // For a symmetrical “traveling salesman” type problem
  // with a single solution
  // D = number of cities
  // Probability of finding this solution by drawing at random:
  // noRepeat 0 => with repeat: the same permutation may be drawn more than once
  // noRepeat 1 => no repeat
  // t = number of draws
  select noRepeat
  case 0
    pr=1-(1-2/factorial(D-1))^t; // 2*D/factorial(D). 2 for symmetry
    //and D for cyclic permutations
  case 1
    //	if t>=factorial(D)-2*D then pr=1;
    //else
    if t==1 then 
      pr=2/factorial(D-1);  
    else
      z=2*D/(factorial(D)-t+1);
      pr=probaTsp(D,1,t-1)*(1-z)+z;
    end
    //end
  end
endfunction
//---------------------------------------------------------
function probaTspPlot(D,tMax)
  // Comparison curves for prbability of success with and without repeat
 for t=1:tMax
    pr0(t)=probaTsp(D,0,t);
    pr1(t)=probaTsp(D,1,t);
  end
  scf();
  plot2d([1:tMax],[pr0,pr1],style=[2,3]);
  xtitle(" ","Number of draws","Probability of success"); 
  legends(["with repeat";"without repeat"],[2,3],1); 
endfunction