14
$\begingroup$

Recently I asked a question about the morphing between two functions and got two excellent answers. The accepted answer is using the continuous optimal transport theory. This seems to be very suitable for this problem.

enter image description here This figure is made with Adobe Illustrator, MA solution is desirable.

However, I have difficulties to convert the code from symbolic to pure numeric one. In particular, I do not know how to numerically construct inverse functions and derivatives efficiently. My intention is to apply the code of Federico to the pair of two functions such as shown below

f[x_]:=UnitBox[x+3]
g[x_]:=UnitTriangle[x-3]

I take the liberty to copy the symbolic code here:

F[x_] = Integrate[f[x], {x, -∞, x}];
G[x_] = Integrate[g[x], {x, -∞, x}];
Ginv[q_] = InverseFunction[G][q];
T[t_, x_] = (1 - t) x + t Ginv[F[x]] // Simplify;
dT[t_, x_] = D[T[t, x], x] // Simplify;
ParametricPlot[Evaluate@Table[
   {T[t, x], f[x]/dT[t, x]}, {t, 0, 1, .1}],
   {x, -10, 5}, PlotRange -> All, AspectRatio -> 1/2]

I am seeking a pure numeric solution that can be further applied to any pair of interpolation functions. f[x] and g[x] presented above is just a simple example that cannot be integrated symbolically and because piecewise functions are hard to invert symbolically too. I've selected them because it is known that MA is not able to integrate UnitBox and UnitTriangle symbolically.

Edit

MichaelE2 suggested to provide interpolation functions. Below are two strongly truncated realistic data to work with

