Skip to main content

equation solving - Convert symbolic to numeric code: speed up morphing



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.



Answer





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




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]

Comments

Popular posts from this blog

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...

plotting - How to draw lines between specified dots on ListPlot?

I would like to create a plot where I have unconnected dots and some connected. So far, I have figured out how to draw the dots. My code is the following: ListPlot[{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4,13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full] I have thought using ListLinePlot command, but I don't know how to specify to the command to draw only selected lines between the dots. Do have any suggestions/hints on how to do that? Thank you. Answer One possibility would be to use Epilog with Line : ListPlot[ {{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4, 13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full, Epilog -> { Line[ ...