ucascnic / CudaOT

Repository for solving discrete optimal transport problems via Cplex, Sinkhorn, FastEMD and the GPU implementation of SparseSinkhorn,
53 stars 2 forks source link

What is the from_matrix function in your python example #5

Open bearinsuke opened 4 months ago

bearinsuke commented 4 months ago

In your example 'Color transfer between images', you called an undefined function from_matrix. This function only showed up once (without define or import). Could you please give some references or corrections?

ucascnic commented 4 months ago

Thanks a lot for your pointing out, the function from_matrix could be imported from:

from pykeops.torch.cluster import cluster_ranges_centroids from pykeops.torch.cluster import sort_clusters from pykeops.torch.cluster import from_matrix

Example code:

from pykeops.torch.cluster import cluster_ranges_centroids
from pykeops.torch.cluster import sort_clusters
from pykeops.torch.cluster import from_matrix

torch.cuda.empty_cache()
D =  -torch.sum((x_centroids[:, None, :] - y_centroids[None, :, :]) ** 2, 2) + f + g.t()
#D =   -(D - f - g.t())   
thetalist = torch.logspace(-1,-12,steps=500)
for theta in  thetalist:
    #print(theta)
    #theta  = 1e-9
    theta = torch.as_tensor(theta).type(dtype)

    keep = D > torch.log(theta)* regk
    ranges_ij   = from_matrix(x_ranges, y_ranges, keep)

    a = D[keep]

    if torch.all(torch.sum(keep,dim=0)>0) and torch.all(torch.sum(keep,dim=1)>0):
        break
areas = (x_ranges[:, 1] - x_ranges[:, 0])[:, None] * (y_ranges[:, 1] - y_ranges[:, 0])[
    None, :
]
total_area = torch.sum(areas)  # should be equal to N*M
sparse_area = torch.sum(areas[keep])
print(
    "We keep {:.3e}/{:.3e} = {:.3f}% of the original kernel matrix.".format(
        sparse_area, total_area, float(100 * sparse_area / total_area)
    )
)
bearinsuke commented 4 months ago

Thank you very much! Happy spring festival!

bearinsuke commented 4 months ago

I am afraid that there is another undefined function: 'compute_distance' in the line "distance_sparse = compute_distance(d,X_i,Y_j,fall,gall,regk,ranges_ij=subset)", could you please give a reference for that?

ucascnic commented 4 months ago

Hi, we are sorry that the file "example estimating subset.ipynb" is duplicated, please refer to the Python example in https://github.com/ucascnic/CudaOT/tree/main/pythonExample/sparseOT. Thanks a lot!

bearinsuke commented 4 months ago

Thanks a lot! Here I have a further question: In your Python example "sparse transport", you implemented a robust example to solve OT via M3S. What are the meanings of your input? In other words, what is the meaning of sp_matrix, indices, a, b? It seems like the 'sp_matrix'; indicates your cost matrix while the vectors 'a' and 'b' are the two initialized dual variables. So, what are your original margins, i.e., what are the two vectors you want to solve OT? Regards

ucascnic commented 4 months ago

Hi, the two vectors are uniform distribution with all elements to be 1/n. The problem we want to solve is to map the 3-D point clouds to a 2-D Delaunay Triangulation, which means that we have some 3 D points with axis like (x, y, z) and we want to map the point to the plane (i, j). Generally, if we do not have any restraint, the cost matrix is a dense matrix. However, we may assum that any 3-D point will be mapped into some piece of Delaunay that close enought to the 3 D point. Under such condition we could compute a sparse cost matrix and we solve a sparse OT problem.

bearinsuke commented 4 months ago

Hi, thanks for your answer. What is the relationship between your mapping (3D point cloud to 2D Delaunay) and these two uniform distributions? In other words, were you trying to explain that each element in your discrete uniform distribution corresponds to a point in your point cloud or a 2D point in your Delaunay Triangulation? Can I understand your previous statement in this way: your cost matrix, which is the 'dist_matrix', defines the cost of moving one point from a 3D point cloud to one point in your Triangulation?