dataA= "1:eJx9m0mOpclShUuIPVAMWQAh75s9MKoVIAZPYsTgsUiWxXd+99/dLlyLLKmyKiPSrlt37FgT//If//XX3/7njz/++Ps/8q9/+8+///ff/sH+31//zH/8a/mZc+Y0z68Wev339vPXn3w1/OSYcu2xhf0rpcYXHUGtxDETv41UYop5ZiMo5dRC60jPvY1R0/AFpRnKqFtU7uP53ldQLrXEOCtPiTm01P0XjZln6GPMUnrKLcfSjKDSIn/+vDW0wBerL6iO3nrcb+qzjmQE1Zk7fxxCDanytJR9QTwjbQPxoWXMaQT11EsNocxUWtJruysIrTBRf5WLOT9e+6claLQRI/LXizC9L6fKoK9XSmvVOK2EODuPLDkRAiPM4GuGdVIo+TVSDsNoVuIMtc4Wz4t877eBV3PjtxIJFxxdjKCMUj2lHkYticBMvtPkqHTCUWoYp5VaeykNBcOspXde+P8F/bkFxUSYrHDDXCGaKCq9zMxTntfk2sYvzsddpEjdXsu9R6sZ1ilhtJ5Gk715si+IzwvbAnxmJnivoIpKgTf1vIyICX1BoWPQuP0/SRIrCPukVFJLg5f1GH/JfUIwxZ1mcUQ+2dio6osBcStnI7J8QRmfj74DCVsH4zT5noxpxPckffucvo0Uu3mEEhIyQsFxBkQaBq6hjdcXpJIrCF36fI0Z9SyjGnrmB9pyjg1UGL+k7PMt7f3O2WM3qdZDlqXzC8O83BVETNf4wmMJGL0bQUVeAF5HKZg9huLbSN7o9UQ2ZjLQL0uH3PMbZeSvL0hPL0dSBeSuoJHniGkmfge2Z0eUKyjOSOyONytjn8ZrwDlRVBeCEPY1+gFJms2a4862QVqWi4+zBn0V0H9iP/dfHpQCj98lgvJX8nFa/CG26hxA1EraHKvvNNKVCvLiFgV3dCNoElg1lJOIv7wIBEWxF0ixdahXECBM7SAia+zULKq5D0foTnnMO/0HpdSolppSjNK/Fa/dQ+z8Q7CmIs9OQqaGWG5Riz9yZlRheTKWMhW8DEFQqTGcch1zmeUKAiBTm4ovvbkg1LMRggIocaIojxHTFdRa1J8RQQ+wUfY8G+UfkGhQlA8g8c1XUG+ZAjPyY+k8W3MBMisJGuCwc2CElQRbkHQZGX13DSnzlxcRJUXFWMgWKM/hQHYijgYm7C+dGz6uZV4/S9hlH6Ro/cnvJ0PST8IB0JE2Z5moD3Pz5Qi4HtcKIfPotdwHFT4EExGRD/DV7IJI/uHDIAaz8WwRIEpOuoJkPAoINfJJthmaV0MQhF/G/kSglkwZV9CktIKZ/APk1kho+PHYhKJvnPQAah+n5R+Rx1wFbo+ROpzLFSSGGOKhtMLJKwi7ATCg70pE0tB/Ua3ASDwcotc+r6AWlWM4czGIPH4JowogwoleyKZaGtVGVJzWQdGu8OuS3DqbfwrY2udJkZnz8RoQAxqRslFslkdXv31AEKgGHr1Pqr2NK0icBgM1Uk5pTbL47qdIE/oPexTYtnarGkUYxKxiXs+TMKXL13AM9bzuTxSQ1BSNoNlxhYj3QwxKdJEWQSplLwAqOmu/gvBXEfOGEhDc9DO/BCTtGA9Kb12bqRvVSiQNIQYLHCLf4L8owXjz49xHO3Deqhbh3BSXXR9JXq+sIYj0nqefSaSbUS1RPPlyX+1DxKt+9idYYX3rbBSJNapBYpI40wps8rH6uAb8hXi830XdX1xTYZDLgZfFjRZNdeTkNvoCZMF2RrWbIcDj0wbu16rz8QUFmN7DQp+STZybVCOmQlfMPw0W6Fd97/PcOhdVE7gB/hZF1MfxzP0xYZFCR5Ca3fDSH+GbUQ2FwaLQN641H/sBUzV6+xtBAIL8qMYXwaIkWvuEKxzBpf18L25pix5KnJqSKwi50BOc/1Qq6ET1jC14JwrH20E/Db4pIrgJqt7UyVInh8/XEmQTJnL6a2zyUdaAOrUQu5zHETyARJA6zX67NdNjUWf5u2Rr3jyF2u+rNpQbzZDjVdZWwY6wmqQcWuHapu80cAucOFwEgxvNgqAsUsxWnGEv/0G9BqDhYv+4VS1qiEMlp5at8CCMvHhEEPk+ys21ng3LAoSzMGSHPZTbtzWexZhHN4ieIZA9R0g3edFX/0hX6gsqZYZdRwVdVDFDIImIquKwIT26+KiQI4VWcouwQWuNakQGFF6cTTUP1tZ9GxFEyCm71XiixnBjAOapB6ua59WlOIIyeH0GPkWTOkP7yRxCIL+zt/ZLqtGv5HL6kAbrN/FIAFKFNRAUNkDr/XgsSvxTHLF1LheNoH2P2jsNVbN9QTynHp8FGK0parB4yiXosWGtZK+EpB9xuvhQ+oeLYu1+XwSbDDFiw/WiGdxBDYJa7n3NRAQ5NBsWRODv6JV2uMK/fXzM4tAb96FTNEY2HlvQaO79lBp/8X7C0BpfzhVKddHf0z3i+cdn61f4JfmpxxCIeLqoNs0wazT15zGfIjK+eO3PLShp2rXqMYyEv2XnNE3JGnewZmLSfZBSop/5YyG0zVSkdBg2qr3BOtx2FkEVI+U7Fct2AEkDAzuGCm7CM93CjyD6zHnJmiZRV1DSmBMIP23KL7mv76oHQ7BTN6opSWN4podPR1O6j0ahKBUPyQYijWrow3MhbPvF7VtDs51G0NPRvKSmEL9XTqQhCJpvrgSCkXq2pkqQHWm+NaRRvY1mxHkQ1XuhKlQvQRBENo1hcrYZzaLqT1a7tlJtuB1WFEMUPr42olQap4kTany5swds8mxN3SIUQ39HlVlV1QgixJ4VyTsy/UW1oWnqOw+nds9pRlkUfToaAO9FPV+zIaK6CZYweXSrGbSZEBvnc34RRGRAVd7eAXFtmAeRrDyypo0N7oD+KcgQ8fpOvGsJ1meQAc3wjktdCIkqNrRf43RYdl8UnyUR1t+dRXMntA+QQmjXeFYejrV/CIJOpKO2T2gQpG5lU1BSCWL/oVnVvGe8Y4WUvdUDDIFwHPEdMBHF88NleYAie8NB8+/HEKakuQqn4+uror0xFLS/Si/bjR6kIYdOJ4R3uijYt4oBBFWTtwVoOXzb8mzF+E5i8WzCogUiPEaZJqzLomrFJcXIqeDhmRiCt719pCv9snZY68XDnTwjCJC4y0CAzU7n4zMgg+CcCY6vmYpePvQKwmEWc5FOHw4SzkTBfw88HH/fMcawQ34FVdFWaoN+SV5jjSA+qbbbNvZgFWvaR5JpOwnzt0ZmKwaDooAemo69jZwegF+w6fWZS0GQA7vs5bCi/H8EqSJi7HcD5HJrBAWKerlVKFrqELX6om94mUX7hrBbM3APAnoSv8yPaqbpFjbcq6TYfFNDvAjrY2u4aDO+n9rukBWbgww/W8HUHu/0WsYwDyIqAJBUFkbl+C1dt2K0AmeELz5j1zIiEjR4VMUdqdVj1lp0QL4O01d/bQlIKDTsavdWvQ8usw4/ujUYpwSRSC1bQZBcEviF/DG/+P7PLaiqKTsWGiuol6XxN0SmzbCDzO3zEEPEh9MIz/Cx21XfAdEhRXZsuFkG44VRtneWDsXqwyqm7jzMZxyhA4DwrTvfiqmwpncGzuvax4MGdUzFfq/kXBaLHAHeW2Iy0Zs+BWm3B6jsOWjzTSRQLPWeGpT4odnoz5pw5+v4tnDYmpHzOdxVSk7TyhExo4d5maB7soCcKAucVRpx+aFZ1d6a0v1UOkGEK6h1Fd+7uligdYJIJwPxpEfwFdPdy0h3UDwszQ+inJoArRo+3XMVrbfJ9GTmcsMC49SVie4MFs3vLrHSwh0qGM8cpNX2AYyzqCnFgguC3SKEIGh0C7fv/PAZAYaUQe6/5ddXjTyi2zhEv6dgkVFjr6KGcXVvwyNW9HFa3s3DrEYYxmURaIEQgQfLgO68UXckIoy3mH1QT9iQGkWyYr23uqN9+kGAaN5RGo83zpcM9cl8yAp5d1SAIB1edeM0W6YDJH8C36JXqzHxEz+D0Lcd0HnAtIKyhlDPMdYa7nvlldAVgB523vS0D0E9FJ0svAMp3/laUudbPGLr40PQCAr2N/XLL4KoWfNqRkB9qsbf1Vv7quMu2VMOaIn+zktk3PFpbJGrFuPyxLfToFdQ0dJy3rncpyBIJGjAexZc/ZJp1Aa1wvHl6N1W2ACxrbijaSX77FH9DBGVz/FedJltgwStMaG2RA9NaX5k410K7LmgqcFS/QBKzaHeBI71pLUfR881UrgXC82CSNC5WwNH4AIUNyXJd5yVOWkGLzpSc9LHg6qWx0qR5wTJY0V8K0/Plw+TENZngLU8AWRP6Z2+TZz2t6oFvHsdMsreqmibjQx+E+z3Ml05UZOAO/zWVl9yZD36jh+KVMTQpav2y7VfUG3jH9rUu/mKxY4/H8zUmUl+cjB9258fYqj54ilEI5qFRSXTgSnyZi1Aonc39wymWi63pUrR7PM1qA2o1KkvInEjuPbRRrTfXSX5eDA/YmeKQNAMdY1Hv0H1ntyFNqPdna6R/hLD/1fRKoW0x4W3GApwvbvcUaJdMUztQ0LVNYX64/KtJd9yuhpTs1zOw2gF+dXSmMBZ/N99Dsal6TDrzktgEx7QKDfp0O0Rlr+dAm45BYpyt++khbnhocSTt9TL+fRH7WsMbjlQ+nhXQlBU4yzq6Zg6DysrZbwBo44Ye2/3NTXaeTdU4ZnBou4D0PHbiduWI9i89xK6VLlyBAJa5Nay+uw0nH4DOdq33AOOUqZZ45CxQC3olNbpa/s2yt0DbS0T7kUJvYDRS+DfNUgfzyGVdvCuHFhFuicuuhu+cqCH5EHR1OuhU+HbRcGWo0C+Nzdx2FHepCBpx6NS+DjBFaPiFw6RLjRAhgUBOk13UjPsRYefFCp85d4k4bmPpYIGK5E437XNG8IgBx447pHUTNmohbVoNRSAmyC50fOMEOzRliV3/YknIGUuFhm/dYfn5ZCe2/r0YO8siVF1Kik8rAUkc6OQ/MzBnrUloxaR3DV2FD4Blt3P9amx8z2zE4bd5yj8NFPOO/W61/goPrSXOxPlsS7yXjnasjf1mU+BjOvS/YucKIJAZ3g3LnEatR5QrrqZf/SuwQmeKMIy+qUZwgfzHP52UhDu4IrfDtG2HApKvO1la81M8Krwi9qXX57m5ER8xrz1zhPp1Q09UEPY1KbtxqB+63m2HPSZ55JVcw17o63JIFX2bayGRzO03Ksl3+uBEiwxrOJfQFt+OysPeqJIeO0XmXOsVq+olX/XHHgNub4t67ecqgH7EVNqSNfM69izz7qxySvrUV0KxPr82EksNrfW3RXAtQdK/dvwbhcDwHve6+xQ7RJBLBrpulBb3eJ0zUNC0e295+JttmrcVbSkDWG85Tp9a7+2HCzQ7/36WCfBr3kUnzCPsFsPP7dgIFSBF3iQaBsUXQhqb5ZfN3gDt6jGGyu/2yptLYxauWlqFwmeVUdGduWQOz2E25/kZtyl4+BetKJeWVM8nhE1mcANp4MHUu0P0hCUAFiPL2f0QUPrgXH5Qck2u6jVVFkqzp7DN99dkNmY7M+I2D1EDlR6+sS6Of7X1nTLUVDYH1qx2aXlrH46Ymw2/HU/v+U8GX62eS3Zn6OgpsW1AFyn39HPdrQCG85aHasaf2lFwxcx726kv53TbTnlmV7cnzOymywdEabn4Hj1h92j3yLqWimauyzbBCaqiAjvYY7fthBLDp0H2W4OF6fVC320fukbDLLXLCGnPj/P8jJDGs+PqSY1nayp78zi2955y6GYF/ujanZUrz091GfMd1717VJkU3W1A5enqi02coZ+AE93ZHtO7cfhFHRfVNWZmpEDjuoeqf22E9lyMF4P91fNVi/Cs1OKx/vThY99/hcRQ4hM";
dataB= "1:eJx9m0uuNLeRhYWG12DJPfPY6AsG39xDj7SChgcGPOqBe5FeVn+HjyTLTpYA6er/q24kGY8TJx7557/+7+9/++cvv/zyjz/wn//++z/+72//cf7p9z/xP/8Vf1procS2/snVhf/JP7//xqfuJzpfXA5lfcoX+fBfBP02BeXoQ/T8qD6a9+YPOVajpVZ9tuByCrE/5HIg31wLqQ1JoaaUJejXISikFFJwrcbmY/Y5X+XUFlrytbYYiw851BKPA6UQYi1myZUUCp/Gu6BUq3mr80jNmh2Cik85mQtlfCG2+824TS421OOyS9mXQ1DjrCnVcZ5ai0tXQYVvBlfKvBzCDpsliy1ZspxDsJxKdHa1Gc8Kza1HxhLqISdU45ddtOJqih4t3A9kXKotHaBNy4egHIOVllMrrbaWqn9xoiko12CuZn70xzYLh9FSs4g/1JCtuubwtbs35sgTy3JGl9Kh6uzNoveeSxXDoVr+ciDz1Zflt8FhmENQkhj8HmsGZxHVXwWlgooz5+hWw9/i4UUon5tGpCQrqZbwxa9lipjqPDxHG4JGfOCNiaB1y8m6h13kOKybVgDU5sphs5Iz4RODt+yqiy5/Mb48MMcRZVYxUTtsVgkNCyHIvVCWtfISH9NmEQfygNHwI3R9eHXNhF3TbfvNzJW7qvH4VoOTHgz/LQMlpqBmUlDKuXCoDCKlu/H1BUAtDH8EB91hs1aSywqLLDUlwPEesEQivzsCsnp0bW0Jsh8cmVhLmfv1EMntLojgStlNbMR8MaUtyCz4kHzpYgjr/AXU8Gjc1i/HRmw8BNWaQZUoXETtePY90BChJy5JHG7L8RgRxbgBIMVKuLsjz7FAGAD+KRbcyNUtCPUUQiMtp89frM+3Wqsgew81gi37LSi2YgBtySOBpG+BRvzw2IHVim6uuQVlWdTlIjQj8Lnj3Y0kxvJSQpB5VsTaTyH8sANolD23ti8Q4pLj6o4fPhnfdyFsOWhNkCiP1+cufpFDCCiJztgnjz6q9j+on1Alj3Gr4skzV1WHn1YdGW/iOqotNW9BoCbYKi8ireDT8XozBAGdxEBbIJJ3QvM/Mbds2B9EikQtKekloS1BXN/CdKKCqkPbgki+EuOXQXPxV0GwhYpm6kCjTBY8rtYIPEzBfYOCFte8eRGCApBn85v4pW/P1QIuZimT9afVqk+3AEGQQ9PW0VhWBpnbFsQJsBu0yYYfJbvhdQDbUbcbGR9IBijzFoTxHXm3dVTwaCh8ERQiTGwCKUQCSNqCSoWeBRJePxEB/ZZkfx2CcmtGBppBCbRtMfxVdvCQafziv9gs5ygdTuKDqsMDjkQzRxUH6X4Bu8rlfjE0SJ4fAQl4QRf8FkT2xcUgO17IbzCbG1yjTXHrvJij98TJFhSQRKDFEYjig/cAgRhmlyemKTP742pcN3QOOqla+RJpCY4O1i+8jiWeV8vilQDLogRvBGsKirUGH8uKEPRQD0GhEh9w0O5ixonvAQIIUx24eSRz47vLaq4qxS4egi3uAUKR4WJYSkCohdOPiC9ONeDcv1cyU1AA053vuU84Imq7BXWXAOtEZmxc8i4IX4O3dM8lRuQJC7DlZJyXf3tGI5Fc6QNyjCSLO8+kRkVkWw7Mi991boaPveH1b0OOL9jLwRBGXifA9708ZkyCqCFGlPd6Hi9S1FYdA7qlA0FwU24TV3pBl3dv9I503J4My6GeE3kgGGpLvTATHme7lVYAaYHj11WqAMjuwP0e7FQV0xA43F3VRFKzMkxSxYDPlEY2gaNgR0EjdDB+wX1YEwlipWNL/kxpQSUuwTEeYl8DzWUelSceQzuwz5FkkUuIlO6LVLv2xYswFJVh17Qs44fHTUGoN2cjwmaiSldlk/6o8AiLnvKJfxR6MKymTA+uLIblvwmiuoLytMUggj84H9ARVDjM+OB0N4dEEFrgHL0lEHSVndJgWCSzRlZZlLDUm7LxORnGr8IaZe8qDc7XitC1LYSp+QZGCOJZzS/Kp+x8kkfqoADxW4U3UXsXRHAk21jky9MOMWpBbI9/rWZJux9IsGEuP0Sk2EmL0T9+bbMPoPLoLgjukFtbyC/acQiCVSltjqoC/nRFEQQRFJx8xRqofJQgGJ4whZz0hA3nc/cTZVGWJ5ZIIvG8GswkVQrrQdShB7e8jyDl4zIKOsJK5eoBR3AQ9aOGZ0DFyv1qVEEVaMuTrekWBzOijIU0lFHq4/f5RmlAHOgDiSNNsgEE2hbUxAF15mGLeEfaDl094iftc0d5JYiBQXjhZgdJQvaLIKcezFODqAFxFA+h9kPNeoAy7a5s1OkhP3EV4paPUAvqpMkPRyMIynY/ETLgPstoJOddp1GlOyozU+OtPyVCfa6CSGnYAnY1sAtnSluQUqYpAchoYmDX4hpBfJmKbvWpONDZMSSciZm2oKHmLyfC7NUG8IvT5nw06IL4BCSyc1RVx9c2llclDm0dISBfGmqYgpAORRVdHpf7djXImAqwp1qt+WjzwlGci62uVpC9NXuWIM/R22RRSX3eo7GmNlKlWFq1hV3bT17NCo6xWo+xWTmuZgp4oe+C2nCPNeAIgF9VMdV9Oq4mxq/EsLzDXQss5VKQqCy2Bm77w2rmsQPuMVt99V7NEk4cnrS8aj686byalE/+DKtb6u9WI6+nlNok2V7M7RCEOdXkecrvLxmbr5IlVpcBZdtxNaVOYf7ockGVwu1q0IQCCD7QVXCD42pwNZQDGZzOUa5QiyDqixDKsgtQa6cgyBG8Z1cXt6shyGWVvUtJNrBr9lWpzzgsAPHA/1VOd0UfFt2EuR/tUFe8OuC2SFhIt+BHkKYFtlrhfPXsPhLypNq6Wyb+lh8RxLOKe7o0qPs0Wm7du1bnrVzrGVNXGF4/ywenhHJeLZMMgP/VWHD1RkYQREVT28xqQv5wXk0By+er83Z3bFGyDFI81J6g/bgaKaZRXK7y8YprCIL8ooM8aYImLKcgcnYTfxp5Dze/C1JHuAweSvxTxNinssGj6dVZmHIVpJ6izzPbQD3xgtP8xAc14yrBXLvBEYJS7x6kWfdlVz4dEpRGd9P+KOAuyJPL0soSYTKgRxDsIsMNlzHy/WpRbYb4zGZK9R9XayqMw9N9uTaNEZTUVV8xgLLdeTXyAmBW1xygXRs14q1qDzza5ADHUMWNeM6L8l47R/AWEWmfnuZa/gi1quEDtdXK5/UOa4E7RVuxj0OeWQ2gVbFDnAx3zV9UpA6dW3wcbu7tNBqOBcmFqU2mci0g7EdwnV18+Lo7J2EUsKrUol/tyWv7GUG4vUY5kx7h3IfRzBH4TeE6rp6/uJHXbLK5JxvVz6wGv8oImlwlfLEajBW/DTv6z1CjaMr4NIqbOvziRmSiKmUubaf2cbUi9G/5KTDuJyKeiabVxsbtyuGPaooV9L/w6FoamZg9OphOnaTU8iGH8guYmSH7BWdlFSjEJP1Ubf7jYlV6s3nYXN+8qLeg4OO4bohPKiq1fsiBmvTJ1Ej8V00jhyTaonseWT7oI7ChDkeYqJb8TdMIosbUFHx1IXL4uBnlDl8oq6d87WO6n1pVaKw5J79pH97o7eTq+drwQVAnPGvDgcrsM9A0PzRKtTlLvGKIZq8Avk0Y5aH2ka7hOjE0QmxwiJFe3gUVtUBdWgsHweIHxyKjwAMfMpduFBtBUV2auka8RO8H8KsO0VRl+mO6X62Y9zC9NVPz+SOnZVWELVEVdHPYW108/TFD9CkLnyFGOmGW+kYJIcUJVVeqjhyegtmepnH0R6CBv07ZnLQwx9d3OWrRlBWv6oCeF8OEHElcpDc7vHtR9byYRo4kht0VqR9yULNmKStT3QMtBZJ6yitnFXcGmlOvTaPEMqHK321PqqZIb8tLqPg/buYVO1G0cO7r/LugX4cgdZh89A/FKp9isrKuX/hygzStEoE3IT0JzXn/KYjjQgOWs941pIEuie+BkBDDhyCUo32NaYc3sJ4mCyKp3h6bfR7INLdrBNvMd3cACerf2V7ecO3TZBrOwC7nCOjaoaOeJ+g5/mLpyPkwWSdYAVuszZ67oJR46MoeGoKeYdYHc9xIc77uGnejEfWiD+tmqllPQbB2wr4sdLzSInUY0GCdtIeCr3xwWSgyVhXTW2nhajQTryyrrxJDOHXtEvLxZq3V9efcjW8af6w+vnDkdCIhkdMCUZrE4i7Hae8pP/2O+OGMGoH19SG8cqTOq4IG1V90TjyjnYIc4NNr6zx53F2Q0yKe7bWF0Q35yxLkNAEPou9hWIO85InfmjT6VZK4oDe/7ZNmME8xPe6zz6gFnaQ9J+tRflGaEgaes0cswT78nDuSIrykzFn2VQ7h7ffM55ywLIiDbVO2D6Z5PQ/lsOU9hLL64VR4P6hEiilz2+6SBJDDZdyeipk/9/ucRDetSObBxuL1PHgA2d0eGplPOeCUakGnzqlKw9fFxSmHXFL33FCTwy1HQUiw+T5/E6WLL9E78mlfJXrmmK7Fx1xes0mKGt08i/dhjEvsgpOYteyxak1jRPeffAgwUslHl5WVcHFtW7Zb1KkVq/Zc283sMYFQBJg6SxoDRy6HO/v6jkvjZrhwHxWubkw7pvPEg9cgCOXIXKm9TTJmjxl61J5RuC+pHRNsgJ8TWa0DQuGBlwVR5ETNFVb3jEywk63/Ifi12YJeRq3m/NXwUZste1mAmxyjR6erhJ6yxXmRenXopNHs3l7AwMfIiLQApJP3Ys+e87DvcvCM9KxTUNfGsuWgMRgGxbfvH+OyL3Yf5lIfpO39jlA+xDjtLKnrOjz+jj+d5R/7JjOrzcWuqI6thsA9hl8nvFMMoFL2/gsueVqraJoZxEQ6Nr1yvikHR7a9jiNs3nIEg1qOSuMxFM5X79H6WHr2g0gko9D/4zgPJiIBaz2EjBdutEirqGrrPR0V9a+OaAeao/YgfXctKN3LaeZqbEPA2p9y7VgL9NpZ0JZJjWNogwGux9GiU332uYrGOIfNtX+mBsagg/EVw6Yc9TiPBbPZTplLgQSeVguD6zzH4kue/uMvY74VgYtn4a2e42pOSQINbazv19dZ7JySBWrpY/8uhONS+LGpX4elxkrgHXcgiOSlvRAY0xEQUZyhhr4XqBZifus0j69qm6XuDUWY78dsEGPje72/Rmi4W8LRai2xd2xMVjvuRQigfeJy2eCtcJ1yFOZ7g1NrXVuOttNgR9RIo1cz3g94laP6ze2VUjLpcS8oHuyxwn+kn6wl5quc3u6yB8CsHvdyFL6KEpuglG6ERYu+8PJn6bYoze2ErPCEusLNR/CNBfFXOdDy5PcSsI9nP7/B7dWmj+PONroH73ISjGdvJWP5g4g1AWPV8u5o07V6tTsJrtqzJq2i8iAsGFEHAU+n8t5GecvVVE/sve16jmCKxkWxPM1Vu8UFcrJO8yyS13bWviWrr07ROkfUzq5+SCEhLrI5+TkP7vvVOFcYxVpOt9cjTLUvYLnaMJ2/bTkaWsW+INLTkUr8q5yckz2r/+pIPolLE3NOondsRqbIb9syC2FMi6LPqwg1ne+hNM3IQ5nLC+112XbKCRkK+ezcaIJ8yNGSXVaV2dth72XmlEMMh+ddDQLtnLsmLROgMT/2mrWXd5NDEKrnu/rBxZ+DYOCxt01WIzy8baVNOcp7+2UWPPdjnwCCYGIsowL94s5ysbTfrknx45UoVRtQxTrroPq2RTrk4Bqx7dd9APzzXqZhpIw29oPsbcY15cSS4379COJ33IuUhN/ApJ8XOe7nEe853oeq8bhX0Kseaj5MQpPfNppn+kbJfr+gZe1z3YLIHTHWn9NudamoHzrcb4zh/ee9wJTcB64dvvGAqz8Dzd7tV9icP5uBQW3Oov3Y1TO4wqqgJz6v1EWee9xLG/y4jcCkx6m/ETq9Y1ChNPsVv+bO5Q81K2XFqb9yh58+LNrvHLaPvQbU0ymLrTcSO1H7f9HejnQ=";
ListLinePlot[{Uncompress[dataA],Uncompress[dataB]},PlotRange->{0,10},PlotTheme->{"VibrantColor","Frame"}]

