Correction TP transport optimal


/// Transport 1D

function [iX,iY] = transport1D(X,Y)
 [Xs,iX] = gsort(X,'g','i'); // on a Xs = X(iX)
 [Ys,iY] = gsort(Y,'g','i');
endfunction


// Calcule l'assignent en 2D et evalue le coût total du transport

function [Z,c]=transport2D(X,Y,N,e) 
// X et Y deux nuages de points en 2D ; matrice 2 colonnes et n lignes. La fonction renvoie un nuage de points Z qui est obtenu par déplacement de X vers les points de Y.
 
  theta = %pi/2*(2*rand(N,1)-1);
  Z=X;
  for k=1:N
    utheta = [cos(theta(k));sin(theta(k))];
    Zt = Z*utheta;
    Yt = Y*utheta;
    [iZ,iY] = transport1D(Zt,Yt);
    Z(iZ,:) = Z(iZ,:)+e*(Yt(iY)-Zt(iZ))*utheta';
    D=sqrt(sum((Z(iZ,:)-Y(iY,:)).^2,2));
   if (max(D)<1D-05) then
      mprintf('stop a k=%d\n',k);
      break;
    end
  end
 c = sum(sum((X(iZ,:)-Y(iY,:)).^2,2))
endfunction