add mapping estimation (still debugging)

Author: rflamary

ot/da.py

@@ -6,13 +6,15 @@
 import numpy as np
 from .bregman import sinkhorn
 from .lp import emd
-from .utils import unif,dist
+from .utils import unif,dist,kernel
 from .optim import cg
 
 
 def indices(a, func):
     return [i for (i, val) in enumerate(a) if func(val)]
 
+
+
 def sinkhorn_lpl1_mm(a,labels_a, b, M, reg, eta=0.1,numItermax = 10,numInnerItermax = 200,stopInnerThr=1e-9,verbose=False,log=False):
     """
     Solve the entropic regularization optimal transport problem with nonconvex group lasso regularization
@@ -129,34 +131,38 @@ def joint_OT_mapping_linear(xs,xt,mu=1,eta=0.001,bias=False,verbose=False,verbos
 
     if bias:
         xs1=np.hstack((xs,np.ones((ns,1))))
-        I=eta*np.eye(d+1)
+        xstxs=xs1.T.dot(xs1)
+        I=np.eye(d+1)
         I[-1]=0
         I0=I[:,:-1]
         sel=lambda x : x[:-1,:]
     else:
         xs1=xs
-        I=eta*np.eye(d)
+        xstxs=xs1.T.dot(xs1)
+        I=np.eye(d)
         I0=I
         sel=lambda x : x
 
     if log:
         log={'err':[]}
 
     a,b=unif(ns),unif(nt)
-    M=dist(xs,xt)
+    M=dist(xs,xt)*ns
     G=emd(a,b,M)
 
     vloss=[]
 
     def loss(L,G):
+        """Compute full loss"""
         return np.sum((xs1.dot(L)-ns*G.dot(xt))**2)+mu*np.sum(G*M)+eta*np.sum(sel(L-I0)**2)
 
     def solve_L(G):
-        """ solve problem with fixed G"""
+        """ solve L problem with fixed G (least square)"""
         xst=ns*G.dot(xt)
-        return np.linalg.solve(xs1.T.dot(xs1)+I,xs1.T.dot(xst)+I0)
+        return np.linalg.solve(xstxs+eta*I,xs1.T.dot(xst)+eta*I0)
 
     def solve_G(L,G0):
+        """Update G with CG algorithm"""
         xsi=xs1.dot(L)
         def f(G):
             return np.sum((xsi-ns*G.dot(xt))**2)
@@ -175,8 +181,11 @@ def df(G):
         print('{:5d}|{:8e}|{:8e}'.format(0,vloss[-1],0))
 
 
-    # regul matrix
-    loop=1
+    # init loop
+    if numItermax>0:
+        loop=1
+    else:
+        loop=0
     it=0
 
     while loop:
@@ -191,18 +200,116 @@ def df(G):
 
         vloss.append(loss(L,G))
 
+        if it>=numItermax:
+            loop=0
+
         if abs(vloss[-1]-vloss[-2])<stopThr:
             loop=0
 
         if verbose:
             if it%20==0:
                 print('{:5s}|{:12s}|{:8s}'.format('It.','Loss','Delta loss')+'\n'+'-'*32)
             print('{:5d}|{:8e}|{:8e}'.format(it,vloss[-1],abs(vloss[-1]-vloss[-2])/abs(vloss[-2])))
+    if log:
+        log['loss']=vloss
+        return G,L,log
+    else:
+        return G,L
 
-    return G,L
 
+def joint_OT_mapping_kernel(xs,xt,mu=1,eta=0.001,kernel='gaussian',sigma=1,bias=False,verbose=False,verbose2=False,numItermax = 100,numInnerItermax = 20,stopInnerThr=1e-9,stopThr=1e-6,log=False,**kwargs):
+    """Joint Ot and mapping estimation (uniform weights and )
+    """
 
+    ns,nt,d=xs.shape[0],xt.shape[0],xt.shape[1]
 
+    if bias:
+        K=
+        xs1=np.hstack((xs,np.ones((ns,1))))
+        xstxs=xs1.T.dot(xs1)
+        I=np.eye(d+1)
+        I[-1]=0
+        I0=I[:,:-1]
+        sel=lambda x : x[:-1,:]
+    else:
+        xs1=xs
+        xstxs=xs1.T.dot(xs1)
+        I=np.eye(d)
+        I0=I
+        sel=lambda x : x
+
+    if log:
+        log={'err':[]}
+
+    a,b=unif(ns),unif(nt)
+    M=dist(xs,xt)*ns
+    G=emd(a,b,M)
+
+    vloss=[]
+
+    def loss(L,G):
+        """Compute full loss"""
+        return np.sum((xs1.dot(L)-ns*G.dot(xt))**2)+mu*np.sum(G*M)+eta*np.sum(sel(L-I0)**2)
+
+    def solve_L(G):
+        """ solve L problem with fixed G (least square)"""
+        xst=ns*G.dot(xt)
+        return np.linalg.solve(xstxs+eta*I,xs1.T.dot(xst)+eta*I0)
+
+    def solve_G(L,G0):
+        """Update G with CG algorithm"""
+        xsi=xs1.dot(L)
+        def f(G):
+            return np.sum((xsi-ns*G.dot(xt))**2)
+        def df(G):
+            return -2*ns*(xsi-ns*G.dot(xt)).dot(xt.T)
+        G=cg(a,b,M,1.0/mu,f,df,G0=G0,numItermax=numInnerItermax,stopThr=stopInnerThr)
+        return G
+
+
+    L=solve_L(G)
+
+    vloss.append(loss(L,G))
+
+    if verbose:
+        print('{:5s}|{:12s}|{:8s}'.format('It.','Loss','Delta loss')+'\n'+'-'*32)
+        print('{:5d}|{:8e}|{:8e}'.format(0,vloss[-1],0))
+
+
+    # init loop
+    if numItermax>0:
+        loop=1
+    else:
+        loop=0
+    it=0
+
+    while loop:
+
+        it+=1
+
+        # update G
+        G=solve_G(L,G)
+
+        #update L
+        L=solve_L(G)
+
+        vloss.append(loss(L,G))
+
+        if it>=numItermax:
+            loop=0
+
+        if abs(vloss[-1]-vloss[-2])<stopThr:
+            loop=0
+
+        if verbose:
+            if it%20==0:
+                print('{:5s}|{:12s}|{:8s}'.format('It.','Loss','Delta loss')+'\n'+'-'*32)
+            print('{:5d}|{:8e}|{:8e}'.format(it,vloss[-1],abs(vloss[-1]-vloss[-2])/abs(vloss[-2])))
+    if log:
+        log['loss']=vloss
+        return G,L,log
+    else:
+        return G,L
 
 
 class OTDA(object):
@@ -294,6 +401,7 @@ def predict(self,x,direction=1):
 
 class OTDA_sinkhorn(OTDA):
     """Class for domain adaptation with optimal transport with entropic regularization"""
+
     def fit(self,xs,xt,reg=1,ws=None,wt=None,**kwargs):
         """ Fit domain adaptation between samples is xs and xt (with optional
             weights)"""
@@ -335,3 +443,51 @@ def fit(self,xs,ys,xt,reg=1,eta=1,ws=None,wt=None,**kwargs):
         self.G=sinkhorn_lpl1_mm(ws,ys,wt,self.M,reg,eta,**kwargs)
         self.computed=True
 
+class OTDA_mapping(OTDA):
+    """Class for optimal transport with joint linear mapping estimation"""
+
+
+    def __init__(self,metric='sqeuclidean'):
+        """ Class initialization"""
+
+
+        self.xs=0
+        self.xt=0
+        self.G=0
+        self.L=0
+        self.bias=False
+        self.metric=metric
+        self.computed=False
+
+    def fit(self,xs,xt,mu=1,eta=1,bias=False,**kwargs):
+        """ Fit domain adaptation between samples is xs and xt (with optional
+            weights)"""
+        self.xs=xs
+        self.xt=xt
+        self.bias=bias
+
+        self.ws=unif(xs.shape[0])
+        self.wt=unif(xt.shape[0])
+
+        self.G,self.L=joint_OT_mapping_linear(xs,xt,mu=mu,eta=eta,bias=bias,**kwargs)
+        self.computed=True
+
+    def mapping(self):
+        return lambda x: self.predict(x)
+
+
+    def predict(self,x):
+        """ Out of sample mapping using the formulation from Ferradans
+
+        It basically find the source sample the nearset to the nex sample and
+        apply the difference to the displaced source sample.
+
+        """
+        if self.computed:
+            if self.bias:
+                x=np.hstack((x,np.ones((x.shape[0],1))))
+            return x.dot(self.L) # aply the delta to the interpolation
+        else:
+            print("Warning, model not fitted yet, returning None")
+            return None
+

ot/datasets.py

@@ -8,8 +8,8 @@
 
 
 def get_1D_gauss(n,m,s):
-    """return a 1D histogram for a gaussian distribution (n bins, mean m and std s) 
-    
+    """return a 1D histogram for a gaussian distribution (n bins, mean m and std s)
+
     Parameters
     ----------
 
@@ -20,21 +20,21 @@ def get_1D_gauss(n,m,s):
     s : float
         standard deviaton of the gaussian distribution
 
-  
+
     Returns
     -------
     h : np.array (n,)
-          1D histogram for a gaussian distribution      
-    
+          1D histogram for a gaussian distribution
+
     """
     x=np.arange(n,dtype=np.float64)
     h=np.exp(-(x-m)**2/(2*s^2))
     return h/h.sum()