enter image description here

I need 5 curves in between.

Solution of Federico is very nice, however it takes 52s to compute InverseCDFon 61 point. I have at least 200 points and many function-pairs. Therefore, speed is an issue. I still have to see how the solution of Carl Woll performs.

Context

I need 9 min to generate 1 curve by doing calculations on 24-threads. My hope is to generate intermediate curves by morphing at least an order of magnitude faster then it takes to generate the original ones.

$\endgroup$
1
  • $\begingroup$ Why not give "two interpolation functions" as an example? Then you might avoid getting answers that focus on the peculiarities of UnitTriangle and so forth. $\endgroup$
    – Michael E2
    Nov 5, 2019 at 16:49

2 Answers 2

10
+100
$\begingroup$

Numeric solution

In this example I compute the $W_2$ geodesic (Wasserstein distance) between two densities defined as InterpolatingFunction.

(* unnormalized density functions *)
uf = Interpolation[{{-2, .5}, {0, 2}, {.5, 1}, {1, .5}}];
ug = Interpolation[{{-1, 1}, {0, .5}, {1, 2}, {2, .5}}];

(* normalized density functions *)
f[x_] = uf[x]/NIntegrate[uf[x], {x, -2, 1}];
g[x_] = ug[x]/NIntegrate[ug[x], {x, -1, 2}];
ℱ = ProbabilityDistribution[f[x], {x, -2, 1}];
\[ScriptCapitalG] = ProbabilityDistribution[g[x], {x, -1, 2}];

