GITHUB_abc_position_reconstruction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ABC applied to position reconstruction\n",
    "\n",
    "*Bart Pelssers, 26-02-2018*\n",
    "\n",
    "This notebook provides some ingredients for applying the ABC algorithm to position reconstruction, the most basic case. Just to test the framework.\n",
    "\n",
    "\n",
    "* Provides:\n",
    "  * prior mean\n",
    "  * forward model\n",
    "  * summary statistic\n",
    "  \n",
    "*Umberto Simola, 22-03-2018*\n",
    "\n",
    "Provided the ABC-PMC algorithm for running the analyses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from abc_reconstruction.simpleModel import SimpleModel as Model\n",
    "from abc_reconstruction.utils import Generator, PriorPosition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup the Models\n",
    "model = Model()\n",
    "prior_mean = PriorPosition()\n",
    "generator = Generator(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Length of pattern: 127, Sum of pattern: 500.00\n",
      "[  0.14841595   0.14397198   0.13988137   0.24159058   0.18933869\n",
      "   0.22591982   0.20805915   0.23506457   0.31582144   0.28612889\n",
      "   0.33985588   0.38021411   0.48340495   0.56819242   0.67613899\n",
      "   0.8470959    1.0991352    1.25354921   1.30006622   1.43827737\n",
      "   1.17062482   1.07080343   0.89733495   0.78772861   0.62917211\n",
      "   0.40026487   0.44570271   0.32362463   0.301934     0.27624214\n",
      "   0.28848058   0.21737252   0.28141832   0.23859565   0.16602611\n",
      "   0.21634046   0.36649959   0.47636725   0.44341384   0.45075319\n",
      "   0.5022231    0.53645639   0.5805464    0.68153719   0.71537357\n",
      "   0.79270054   0.91853777   1.15995621   1.30145776   1.76802897\n",
      "   2.17324897   2.91035191   2.87717576   2.39581913   1.79414895\n",
      "   1.44549592   1.05471978   0.81379478   0.74454349   0.69864421\n",
      "   0.590455     0.62706968   0.5057492    0.42780855   0.41207846\n",
      "   0.43831424   0.58574815   0.56964808   0.68396244   0.58466108\n",
      "   0.61196534   0.84590827   0.83313924   0.99069221   1.20798275\n",
      "   1.83992028   3.19752242   5.82342523  10.2681775    9.54624032\n",
      "   5.59288128   2.95237534   1.79524267   1.28391557   0.929167\n",
      "   0.83422188   0.70135432   0.63948979   0.5272279    0.47227546\n",
      "   0.69310408   0.60975179   0.68438211   0.8536707    0.97884602\n",
      "   1.18834598   1.91369186   4.96917712  18.36824766  62.15840859\n",
      "  46.47265287  12.05379726   3.49258965   1.68795161   0.99465714\n",
      "   0.96993554   0.8280701    0.70821068   0.93718629   0.91659382\n",
      "   1.1246014    1.4322029    3.6739197   19.57090981 103.9310611\n",
      "  48.17341292   7.18703857   2.04794571   1.44048082   0.85647978\n",
      "   1.535153     1.93117414   5.56232536  22.95846181   8.18196078\n",
      "   2.41219596   4.34920249]\n"
     ]
    }
   ],
   "source": [
    "# Example pattern from x=2.6264, y=-17.96082\n",
    "\n",
    "# The range of x and y is [-47.884375 cm, -47.884375 cm]\n",
    "# But x**2 + y**2 < 47.884375**2\n",
    "# Otherwise m() will raise and exception\n",
    "\n",
    "pattern = model(2.6264, -17.96082)\n",
    "\n",
    "print(\"Length of pattern: %d, Sum of pattern: %.2f\" % (len(pattern), pattern.sum()))\n",
    "print(pattern)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4.115222817130077, -15.358220637996187)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# If using a 2D Normal prior this would be a good guess for the mean (x,y) of that prior.\n",
    "prior_mean(pattern)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "###Function for sampling proposal coordinates from the prior: Normal Centered on the naiveCoord, sd=15\n",
    "def priorFunction(coord):\n",
    "    coordX = coord[0]\n",
    "    coordY = coord[1]\n",
    "    xProp=np.random.normal(coordX,15,1)[0]\n",
    "    yProp=np.random.normal(coordY,15,1)[0]\n",
    "    while xProp**2+yProp**2>47.884375**2:\n",
    "        xProp=np.random.normal(coordX,15,1)[0]\n",
    "        yProp=np.random.normal(coordY,15,1)[0]\n",
    "    return xProp,yProp\n",
    "#priorFunction(naiveCoord)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "###Distance function for comparing the real and the simulated dataset\n",
    "def rho(x,y):\n",
    "    return sum((x-y)**2)/np.shape(x)[0]\n",
    "    #return sum(np.abs(x-y))/np.shape(x)[0]\n",
    "#print(rho(event,SimulatedModel))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "###Load the true position matrix for generating the Observed dataset\n",
    "\n",
    "###For sure you can load only the positions from the original files but I got errors when trying to call the data\n",
    "\n",
    "trueCoordMatrix=np.loadtxt('data/truepos')\n",
    "#print(trueCoordMatrix[0])\n",
    "#np.shape(trueCoordMatrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "###ABC-SMC Definition of the necessary quantities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "###Transformation kernel for resampling for t>0 (i.e. rather the using the prior we use this kernel)\n",
    "def transfKernel(propx0,propy0,varx,vary):\n",
    "    propx=np.random.normal(propx0,2*varx,1)[0]\n",
    "    propy=np.random.normal(propy0,2*vary,1)[0]\n",
    "    while propx**2+propy**2 > 47.884375**2:\n",
    "        propx=np.random.normal(propx0,2*varx,1)[0]\n",
    "        propy=np.random.normal(propy0,2*vary,1)[0]\n",
    "    return propx,propy "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "###ABC-SMC on a single selected event event j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-16-5ebc29dcee26>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     56\u001b[0m                 \u001b[0msimulatedData\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpropCoord\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpropCoord\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m500\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnormal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mTrueModel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     57\u001b[0m                 \u001b[0mtotDraws\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtotDraws\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 58\u001b[0;31m                 \u001b[0md\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrho\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mTrueModel\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0msimulatedData\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     59\u001b[0m             \u001b[0mabcCoordsx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mpropCoord\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     60\u001b[0m             \u001b[0mabcCoordsy\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mpropCoord\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-11-2423979e929e>\u001b[0m in \u001b[0;36mrho\u001b[0;34m(x, y)\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m###Distance function for comparing the real and the simulated dataset\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mrho\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m     \u001b[0;31m#return sum(np.abs(x-y))/np.shape(x)[0]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;31m#print(rho(event,SimulatedModel))\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "import scipy.stats\n",
    "\n",
    "###Number of particles for recunstructing the posterior for each coordinate\n",
    "N=5000\n",
    "\n",
    "###Number of iterations before stopping the algorithm\n",
    "nIter=40\n",
    "\n",
    "###Quantile used for shrinking the tolerances through the iterations\n",
    "quantile=0.85\n",
    "\n",
    "###Importance weigths for the ABC-PMC: because for t>0 we using a kernel and not the prior as proposal distribution\n",
    "weights=np.zeros((nIter,N))\n",
    "weights[0,:]=1/N\n",
    "\n",
    "###Accepted elements are stored here\n",
    "###abcCooord x\n",
    "abcCoordsx=np.zeros((nIter,N))\n",
    "###abcCooord y\n",
    "abcCoordsy=np.zeros((nIter,N))\n",
    "\n",
    "###Distance of accepted elements are stored here\n",
    "d=np.zeros((nIter,N))\n",
    "\n",
    "###First tolerance epsilon1\n",
    "epsilon=4\n",
    "\n",
    "###Total number of draws required for covering the entire analyses\n",
    "totDraws=0\n",
    "\n",
    "#########################################################################\n",
    "#########################################################################\n",
    "###pick the event j\n",
    "j=0\n",
    "\n",
    "###given an event, we need the true coordinates of this event and getting the 'true' dataset. The trueDataset called here\n",
    "###TrueModel\n",
    "trueCoord = trueCoordMatrix[j] \n",
    "TrueModel=model(trueCoord[0], trueCoord[1], 500)\n",
    "###means of the priors\n",
    "naiveCoord=prior_mean(TrueModel)\n",
    "\n",
    "###Simulated accepted dataset from the last iteration\n",
    "simAccData=np.zeros((N,np.shape(TrueModel)[0]))\n",
    "    \n",
    "for t in range(0,nIter):\n",
    "    print(t)\n",
    "    if(t==0):\n",
    "        for i in range(0,N):\n",
    "            d[t,i]=epsilon+1\n",
    "            while d[t,i] > epsilon: \n",
    "                propCoord = priorFunction(naiveCoord)\n",
    "                ###Deterministic FM\n",
    "                #simulatedData = model(propCoord[0], propCoord[1], 500)\n",
    "                ###Stochastic FM, for now adding gaussian noise (mean=0,sd=1)\n",
    "                simulatedData = model(propCoord[0], propCoord[1], 500)+np.random.normal(0,1,np.shape(TrueModel)[0])\n",
    "                totDraws = totDraws+1\n",
    "                d[t,i] = rho(TrueModel,simulatedData) \n",
    "            abcCoordsx[t,i]=propCoord[0]\n",
    "            abcCoordsy[t,i]=propCoord[1]\n",
    "    else:\n",
    "        epsilon= np.percentile(d[t-1,],quantile*100)\n",
    "        print(epsilon)\n",
    "        meanx=np.sum(abcCoordsx[t-1,:]*weights[t-1,:])\n",
    "        varx=np.sum((abcCoordsx[t-1,:]-meanx)**2*weights[t-1,:])\n",
    "        meany=np.sum(abcCoordsy[t-1,:]*weights[t-1,:])\n",
    "        vary=np.sum((abcCoordsy[t-1,:]-meany)**2*weights[t-1,:])\n",
    "        for i in range(0,N):\n",
    "            d[t,i]=epsilon+1\n",
    "            while d[t,i] > epsilon: \n",
    "                sample=np.random.choice(N,1,p=weights[t-1,:])\n",
    "                propx0=abcCoordsx[t-1,sample]\n",
    "                propy0=abcCoordsy[t-1,sample]\n",
    "                prop=transfKernel(propx0,propy0,varx,vary)\n",
    "                propx=prop[0]\n",
    "                propy=prop[1]\n",
    "                ###Deterministic FM\n",
    "                #simulatedData = model(propx, propy, 500)\n",
    "                ###Stochastic FM, for now adding gaussian noise (mean=0,sd=1)\n",
    "                simulatedData = model(propx, propy, 500)+np.random.normal(0,1,np.shape(TrueModel)[0])\n",
    "                totDraws = totDraws+1\n",
    "                d[t,i] = rho(TrueModel,simulatedData)\n",
    "            simAccData[i,:]=simulatedData\n",
    "            abcCoordsx[t,i]=propx\n",
    "            abcCoordsy[t,i]=propy \n",
    "            weightsDen=np.sum(weights[t-1,:]*scipy.stats.norm.pdf(propx,abcCoordsx[t-1,:],np.sqrt(2*varx))*scipy.stats.norm.pdf(propy,abcCoordsy[t-1,:],np.sqrt(2*vary)))\n",
    "            weightsNum=scipy.stats.norm.pdf(propx,propCoord[0],10)*scipy.stats.norm.pdf(propy,propCoord[1],10)     \n",
    "            weights[t,i]=weightsNum/weightsDen\n",
    "    weights[t,:]=weights[t,:]/sum(weights[t,:])\n",
    "print(totDraws)\n",
    "abcOutput=np.column_stack((abcCoordsx[nIter-1,:],abcCoordsy[nIter-1,:]))\n",
    "np.savetxt('abcPosteriorCoords_%d.dat'%(j), abcOutput)\n",
    "importanceWeigths=weights[nIter-1,:]\n",
    "np.savetxt('WeightsabcPosteriorCoords_%d.dat'%(j), importanceWeigths)\n",
    "np.savetxt('simAccData_%d.dat'%(j), simAccData)\n",
    "np.savetxt('TrueData_%d.dat'%(j), TrueModel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3246]\n"
     ]
    }
   ],
   "source": [
    "###I pick randomly 10 elements for the analyses\n",
    "elements=np.random.choice(np.shape(trueCoordMatrix)[0],1)\n",
    "print(elements)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3246\n",
      "0.875631205344\n"
     ]
    }
   ],
   "source": [
    "import scipy.stats\n",
    "###ABCPMC for all the events or a random selection among the possbilities\n",
    "#for j in np.random.choice(np.shape(trueCoordMatrix)[0],10):\n",
    "#for j in range(0,1): \n",
    "for j in elements:\n",
    "    print(j)\n",
    "    N=5000\n",
    "    nIter=40\n",
    "    quantile=0.85\n",
    "    epsilon=4\n",
    "    ###Importance weigths for the ABC-PMC\n",
    "    weights=np.zeros((nIter,N))\n",
    "    weights[0,:]=1/N\n",
    "    ###abcCooord x\n",
    "    abcCoordsx=np.zeros((nIter,N))\n",
    "    ###abcCooord y\n",
    "    abcCoordsy=np.zeros((nIter,N))\n",
    "    ###Distance\n",
    "    d=np.zeros((nIter,N)) \n",
    "    \n",
    "    trueCoord = trueCoordMatrix[j]\n",
    "    TrueModel=model(trueCoord[0], trueCoord[1], 500)\n",
    "    naiveCoord=prior_mean(TrueModel)\n",
    "    totDraws=0\n",
    "    t=0\n",
    "    for t in range(0,nIter):\n",
    "        #print(t)\n",
    "        if(t==0):\n",
    "            for i in range(0,N):\n",
    "                d[t,i]=epsilon+1\n",
    "                while d[t,i] > epsilon: \n",
    "                    propCoord = priorFunction(naiveCoord)\n",
    "                    ###Deterministic FM\n",
    "                    #simulatedData = model(propCoord[0], propCoord[1], 500)\n",
    "                    ###Stochastic FM, for now adding gaussian noise (mean=0,sd=1)\n",
    "                    simulatedData = model(propCoord[0], propCoord[1], 500)+np.random.normal(0,1,np.shape(TrueModel)[0])\n",
    "                    totDraws = totDraws+1\n",
    "                    d[t,i] = rho(TrueModel,simulatedData)\n",
    "                abcCoordsx[t,i]=propCoord[0]\n",
    "                abcCoordsy[t,i]=propCoord[1]\n",
    "        else:\n",
    "            epsilon= np.percentile(d[t-1,],quantile*100)\n",
    "            #print(epsilon)\n",
    "            meanx=np.sum(abcCoordsx[t-1,:]*weights[t-1,:])\n",
    "            varx=np.sum((abcCoordsx[t-1,:]-meanx)**2*weights[t-1,:])\n",
    "            meany=np.sum(abcCoordsy[t-1,:]*weights[t-1,:])\n",
    "            vary=np.sum((abcCoordsy[t-1,:]-meany)**2*weights[t-1,:])\n",
    "            for i in range(0,N):\n",
    "                d[t,i]=epsilon+1\n",
    "                while d[t,i] > epsilon: \n",
    "                    sample=np.random.choice(N,1,p=weights[t-1,:])\n",
    "                    propx0=abcCoordsx[t-1,sample]\n",
    "                    propy0=abcCoordsy[t-1,sample]\n",
    "                    prop=transfKernel(propx0,propy0,varx,vary)\n",
    "                    propx=prop[0]\n",
    "                    propy=prop[1]\n",
    "                    ###Deterministic FM\n",
    "                    #simulatedData = model(propx, propy, 500)\n",
    "                    ###Stochastic FM, for now adding gaussian noise (mean=0,sd=1)\n",
    "                    simulatedData = model(propx, propy, 500)+np.random.normal(0,1,np.shape(TrueModel)[0])\n",
    "                    totDraws = totDraws+1\n",
    "                    d[t,i] = rho(TrueModel,simulatedData)\n",
    "                abcCoordsx[t,i]=propx\n",
    "                abcCoordsy[t,i]=propy\n",
    "                weightsDen=np.sum(weights[t-1,:]*scipy.stats.norm.pdf(propx,abcCoordsx[t-1,:],np.sqrt(2*varx))*scipy.stats.norm.pdf(propy,abcCoordsy[t-1,:],np.sqrt(2*vary)))\n",
    "                weightsNum=scipy.stats.norm.pdf(propx,propCoord[0],10)*scipy.stats.norm.pdf(propy,propCoord[1],10)     \n",
    "                weights[t,i]=weightsNum/weightsDen\n",
    "        weights[t,:]=weights[t,:]/sum(weights[t,:])\n",
    "    ###here we save the coordinated saved in the last iteration, and corresponding weights.\n",
    "    abcOutput=np.column_stack((abcCoordsx[nIter-1,:],abcCoordsy[nIter-1,:]))\n",
    "    np.savetxt('abcPosteriorCoords_%d.dat'%(j), abcOutput)\n",
    "    importanceWeigths=weights[nIter-1,:]\n",
    "    np.savetxt('WeightsabcPosteriorCoords_%d.dat'%(j), importanceWeigths)\n",
    "    print(epsilon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}