-    
-    
+
+
 def get_2D_samples_gauss(n,m,sigma):
-    """return n samples drawn from 2D gaussian N(m,sigma) 
-    
+    """return n samples drawn from 2D gaussian N(m,sigma)
+
     Parameters
     ----------
 
@@ -45,12 +45,12 @@ def get_2D_samples_gauss(n,m,sigma):
     sigma : np.array (2,2)
         covariance matrix of the gaussian distribution
 
-  
+
     Returns
     -------
     X : np.array (n,2)
-          n samples drawn from  N(m,sigma)       
-    
+          n samples drawn from  N(m,sigma)
+
     """
     if  np.isscalar(sigma):
         sigma=np.array([sigma,])
@@ -61,9 +61,10 @@ def get_2D_samples_gauss(n,m,sigma):
         res= np.random.randn(n,2)*np.sqrt(sigma)+m
     return res
 
-def get_data_classif(dataset,n,nz=.5,**kwargs):
+
+def get_data_classif(dataset,n,nz=.5,theta=0,**kwargs):
     """ dataset generation for classification problems
-    
+
     Parameters
     ----------
 
@@ -74,13 +75,13 @@ def get_data_classif(dataset,n,nz=.5,**kwargs):
     nz : float
         noise level (>0)
 
-  
+
     Returns
     -------
     X : np.array (n,d)