Show[
 Plot[f[x], {x, -2, 1}, PlotStyle -> Blue, Filling -> 0],
 Plot[g[x], {x, -1, 2}, PlotStyle -> Red, Filling -> 0],
 PlotRange -> {All, {0, All}}, AxesOrigin -> {0, 0}]

pdf

The points xF are a linear sampling of the domain of f. The points qF are the quantiles associated to the points xG. The points xℱ are the union of the two, in order to ensure that both densities are discretized sufficiently well.

xF = Range[-2, 1, .05];
xG = Range[-1, 2, .05];
qF = InverseCDF[ℱ, CDF[\[ScriptCapitalG], xG]];
qG = InverseCDF[\[ScriptCapitalG], CDF[ℱ, xF]];
xℱ = Union[xF, qF];
x\[ScriptCapitalG] = Union[xG, qG];

X[t] is the interpolation between the starting and final points, whereas dens[t] is the intermediate density at those points.

X[t_] := (1 - t) xℱ + t x\[ScriptCapitalG]
dens[t_] := 1/((1 - t)/f /@ xℱ + t/g /@ x\[ScriptCapitalG])

The resulting density can be visualized as

ListLinePlot[Evaluate@Table[{X[t], dens[t]}\[Transpose], {t, 0, 1, .1}]]

