{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3f591549-cf33-4816-9a2c-7da947efafd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as ssp\n",
    "import scipy.sparse.linalg as sspl\n",
    "from scipy.special import gamma\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "353a0a02-6665-455a-a16d-a910df6a682a",
   "metadata": {},
   "outputs": [],
   "source": [
    "steps = 4\n",
    "Nx = 2**steps\n",
    "Ny = Nx\n",
    "Nz = 2*(steps+2)-2\n",
    "NxNy = Nx*Ny\n",
    "N = NxNy*Nz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "a01fe040-251a-435c-a60c-7cd5fbb855c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "el0 = np.loadtxt('./edgelist_0_seed1_0', dtype=int)\n",
    "nedges0 = len(el0)\n",
    "indices = np.lexsort((el0[:,1],el0[:,0])) #needed for CSR conversion below\n",
    "el0 = el0[indices]\n",
    "###CSR static, do not modify\n",
    "dat = np.tile([-1,1],nedges0) \n",
    "ind = el0.flatten()\n",
    "ptr = np.arange(0,2*(1+nedges0),2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "7d2c781c-ba83-4fd1-8446-ac531782fcd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rfm4(tt, boundary, niter): \n",
    "    conductances = np.ones(nedges0, dtype=int)\n",
    "    inc = ssp.csr_array((np.copy(dat),np.copy(ind),np.copy(ptr)), shape=(nedges0, N))\n",
    "    u = np.zeros(N, dtype=float) #voltages\n",
    "    Ntop = N-boundary\n",
    "    Nbot = boundary\n",
    "    ebot = boundary\n",
    "    u[Ntop:N] = 1 #top boundary condition\n",
    "    I,V = [0],[0]\n",
    "    ns = 0\n",
    "    nr = 0\n",
    "    #while(True):\n",
    "    for _ in range(niter):\n",
    "        ci = (conductances*inc.T).T\n",
    "        lap = inc.T@ci\n",
    "        rhs = -lap[Nbot:Ntop,Ntop:N]@u[Ntop:N]\n",
    "        x = sspl.spsolve(lap[Nbot:Ntop,Nbot:Ntop], rhs )\n",
    "        ns += 1\n",
    "        u[Nbot:Ntop] = x[:]\n",
    "        loc_i = ci@u\n",
    "        if np.mean(loc_i[0:ebot])<1e-6:\n",
    "            print('Failure after '+str(ns)+' steps')\n",
    "            break\n",
    "        ratios = np.abs(loc_i/tt)\n",
    "        weakid = np.argmax(ratios)\n",
    "        ivscale = ratios[weakid]\n",
    "        topb_i = ( \n",
    "                lap[Ntop:, :Nbot]@u[:Nbot] + \n",
    "                lap[Ntop:, Nbot:Ntop]@u[Nbot:Ntop]+\n",
    "                lap[Ntop:, Ntop:]@u[Ntop:]\n",
    "        )\n",
    "        glo_i = np.sum(topb_i)/ivscale\n",
    "        I.append(glo_i)\n",
    "        V.append(1/ivscale)\n",
    "        conductances[weakid] = 0\n",
    "        nr += 1\n",
    "    return np.array(V),np.array(I)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a3e98945-3497-48a2-9685-d303fd83bf2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 1\n",
    "weik = 4\n",
    "ra = np.random.default_rng(seed=seed)\n",
    "ths0 = ra.weibull(weik,nedges0)/gamma(1+1/weik)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "7c1a20b7-0aa2-4b39-a82a-c80fde91f450",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Failure after 480 steps\n"
     ]
    }
   ],
   "source": [
    "V, I = rfm4(ths0, NxNy,10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "8509fa6c-e1db-4f7c-8698-d90204fcd4a2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'current')"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 71892 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots()\n",
    "ax.plot(V,I)\n",
    "ax.set_xlabel('voltage')\n",
    "ax.set_ylabel('current')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "5fadb170-3251-4fec-8f85-a6285b2489da",
   "metadata": {},
   "outputs": [],
   "source": [
    "def envelope_from(x,y):\n",
    "    #x,y = data[:,0],data[:,1]\n",
    "    xx, yy = [], []\n",
    "    xx.append(0)\n",
    "    yy.append(0)\n",
    "    xx.append(x[0])\n",
    "    yy.append(y[0])\n",
    "    for t in range(1,len(x)):\n",
    "        if x[t]>xx[-1]:\n",
    "            xx.append(x[t])\n",
    "            yy.append(y[t])\n",
    "    xx.append(xx[-1])\n",
    "    yy.append(0)\n",
    "    return np.array(xx),np.array(yy)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "c6a00a5d-6f85-4a06-b35d-98d28b974626",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'current')"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 74304 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots()\n",
    "ax.plot(V,I)\n",
    "Ve,Ie = envelope_from(V,I)\n",
    "ax.plot(Ve,Ie)\n",
    "ax.set_xlabel('voltage')\n",
    "ax.set_ylabel('current')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfea5c4c-2f6d-4060-b3e5-cfb0c69d9505",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