-          n observation of size d      
+          n observation of size d
     y : np.array (n,)
-          labels of the samples         
+          labels of the samples
 
     """
     if dataset.lower()=='3gauss':
@@ -90,10 +91,10 @@ def get_data_classif(dataset,n,nz=.5,**kwargs):
         x[y==1,0]=-1.; x[y==1,1]=-1.
         x[y==2,0]=-1.; x[y==2,1]=1.
         x[y==3,0]=1. ; x[y==3,1]=0
-        
+
         x[y!=3,:]+=1.5*nz*np.random.randn(sum(y!=3),2)
         x[y==3,:]+=2*nz*np.random.randn(sum(y==3),2)
-        
+
     elif dataset.lower()=='3gauss2':
         y=np.floor((np.arange(n)*1.0/n*3))+1
         x=np.zeros((n,2))
@@ -102,12 +103,29 @@ def get_data_classif(dataset,n,nz=.5,**kwargs):
         x[y==1,0]=-2.; x[y==1,1]=-2.
         x[y==2,0]=-2.; x[y==2,1]=2.
         x[y==3,0]=2. ; x[y==3,1]=0
-        
+
         x[y!=3,:]+=nz*np.random.randn(sum(y!=3),2)
-        x[y==3,:]+=2*nz*np.random.randn(sum(y==3),2)          
+        x[y==3,:]+=2*nz*np.random.randn(sum(y==3),2)
+
+    elif  dataset.lower()=='gaussrot'      :
+        rot=np.array([[np.cos(theta),-np.sin(theta)],[np.sin(theta),np.cos(theta)]])
+        m1=np.array([-1,-1])
+        m2=np.array([1,1])
+        y=np.floor((np.arange(n)*1.0/n*2))+1
+        n1=np.sum(y==1)
+        n2=np.sum(y==2)
+        x=np.zeros((n,2))
+
+        x[y==1,:]=get_2D_samples_gauss(n1,m1,nz)
+        x[y==2,:]=get_2D_samples_gauss(n2,m2,nz)
+
+        x=x.dot(rot)
+
+
+
     else:
         x=0
         y=0
         print("unknown dataset")
-    
+
     return x,y.astype(int)