interpolation

The transport map can also be computed and plotted with

dT = f /@ xℱ/g /@ x\[ScriptCapitalG];
T = Interpolation[{{xℱ}\[Transpose], x\[ScriptCapitalG], dT}\[Transpose]];
Plot[T[x], {x, xℱ[[1]], xℱ[[-1]]}]

transport map

Symbolic solution

Mathematica appears to be able to deal with distributions, CDF, inverse CDF and pushforwards of distributions:

    ℱ = UniformDistribution[-1 + {-1, 1}/2];
    \[ScriptCapitalG] = TriangularDistribution[1 + {-1, 1}];
    T[x_] = InverseCDF[\[ScriptCapitalG], CDF[ℱ, x]] // Simplify;
    \[ScriptCapitalD][t_] := TransformedDistribution[(1 - t) x + t T[x], x \[Distributed] ℱ]
    Plot[{PDF[ℱ, x], PDF[\[ScriptCapitalG], x]}, {x, -2, 2}]
    Plot[Evaluate@Table[PDF[\[ScriptCapitalD][t], x], {t, 0., 1., .1}], {x, -2, 3}]

morphing between two densities

Symbolic integration of UnitBox and UnitTriangle

While it's true that

    Integrate[UnitBox[y], {y, -∞, x}]

and

    Integrate[UnitTriangle[y], {y, -∞, x}]

