.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_0010_multiscale-setup.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_plot_0010_multiscale-setup.py: 1d Sparse Sinkhorn =============================================================================== This example demonstrates on a simple 1-dimensional example the basic usage of the TMultiScaleSetupGrid class for representing a point cloud with a measure on multiple resolution levels and how to use the SparseSinkhorn solver. .. GENERATED FROM PYTHON SOURCE LINES 7-24 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np import scipy.sparse import MultiScaleOT # Generate a 1D Gaussian measure over a 1D list of points pos=np.arange(32,dtype=np.double).reshape((-1,1)) mu=np.exp(-0.5*((pos-16.)/4.)**2).ravel() mu=mu/np.sum(mu) # Simple visualization plt.plot(mu) plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_0010_multiscale-setup_001.png :alt: plot 0010 multiscale setup :srcset: /auto_examples/images/sphx_glr_plot_0010_multiscale-setup_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 25-26 Now we generate the TMultiScaleSetup object .. GENERATED FROM PYTHON SOURCE LINES 26-33 .. code-block:: Python # determines how many layers the multiscale representation will have hierarchyDepth=5 # generate object MultiScaleSetup=MultiScaleOT.TMultiScaleSetupGrid(mu,hierarchyDepth) .. GENERATED FROM PYTHON SOURCE LINES 34-35 How many layers are there? .. GENERATED FROM PYTHON SOURCE LINES 35-39 .. code-block:: Python nLayers=MultiScaleSetup.getNLayers() print(nLayers) .. rst-class:: sphx-glr-script-out .. code-block:: none 6 .. GENERATED FROM PYTHON SOURCE LINES 40-41 How many points are on each layer? .. GENERATED FROM PYTHON SOURCE LINES 41-44 .. code-block:: Python print([MultiScaleSetup.getNPoints(l) for l in range(nLayers)]) .. rst-class:: sphx-glr-script-out .. code-block:: none [1, 2, 4, 8, 16, 32] .. GENERATED FROM PYTHON SOURCE LINES 45-48 Plot all versions of the measure at all layers. At the coarsest layer it is only a single point with mass 1. At each subsequent finer layer, the mass is split over more points. .. GENERATED FROM PYTHON SOURCE LINES 48-57 .. code-block:: Python for l in range(nLayers): posL=MultiScaleSetup.getPoints(l) muL=MultiScaleSetup.getMeasure(l) plt.plot(posL,muL,marker="x",label=l) plt.legend() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_0010_multiscale-setup_002.png :alt: plot 0010 multiscale setup :srcset: /auto_examples/images/sphx_glr_plot_0010_multiscale-setup_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 58-59 Create a second measure, a sum of two gaussians. Create a corresponding multiscale object. Plot both measures for comparison. .. GENERATED FROM PYTHON SOURCE LINES 59-69 .. code-block:: Python nu=np.exp(-0.5*((pos-8.)/2.)**2).ravel()+np.exp(-0.5*((pos-24.)/2.)**2).ravel() nu=nu/np.sum(nu) MultiScaleSetup2=MultiScaleOT.TMultiScaleSetupGrid(nu,hierarchyDepth) plt.plot(mu) plt.plot(nu) plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_0010_multiscale-setup_003.png :alt: plot 0010 multiscale setup :srcset: /auto_examples/images/sphx_glr_plot_0010_multiscale-setup_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 70-71 Create a cost function object for the two measures. Evaluate all pairwise costs and display as image. .. GENERATED FROM PYTHON SOURCE LINES 71-83 .. code-block:: Python costFunction=MultiScaleOT.THierarchicalCostFunctionProvider_SquaredEuclidean( MultiScaleSetup,MultiScaleSetup2) # number of points in the two measures: xres=mu.shape[0] yres=nu.shape[0] c=np.array([[costFunction.getCost(hierarchyDepth,x,y) for y in range(yres)] for x in range(xres)]) plt.imshow(c) plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_0010_multiscale-setup_004.png :alt: plot 0010 multiscale setup :srcset: /auto_examples/images/sphx_glr_plot_0010_multiscale-setup_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 84-85 Create an epsilon scaling object. Choosing the proper values for epsilon scaling and the scheduling over the multiple layers is not trivial. The following parameters should work well on most Wasserstein-2-type problems. .. GENERATED FROM PYTHON SOURCE LINES 85-91 .. code-block:: Python epsScalingHandler=MultiScaleOT.TEpsScalingHandler() epsScalingHandler.setupGeometricMultiLayerB(nLayers,1.,4.,2,2) # Check which values for epsilon scaling have been generated. This returns a list of eps values to be used on each layer. print(epsScalingHandler.get()) .. rst-class:: sphx-glr-script-out .. code-block:: none [array([4096., 2048., 1024.]), array([1024., 512., 256.]), array([256., 128., 64.]), array([64., 32., 16.]), array([16., 8., 4.]), array([4. , 2. , 1. , 0.5 , 0.25])] .. GENERATED FROM PYTHON SOURCE LINES 92-94 Now generate Sinkhorn solver object, initialize, solve, extract optimal coupling and convert it to scipy.sparse.csr_matrix. Visualize optimal coupling as image. .. GENERATED FROM PYTHON SOURCE LINES 94-115 .. code-block:: Python # error goal errorGoal=1E-3 # Sinkhorn solver object SinkhornSolver=MultiScaleOT.TSinkhornSolverStandard(epsScalingHandler, 0,hierarchyDepth,errorGoal, MultiScaleSetup,MultiScaleSetup2,costFunction ) # initialize and solve SinkhornSolver.initialize() SinkhornSolver.solve() # extract optimal coupling kernelData=SinkhornSolver.getKernelCSRDataTuple() kernel=scipy.sparse.csr_matrix(kernelData,shape=(xres,yres)) plt.imshow(kernel.toarray()) plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_0010_multiscale-setup_005.png :alt: plot 0010 multiscale setup :srcset: /auto_examples/images/sphx_glr_plot_0010_multiscale-setup_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 116-117 Print the optimal transport cost part of the primal objective (cost function integrated against optimal coupling) and compare it with manually computed value. .. GENERATED FROM PYTHON SOURCE LINES 117-119 .. code-block:: Python print(SinkhornSolver.getScoreTransportCost()) print(np.sum(kernel.toarray()*c)) .. rst-class:: sphx-glr-script-out .. code-block:: none 23.965531931307662 23.96553193130768 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.425 seconds) .. _sphx_glr_download_auto_examples_plot_0010_multiscale-setup.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_0010_multiscale-setup.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_0010_multiscale-setup.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_0010_multiscale-setup.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_