From 6239dee8b39dfa09c7199d8288d64b61940991f2 Mon Sep 17 00:00:00 2001 From: Laura Messio Date: Tue, 30 Jul 2013 21:46:03 +0200 Subject: [PATCH] test/triqs/mc_tools directory added, with an example of simple random walker --- different_moves.ipynb | 63 ++++++++++++ test/triqs/mc_tools/CMakeLists.txt | 3 + test/triqs/mc_tools/different_moves_mc.cpp | 110 +++++++++++++++++++++ 3 files changed, 176 insertions(+) create mode 100644 different_moves.ipynb create mode 100644 test/triqs/mc_tools/CMakeLists.txt create mode 100644 test/triqs/mc_tools/different_moves_mc.cpp diff --git a/different_moves.ipynb b/different_moves.ipynb new file mode 100644 index 00000000..02de6a0e --- /dev/null +++ b/different_moves.ipynb @@ -0,0 +1,63 @@ +{ + "metadata": { + "name": "different_moves.ipynb" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from pytriqs.archive import HDFArchive\n", + "R = HDFArchive('histo.h5', 'r') # Opens the file myfile.h5 in readonly mode \n", + "H = R['H'] \n", + "P = HDFArchive('params.h5', 'r') # Opens the file myfile.h5 in readonly mode \n", + "pl = P['pl'].real[0]\n", + "pr = P['pr'].real[0]\n", + "xmax = P['xmax'].real[0]\n", + "N_Cycles = P['N_Cycles'].real[0]\n", + "Length_Cycle = P['Length_Cycle'].real[0]\n", + "\n", + "temp=(pr+pl)/2\n", + "pr/=temp\n", + "pl/=temp\n", + "sigma=sqrt(Length_Cycle*min(pr,pl))\n", + "\n", + "plot(H, 'b')\n", + "x=[i for i in range(100)]\n", + "y=[ 1/ (sqrt(2*pi)*sigma) * exp( - (x[i]-xmax)**2 / (2*sigma**2) ) for i in range(100)]\n", + "plot(x, y, 'g')" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "pyout", + "prompt_number": 8, + "text": [ + "[]" + ] + }, + { + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD9CAYAAAC7iRw+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVPX+x/HXsKiIJuJVTMBQQQEXRFHUFgcVKRcyrcRb\nZmbpzzKrW7dbv7uk3V+Wv1ZvWldbvKml3BY3RBSUUUuBClwSSkjIAQVXUFAChvP74/6icBlRBr4z\ncz7Px8PHw4HvnPM+p5m3p++cOcegaZqGEEII3XBRHUAIIUTzkuIXQgidkeIXQgidkeIXQgidkeIX\nQgidkeIXQgiduWrxJyUlERwcTFBQEAsXLrzsmLlz5xIUFERYWBhZWVl1Pw8ICKBfv36Eh4czePBg\n26UWQghx3dys/dJisTBnzhxSUlLw9fVl0KBBxMbGEhISUjcmMTGRvLw8cnNzSU9PZ/bs2aSlpQFg\nMBgwmUx4e3s37VYIIYRoMKtH/BkZGQQGBhIQEIC7uztxcXGsX7++3pgNGzYwbdo0ACIjIyktLaWk\npKTu9/L9MCGEsC9Wj/iLiorw9/eve+zn50d6evpVxxQVFeHj44PBYGDUqFG4uroya9YsHnnkkXrP\nNRgMttgGIYTQncYcVFs94m9oMV8pwJdffklWVhabN29myZIl7Nq167LPlT8aL7zwgvIM9vJH9oXs\nC9kX1v80ltXi9/X1xWw21z02m834+flZHVNYWIivry8AXbp0AaBjx47cddddZGRkNDqwEEKIxrFa\n/BEREeTm5lJQUEBVVRXx8fHExsbWGxMbG8uKFSsASEtLw8vLCx8fH86fP8+5c+cAqKioYOvWrfTt\n27eJNkMIIURDWZ3jd3NzY/HixcTExGCxWJgxYwYhISEsXboUgFmzZjFmzBgSExMJDAzE09OT5cuX\nA1BcXMzEiRMBqKmp4b777mP06NFNvDmOy2g0qo5gN2Rf/Er2xa9kX9iOQbPFhNH1rtxgsMl8lRBC\n6Elju1O+uSuEEDojxS+EEDojxS+EEDojxS+EEDojxS+EEDojxS+EEDojxS+EEDojxS+EEDojxS+E\nEDojxS+EEDojxS+EEDojxS+EEDojxS+EEDojxS/szpEj8P9X9xZCNAEpfmF35s+Hhx8GuWGbEE1D\nil/YlcJCWLsW3noLZs6EmhrViYRwPlL8wq68+SY8+CDMmQOdOsGiRaoTCeF85A5cwm6cPg2BgbB/\nP/j5QV4eDBkC334LN930nzFlZdC6Nbi7q80qhEqN7U4pfmE3XnxR45uj33LTnR+RdzqPyb0nczjh\nbr5KbUNkJKSk/Ocfgfnz4S9/UZ1WCHXk1ovCKXz07RpePNWHvYGT+Z3H73ig3wN8kfMF/3Dzo3Dg\ndGpcz/Laa5CYCBs2qE4rhGOTI36h3GfZn/HwZ0/QL/cTdqy4DYPBUPe74xXHeX7b8+SfySfxvkRc\ntVZ06gQ5OdC5s8LQQigkR/zCoW3P387M9Y/ismYT7z43vF7pA3Ty7MSyccvo0LoD931xHy6uFkaP\n/s+RvxDi+kjxC2Uyj2VyT3wchs/+zcev96d378uPc3VxZdVdqyitLOWxxMcYO1YjIaF5swrhTKT4\nhRKllaWM/Xg8rVL+yf88bOSOO6yPb+nWkrWT15JRlEGR32K2bYPKyubJKoSzkeIXSry08yVc8sYw\npf9EZs9u2HNuaHkDK+9ayRvfvkhw+Gl27GjajEI4K/lwVzS7w2cOM+DdQbRffZAf93XG5RoPP2Zv\nms2BrFaEl7zJ2283TUYh7Jl8uCscznMpz9H73FNMvevaSx9gvnE+B11XsnZHLnLcIMS1k+IXzWq3\neTd7CveQ/8kfiIu7vmV08uzEs7c8w+mBfyI727b5hNADKX7RbDRN4w9b/sADXRbQ0as1oaHXv6yn\nhj6Jm38mb2/YabuAQuiEFL9oNuu+X0d1bTUnU+9jypTGLauVWyse6/UyH5c8b5twQuiIFL9oNm9n\nvM0fIv/IF5+7XPc0z2/Nu+dezrsV8sWezMYvTAgdkeIXzSLnRA45J3NoY55IUBAEBDR+mS1buHJz\ni/9ifuI7jV+YEDoixS+axbvfvMvDAx7ms/gWjZ7m+a2X7p7BAcvnnCw/Y7uFCuHkpPhFkyuvKmfV\n/lVMDZnJxo1w7722W/atAzrhfXIs//3pv2y3UCGcnBS/aHIf7/+YW7sO560X/RkyBHx8bLv86X0e\nZXXuO9RqtbZdsBBO6qrFn5SURHBwMEFBQSxcuPCyY+bOnUtQUBBhYWFkZWXV+53FYiE8PJzx48fb\nJrFwKJqmsWj3O+R+/BiFhbB6te3X8fz9Qzlf5snafSm2X7gQTshq8VssFubMmUNSUhLZ2dmsXr2a\nnJycemMSExPJy8sjNzeXZcuWMfuiC68sWrSI0NDQSy63K/Rh8frdHMqvZNptI1i3Dtq3t/06vL0N\nDLA8Kh/yCtFAVos/IyODwMBAAgICcHd3Jy4ujvXr19cbs2HDBqZNmwZAZGQkpaWllJSUAFBYWEhi\nYiIPP/ywXJNHp+Zv+if395rNn551ua7LMzTUX++8j+yKXRSeLWy6lQjhJNys/bKoqAh/f/+6x35+\nfqSnp191TFFRET4+Pjz11FO8+uqrnD179orrmDdvXt3fjUYjRqPxGjdB2KtjJy9w6ncbefHe15p8\nXWNHe9JyVSyLUj7l1YlPNfn6hGhOJpMJk8lks+VZLf6GTs9cfDSvaRoJCQl06tSJ8PBwq4F/W/zC\nuby1KRHvygi6etv409zLcHWFcQGTWZU1T4pfOJ2LD4rnz5/fqOVZ/Z9vX19fzGZz3WOz2Yyfn5/V\nMYWFhfj6+rJ79242bNhAt27dmDJlCtu3b+eBBx5oVFjhWL74IZ7hv5vcbOv76/0jKan+kbyTBc22\nTiEckdXij4iIIDc3l4KCAqqqqoiPjyc2NrbemNjYWFasWAFAWloaXl5edO7cmQULFmA2m8nPz2fN\nmjWMGDGibpxwfuVV5Rw2bOFR48RmW2efEHc6nryLBev/3WzrFMIRWZ3qcXNzY/HixcTExGCxWJgx\nYwYhISEsXboUgFmzZjFmzBgSExMJDAzE09OT5cuXX3ZZclaPvqzKSMDl6DCihnRo1vVO6TuZj3Kf\n5UOebdb1CuFI5A5coklEvDYBLWcC337wYLOu98SpGnwW+vL17K8Y2C2wWdctRHORO3AJu1NWWcb+\nc9u5t9+EZl93xw5udP/5bl78PL7Z1y2Eo5DiFza3/of1uBUaGR/tpWT9jwydTEqxFL8QVyLFL2xu\n+dfxtDg0mZAQNet/atItVBpOkZiRc/XBQuiQFL+wqbLKMtKO7mL0TbGo+jy/hbsLYe73sDDhMzUB\nhLBzUvzCprb+uJX2527h9hFtleaYZYwl/XQCcu6AEJeS4hc2tfFQAue+HcfIkWpzPDjiFqrbHWKT\nqURtECHskBS/sBlLrYWE7zfjUzaWm25Sm6WlWwtCW0bzRkKi2iBC2CEpfmEzXx/9GrdKH+6NUdz6\n/2/asLF8dXwT1dWqkwhhX6T4hc0kHErA8v04JjT/6fuX9cDQO7DclMKmpCrVUYSwK1L8wmY+P5CA\nS944IiJUJ/mPTp6d8PcI4R/rd6qOIoRdkeIXNmEuM3OktJB7hg5p0huuXKspA8fx1YkEystVJxHC\nftjRW1Q4ssTcRDyKbmfiBFfVUeq5t/9Y3EITWLtWzusU4hdS/MImPjuQQOX+cQwfrjpJfWE+YbRq\nU8l7aw+pjiKE3ZDiF412vvo8X5p3MDY4Bnd31WnqMxgM3BkyjvQzCZw/rzqNEPZBil80mqnARKvS\ncCbHtlcd5bIm9hmHR78Edu1SnUQI+yDFLxptw8EtXNh/OzExqpNcXlRAFBe8vyFhq3zCKwRI8YtG\n0jT4NDOZYZ2j8fRUnebyPFt40tc7go375bROIUCKXzRCdTXcM6OQs5bjrH4jXHUcqyb0i6bYM5mi\nItVJhFBPil9cl3PnYPx4yKtNYXyfEfh0sq/TOC8WExhNy5BkkpNVJxFCPSl+cV1mzoQbb4SQcVu5\no2e06jhXNeDGAVhaH2PddjnkF0KKX1yXHTvgr3+rZVt+CtE97L/4XV1cGd51BNvyU6itVZ1GCLWk\n+MU1Ky6Gykooa7Ufr1ZeBHgFqI7UIHf2icYlMJm9e1UnEUItKX5xzbKyIDwcUg4nO8TR/i+iu0dT\n0zWFLVvk8g1C36T4xTX7pfiTDycT3d1xir9b+26082jD2t0HVEcRQikpfnHNsrKgd1glewr3EBUQ\npTrONRkTHM2+8mS5WqfQNSl+cc0yM6HW70v6dupLu1btVMe5JmN7jaZ1n2QyMlQnEUIdKX5xTcrK\noKQEfqh2rPn9X0R1i+J8h93szqhUHUUIZaT4xTXZuxf69oVtBSmM6jZKdZxr5tXKi64eoWzN2aM6\nihDKSPGLa5KVBaEDznDo1CEi/SJVx7kuI7uPYF9ZquoYQigjxS+uSVYWtOq1k6F+Q2nh2kJ1nOsy\ncUAU5zulUlKiOokQakjxi2uSmQmn2213uLN5fuuWrjejdc5iV1qF6ihCKCHFLxrswgXIy4PvKlIZ\n0W2E6jjXrbV7a7oYwlmX+ZXqKEIoIcUvGuy776B7nxMcOfsTA7sMVB2nUYbeOII9xdtVxxBCCSl+\n0WBZWdBxkIlbu96Km4ub6jiNck9EFEdcUtHk6g1Ch6T4RYNlZYGla6pDz+//Ylz/IVg6HGT/oTLV\nUYRodlL8osEyM8HsmkpUN8cv/lZurehQGcknX8od2IX+XLX4k5KSCA4OJigoiIULF152zNy5cwkK\nCiIsLIysrCwAKisriYyMpH///oSGhvL888/bNrloVjU1cCD/GGWWEsJ8wlTHsYn+7aLYfljO5xf6\nY7X4LRYLc+bMISkpiezsbFavXk1OTk69MYmJieTl5ZGbm8uyZcuYPXs2AK1atSI1NZW9e/eyf/9+\nUlNT+fLLL5tuS0STysuDtv1SMXYbjquLfd9msaHGhETx/c/yAa/QH6vFn5GRQWBgIAEBAbi7uxMX\nF8f69evrjdmwYQPTpk0DIDIyktLSUkr+/5sxrVu3BqCqqgqLxYK3t3dTbINoBgcOgEeIc8zv/2LK\n8EGUt8zjRPlp1VGEaFZWT80oKirC39+/7rGfnx/p6elXHVNYWIiPjw8Wi4WBAwfy448/Mnv2bEJD\nQy9Zx7x58+r+bjQaMRqN17kpoikdOABnO6QyotsTqqPYTOeOLfA4cTOf7N7BE6PvUh1HiCsymUyY\nTCabLc9q8RsMhgYtRLvonLhfnufq6srevXspKysjJiYGk8l0SbH/tviF/Ur/3kxNvzJ6d+ytOopN\n9WwRxcbvtkvxC7t28UHx/PnzG7U8q1M9vr6+mM3musdmsxk/Pz+rYwoLC/H19a03pl27dowdO5Zv\nvvmmUWGFOnvPmBh2o7HBBwOOYkR3I5mndqiOIUSzslr8ERER5ObmUlBQQFVVFfHx8cTGxtYbExsb\ny4oVKwBIS0vDy8sLHx8fTp48SWlpKQAXLlwgOTmZ8PDwJtoM0ZQqKuBU2x2M7W1UHcXmHogeQCkF\nnDp/SnUUIZqN1akeNzc3Fi9eTExMDBaLhRkzZhASEsLSpUsBmDVrFmPGjCExMZHAwEA8PT1Zvnw5\nAMeOHWPatGnU1tZSW1vL1KlTGTlyZNNvkbC57Gxw7WEiqvuTqqPYXFgfd9yLh/Fpxk7+yyjTPUIf\nDNrFE/TNuXKD4ZLPB4T9eW1ZIX8pDOf8/BJcDM73nb/wx1/Bu+sxtv1xkeooQjRIY7vT+d7Fwua2\n5e0gyH24U5Y+wPi+Rr6VeX6hI875ThY2tf+siVv8h6uO0WSm3z6Qs66HOVkh5/MLfZDiF1dV4mFi\n0gCj6hhNpltXd1qfGsqKHTtVRxGiWUjxC6v2HS6ituUZovo41/n7F+t7g5G1WSbVMYRoFlL8wqr4\n9B14nxuOq4tzv1Qm9Deyr1Tm+YU+OPe7WTSaqcBEiIfzzu//4qGYCM61yOPoGZnnF85Pil9YdbDC\nxPCbjKpjNLmOHdxpd3YoH6bI9fmF85PiF1d09NxRzmunGN2/j+oozSK8vZGNB2S6Rzg/KX5xRab8\nHWgFt9Gvrz5eJhMHDue7crkxi3B++nhHi+uy8TsTbU8Z8fJSnaR5TBs1iPMeefx0/IzqKEI0KSl+\ncUW7juygTxuj6hjN5gbPFrSvGMKHKXKnOOHcpPjFZR07d4xTlScY0aev6ijNaoC3kcSDJtUxhGhS\nUvzisnb8tIOWxbcxOlpfL5GJ4Uayz5tUxxCiSenrXS0abHOOiZ9/GM7gwaqTNK+pIwdxvvUhjhwv\nVR1FiCYjxS8uKznXxJDORtzdVSdpXm1bt6B9RSTvb5Xz+YXzkuIXlyguL+ZU5XEm3txPdRQlBnYw\nkpgt5/ML5yXFLy5hKtiBi/k2Ykbr8+UxcYCRnAsm1TGEaDL6fGcLq9btNeFRMpygINVJ1Lg/ahDn\nW39PQbHM8wvnJMUvLmHKNzE8wIjBoDqJGm1bt8T7fCQfJMv5/MI5SfGLeorLizlVVczk4fqc3//F\nAG8jidkm1TGEaBJS/KKe7T/uhILbiB7lqjqKUhMHGsm5IB/wCuckxS/q+XeGic4/D6dDB9VJ1Joa\nNZgLnjkUHCtTHUUIm5PiF/V8WZRKdJBRdQzl2ni0pEPlYJYkyPn8wvlI8Ys6x84do7S6hN+PDFMd\nxS7c6jeC9fvkMs3C+Ujxizpr96bicmQ4w2/V9/z+Lx4eFcWP2nbOn1edRAjbkuIXddakp9KnTZTu\nLtNwJdGhgzB45/HZJrkPr3AuUvyiTubpVO6JGKE6ht1o4dqCnh4380GKnN0jnIsUvwAg9/gRzlvO\n8nBsb9VR7MrE8CjSjm+nqkp1EiFsR4pfALAsORXvs0Y6dtTp13Wv4M5+I3Dtkcq2baqTCGE7UvwC\ngISD27m5i0zzXCz8xnC0toV8vK5EdRQhbEaKX1Bbq5FbncqMEVGqo9gdNxc3bva7jQ0HTNTUqE4j\nhG1I8QuSvz2MZqhh3NCeqqPYpbGhUbTslcou+S6XcBJS/IL3UrbTw2UELi4yv385Ud2i0AK2k5Cg\nOokQtiHFL9h5JJUxITLNcyX9fPpR5XaKtOxC1VGEsAkpfp07c0bjZNvtzBwtxX8lLgYXbvUzsv9c\nKpqmOo0QjSfFr3Mfb82hpWsrQjp3Ux3Frt0RMoIq3+2UyMk9wglctfiTkpIIDg4mKCiIhQsXXnbM\n3LlzCQoKIiwsjKysLADMZjNRUVH07t2bPn368I9//MO2yYVNfJ6VQm+PURj0erutBoruPgpD9xT2\n75dDfuH4rBa/xWJhzpw5JCUlkZ2dzerVq8nJyak3JjExkby8PHJzc1m2bBmzZ88GwN3dnTfffJOD\nBw+SlpbGkiVLLnmuUC+zLJk7+0SrjmH3enboibu7gW37flAdRYhGs1r8GRkZBAYGEhAQgLu7O3Fx\ncaxfv77emA0bNjBt2jQAIiMjKS0tpaSkhM6dO9O/f38A2rRpQ0hICEePHm2izRDX4/jJas6238mM\nESNVR7F7BoOBfm1GscOcojqKEI3mZu2XRUVF+Pv71z328/MjPT39qmMKCwvx8fGp+1lBQQFZWVlE\nRkZeso558+bV/d1oNGI0Gq91G8R1Wp6cTtvqHnTx+p3qKA4hJiiaN/LXAHNURxE6YzKZMJlMNlue\n1eJv6LyvdtGpDr99Xnl5OXfffTeLFi2iTZs2lzz3t8Uvmtf6Ayn0byvTPA11/7CRvJAxm5+ra2jp\nbvWtI4RNXXxQPH/+/EYtz+pUj6+vL2azue6x2WzGz8/P6pjCwkJ8fX0BqK6uZtKkSdx///1MmDCh\nUUGF7e2vSGZi/1GqYziM7j6daHE+gLUZX6uOIkSjWC3+iIgIcnNzKSgooKqqivj4eGJjY+uNiY2N\nZcWKFQCkpaXh5eWFj48PmqYxY8YMQkNDefLJJ5tuC8R1OVJyloo2+5k+8hbVURyKf/Uo1u5LVh1D\niEaxWvxubm4sXryYmJgYQkNDmTx5MiEhISxdupSlS5cCMGbMGLp3705gYCCzZs3inXfeAeCrr75i\n1apVpKamEh4eTnh4OElJSU2/RaJBlm010b4iknaeHqqjOJRBHaJJOy4f8ArHZtAunqBvzpUbDJd8\nPiCax4A/z6VjS1+2/O1PqqM4lBVrKngo24czfz5G25ZtVccROtXY7pRv7upUzs/J3Bsh8/vXalCY\nJy1ODGLnTztVRxHiuknx69BBcyGVrif4/Yhw1VEcTlAQVP8QTeIPMs8vHJcUvw69ty2FThUj8Ggl\n//mvlZsbdNNGsfmQFL9wXPLO16HEQ0nc3DlGdQyHNdh/ICfOl2AuM199sBB2SIpfZ6otNfzIVh6/\n43bVURxW/36u+FaOZnPeZtVRhLguUvw6szI1HbeKrhgH+KqO4rD69gXXw2Ok+IXDkuLXmeW7NtO/\nzR3IVZivX//+ULQzhtT8VKosVarjCHHNpPh15tuyzTww9A7VMRyajw9ED+vIDdU9+erIV6rjCHHN\npPh1JP1gMZUeh5kxeqjqKA7v+eeh9Os72Ph9ouooQlwzKX4dWZSwhZtqR9KqhbvqKA5vwAAIbXEH\n8Zkyzy8cjxS/jqT8lMj4EJnmsZUFcwZRXF5MwWk5rVM4Fil+nSg+XsOJG5J5cqwUv61EDXfF+0wM\nCz6Vo37hWKT4dWLx2nTa0ZXuHbuojuI0DAZ4YOgdxGcmItcaFI5Eil8nPt27mVtvlKN9W/vTpBjK\nO6ayfYec1ikchxS/DlRVQZ4hkZlRUvy21qlNRzq5BLNyp1ytUzgOKX4d+GLbEVy8zNzRZ5jqKE5p\nhG8sqUc3qI4hRINJ8evAB1+up7f7ONxc5AbhTWHWbRMwe67DYpGJfuEYpPh1IK10HXHhd6qO4bRu\nDQ7FzdCCz77KUh1FiAaR4ndyOQVnKG/3NY/GjFYdxWkZDAZ61k5gRcZ61VGEaBApfif31qZNdPl5\nBDd4tFYdxamN6T6BPafXqY4hRINI8Tu5xMPriLlpguoYTu9+41DKao+RfyZfdRQhrkqK34lV/FxJ\nUcsUnhwzTnUUp9cn1BW3w+NZ+bVM9wj7J8XvxJZu3YbH2TD6Bf5OdRSn5+ICfd0n8O/9Mt0j7J8U\nvxP7OHMdg9vKNE9zGRcyitzyLE6eP6k6ihBWSfE7KUuthQM/b2T6MDmNs7kYb/HAs2Qkmw5tUh1F\nCKuk+J3UluzdWMp8uDe6u+ooujF4MFR8cxfx332mOooQVknxO6m3tq0m6Oc4WrVSnUQ/WreGPu53\nsrNgJ2cunFEdR4grkuJ3QjW1New69RlxfSerjqI7tw6+gQBtFF/kfKE6ihBXJMXvhFJ+3E7NyW5M\nHSfTPM0tNhbO7JzC6u/WqI4ixBVJ8Tuhd3etxss8hR49VCfRn6go6FE7lj0F31BSXqI6jhCXJcXv\nZH6u+Zlk83pie9yjOoouGQzwvy95wKFxfLLvU9VxhLgsKX4nszlvMy3O9OOe231VR9GtIUOgr0sc\nb6fKdI+wT1L8TmZl1hoqv5nC8OGqk+jbu89EU1D+PQcLj6iOIsQlpPidSHlVOZvzNnOL9yQ8PFSn\n0bfwfi0IrJ7Ik+/Hq44ixCWk+J3Ixh824l0xjDuj5do89uCFSXGknlhNRYXqJELUd9XiT0pKIjg4\nmKCgIBYuXHjZMXPnziUoKIiwsDCysn69C9FDDz2Ej48Pffv2tV1icUX/2vsvKnbfzx1yT3W7EDdk\nOO7tTvCP+P2qowhRj9Xit1gszJkzh6SkJLKzs1m9ejU5OTn1xiQmJpKXl0dubi7Lli1j9uzZdb+b\nPn06SUlJTZNc1PNT6U+km7/B+/hEAgNVpxEAri6ujOkynX+mf6A6ihD1WC3+jIwMAgMDCQgIwN3d\nnbi4ONavr3+98Q0bNjBt2jQAIiMjKS0tpbi4GIBbb72V9u3bN1F08VvL9y6nV/XvGRcjk/v25KW7\nH8Ls9QmHj1SqjiJEHTdrvywqKsLf37/usZ+fH+np6VcdU1RUROfOnRsUYN68eXV/NxqNGI3GBj1P\n/MpSa+GDzA+58OlGPpAzCO1KcOcAbnTpz98+Wceq5+JUxxEOymQyYTKZbLY8q8VvMBgatBBN067r\neVC/+MX1ST6cjFulD4O7htGnj+o04mIPD3iY11PfQ9PiuIa3hhB1Lj4onj9/fqOWZ3Wqx9fXF7PZ\nXPfYbDbj5+dndUxhYSG+vvLloeb0fub7XNj9MM88ozqJuJzn7pzAhXb72LT7sOooQgBXKf6IiAhy\nc3MpKCigqqqK+Ph4YmNj642JjY1lxYoVAKSlpeHl5YWPj0/TJRb1HK84TtKhFHyOxxEVpTqNuByP\nFi2JaHE//7NpueooQgBXKX43NzcWL15MTEwMoaGhTJ48mZCQEJYuXcrSpUsBGDNmDN27dycwMJBZ\ns2bxzjvv1D1/ypQpDBs2jEOHDuHv78/y5fLCt7WV+1bSpnACzz3VTqYR7Nhfx87g6+rlVFbVqI4i\nBAbt4gn65ly5wXDJ5wOi4Wq1Wrq/HsqF+Pcp3H0L7u6qEwlr2j41hKci/psX74u9+mAhrGhsd8o3\ndx3Y5tzNlJ1szZ+m3Cyl7wAmB8zl3aw3VccQQorfkT234Q3cv/kDM2fKHI8j+N9p93CKPLYeyFQd\nReicFL+DSszcS3bJD6x76V7atFGdRjSEt5c7g2rn8uzaN1RHETonxe+Aqqth2j/f4A7vxxkW2UJ1\nHHENXrnnEQ5UJnKktFB1FKFjUvwO6Mm/FlHmk8BHc2eqjiKukXGIFx0KH+CPn76tOorQMSl+B7Nv\nH3yUs4Rp4ffTwVOug+RoDAZ4fPATrDvyAeVV5arjCJ2S4ncw739UgRb+Hs8Zn1AdRVynuVO7oR2O\nYtHOD1VHETolxe9AamrgXzmLua1rFD28e6iOI65Tu3YQ0/ZZXv3qVSpr5KqdovlJ8TuQjVvLqAx/\nnTdjX1QdRTTSvEcGUfXTQBanv6s6itAhKX4H8sKW1xnQdgzBvwtWHUU00sCB0PfE3/n79lc49/M5\n1XGEzkgbMWN7AAAL70lEQVTxO4iCEyf4zmMJi+9+QXUUYSN/f7wvhvxRvJW2SHUUoTNyrR4HMXbR\nM3z3wwV+emeJ6ijCRjQN+tyWhzlmCAVPH8Lbw1t1JOEg5Fo9OlB0tojkkx/y7JA/q44ibMhggL8+\nFkjrn+7i1d2vqo4jdESO+B3AffEz+GxlB07H/y+enqrTCFuqqYHu/c2U/b4/Bx7Lomu7rqojCQcg\nR/xO7ssjX7Lphy1M8P6LlL4TcnOD5x/zp3PBEzyRJN/NEM1Dit+OVVuqeXjdf6Ftfos/P3OD6jii\niTz4INTu+hO7vs9m/fcbVMcROiDFb8fe2PMmpw535ek7JtGvn+o0oql4eMCXppa0+/Idpn4yl3OV\nFaojCScnc/x2qqC0gL5vR3DT1gyytneXG63oQGkp9Hzufnw8fMl6dSFubqoTCXslc/xOSNM0Znw+\nh9qv/sDqd6T09cLLC9L//jo/eH7IjP/erzqOcGJS/HZoScY7fP39MZ695Rn69lWdRjSnbh19ePP2\n1/ikKo71iTLlI5qGTPXYmaxjWdy6bDRBO/eQkRQoR/s6FfPPaeza4UrBog/p1El1GmFvZKrHiZz7\n+Rx3rroXly1vs+5DKX09+/yhJXj03E300yuRYyNha1L8dkLTNB5aO4vSvVH865k4brpJdSKhUpsW\nbUie+W+yu/6B5177QXUc4WSk+O3Ea7tfJ3nfd8R5vcXEiarTCHswwLcfLw5/iTeOTuD1d0+pjiOc\niMzx24GV+1Yyd/2fuXHzV2Sm+tOqlepEwp488ukfWWH6koWh23jysdaq4wg70NjulOJXLCkvicmr\np9H6s+18vak3fn6qEwl7U6vVcveqB0nacZoF/dby5Fz58Efv5MNdB5ZemM7k+Km4fraW7Wuk9MXl\nuRhciP/9B0QOqeUvGTN57/1a1ZGEg5PiV2RL3hZuXzEe1n7Epn8OIyREdSJhz9xd3Ul44FNChx1m\n7o7f8+kXP6uOJByYFL8CH3z7EZNWTYM1a4n/nzEMHao6kXAEni082TlzCzffYuG+zbezaVup6kjC\nQckcfzOq1WqZE7+A9zM/YFh+Ip8sCqFLF9WphKOx1Fq4+/2n2PhdKp9O2MRdI+Qa/nojH+46AE2D\nz7Yc47GtD3LmQhmvD/6Cxx/sgsGgOplwVJqmMXP5W3zww8tMafc2K5+bjIv8/7tuSPHbuZ074YEF\n6yjs/1/c3nEWH8/8C+3aylkZwjY2fvsN98bfR/uKSHY+t5hAf7lvgx5I8dup2lr408I83s55jnbB\nmXx+3ypuuWmY6ljCCZWdryDqlafZd34TD3Z9iaWP3Y+bqxz+OzMpfjuUV3SK0S/9HXP7VTw99Gle\niHkSD3cP1bGEk/vXtt08vvEZal0u8Oro15gdMwKDzCc6JSl+O3Kg5Due+ORtTCf+TT/DFDb9cR6+\nXnJpRdF8LBaNuf/8nKV5z9NC82JU27n8IeZebh3aEldX1emErcgXuBQrKT/O27uXEbE4iohFo/k+\nw5eE23PY+9I711T6JpOp6UI6GNkXv7rWfeHqamDJY3dTsfB75hlf4IBhJaMSbqLLg8/w3uY0ajXH\n/fKXvC5s56rFn5SURHBwMEFBQSxcuPCyY+bOnUtQUBBhYWFkZWVd03Mdzakz1bzwXhrBMxbi/nAU\nnV/qyVP/2M5Pnz7G6/4FmFf9jTG3db7m5cqL+leyL351vfuiZQtXnp0wjvwXt3LgaRM3D/bk0aQZ\ntP1bV6Z/9hifHvyUkvIS24ZtYvK6sB2rd/W0WCzMmTOHlJQUfH19GTRoELGxsYT85mumiYmJ5OXl\nkZubS3p6OrNnzyYtLa1Bz7VnJSWwfHUZ63bkcablfio891PeZh9lnt9wg6UbQ/oO58/DnuLOPtHc\n0Frm74X9CukYzBePz+fs2fnM/XsOHy9O4OOAj6j1ewQPS2falg+kxekwDMf70cs7lJlTfBk/1lXu\nB+HErBZ/RkYGgYGBBAQEABAXF8f69evrlfeGDRuYNm0aAJGRkZSWllJcXEx+fv5Vn9ucamprOF99\nnhNl5RwpLufoqbOcrDjD6fNnOFN5hmNnSyguL+ZE5TGOVRRyzvUwri1/xm9ID4La9aW7Zz96tIlh\nYuRgenTxVrINQjTGDTfAv14NYbkWwsmTf+RwgYXdeQcwV++lsGYfBReSSD99iK2ZJ3HZ0RWflgHc\n4HIjNxg6086tM/4dOhDo157Qbt5079KOti3b0KbFf/60dG0pHyQ7EKvFX1RUhL+/f91jPz8/0tPT\nrzqmqKiIo0ePXvW5gF2/WGqAAvZTwH7gYwCebcL1zZ8/vwmX7lhkX/xKxb6wkMtRcjna7Gu2Tl4X\ntmG1+Btaytf76bIzndEjhBCOwmrx+/r6Yjab6x6bzWb8Lrp28MVjCgsL8fPzo7q6+qrPFUII0fys\nntUTERFBbm4uBQUFVFVVER8fT2xsbL0xsbGxrFixAoC0tDS8vLzw8fFp0HOFEEI0P6tH/G5ubixe\nvJiYmBgsFgszZswgJCSEpUuXAjBr1izGjBlDYmIigYGBeHp6snz5cqvPFUIIoZimyObNm7VevXpp\ngYGB2iuvvKIqhhJHjhzRjEajFhoaqvXu3VtbtGiRpmmadurUKW3UqFFaUFCQFh0drZ05c0Zx0uZT\nU1Oj9e/fXxs3bpymafrdF2fOnNEmTZqkBQcHayEhIVpaWppu98WCBQu00NBQrU+fPtqUKVO0yspK\n3eyL6dOna506ddL69OlT9zNr275gwQItMDBQ69Wrl7Zly5arLl/JN3d/Occ/KSmJ7OxsVq9eTU5O\njoooSri7u/Pmm29y8OBB0tLSWLJkCTk5ObzyyitER0dz6NAhRo4cySuvvKI6arNZtGgRoaGhdScU\n6HVfPPHEE4wZM4acnBz2799PcHCwLvdFQUEB7733HpmZmRw4cACLxcKaNWt0sy+mT59OUlJSvZ9d\naduzs7OJj48nOzubpKQkHn30UWprr/IN7Sb55+oqdu/ercXExNQ9fvnll7WXX35ZRRS7cOedd2rJ\nyclar169tOLiYk3TNO3YsWNar169FCdrHmazWRs5cqS2ffv2uiN+Pe6L0tJSrVu3bpf8XI/74tSp\nU1rPnj2106dPa9XV1dq4ceO0rVu36mpf5Ofn1zviv9K2L1iwoN6sSUxMjLZnzx6ry1ZyxH+lc//1\nqKCggKysLCIjIykpKcHHxwcAHx8fSkoc6yv11+upp57i1VdfxeU3dxLR477Iz8+nY8eOTJ8+nQED\nBvDII49QUVGhy33h7e3N008/TdeuXenSpQteXl5ER0frcl/84krbfvTo0XpnTDakT5UUvz1/aas5\nlZeXM2nSJBYtWkTbtm3r/c5gMOhiPyUkJNCpUyfCw8Ov+L0OveyLmpoaMjMzefTRR8nMzMTT0/OS\nqQy97Isff/yRt956i4KCAo4ePUp5eTmrVq2qN0Yv++JyrrbtV9svSoq/Id8PcHbV1dVMmjSJqVOn\nMmHCBOA//4oXFxcDcOzYMTp1cv5LOu/evZsNGzbQrVs3pkyZwvbt25k6daou94Wfnx9+fn4MGjQI\ngLvvvpvMzEw6d+6su33xzTffMGzYMDp06ICbmxsTJ05kz549utwXv7jSe+Jy36Xy9fW1uiwlxa/3\nc/w1TWPGjBmEhoby5JNP1v08NjaWjz76CICPPvqo7h8EZ7ZgwQLMZjP5+fmsWbOGESNGsHLlSl3u\ni86dO+Pv78+hQ4cASElJoXfv3owfP153+yI4OJi0tDQuXLiApmmkpKQQGhqqy33xiyu9J2JjY1mz\nZg1VVVXk5+eTm5vL4MGDrS/M1h9INFRiYqLWs2dPrUePHtqCBQtUxVBi165dmsFg0MLCwrT+/ftr\n/fv31zZv3qydOnVKGzlypNOfqnYlJpNJGz9+vKZpmm73xd69e7WIiAitX79+2l133aWVlpbqdl8s\nXLiw7nTOBx54QKuqqtLNvoiLi9NuvPFGzd3dXfPz89M+/PBDq9v+0ksvaT169NB69eqlJSUlXXX5\nSu/AJYQQovnJHbiEEEJnpPiFEEJnpPiFEEJnpPiFEEJnpPiFEEJnpPiFEEJn/g9QoTIBnbU/lAAA\nAABJRU5ErkJggg==\n" + } + ], + "prompt_number": 8 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [], + "language": "python", + "metadata": {}, + "outputs": [] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file diff --git a/test/triqs/mc_tools/CMakeLists.txt b/test/triqs/mc_tools/CMakeLists.txt new file mode 100644 index 00000000..2b4bb98a --- /dev/null +++ b/test/triqs/mc_tools/CMakeLists.txt @@ -0,0 +1,3 @@ +link_libraries (triqs) +all_tests() + diff --git a/test/triqs/mc_tools/different_moves_mc.cpp b/test/triqs/mc_tools/different_moves_mc.cpp new file mode 100644 index 00000000..87d0bd72 --- /dev/null +++ b/test/triqs/mc_tools/different_moves_mc.cpp @@ -0,0 +1,110 @@ +// The goal of this program is to give an example of the use of the probability for proposing a move (here, pl and pr) +// The walker try to go towards its left with probability 2*pl/(pl+pr) (pl is given to the move_left constructor), +// or towards the right with probability 2*pr/(pl+pr) (pr is given to the move_right constructor). +// The left move is accepted with probability proba=pr/pl and the right move with probability proba=pl/pr (in reality, each time with probability min(1,proba) ). +// A python notebook is associated to this example, with the comparison to the theoretical gaussian histogram. + +#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK +#include +#include + +triqs::arrays::array,1> make_array( std::complex c){return {c}; }; + +struct configuration { + int x; //position of the walker + configuration():x(0){} //constructor +}; + +//------------------------------------------------------------ +// MOVES +//------------------------------------------------------------ +struct move_left { + configuration *config; + double proba; + triqs::mc_tools::random_generator &RND; + move_left(configuration &config_, double pl, double pr, triqs::mc_tools::random_generator &RND_):config(&config_), proba(pr/pl), RND(RND_){} //constructor + double attempt() {return proba; } + double accept() {config->x += -1; return 1; } + void reject() {} +}; +struct move_right { + configuration *config; + double proba; + triqs::mc_tools::random_generator &RND; + move_right( configuration &config_, double pl, double pr, triqs::mc_tools::random_generator &RND_):config(&config_), proba(pl/pr), RND(RND_){} //constructor + double attempt() {return proba; } + double accept() {config->x += 1; return 1; } + void reject() {} +}; + +//------------------------------------------------------------ +// MEASUREMENTS +//------------------------------------------------------------ +struct compute_histo{ + configuration *config; + triqs::arrays::array H; // an histogram + int tot, xmax; + compute_histo(configuration &config_, triqs::arrays::array &H_, int xmax_) : config(&config_), H(H_), xmax(xmax_), tot(0) {} + void accumulate(double sign) { + if( config->x+xmax>=0 && config->x-xmax<0){ + H( floor( (config->x+xmax) ) ) += 1; + tot += 1; + } + config->x=0; + } + void collect_results(boost::mpi::communicator const &c) { + H/=tot; + H5::H5File file("histo.h5",H5F_ACC_TRUNC); + h5_write(file,"H",H); + } +}; + +//------------------------------------------------------------ +// MAIN +//------------------------------------------------------------ +int main(int argc, char* argv[]) { + + // initialize mpi + boost::mpi::environment env(argc, argv); + boost::mpi::communicator world; + + // greeting + if (world.rank() == 0) std::cout << "Random walk calculation" << std::endl; + + // prepare the MC parameters + int N_Cycles = 100000; + int Length_Cycle = 100; + int N_Warmup_Cycles = 0; + std::string Random_Name = ""; + int Random_Seed = 374982 + world.rank() * 273894; + int Verbosity = (world.rank() == 0 ? 3 : 0); + int xmax=floor(4*sqrt(Length_Cycle) ); //max of the position registered in the histogram + double pl=1.5, pr=1; //non normalized probabilities for proposing a left or right move + + H5::H5File file("params.h5",H5F_ACC_TRUNC); + h5_write(file,"pr",make_array(pr)); + h5_write(file,"pl",make_array(pl)); + h5_write(file,"xmax",make_array(xmax)); + h5_write(file,"N_Cycles",make_array(N_Cycles)); + h5_write(file,"Length_Cycle",make_array(Length_Cycle)); + + // construct a Monte Carlo loop + triqs::mc_tools::mc_generic IntMC(N_Cycles, Length_Cycle, N_Warmup_Cycles, Random_Name, Random_Seed, Verbosity); + + // construct configuration + configuration config; + triqs::arrays::array H(2*xmax); + H()=0; + + // add moves and measures + IntMC.add_move(move_left (config, pl, pr, IntMC.rng() ), "left move", pl); + IntMC.add_move(move_right(config, pl, pr, IntMC.rng() ), "right move", pr); + IntMC.add_measure(compute_histo(config, H, xmax), "histogramn measure"); + + // Run and collect results + IntMC.start(1.0, triqs::utility::clock_callback(600)); + IntMC.collect_results(world); + + return 0; + +}