do not work as intended, giving a slight hint regarding the domain of x helps in both cases

    Integrate[UnitBox[y], {y, -∞, x}, Assumptions -> x ∈ Reals]
    Integrate[UnitTriangle[y], {y, -∞, x}, Assumptions -> x ∈ Reals]

and the returned result are piecewise functions. An antiderivative can also be found with

    Derivative[-1][UnitBox][x]
    Derivative[-1][UnitTriangle][x]
$\endgroup$
4
  • 1
    $\begingroup$ Thank you for providing a reference solution. However, it is explicitly asked for a numeric approach. $\endgroup$
    – yarchik
    Nov 5, 2019 at 15:34
  • $\begingroup$ @yarchik You are right, I modified the answer to include a numeric approach. $\endgroup$
    – Federico
    Nov 5, 2019 at 17:33
  • $\begingroup$ Thanks a lot! Numerical solution works fine, but it is rather slow: 52s to compute InverseCDF for 61 point. This seems to be a well known problem of this function as shown here and here. In fact, efficiency is my biggest concern. I have 40 function pairs and some of them with sharp peaks that require many points. I hope someone can come up with a compiled solution. $\endgroup$
    – yarchik
    Nov 5, 2019 at 22:13
  • $\begingroup$ This method is very lovely! Do you happen to have a good reference discussing this and related topics in detail? $\endgroup$ Nov 16, 2019 at 0:23
