sphere_stacking.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\"\"\"IMPORT PACKAGES\"\"\"\n",
    "import numpy as np #array magic from numpy \n",
    "from scipy.stats import lognorm #import the log-normal distribution for choosing sphere sizes from it \n",
    "import copy #to make deep copies which do not change the original array they are copied from\n",
    "import time # to time stuff \n",
    "from stack_tools import *\n",
    "from move_tools import *\n",
    "import mayavi\n",
    "from mayavi import mlab\n",
    "%gui qt\n",
    "from IPython.display import clear_output\n",
    "import itertools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sequential_addition(particles,mgrid,dims,delta,grid_size,debug): \n",
    "    \"\"\"Take in an array of spherical particles [r,x,y,z,t] where r is radius, x,y,z are position coordinates-- z needs to be\n",
    "    higher tha any point in the mgrid, mgrid is a grid containing the particles, save==True generates folder named folder\"\"\"\n",
    "    for i,p in enumerate(particles):\n",
    "        clear_output(wait=True)\n",
    "        print('iteration ' +str(i))\n",
    "        p,debug = rearrange_sa(p,grid_size,dims,delta,mgrid)\n",
    "        \n",
    "def collective_rearrangement(particles,mgrid,dims,delta,grid_size):\n",
    "    \"\"\"take in a spatial mgrid filled with potentially overlapping particles [r,x,y,z,t] (radius,position,tracer ind)\n",
    "    and rearrange them until they no longer overlap\"\"\"\n",
    "    for i,p in enumerate(particles):\n",
    "        clear_output(wait=True)\n",
    "        print('iteration ' +str(i))\n",
    "        s = rearrange_cr(p,grid_size,dims,delta,mgrid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example A: Generate 500 uniform spheres with uniformly distributed radius within a 5x5x5 grid, then sequentially add them into a dense packing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I think find Odir is broken when N_cons >=2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dims = [10,10,10]\n",
    "mdims = list(np.array(dims)+np.array([2,2,0]))\n",
    "NX,NY,NZ = dims\n",
    "delta = 0.1\n",
    "radii_array = np.ones(5000,dtype=float)\n",
    "#radii_array = np.random.uniform(size=5000)*0.8+0.2\n",
    "#radii_array = np.random.lognormal(mean=0.5,sigma=0.2,size=5000)\n",
    "#delta = np.amin(radii_array)/4\n",
    "grid_size = np.amax(radii_array)*2+delta\n",
    "spheres = initialize_spheres(radii_array,dims,grid_size,delta,seq_add=True)\n",
    "#for s in spheres:\n",
    "#    s.position = np.array([NX/2,NY/2,NZ/2],dtype=float)*grid_size+np.random.random(size=3)*delta\n",
    "mgrid = initialize_grid(mdims)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iteration 280\n",
      "(array([6.52126196e-01, 8.84789856e-03, 2.31000000e+01]), 1.0, 0)\n",
      "(array([6.52126196e-01, 8.84789856e-03, 2.31000000e+01]), 1.0, 1)\n",
      "0 before D\n",
      "0 after D\n",
      "(array([0.6521262, 0.0088479, 5.379    ]), 1.0, 3)\n",
      "0 before O\n",
      "[5.338698177762397]: zees\n",
      "[-0.04897815402812743]: thes\n",
      "[array([-0.45938481, -0.88823735,  0.        ])]: axes\n",
      "[array([0.6521262, 0.0088479, 5.379    ])] p position\n",
      "[array([-0.03846166,  0.36601096,  3.51319873])]: contacts\n",
      "attempting drop out of overlap generating O\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "'numpy.ndarray' object has no attribute 'position'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-69-d04140529b0a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     22\u001b[0m             \u001b[0;32mcontinue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     23\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0ms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstate\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 24\u001b[0;31m             \u001b[0ms\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mO\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcons\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mgrid_size\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdims\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdelta\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mmgrid\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     25\u001b[0m             \u001b[0;32mcontinue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     26\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0ms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstate\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/home/kvnp/Desktop/sphere_stacking/stack_tools.py\u001b[0m in \u001b[0;36mO\u001b[0;34m(p, cons, grid_size, dims, delta, mgrid)\u001b[0m\n\u001b[1;32m    512\u001b[0m                     \u001b[0mo\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstate\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    513\u001b[0m                     \u001b[0;32mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'attempting drop out of overlap generating O'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 514\u001b[0;31m                     \u001b[0;32mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'generated overlap is '\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mo\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mposition\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mo\u001b[0m \u001b[0;32min\u001b[0m \u001b[0moverlaps\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdo\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mgrid_size\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdelta\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mmgrid\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mposition\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    515\u001b[0m                     \u001b[0;32mbreak\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    516\u001b[0m     \u001b[0;32mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moverlaps\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mo\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mgrid_size\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdelta\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mmgrid\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m' after O'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'numpy.ndarray' object has no attribute 'position'"
     ]
    }
   ],
   "source": [
    "debug = []\n",
    "for i,p in enumerate(spheres):\n",
    "    clear_output(wait=True)\n",
    "    print('iteration ' +str(i))\n",
    "    \n",
    "    NX,NY,NZ = dims\n",
    "    s = copy.deepcopy(p)\n",
    "    s.state=0    \n",
    "    \n",
    "    while True: #this loop just passes the particle through G,D,W,O,OO,R as necessary\n",
    "        print(s.position,s.radius,s.state)\n",
    "        \n",
    "        cons = c_sort(s,grid_size,delta,mgrid)\n",
    "        if s.state==0:\n",
    "            s = G_sa(s,cons,grid_size,dims,delta,mgrid)\n",
    "            continue \n",
    "        elif s.state==1:\n",
    "            s = D(s,grid_size,delta,mgrid)\n",
    "            continue\n",
    "        elif s.state==2:\n",
    "            s = W(s,grid_size,dims)\n",
    "            continue\n",
    "        elif s.state==3:\n",
    "            s = O(s,cons,grid_size,dims,delta,mgrid)\n",
    "            continue \n",
    "        elif s.state==4: \n",
    "            s = OO(s,cons,grid_size,dims,delta,mgrid)\n",
    "            continue\n",
    "        elif s.state==5:\n",
    "            if indices(s,grid_size)[-1] < NZ: \n",
    "                s = R(s,grid_size,dims,mgrid)\n",
    "            else: \n",
    "                print('Particle past max height')\n",
    "                s = None \n",
    "            break    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-43.03943838, -17.16561086,   1.0014588 ])"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s.position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "view(mgrid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "s.tag [q.tag for q in itr(mgrid)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view_condition(debug,lambda x: x.tag == s.tag)\n",
    "view(debug)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p = None\n",
    "q = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p.position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "\n",
    "def rotation_matrix(axis, theta):\n",
    "    \"\"\"\n",
    "    Return the rotation matrix associated with counterclockwise rotation about\n",
    "    the given axis by theta radians.\n",
    "    \"\"\"\n",
    "    axis = np.asarray(axis)\n",
    "    axis = axis/math.sqrt(np.dot(axis, axis))\n",
    "    a = math.cos(theta/2.0)\n",
    "    b, c, d = -axis*math.sin(theta/2.0)\n",
    "    aa, bb, cc, dd = a*a, b*b, c*c, d*d\n",
    "    bc, ad, ac, ab, bd, cd = b*c, a*d, a*c, a*b, b*d, c*d\n",
    "    return np.array([[aa+bb-cc-dd, 2*(bc+ad), 2*(bd-ac)],\n",
    "                     [2*(bc-ad), aa+cc-bb-dd, 2*(cd+ab)],\n",
    "                     [2*(bd+ac), 2*(cd-ab), aa+dd-bb-cc]])\n",
    "\n",
    "v = [3, 5, 0]\n",
    "axis = [4, 4, 1]\n",
    "theta = 1.2 \n",
    "\n",
    "print(np.dot(rotation_matrix(axis,theta), v)) \n",
    "# [ 2.74911638  4.77180932  1.91629719]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#need to test the rotator \n",
    "x = 2.1*5/2\n",
    "dims = (5,5,5)\n",
    "delta = 0.5\n",
    "grid_size = 2.1\n",
    "mgrid = initialize_grid((7,7,5))\n",
    "\n",
    "p = sphere(1.0,1,dims,grid_size,delta)\n",
    "p.position = np.array([x,x,x],dtype=float)\n",
    "\n",
    "q = sphere(1.0,2,dims,grid_size,delta)\n",
    "q.position = np.array([x-delta/10,x,x-1.85],dtype=float)\n",
    "\n",
    "q = R(q,grid_size,dims,mgrid)\n",
    "p  = R(p,grid_size,dims,mgrid)\n",
    "\n",
    "#m = sphere(1.0,3,dims,grid_size,delta)\n",
    "#m.position = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "[p.position for p in itr(mgrid)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "contacts(p,grid_size,delta,mgrid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cons = contacts(p,grid_size,delta,mgrid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "dthe,A,P = find_Odir(p,cons,grid_size,(5,5,5),delta,mgrid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parts = []\n",
    "for i in range(200):\n",
    "    p.position = rot(dthe,A,P,p)\n",
    "    parts.append(p.position)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for p in parts:\n",
    "    R(p,grid_size,dims,mgrid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "view(mgrid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[p.position for p in parts]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "h = sphere(1.0,3,dims,grid_size,delta)\n",
    "h.position = np.array([5.70,5.25,5.061])\n",
    "h = R(h,grid_size,dims,mgrid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "view(mgrid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q.position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def find_Odir(p,cons,grid_size,dims,delta,mgrid):\n",
    "    NX,NY,NZ = dims \n",
    "    zees = [] #z value after rotation\n",
    "    thes = [] #rotation angle increments to sweep through delta\n",
    "    axes = [] #rotation axes\n",
    "    pivs = [] #pivot points \n",
    "    x,y,z = p.position #unpack particle position\n",
    "    for c in cons:\n",
    "        #print(p.position)\n",
    "        #print(c.position)\n",
    "        c1,c2,c3 = P = c.position #pivot point\n",
    "        dthe = 0.99*delta/np.linalg.norm(p.position-c.position) #unsigned rotation increment \n",
    "        A = np.array([y-c2,c1-x,0]) #rotation axis \n",
    "        if np.linalg.norm(A)>0:\n",
    "            A = A/np.linalg.norm(A)\n",
    "        for sign in [-1,1]:\n",
    "            o = copy.deepcopy(p) #a copy of p to test the rotation \n",
    "            dthe = sign*dthe #signed rotation increment\n",
    "            print('before: '+str(o.position))\n",
    "            o.position = rot(dthe,A,P,o) #perform the test rotation\n",
    "            print('med' + str(o.position))\n",
    "            if not within_grid(o,grid_size,dims): #if test rotation left boundary, wrap it back\n",
    "                o = W(o,grid_size,dims)\n",
    "            if np.isnan(o.position[0]):\n",
    "                return (dthe,A,P,o)\n",
    "        \n",
    "            print('after: '+str(o.position))\n",
    "\n",
    "            Z = o.position[-1] #what is the new value of Z induced by the test rotation?\n",
    "            Z0 = p.position[-1]\n",
    "            N_ov = len(overlaps(o,grid_size,delta,mgrid)) #how many overlaps were generated by the test rotation?\n",
    "            if  Z<Z0 and N_ov==0:#it o is lower than p and has no overlaps\n",
    "                zees.append(Z) #add potential movement parameters \n",
    "                thes.append(dthe)\n",
    "                axes.append(A)\n",
    "                pivs.append(P)\n",
    "    #now among the possibilities find the optimal set of parameters \n",
    "    if len(zees)>0: #find largest decrease in z \n",
    "        min_z, idx = min((val, idx) for (idx, val) in enumerate(zees)) #smallest z coordinate after rotation\n",
    "        A = axes[idx] #largest decrease in z through rot for particle p induced by these parameters\n",
    "        P = pivs[idx]\n",
    "        dthe = thes[idx]\n",
    "        return (dthe,A,P) #return the tuple of optimal rotation parameters\n",
    "    else:\n",
    "        return None #return statement that O fails "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rot(,A,P,o) #perform the test rotation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example B: Generate 300 uniform size spheres within a 5x5x5 grid, then collectively rearrange them into a dense packing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s.tag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "NX,NY,NZ = dims = (5,5,5)\n",
    "mgrid = initialize_grid(dims = (NX+2,NY+2,NZ))\n",
    "delta = 0.05\n",
    "radius = 0.5  \n",
    "grid_size = 2*radius+delta #2*max_radius + delta\n",
    "N = 300\n",
    "radii_array = np.ones(N)*radius\n",
    "spheres = initialize_spheres(radii_array,dims,grid_size,delta)\n",
    "collective_rearrangement(spheres,mgrid,dims,delta,grid_size)\n",
    "view(mgrid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "view(mgrid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# transport model development"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#just make a bed for testing the transport model \n",
    "NX,NY,NZ = dims = [20,20,4]\n",
    "mgrid0 = initialize_grid((NX+2,NY+2,NZ))\n",
    "delta = 0.05\n",
    "radius = 0.5  \n",
    "grid_size = 2*radius+delta #2*max_radius + delta\n",
    "N = 2000\n",
    "radii_array = np.ones(N)*radius\n",
    "spheres = initialize_spheres(radii_array,dims,grid_size,delta)\n",
    "collective_rearrangement(spheres,mgrid0,dims,delta,grid_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#now make z bigger \n",
    "NZ = 2*NZ\n",
    "dims[-1]=NZ\n",
    "mgrid = extend_mgrid(mgrid0,NZ)\n",
    "mgrid0 = copy.deepcopy(mgrid)\n",
    "view(mgrid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mgrid = copy.deepcopy(mgrid0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## verbose globally updating transport model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#reset\n",
    "mgrid = copy.deepcopy(mgrid0)\n",
    "#initialize the model only once\n",
    "clusters = determine_clusters(mgrid,grid_size,delta)  \n",
    "#generate transport array\n",
    "transport = initialize_grid(dims[:-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "view_clusters(mgrid,clusters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i=0\n",
    "transport_rate=[]\n",
    "travel_times = []\n",
    "entrainment_rate = []\n",
    "deposition_rate = []\n",
    "N=len(itr(mgrid))\n",
    "print(str(len(itr(mgrid)))+' particles initially in system.')\n",
    "while True:\n",
    "    i+=1\n",
    "    \n",
    "    clear_output(wait=True)\n",
    "    #determine the clusters to be entrained\n",
    "    print('-----------------------------------------------------------------------')\n",
    "    print('ITERATION '+str(i))\n",
    "    print('-----------------------------------------------------------------------')\n",
    "    entrained_clusters = determine_entrained(clusters)\n",
    "    #print(str(len(entrained_clusters))+' clusters to be entrained')\n",
    "    #print(str(len([cp for ec in entrained_clusters for cp in ec.particles]))+ ' particles to be entrained')\n",
    "    #print(str(n_unstable(mgrid,delta,grid_size))+ ' particles unstable total')\n",
    "    #print(str(len(fullitr(transport)))+ ' particles in transport')\n",
    "    \n",
    "    #actually entrain the clusters into the transport array \n",
    "    #n1 = len(itr(mgrid))\n",
    "    #m1 = len(fullitr(transport))\n",
    "    #l1 = n_unstable(mgrid,delta,grid_size)\n",
    "    #####################\n",
    "    clusters = entrain_into_transport(transport,clusters, mgrid,grid_size,delta,*entrained_clusters) #modifies transport in place\n",
    "    #####################\n",
    "    #n2 = len(itr(mgrid))\n",
    "    #m2 = len(fullitr(transport))\n",
    "    #entrainment_rate.append(m2-m1)\n",
    "    #l2 = n_unstable(mgrid,delta,grid_size)\n",
    "    #if not -(n2-n1)==m2-m1:\n",
    "    #    print('number conservation fail in entrain_into_transport')\n",
    "    #    break\n",
    "    #if l1==0 and l2!=0:\n",
    "    #    print('unstable particles generated from entrain_into_transport')\n",
    "    #    break\n",
    "\n",
    "    #print(str(n_unstable(mgrid,delta,grid_size))+ ' particles unstable total')\n",
    "\n",
    "    \n",
    "    #reconstruct clusters array following modification of bed surface\n",
    "    #print(len(clusters), ' clusters initially')\n",
    "    #print(len([p for cl in clusters for p in cl.particles]), ' particles in clusters initially')\n",
    "    ###############################\n",
    "    #clusters = determine_clusters(mgrid,grid_size,delta)\n",
    "    ####################################\n",
    "    #print(len(clusters), ' clusters after redeterminiation of clusters')\n",
    "    #print(len([p for cl in clusters for p in cl.particles]), ' particles in clusters after redetermination of clusters')\n",
    "    #print(str(n_unstable(mgrid,delta,grid_size))+ ' particles unstable total')\n",
    "\n",
    "    \n",
    "    #transport the active particles \n",
    "    #m1 = len(fullitr(transport))\n",
    "    ###########################\n",
    "    transport = translate_by_timestep(transport,grid_size,delta,dims)\n",
    "    ############################\n",
    "    #m2 = len(fullitr(transport))\n",
    "    #print(str(n_unstable(mgrid,delta,grid_size))+ ' particles unstable total')\n",
    "    \n",
    "    \n",
    "    #determine particles to be disentrained\n",
    "    particles_to_deposit = determine_deposited(transport)\n",
    "    #print(len(particles_to_deposit), 'particles slated for deposition')\n",
    "    #print(str(n_unstable(mgrid,delta,grid_size))+ ' particles unstable total')\n",
    "    \n",
    "    \n",
    "    #disentrain the particles back onto the bed \n",
    "    #n1 = len(fullitr(transport))\n",
    "    #m1 = len(itr(mgrid))\n",
    "    #l1 = n_unstable(mgrid,delta,grid_size)\n",
    "    transport,travel_iterations = deposit_onto_bed(transport,mgrid,grid_size,delta,dims,*particles_to_deposit)\n",
    "    travel_times.extend(travel_iterations)\n",
    "    #l2 = n_unstable(mgrid,delta,grid_size)\n",
    "    #n2 = len(fullitr(transport))\n",
    "    #deposition_rate.append(n1-n2)\n",
    "    #m2 = len(itr(mgrid))\n",
    "    #if not -(n2-n1)==m2-m1:\n",
    "    #    print('number conservation fail in deposit_onto_bed')\n",
    "    #    break\n",
    "    #if l1==0 and l2!=0:\n",
    "    #    print('deposit_onto_bed generated an instability')\n",
    "    #    break\n",
    "    #print(str(n_unstable(mgrid,delta,grid_size))+ ' particles unstable total')\n",
    "    \n",
    "    \n",
    "    #redetermine clusters array following modification of bed surface\n",
    "\n",
    "    #print(len(clusters), ' clusters initially')\n",
    "    #print(len([p for cl in clusters for p in cl.particles]), ' particles in clusters initially')\n",
    "    clusters = determine_clusters(mgrid,grid_size,delta)\n",
    "    #clus3inclus2(clusters)\n",
    "    #print(len(clusters), ' clusters after redeterminiation of clusters')\n",
    "    #print(len([p for cl in clusters for p in cl.particles]), ' particles in clusters after redetermination of clusters')\n",
    "    #print(str(n_unstable(mgrid,delta,grid_size))+ ' particles unstable total')\n",
    "    #print(str(len(itr(mgrid))+len(fullitr(transport)))+' particles total')\n",
    "    #print(str(len(fullitr(transport)))+ ' particles in transport')\n",
    "    transport_rate.append(len(fullitr(transport)))\n",
    "\n",
    "    if n_unstable(mgrid,delta,grid_size)!=0:\n",
    "        print('unstable break')\n",
    "        break\n",
    "    #if len(itr(mgrid))+len(fullitr(transport))!=N:\n",
    "    #    print('particle conservation break')\n",
    "    #    break\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "?T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.save('travel_times.npy',travel_times)\n",
    "np.save('transport_rate.npy',transport_rate)\n",
    "np.save('entrainment_rate.npy',entrainment_rate)\n",
    "np.save('deposition_rate.npy',deposition_rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "clus3inclus2(clusters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "view_clusters(mgrid,clusters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "question: \n",
    "can a 3 cluster contain a 2 cluster?\n",
    "answer: \n",
    "yes. \n",
    "idea:\n",
    "incorporate equality between particles... p1==p2 if their tags are the same\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def clus3inclus2(clusters):\n",
    "    clusters2 = [clus for clus in clusters if len(clus.particles)==2]\n",
    "    clusters3 = [clus for clus in clusters if len(clus.particles)==3]\n",
    "    tags2 = [p.tag for cl in clusters2 for p in cl.particles]\n",
    "    tags3 = [p.tag for cl in clusters3 for p in cl.particles]\n",
    "    for t2 in tags2:\n",
    "        for t3 in tags3:\n",
    "            for t in t2.split('_'):\n",
    "                for w in t3.split('_'):\n",
    "                    if t==w:\n",
    "                        print('a 2 in a 3')\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "clusters[0].__dict__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def degenerate_clusters(clusters):\n",
    "    cparticles = [clu.particles for clu in clusters]\n",
    "    parts = [p for clup in cparticles for p in clup]\n",
    "    parts2=[]\n",
    "    for clus in cparticles:\n",
    "        for dlus in cparticles:\n",
    "            for c in clus:\n",
    "                for d in dlus:\n",
    "                    if c.tag==d.tag:\n",
    "                        parts2.append(c)\n",
    "    if len(parts)!=len(parts2):\n",
    "        print(str(len(parts2)-len(parts))+' degenerate clusters in system.')             "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## local updating transport model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#here is the model\n",
    "for i in range(100):\n",
    "    print('iteration '+str(i))\n",
    "    print(str(n_unstable(mgrid,delta,grid_size))+ ' particles unstable total')\n",
    "    \n",
    "    entrained_clusters = determine_entrained(clusters)\n",
    "    \n",
    "    print(str(len(entrained_clusters))+ ' entrained clusters')\n",
    "    print(' len clusters array before entrainment: ' + str(len(clusters))+ '; '+ str(n_unstable(mgrid,delta,grid_size))+ ' unstable')\n",
    "    \n",
    "    clusters = entrain_clusters(clusters,mgrid,grid_size,delta,*entrained_clusters)\n",
    "    \n",
    "    print(' len clusters array after entrainment: ' + str(len(clusters))+ '; '+ str(n_unstable(mgrid,delta,grid_size))+ ' unstable')\n",
    "    print(' number of particles in transport array before entrainment: ' + str(len(fullitr(transport)))+ '; '+ str(n_unstable(mgrid,delta,grid_size))+ ' unstable')\n",
    "    \n",
    "    transport = entrain_into_transport(transport,mgrid,grid_size,delta,*entrained_clusters)\n",
    "    \n",
    "    print(' number of particles in transport array after entrainment: ' + str(len(fullitr(transport)))+ '; '+ str(n_unstable(mgrid,delta,grid_size))+ ' unstable')\n",
    "    \n",
    "    transport = translate_by_timestep(transport,grid_size,delta,dims)\n",
    "    \n",
    "    print(' number of particles in transport array after displacement: ' + str(len(fullitr(transport)))+ '; '+ str(n_unstable(mgrid,delta,grid_size))+ ' unstable')\n",
    "    \n",
    "    deposited_clusters = determine_deposited(transport)\n",
    "    \n",
    "    print(str(len(deposited_clusters))+ ' particled slated for deposition from transport array')\n",
    "    print(str(len(itr(mgrid)))+' particles in mgrid before deposition.'+ '; '+ str(n_unstable(mgrid,delta,grid_size))+ ' unstable')\n",
    "    \n",
    "    transport,clusters = deposit_onto_bed(transport,clusters,mgrid,grid_size,delta,dims,*deposited_clusters)     \n",
    "    \n",
    "    print(str(len(itr(mgrid)))+' particles in mgrid after deposition.'+ '; '+ str(n_unstable(mgrid,delta,grid_size))+ ' unstable') \n",
    "    print('--------------------------------------------------------------')    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for i in range(100):\n",
    "    print('iteration '+str(i))\n",
    "    print(str(n_unstable(mgrid,delta,grid_size))+ ' particles unstable total')\n",
    "    #print(len(itr(mgrid))+len(fullitr(transport)),' total particles at iteration ', i)\n",
    "    entrained_clusters = determine_entrained(clusters)[:1]\n",
    "    clusters = determine_clusters(mgrid,grid_size,delta)\n",
    "    print(str(len(entrained_clusters))+ ' entrained clusters')\n",
    "    print(' len clusters array before entrainment: ' + str(len(clusters))+ '; '+ str(n_unstable(mgrid,delta,grid_size))+ ' unstable')\n",
    "    mgrid_bug = copy.deepcopy(mgrid)\n",
    "    clusters_bug = copy.deepcopy(clusters)\n",
    "    entrained_clusters_bug = copy.deepcopy(entrained_clusters)\n",
    "    \n",
    "    clusters = entrain_clusters(clusters,mgrid,grid_size,delta,*entrained_clusters)\n",
    "    \n",
    "    if n_unstable(mgrid,delta,grid_size)!=0:\n",
    "        print('broken')\n",
    "        break\n",
    "        \n",
    "    print(' len clusters arrmay after entrainment: ' + str(len(clusters))+ '; '+ str(n_unstable(mgrid,delta,grid_size))+ ' unstable')\n",
    "    print(' number of particles in transport array before entrainment: ' + str(len(fullitr(transport)))+ '; '+ str(n_unstable(mgrid,delta,grid_size))+ ' unstable')\n",
    "    transport = entrain_into_transport(transport,mgrid,grid_size,delta,*entrained_clusters)\n",
    "    print(' number of particles in transport array after entrainment: ' + str(len(fullitr(transport)))+ '; '+ str(n_unstable(mgrid,delta,grid_size))+ ' unstable')\n",
    "    transport = translate_by_timestep(transport,grid_size,delta,dims)\n",
    "    print(' number of particles in transport array after displacement: ' + str(len(fullitr(transport)))+ '; '+ str(n_unstable(mgrid,delta,grid_size))+ ' unstable')\n",
    "    deposited_clusters = determine_deposited(transport)\n",
    "    print(str(len(deposited_clusters))+ ' particled slated for deposition from transport array')\n",
    "    print(str(len(itr(mgrid)))+' particles in mgrid before deposition.'+ '; '+ str(n_unstable(mgrid,delta,grid_size))+ ' unstable')\n",
    "    contacts_above_in_clusters(clusters,mgrid,flag='before deposition')\n",
    "\n",
    "    transport,clusters = deposit_onto_bed(transport,clusters,mgrid,grid_size,delta,dims,*deposited_clusters)\n",
    "    contacts_above_in_clusters(clusters,mgrid,'after depsosition')\n",
    "    print(str(len(itr(mgrid)))+' particles in mgrid after deposition.'+ '; '+ str(n_unstable(mgrid,delta,grid_size))+ ' unstable') \n",
    "    print('--------------------------------------------------------------')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mgrid =     copy.deepcopy(mgrid_bug) \n",
    "clusters =     copy.deepcopy(clusters_bug)\n",
    "entrained_clusters = copy.deepcopy(entrained_clusters_bug)\n",
    "transport = translate_by_timestep(transport,grid_size,delta,dims)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def filter_invalid_clusters(clusters_to_check,clusters,mgrid,grid_size,delta):\n",
    "    \"\"\"check whether a list of clusters are all valid in terms of contingency\"\"\"\n",
    "    new_clusters = []\n",
    "    for clu in clusters_to_check: \n",
    "        #find the index of clu in clusters\n",
    "        j = clusters.index(clu) #this is the index of clu in clusters\n",
    "        cluclusters = []\n",
    "        cluclusters = update_clusters(clu.particles,cluclusters,mgrid,grid_size,delta)\n",
    "        #does clu remain in cluclusters? if so, ok. if not, must delete clu from clusters\n",
    "        if clu not in cluclusters:\n",
    "            del clusters[j] #remove clu from clusters, because deposition removed it from cluster status \n",
    "\n",
    "        clusters.extend([qlu for qlu in cluclusters if not qlu==clu]) #either way, add any new clusters \n",
    "    return clusters          "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def filter_clusters_from_deposition(settled_particles,clusters,mgrid,grid_size,delta):\n",
    "    \"\"\"depo_cons are the contacts of recently deposited particles. If they are also within clusters,\n",
    "    these clusters need to be reconsidered\"\"\"\n",
    "    clusters_to_reconsider=[]\n",
    "    particles_to_reconsider_for_clustering = []\n",
    "    for clus in clusters: #take a cluster already within clusters. ask if it is affected by deposition\n",
    "        #if settled_particles landed on a cluster, need to reexamine all settled particles AND the settled particle\n",
    "        flag=True# each cluster in need of examination will be added to the examination list only once\n",
    "        for c in clus.particles: #for each particle in the cluster\n",
    "            for s in settled_particles: #for each particle which has settled\n",
    "                if s.radius+c.radius-delta<distance(c,s)<s.radius+c.radius+delta: #if cluster particle in contact with settled particle\n",
    "                #if s in contact with cluster \n",
    "                    if flag: \n",
    "                    #if this is the first such s to contact this cluster \n",
    "                        particles_to_reconsider_for_clustering.extend(clus.particles) #then reconsider its particles\n",
    "                        clusters_to_reconsider.append(clus) # and reconsider its status as a cluster \n",
    "                        flag=False #but only once per cluster\n",
    "                    particles_to_reconsider_for_clustering.append(s) #and also must consider recently deposited s \n",
    "    \n",
    "    #now need to take clusters list and filter it over clusters_to_reconsider. method is filter_invalid_clusters \n",
    "    #need a function that checks if a cluster is valid \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    #finally need to take particles_to_reconsider_for_clustering list and use update_clusters on it\n",
    "    clusters = update_clusters(particles_to_reconsider_for_clustering,clusters,mgrid,grid_size,delta)\n",
    "    return clusters "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:mayavi]",
   "language": "python",
   "name": "conda-env-mayavi-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}