gennet_abc_backup.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"PATH\"] += os.pathsep + 'C:/Users/anura/anaconda3/Library/bin/graphviz'\n",
    "import sys\n",
    "import networkx as nx\n",
    "import numpy as np\n",
    "import random\n",
    "import pickle\n",
    "import sys\n",
    "from numba import jit\n",
    "import elfi\n",
    "import sklearn as sk\n",
    "import scipy\n",
    "from scipy.spatial import distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 20170530  # this will be separately given to ELFI\n",
    "np.random.seed(seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#parameters of the script\n",
    "\n",
    "networktype = 'pref' #pref, smallworld, grid, ER, korea1, korea2, ckm\n",
    "# nAgents = 4**2\n",
    "nAgents = 10**2\n",
    "sidelength = int(nAgents**.5) #10^2 = 100, 32^2 = 1024\n",
    "haltMin = .6 # minimum % of nodes active\n",
    "haltMax = .7 # maximum % of nodes active\n",
    "\n",
    "pRewire = .1 # proportion of edges to rewire\n",
    "numDatasets = 100 # number of data sets to generate\n",
    "saveData = True #save the output to file?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def testThresh(agents, mn, mx):\n",
    "    if np.mean(agents) > mx:\n",
    "        return False\n",
    "    if np.mean(agents) < mn:\n",
    "        return False\n",
    "    if np.mean(agents) <= 0:\n",
    "        return False\n",
    "    if np.mean(agents) >= 1:\n",
    "        return False\n",
    "\n",
    "    return True\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def testThresh(agents, mn, mx):\n",
    "    if np.mean(agents) > mx:\n",
    "        return False\n",
    "    if np.mean(agents) < mn:\n",
    "        return False\n",
    "    if np.mean(agents) <= 0:\n",
    "        return False\n",
    "    if np.mean(agents) >= 1:\n",
    "        return False\n",
    "\n",
    "    return True\n",
    "\n",
    "\n",
    "def genNet(n, k=4, pRewire=.1, type='grid'):\n",
    "  # create net\n",
    "  if type == 'grid': #wrap the grid\n",
    "    net = nx.grid_2d_graph(int(n**.5), int(n**.5), periodic=True)\n",
    "    #net = nx.grid_2d_graph(3, 3, periodic=True)\n",
    "    #net = nx.grid_2d_graph(2, 2, periodic=False)\n",
    "    #plotNet(net)\n",
    "    #print(len(net.edges()))\n",
    "    #assert(0)\n",
    "    # rewire\n",
    "    numRewired = 0\n",
    "    #while numRewired < (pRewire * nx.number_of_nodes(net)):\n",
    "    while numRewired < 1:\n",
    "      tries = 0\n",
    "      while tries < 100:\n",
    "        tries = tries + 1\n",
    "        #print([numRewired, pRewire * nx.number_of_nodes(net), tries])\n",
    "        v1 = random.choice(net.nodes())\n",
    "        v2 = random.choice(net.nodes())\n",
    "        if not( net.has_edge(v1,v2) or v1==v2 or len(net.neighbors(v1)) <= 1): #net.neighbors is sometimes (often?) a blank set, changed so v1 needs 2 nb\n",
    "          #print net.neighbors(v1)\n",
    "          break\n",
    "      v1Neighbors = net.neighbors(v1)\n",
    "      #print v1Neighbors\n",
    "      #print v1\n",
    "      #print v2\n",
    "      #print(len(net.edges()))\n",
    "      tobeDeleted = random.choice(v1Neighbors)\n",
    "      net.remove_edge(v1, tobeDeleted)\n",
    "      #print(len(net.edges()))\n",
    "      #print([v1, tobeDeleted, v2])\n",
    "      net.add_edge(v1, v2)\n",
    "      numRewired = numRewired + 1\n",
    "    #plotNet(net)\n",
    "    #assert(0)\n",
    "    return net, nx.to_numpy_matrix(net, dtype=np.float)\n",
    "\n",
    "  elif type == 'smallworld':\n",
    "    #net = nx.connected_watts_strogatz_graph(n, k, .15)\n",
    "    net = nx.connected_watts_strogatz_graph(n, k, pRewire)\n",
    "    return net, nx.to_numpy_matrix(net, dtype=np.float)  \n",
    "  elif type == 'pref':\n",
    "    net = nx.barabasi_albert_graph(n, 2)\n",
    "    return net, nx.to_numpy_matrix(net, dtype=np.float)\n",
    "  elif type == 'ER':\n",
    "    #net = nx.erdos_renyi_graph(n, .006)\n",
    "    targetDegree = 4.\n",
    "    nEdgesPossible = ((n*n)-n)/2.\n",
    "    pEdge = (n * targetDegree) / (2. * nEdgesPossible)\n",
    "    assert(pEdge <= 1)\n",
    "\n",
    "    # spare networks will likely be disconnected, so try a bunch\n",
    "    tries = 100\n",
    "    while tries > 0:\n",
    "      net = nx.erdos_renyi_graph(n, pEdge)\n",
    "      if nx.number_connected_components(net) > 1:\n",
    "        tries = tries - 1\n",
    "      else:\n",
    "        break\n",
    "    return net, nx.to_numpy_matrix(net, dtype=np.float)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ltProp(agents, adjMatrix, avgDegree=0, haltMin=.49, haltMax=.51, rs = None):\n",
    "    # init some stuff for numba\n",
    "    thresholds = np.zeros_like(nAgents)\n",
    "    inp = np.zeros_like(nAgents)\n",
    "    step = 0\n",
    "    numNeighbors = np.zeros_like(nAgents)\n",
    "    prevMean = 0.\n",
    "    liveEdges = np.zeros_like(adjMatrix)\n",
    "    pInfect = np.zeros_like(adjMatrix)\n",
    "    flips = np.zeros_like(adjMatrix)\n",
    "    \n",
    "    globalThreshold = .5\n",
    "    if rs is None:\n",
    "        thresholds = globalThreshold * np.random.random((1, nAgents))\n",
    "    else:\n",
    "        thresholds = globalThreshold * rs.random((1, nAgents))\n",
    "    \n",
    "    numNeighbors = np.sum(adjMatrix, axis=0)\n",
    "    prevMean = -1\n",
    "    step = 1\n",
    "    \n",
    "    while not testThresh(agents, haltMin, haltMax) and (np.mean(agents) > prevMean):\n",
    "        #while np.mean(agents) > prevMean:\n",
    "        prevMean = np.mean(agents)\n",
    "        inp = np.true_divide(np.dot(agents, adjMatrix), numNeighbors)\n",
    "        agents = np.logical_or(agents, (inp >= thresholds)).astype(int)\n",
    "\n",
    "        step = step + 1\n",
    "        \n",
    "    #print('LT-proportional stopped at step ' + str(step) + ' '+ str(np.mean(agents)))\n",
    "    return agents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ltAbs(agents, adjMatrix, avgDegree=0, haltMin=.49, haltMax=.51, rs = None):\n",
    "    # init some stuff for numba\n",
    "    thresholds = np.zeros_like(nAgents)\n",
    "    inp = np.zeros_like(nAgents)\n",
    "    step = 0\n",
    "    numNeighbors = np.zeros_like(nAgents)\n",
    "    prevMean = 0.\n",
    "    liveEdges = np.zeros_like(adjMatrix)\n",
    "    pInfect = np.zeros_like(adjMatrix)\n",
    "    flips = np.zeros_like(adjMatrix)\n",
    "\n",
    "    # this controls the thresholds/pInfects for all contagion types\n",
    "    globalThreshold = .5\n",
    "    \n",
    "    if rs is None:\n",
    "        thresholds = np.random.randint(low=1, high=round(avgDegree*globalThreshold), size=(1, nAgents))\n",
    "    else:\n",
    "        thresholds = rs.randint(low=1, high=round(avgDegree*globalThreshold), size=(1, nAgents))\n",
    "        \n",
    "    numNeighbors = np.sum(adjMatrix, axis=0)\n",
    "    #prevMean = -1\n",
    "    step = 1\n",
    "    \n",
    "    while not testThresh(agents, haltMin, haltMax) and (np.mean(agents) > prevMean):\n",
    "        prevMean = np.mean(agents)\n",
    "        inp = np.dot(agents, adjMatrix)\n",
    "        agents = np.logical_or(agents, (inp >= thresholds)).astype(int)\n",
    "\n",
    "        step = step + 1\n",
    "        \n",
    "    print('LT-absolute stopped at step ' + str(step) + ' '+ str(np.mean(agents)))\n",
    "\n",
    "    return agents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generateGraph(pp, batch_size = 1, random_state=None):\n",
    "    gg = []\n",
    "    print(pp)\n",
    "    while True:\n",
    "        nAgents = sidelength**2\n",
    "        agents = np.zeros((1, nAgents))\n",
    "        net, adjMatrix = genNet(nAgents, type='pref')\n",
    "\n",
    "        for i in range(1): #make seeding more random, exp with neighbors\n",
    "            #agents[0][1] = 1\n",
    "            agents[0][random.randint(0, nAgents-1)] = 1.\n",
    "\n",
    "        # seed neighbors of seeds\n",
    "        for i in range(1):\n",
    "            agents = np.logical_or(agents, np.dot(agents, adjMatrix)).astype(float) #all neighbors\n",
    "\n",
    "        avgDegree = 2*net.number_of_edges() / float(net.number_of_nodes())\n",
    "\n",
    "        if pp >= 0.5:\n",
    "            agents_type = ltProp(agents, adjMatrix, avgDegree=avgDegree, haltMin=haltMin, haltMax=haltMax, rs = random_state)\n",
    "        else:\n",
    "            agents_type = ltAbs(agents, adjMatrix, avgDegree=avgDegree, haltMin=haltMin, haltMax=haltMax, rs = random_state)\n",
    "\n",
    "        #if not testThresh(agents_type, haltMin, haltMax):\n",
    "         #   print('bad data, LTabs1:\\t'+str(np.mean(agents_type)))\n",
    "       #     continue\n",
    "        gparam = {}\n",
    "        gparam['init'] = net\n",
    "        gparam['graph'] = agents_type\n",
    "        \n",
    "        gg.append(gparam['graph'])\n",
    "        if len(gg) >= batch_size:\n",
    "            break\n",
    "\n",
    "    return gg\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compileG(ll, batch_size = 1, random_state=None):\n",
    "    gall = []\n",
    "    for g in ll:\n",
    "        gall.append(generateGraph(g, batch_size, random_state))\n",
    "    return gall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ga = compileG([0.5])\n",
    "gr_obs = ga[0] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(gr_obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prop_prob = elfi.Prior(scipy.stats.uniform, 0, 1)\n",
    "print(type(prop_prob))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Y = elfi.Simulator(generateGraph, prop_prob, observed = gr_obs)\n",
    "Y = elfi.Simulator(compileG, prop_prob, observed = gr_obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def yMean(x):\n",
    "    return np.mean(x)\n",
    "s = elfi.Summary(yMean, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = elfi.Distance('euclidean', s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "elfi.draw(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rej = elfi.Rejection(d, batch_size=1, seed=seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10\n",
    "# You can give the sample method a `vis` keyword to see an animation how the prior transforms towards the\n",
    "# posterior with a decreasing threshold.\n",
    "%time result = rej.sample(N, quantile=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result.samples['prop_prob'].mean()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}