7
$\begingroup$

You can use NDSolveValue to create an interpolating function representation of the inverse. Basically, suppose you want to invert f. Then:

f[finv[x]] == x

where finv is the inverse function. So, an ODE for the inverse function is:

D[f[finv[x]] == x, x]

f'[finv[x]] finv'[x] == 1

Let's use this for your G function:

g[x_] := UnitTriangle[x-3]
G[x_] := Integrate[g[s], {s, -Infinity, x}]

Then we have:

Ginv = NDSolveValue[{G'[inv[x]] inv'[x] == 1, inv[G[3]] == 3}, inv, {x, 0, 1}]

However, it's easy to see that we can use g instead of G', so it will be quicker to do:

Ginv = Quiet @ NDSolveValue[{g[inv[x]] inv'[x] == 1, inv[G[3]] == 3}, inv, {x, 0, 1}];

The quieted messages are associated with the the fact that g is zero when x is at one of the endpoints, 0 or 1. Let's check:

G[Ginv[0]]
G[Ginv[.5]]
G[Ginv[.75]]
G[Ginv[1]]

0.

0.5

0.75

1.

So, Ginv is an interpolating function representation of the inverse of G, and you can take derivatives of it as desired, e.g.:

D[Ginv[Sin[x]], x] /. x->3

-1.86349

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.