Skip to content

Instantly share code, notes, and snippets.

@kvedala

kvedala/Jacobian to solve linear system.ipynb

Last active June 25, 2020 19:28
Show Gist options
  • Save kvedala/2dbfe7be8a3e2e36ade606152fe5db09 to your computer and use it in GitHub Desktop.
Save kvedala/2dbfe7be8a3e2e36ade606152fe5db09 to your computer and use it in GitHub Desktop.
Download ZIP
Jacobian to solve linear system.ipynb
Display the source blob
Display the rendered blob
Raw
Jacobian to solve linear system.ipynb
{
"cells": [
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "import numpy as np\nimport numba\n%reload_ext cython\n%reload_ext line_profiler\n%reload_ext memory_profiler",
"execution_count": 1,
"outputs": []
},
{
"metadata": {
"scrolled": false,
"trusted": true
},
"cell_type": "code",
"source": "def jacobian(A: np.ndarray, b: np.ndarray) -> np.ndarray:\n MAX_ITER = 1000\n N = A.shape[0]\n epsilon = 1e-10\n x0 = np.random.random(A.shape[1]) - 0.5\n x1 = x0.copy()\n delta, k = 1, 1\n max_delta = 0\n for k in range(MAX_ITER):\n delta = 0\n max_delta = 0\n for i in range(N):\n sigma = 0\n for j in range(N):\n if i == j: continue \n sigma += A[i,j] * x0[j]\n x1[i] = (b[i] - sigma) / A[i,i]\n delta = abs(x1[i] - x0[i])\n max_delta = delta if delta > max_delta else max_delta\n# print(k, sigma, max_delta)\n# print('-----')\n if max_delta < epsilon:\n break\n x0 = x1.copy()\n return x1",
"execution_count": 2,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "@numba.njit\ndef jacobian_n(A, b):\n MAX_ITER = 1000\n N = A.shape[0]\n x0 = np.random.random(A.shape[1]) - 0.5\n x1 = x0.copy()\n epsilon = 1e-10\n delta, k = 1, 1\n max_delta = 0\n for k in range(MAX_ITER):\n max_delta = 0\n for i in range(N):\n sigma = 0\n for j in range(N):\n if i == j: continue \n sigma += A[i,j] * x0[j]\n x1[i] = (b[i] - sigma) / A[i,i]\n delta = abs(x1[i] - x0[i])\n max_delta = delta if delta > max_delta else max_delta\n# print(k, sigma, max_delta)\n# print('-----')\n if max_delta < epsilon:\n break\n x0 = x1.copy()\n return x1",
"execution_count": 3,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "%%cython -a -c=-flto=thin\n#cython: boundscheck=False, wraparound=False, cdivision=True\n#-#cython: profile=True, linetrace=True, binding=True\n#-#distutils: define_macros=CYTHON_TRACE_NOGIL=1\n\ncimport numpy as np\nimport numpy as np\n\ncdef extern from \"math.h\":\n double fabs(double) nogil\n double fmax(double, double) nogil\n \nfrom libc.stdio cimport printf\n\ndef jacobian_c(double[:,:] A, double[:] b):\n cdef:\n int MAX_ITER = 100\n int N = A.shape[0]\n double epsilon = 1e-10\n np.ndarray[ndim=1,dtype=double] x = np.random.random(A.shape[1]) - 0.5\n double[:] x0 = x.copy()\n# double[:] x = _x\n# double[:,:] _A = A\n# double[:] _b = b\n double sigma, delta, max_delta\n int k = 1, i, j\n \n with nogil:\n for k in range(MAX_ITER):\n max_delta = 0\n for i in range(N):\n sigma = 0.0\n for j in range(N):\n if i != j: \n sigma += A[i,j] * x0[j]\n x[i] = (b[i] - sigma) / A[i,i]\n delta = fabs(x[i] - x0[i])\n max_delta = fmax(delta, max_delta)\n# with gil:\n# print(k, sigma, max_delta)\n# print('-----')\n if max_delta < epsilon:\n break\n for i in range(N):\n x0[i] = x[i]\n return x",
"execution_count": 4,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 4,
"data": {
"text/plain": "<IPython.core.display.HTML object>",
"text/html": "<!DOCTYPE html>\n<!-- Generated by Cython 0.29.15 -->\n<html>\n<head>\n <meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n <title>Cython: _cython_magic_89adf65e112d455b90bc5490de5d61a1.pyx</title>\n <style type=\"text/css\">\n \nbody.cython { font-family: courier; font-size: 12; }\n\n.cython.tag { }\n.cython.line { margin: 0em }\n.cython.code { font-size: 9; color: #444444; display: none; margin: 0px 0px 0px 8px; border-left: 8px none; }\n\n.cython.line .run { background-color: #B0FFB0; }\n.cython.line .mis { background-color: #FFB0B0; }\n.cython.code.run { border-left: 8px solid #B0FFB0; }\n.cython.code.mis { border-left: 8px solid #FFB0B0; }\n\n.cython.code .py_c_api { color: red; }\n.cython.code .py_macro_api { color: #FF7000; }\n.cython.code .pyx_c_api { color: #FF3000; }\n.cython.code .pyx_macro_api { color: #FF7000; }\n.cython.code .refnanny { color: #FFA000; }\n.cython.code .trace { color: #FFA000; }\n.cython.code .error_goto { color: #FFA000; }\n\n.cython.code .coerce { color: #008000; border: 1px dotted #008000 }\n.cython.code .py_attr { color: #FF0000; font-weight: bold; }\n.cython.code .c_attr { color: #0000FF; }\n.cython.code .py_call { color: #FF0000; font-weight: bold; }\n.cython.code .c_call { color: #0000FF; }\n\n.cython.score-0 {background-color: #FFFFff;}\n.cython.score-1 {background-color: #FFFFe7;}\n.cython.score-2 {background-color: #FFFFd4;}\n.cython.score-3 {background-color: #FFFFc4;}\n.cython.score-4 {background-color: #FFFFb6;}\n.cython.score-5 {background-color: #FFFFaa;}\n.cython.score-6 {background-color: #FFFF9f;}\n.cython.score-7 {background-color: #FFFF96;}\n.cython.score-8 {background-color: #FFFF8d;}\n.cython.score-9 {background-color: #FFFF86;}\n.cython.score-10 {background-color: #FFFF7f;}\n.cython.score-11 {background-color: #FFFF79;}\n.cython.score-12 {background-color: #FFFF73;}\n.cython.score-13 {background-color: #FFFF6e;}\n.cython.score-14 {background-color: #FFFF6a;}\n.cython.score-15 {background-color: #FFFF66;}\n.cython.score-16 {background-color: #FFFF62;}\n.cython.score-17 {background-color: #FFFF5e;}\n.cython.score-18 {background-color: #FFFF5b;}\n.cython.score-19 {background-color: #FFFF57;}\n.cython.score-20 {background-color: #FFFF55;}\n.cython.score-21 {background-color: #FFFF52;}\n.cython.score-22 {background-color: #FFFF4f;}\n.cython.score-23 {background-color: #FFFF4d;}\n.cython.score-24 {background-color: #FFFF4b;}\n.cython.score-25 {background-color: #FFFF48;}\n.cython.score-26 {background-color: #FFFF46;}\n.cython.score-27 {background-color: #FFFF44;}\n.cython.score-28 {background-color: #FFFF43;}\n.cython.score-29 {background-color: #FFFF41;}\n.cython.score-30 {background-color: #FFFF3f;}\n.cython.score-31 {background-color: #FFFF3e;}\n.cython.score-32 {background-color: #FFFF3c;}\n.cython.score-33 {background-color: #FFFF3b;}\n.cython.score-34 {background-color: #FFFF39;}\n.cython.score-35 {background-color: #FFFF38;}\n.cython.score-36 {background-color: #FFFF37;}\n.cython.score-37 {background-color: #FFFF36;}\n.cython.score-38 {background-color: #FFFF35;}\n.cython.score-39 {background-color: #FFFF34;}\n.cython.score-40 {background-color: #FFFF33;}\n.cython.score-41 {background-color: #FFFF32;}\n.cython.score-42 {background-color: #FFFF31;}\n.cython.score-43 {background-color: #FFFF30;}\n.cython.score-44 {background-color: #FFFF2f;}\n.cython.score-45 {background-color: #FFFF2e;}\n.cython.score-46 {background-color: #FFFF2d;}\n.cython.score-47 {background-color: #FFFF2c;}\n.cython.score-48 {background-color: #FFFF2b;}\n.cython.score-49 {background-color: #FFFF2b;}\n.cython.score-50 {background-color: #FFFF2a;}\n.cython.score-51 {background-color: #FFFF29;}\n.cython.score-52 {background-color: #FFFF29;}\n.cython.score-53 {background-color: #FFFF28;}\n.cython.score-54 {background-color: #FFFF27;}\n.cython.score-55 {background-color: #FFFF27;}\n.cython.score-56 {background-color: #FFFF26;}\n.cython.score-57 {background-color: #FFFF26;}\n.cython.score-58 {background-color: #FFFF25;}\n.cython.score-59 {background-color: #FFFF24;}\n.cython.score-60 {background-color: #FFFF24;}\n.cython.score-61 {background-color: #FFFF23;}\n.cython.score-62 {background-color: #FFFF23;}\n.cython.score-63 {background-color: #FFFF22;}\n.cython.score-64 {background-color: #FFFF22;}\n.cython.score-65 {background-color: #FFFF22;}\n.cython.score-66 {background-color: #FFFF21;}\n.cython.score-67 {background-color: #FFFF21;}\n.cython.score-68 {background-color: #FFFF20;}\n.cython.score-69 {background-color: #FFFF20;}\n.cython.score-70 {background-color: #FFFF1f;}\n.cython.score-71 {background-color: #FFFF1f;}\n.cython.score-72 {background-color: #FFFF1f;}\n.cython.score-73 {background-color: #FFFF1e;}\n.cython.score-74 {background-color: #FFFF1e;}\n.cython.score-75 {background-color: #FFFF1e;}\n.cython.score-76 {background-color: #FFFF1d;}\n.cython.score-77 {background-color: #FFFF1d;}\n.cython.score-78 {background-color: #FFFF1c;}\n.cython.score-79 {background-color: #FFFF1c;}\n.cython.score-80 {background-color: #FFFF1c;}\n.cython.score-81 {background-color: #FFFF1c;}\n.cython.score-82 {background-color: #FFFF1b;}\n.cython.score-83 {background-color: #FFFF1b;}\n.cython.score-84 {background-color: #FFFF1b;}\n.cython.score-85 {background-color: #FFFF1a;}\n.cython.score-86 {background-color: #FFFF1a;}\n.cython.score-87 {background-color: #FFFF1a;}\n.cython.score-88 {background-color: #FFFF1a;}\n.cython.score-89 {background-color: #FFFF19;}\n.cython.score-90 {background-color: #FFFF19;}\n.cython.score-91 {background-color: #FFFF19;}\n.cython.score-92 {background-color: #FFFF19;}\n.cython.score-93 {background-color: #FFFF18;}\n.cython.score-94 {background-color: #FFFF18;}\n.cython.score-95 {background-color: #FFFF18;}\n.cython.score-96 {background-color: #FFFF18;}\n.cython.score-97 {background-color: #FFFF17;}\n.cython.score-98 {background-color: #FFFF17;}\n.cython.score-99 {background-color: #FFFF17;}\n.cython.score-100 {background-color: #FFFF17;}\n.cython.score-101 {background-color: #FFFF16;}\n.cython.score-102 {background-color: #FFFF16;}\n.cython.score-103 {background-color: #FFFF16;}\n.cython.score-104 {background-color: #FFFF16;}\n.cython.score-105 {background-color: #FFFF16;}\n.cython.score-106 {background-color: #FFFF15;}\n.cython.score-107 {background-color: #FFFF15;}\n.cython.score-108 {background-color: #FFFF15;}\n.cython.score-109 {background-color: #FFFF15;}\n.cython.score-110 {background-color: #FFFF15;}\n.cython.score-111 {background-color: #FFFF15;}\n.cython.score-112 {background-color: #FFFF14;}\n.cython.score-113 {background-color: #FFFF14;}\n.cython.score-114 {background-color: #FFFF14;}\n.cython.score-115 {background-color: #FFFF14;}\n.cython.score-116 {background-color: #FFFF14;}\n.cython.score-117 {background-color: #FFFF14;}\n.cython.score-118 {background-color: #FFFF13;}\n.cython.score-119 {background-color: #FFFF13;}\n.cython.score-120 {background-color: #FFFF13;}\n.cython.score-121 {background-color: #FFFF13;}\n.cython.score-122 {background-color: #FFFF13;}\n.cython.score-123 {background-color: #FFFF13;}\n.cython.score-124 {background-color: #FFFF13;}\n.cython.score-125 {background-color: #FFFF12;}\n.cython.score-126 {background-color: #FFFF12;}\n.cython.score-127 {background-color: #FFFF12;}\n.cython.score-128 {background-color: #FFFF12;}\n.cython.score-129 {background-color: #FFFF12;}\n.cython.score-130 {background-color: #FFFF12;}\n.cython.score-131 {background-color: #FFFF12;}\n.cython.score-132 {background-color: #FFFF11;}\n.cython.score-133 {background-color: #FFFF11;}\n.cython.score-134 {background-color: #FFFF11;}\n.cython.score-135 {background-color: #FFFF11;}\n.cython.score-136 {background-color: #FFFF11;}\n.cython.score-137 {background-color: #FFFF11;}\n.cython.score-138 {background-color: #FFFF11;}\n.cython.score-139 {background-color: #FFFF11;}\n.cython.score-140 {background-color: #FFFF11;}\n.cython.score-141 {background-color: #FFFF10;}\n.cython.score-142 {background-color: #FFFF10;}\n.cython.score-143 {background-color: #FFFF10;}\n.cython.score-144 {background-color: #FFFF10;}\n.cython.score-145 {background-color: #FFFF10;}\n.cython.score-146 {background-color: #FFFF10;}\n.cython.score-147 {background-color: #FFFF10;}\n.cython.score-148 {background-color: #FFFF10;}\n.cython.score-149 {background-color: #FFFF10;}\n.cython.score-150 {background-color: #FFFF0f;}\n.cython.score-151 {background-color: #FFFF0f;}\n.cython.score-152 {background-color: #FFFF0f;}\n.cython.score-153 {background-color: #FFFF0f;}\n.cython.score-154 {background-color: #FFFF0f;}\n.cython.score-155 {background-color: #FFFF0f;}\n.cython.score-156 {background-color: #FFFF0f;}\n.cython.score-157 {background-color: #FFFF0f;}\n.cython.score-158 {background-color: #FFFF0f;}\n.cython.score-159 {background-color: #FFFF0f;}\n.cython.score-160 {background-color: #FFFF0f;}\n.cython.score-161 {background-color: #FFFF0e;}\n.cython.score-162 {background-color: #FFFF0e;}\n.cython.score-163 {background-color: #FFFF0e;}\n.cython.score-164 {background-color: #FFFF0e;}\n.cython.score-165 {background-color: #FFFF0e;}\n.cython.score-166 {background-color: #FFFF0e;}\n.cython.score-167 {background-color: #FFFF0e;}\n.cython.score-168 {background-color: #FFFF0e;}\n.cython.score-169 {background-color: #FFFF0e;}\n.cython.score-170 {background-color: #FFFF0e;}\n.cython.score-171 {background-color: #FFFF0e;}\n.cython.score-172 {background-color: #FFFF0e;}\n.cython.score-173 {background-color: #FFFF0d;}\n.cython.score-174 {background-color: #FFFF0d;}\n.cython.score-175 {background-color: #FFFF0d;}\n.cython.score-176 {background-color: #FFFF0d;}\n.cython.score-177 {background-color: #FFFF0d;}\n.cython.score-178 {background-color: #FFFF0d;}\n.cython.score-179 {background-color: #FFFF0d;}\n.cython.score-180 {background-color: #FFFF0d;}\n.cython.score-181 {background-color: #FFFF0d;}\n.cython.score-182 {background-color: #FFFF0d;}\n.cython.score-183 {background-color: #FFFF0d;}\n.cython.score-184 {background-color: #FFFF0d;}\n.cython.score-185 {background-color: #FFFF0d;}\n.cython.score-186 {background-color: #FFFF0d;}\n.cython.score-187 {background-color: #FFFF0c;}\n.cython.score-188 {background-color: #FFFF0c;}\n.cython.score-189 {background-color: #FFFF0c;}\n.cython.score-190 {background-color: #FFFF0c;}\n.cython.score-191 {background-color: #FFFF0c;}\n.cython.score-192 {background-color: #FFFF0c;}\n.cython.score-193 {background-color: #FFFF0c;}\n.cython.score-194 {background-color: #FFFF0c;}\n.cython.score-195 {background-color: #FFFF0c;}\n.cython.score-196 {background-color: #FFFF0c;}\n.cython.score-197 {background-color: #FFFF0c;}\n.cython.score-198 {background-color: #FFFF0c;}\n.cython.score-199 {background-color: #FFFF0c;}\n.cython.score-200 {background-color: #FFFF0c;}\n.cython.score-201 {background-color: #FFFF0c;}\n.cython.score-202 {background-color: #FFFF0c;}\n.cython.score-203 {background-color: #FFFF0b;}\n.cython.score-204 {background-color: #FFFF0b;}\n.cython.score-205 {background-color: #FFFF0b;}\n.cython.score-206 {background-color: #FFFF0b;}\n.cython.score-207 {background-color: #FFFF0b;}\n.cython.score-208 {background-color: #FFFF0b;}\n.cython.score-209 {background-color: #FFFF0b;}\n.cython.score-210 {background-color: #FFFF0b;}\n.cython.score-211 {background-color: #FFFF0b;}\n.cython.score-212 {background-color: #FFFF0b;}\n.cython.score-213 {background-color: #FFFF0b;}\n.cython.score-214 {background-color: #FFFF0b;}\n.cython.score-215 {background-color: #FFFF0b;}\n.cython.score-216 {background-color: #FFFF0b;}\n.cython.score-217 {background-color: #FFFF0b;}\n.cython.score-218 {background-color: #FFFF0b;}\n.cython.score-219 {background-color: #FFFF0b;}\n.cython.score-220 {background-color: #FFFF0b;}\n.cython.score-221 {background-color: #FFFF0b;}\n.cython.score-222 {background-color: #FFFF0a;}\n.cython.score-223 {background-color: #FFFF0a;}\n.cython.score-224 {background-color: #FFFF0a;}\n.cython.score-225 {background-color: #FFFF0a;}\n.cython.score-226 {background-color: #FFFF0a;}\n.cython.score-227 {background-color: #FFFF0a;}\n.cython.score-228 {background-color: #FFFF0a;}\n.cython.score-229 {background-color: #FFFF0a;}\n.cython.score-230 {background-color: #FFFF0a;}\n.cython.score-231 {background-color: #FFFF0a;}\n.cython.score-232 {background-color: #FFFF0a;}\n.cython.score-233 {background-color: #FFFF0a;}\n.cython.score-234 {background-color: #FFFF0a;}\n.cython.score-235 {background-color: #FFFF0a;}\n.cython.score-236 {background-color: #FFFF0a;}\n.cython.score-237 {background-color: #FFFF0a;}\n.cython.score-238 {background-color: #FFFF0a;}\n.cython.score-239 {background-color: #FFFF0a;}\n.cython.score-240 {background-color: #FFFF0a;}\n.cython.score-241 {background-color: #FFFF0a;}\n.cython.score-242 {background-color: #FFFF0a;}\n.cython.score-243 {background-color: #FFFF0a;}\n.cython.score-244 {background-color: #FFFF0a;}\n.cython.score-245 {background-color: #FFFF0a;}\n.cython.score-246 {background-color: #FFFF09;}\n.cython.score-247 {background-color: #FFFF09;}\n.cython.score-248 {background-color: #FFFF09;}\n.cython.score-249 {background-color: #FFFF09;}\n.cython.score-250 {background-color: #FFFF09;}\n.cython.score-251 {background-color: #FFFF09;}\n.cython.score-252 {background-color: #FFFF09;}\n.cython.score-253 {background-color: #FFFF09;}\n.cython.score-254 {background-color: #FFFF09;}\n.cython .hll { background-color: #ffffcc }\n.cython { background: #f8f8f8; }\n.cython .c { color: #408080; font-style: italic } /* Comment */\n.cython .err { border: 1px solid #FF0000 } /* Error */\n.cython .k { color: #008000; font-weight: bold } /* Keyword */\n.cython .o { color: #666666 } /* Operator */\n.cython .ch { color: #408080; font-style: italic } /* Comment.Hashbang */\n.cython .cm { color: #408080; font-style: italic } /* Comment.Multiline */\n.cython .cp { color: #BC7A00 } /* Comment.Preproc */\n.cython .cpf { color: #408080; font-style: italic } /* Comment.PreprocFile */\n.cython .c1 { color: #408080; font-style: italic } /* Comment.Single */\n.cython .cs { color: #408080; font-style: italic } /* Comment.Special */\n.cython .gd { color: #A00000 } /* Generic.Deleted */\n.cython .ge { font-style: italic } /* Generic.Emph */\n.cython .gr { color: #FF0000 } /* Generic.Error */\n.cython .gh { color: #000080; font-weight: bold } /* Generic.Heading */\n.cython .gi { color: #00A000 } /* Generic.Inserted */\n.cython .go { color: #888888 } /* Generic.Output */\n.cython .gp { color: #000080; font-weight: bold } /* Generic.Prompt */\n.cython .gs { font-weight: bold } /* Generic.Strong */\n.cython .gu { color: #800080; font-weight: bold } /* Generic.Subheading */\n.cython .gt { color: #0044DD } /* Generic.Traceback */\n.cython .kc { color: #008000; font-weight: bold } /* Keyword.Constant */\n.cython .kd { color: #008000; font-weight: bold } /* Keyword.Declaration */\n.cython .kn { color: #008000; font-weight: bold } /* Keyword.Namespace */\n.cython .kp { color: #008000 } /* Keyword.Pseudo */\n.cython .kr { color: #008000; font-weight: bold } /* Keyword.Reserved */\n.cython .kt { color: #B00040 } /* Keyword.Type */\n.cython .m { color: #666666 } /* Literal.Number */\n.cython .s { color: #BA2121 } /* Literal.String */\n.cython .na { color: #7D9029 } /* Name.Attribute */\n.cython .nb { color: #008000 } /* Name.Builtin */\n.cython .nc { color: #0000FF; font-weight: bold } /* Name.Class */\n.cython .no { color: #880000 } /* Name.Constant */\n.cython .nd { color: #AA22FF } /* Name.Decorator */\n.cython .ni { color: #999999; font-weight: bold } /* Name.Entity */\n.cython .ne { color: #D2413A; font-weight: bold } /* Name.Exception */\n.cython .nf { color: #0000FF } /* Name.Function */\n.cython .nl { color: #A0A000 } /* Name.Label */\n.cython .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */\n.cython .nt { color: #008000; font-weight: bold } /* Name.Tag */\n.cython .nv { color: #19177C } /* Name.Variable */\n.cython .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */\n.cython .w { color: #bbbbbb } /* Text.Whitespace */\n.cython .mb { color: #666666 } /* Literal.Number.Bin */\n.cython .mf { color: #666666 } /* Literal.Number.Float */\n.cython .mh { color: #666666 } /* Literal.Number.Hex */\n.cython .mi { color: #666666 } /* Literal.Number.Integer */\n.cython .mo { color: #666666 } /* Literal.Number.Oct */\n.cython .sa { color: #BA2121 } /* Literal.String.Affix */\n.cython .sb { color: #BA2121 } /* Literal.String.Backtick */\n.cython .sc { color: #BA2121 } /* Literal.String.Char */\n.cython .dl { color: #BA2121 } /* Literal.String.Delimiter */\n.cython .sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */\n.cython .s2 { color: #BA2121 } /* Literal.String.Double */\n.cython .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */\n.cython .sh { color: #BA2121 } /* Literal.String.Heredoc */\n.cython .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */\n.cython .sx { color: #008000 } /* Literal.String.Other */\n.cython .sr { color: #BB6688 } /* Literal.String.Regex */\n.cython .s1 { color: #BA2121 } /* Literal.String.Single */\n.cython .ss { color: #19177C } /* Literal.String.Symbol */\n.cython .bp { color: #008000 } /* Name.Builtin.Pseudo */\n.cython .fm { color: #0000FF } /* Name.Function.Magic */\n.cython .vc { color: #19177C } /* Name.Variable.Class */\n.cython .vg { color: #19177C } /* Name.Variable.Global */\n.cython .vi { color: #19177C } /* Name.Variable.Instance */\n.cython .vm { color: #19177C } /* Name.Variable.Magic */\n.cython .il { color: #666666 } /* Literal.Number.Integer.Long */\n </style>\n</head>\n<body class=\"cython\">\n<p><span style=\"border-bottom: solid 1px grey;\">Generated by Cython 0.29.15</span></p>\n<p>\n <span style=\"background-color: #FFFF00\">Yellow lines</span> hint at Python interaction.<br />\n Click on a line that starts with a \"<code>+</code>\" to see the C code that Cython generated for it.\n</p>\n<div class=\"cython\"><pre class=\"cython line score-8\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">01</span>: <span class=\"c\">#cython: boundscheck=False, wraparound=False, cdivision=True</span></pre>\n<pre class='cython code score-8 '> __pyx_t_1 = <span class='pyx_c_api'>__Pyx_PyDict_NewPresized</span>(0);<span class='error_goto'> if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 1, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);\n if (<span class='py_c_api'>PyDict_SetItem</span>(__pyx_d, __pyx_n_s_test, __pyx_t_1) &lt; 0) <span class='error_goto'>__PYX_ERR(0, 1, __pyx_L1_error)</span>\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_1); __pyx_t_1 = 0;\n</pre><pre class=\"cython line score-0\">&#xA0;<span class=\"\">02</span>: <span class=\"c\">#-#cython: profile=True, linetrace=True, binding=True</span></pre>\n<pre class=\"cython line score-0\">&#xA0;<span class=\"\">03</span>: <span class=\"c\">#-#distutils: define_macros=CYTHON_TRACE_NOGIL=1</span></pre>\n<pre class=\"cython line score-0\">&#xA0;<span class=\"\">04</span>: </pre>\n<pre class=\"cython line score-0\">&#xA0;<span class=\"\">05</span>: <span class=\"k\">cimport</span> <span class=\"nn\">numpy</span> <span class=\"k\">as</span> <span class=\"nn\">np</span></pre>\n<pre class=\"cython line score-8\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">06</span>: <span class=\"k\">import</span> <span class=\"nn\">numpy</span> <span class=\"k\">as</span> <span class=\"nn\">np</span></pre>\n<pre class='cython code score-8 '> __pyx_t_1 = <span class='pyx_c_api'>__Pyx_Import</span>(__pyx_n_s_numpy, 0, 0);<span class='error_goto'> if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 6, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);\n if (<span class='py_c_api'>PyDict_SetItem</span>(__pyx_d, __pyx_n_s_np, __pyx_t_1) &lt; 0) <span class='error_goto'>__PYX_ERR(0, 6, __pyx_L1_error)</span>\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_1); __pyx_t_1 = 0;\n</pre><pre class=\"cython line score-0\">&#xA0;<span class=\"\">07</span>: </pre>\n<pre class=\"cython line score-0\">&#xA0;<span class=\"\">08</span>: <span class=\"k\">cdef</span> <span class=\"kr\">extern</span> <span class=\"k\">from</span> <span class=\"s\">&quot;math.h&quot;</span><span class=\"p\">:</span></pre>\n<pre class=\"cython line score-0\">&#xA0;<span class=\"\">09</span>: <span class=\"n\">double</span> <span class=\"n\">fabs</span><span class=\"p\">(</span><span class=\"n\">double</span><span class=\"p\">)</span> <span class=\"k\">nogil</span></pre>\n<pre class=\"cython line score-0\">&#xA0;<span class=\"\">10</span>: <span class=\"n\">double</span> <span class=\"n\">fmax</span><span class=\"p\">(</span><span class=\"n\">double</span><span class=\"p\">,</span> <span class=\"n\">double</span><span class=\"p\">)</span> <span class=\"k\">nogil</span></pre>\n<pre class=\"cython line score-0\">&#xA0;<span class=\"\">11</span>: </pre>\n<pre class=\"cython line score-0\">&#xA0;<span class=\"\">12</span>: <span class=\"k\">from</span> <span class=\"nn\">libc.stdio</span> <span class=\"k\">cimport</span> <span class=\"n\">printf</span></pre>\n<pre class=\"cython line score-0\">&#xA0;<span class=\"\">13</span>: </pre>\n<pre class=\"cython line score-56\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">14</span>: <span class=\"k\">def</span> <span class=\"nf\">jacobian_c</span><span class=\"p\">(</span><span class=\"n\">double</span><span class=\"p\">[:,:]</span> <span class=\"n\">A</span><span class=\"p\">,</span> <span class=\"n\">double</span><span class=\"p\">[:]</span> <span class=\"n\">b</span><span class=\"p\">):</span></pre>\n<pre class='cython code score-56 '>/* Python wrapper */\nstatic PyObject *__pyx_pw_46_cython_magic_89adf65e112d455b90bc5490de5d61a1_1jacobian_c(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/\nstatic PyMethodDef __pyx_mdef_46_cython_magic_89adf65e112d455b90bc5490de5d61a1_1jacobian_c = {\"jacobian_c\", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_46_cython_magic_89adf65e112d455b90bc5490de5d61a1_1jacobian_c, METH_VARARGS|METH_KEYWORDS, 0};\nstatic PyObject *__pyx_pw_46_cython_magic_89adf65e112d455b90bc5490de5d61a1_1jacobian_c(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {\n __Pyx_memviewslice __pyx_v_A = { 0, 0, { 0 }, { 0 }, { 0 } };\n __Pyx_memviewslice __pyx_v_b = { 0, 0, { 0 }, { 0 }, { 0 } };\n PyObject *__pyx_r = 0;\n <span class='refnanny'>__Pyx_RefNannyDeclarations</span>\n <span class='refnanny'>__Pyx_RefNannySetupContext</span>(\"jacobian_c (wrapper)\", 0);\n {\n static PyObject **__pyx_pyargnames[] = {&amp;__pyx_n_s_A,&amp;__pyx_n_s_b,0};\n PyObject* values[2] = {0,0};\n if (unlikely(__pyx_kwds)) {\n Py_ssize_t kw_args;\n const Py_ssize_t pos_args = <span class='py_macro_api'>PyTuple_GET_SIZE</span>(__pyx_args);\n switch (pos_args) {\n case 2: values[1] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 1);\n CYTHON_FALLTHROUGH;\n case 1: values[0] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 0);\n CYTHON_FALLTHROUGH;\n case 0: break;\n default: goto __pyx_L5_argtuple_error;\n }\n kw_args = <span class='py_c_api'>PyDict_Size</span>(__pyx_kwds);\n switch (pos_args) {\n case 0:\n if (likely((values[0] = <span class='pyx_c_api'>__Pyx_PyDict_GetItemStr</span>(__pyx_kwds, __pyx_n_s_A)) != 0)) kw_args--;\n else goto __pyx_L5_argtuple_error;\n CYTHON_FALLTHROUGH;\n case 1:\n if (likely((values[1] = <span class='pyx_c_api'>__Pyx_PyDict_GetItemStr</span>(__pyx_kwds, __pyx_n_s_b)) != 0)) kw_args--;\n else {\n <span class='pyx_c_api'>__Pyx_RaiseArgtupleInvalid</span>(\"jacobian_c\", 1, 2, 2, 1); <span class='error_goto'>__PYX_ERR(0, 14, __pyx_L3_error)</span>\n }\n }\n if (unlikely(kw_args &gt; 0)) {\n if (unlikely(<span class='pyx_c_api'>__Pyx_ParseOptionalKeywords</span>(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, \"jacobian_c\") &lt; 0)) <span class='error_goto'>__PYX_ERR(0, 14, __pyx_L3_error)</span>\n }\n } else if (<span class='py_macro_api'>PyTuple_GET_SIZE</span>(__pyx_args) != 2) {\n goto __pyx_L5_argtuple_error;\n } else {\n values[0] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 0);\n values[1] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 1);\n }\n __pyx_v_A = <span class='pyx_c_api'>__Pyx_PyObject_to_MemoryviewSlice_dsds_double</span>(values[0], PyBUF_WRITABLE);<span class='error_goto'> if (unlikely(!__pyx_v_A.memview)) __PYX_ERR(0, 14, __pyx_L3_error)</span>\n __pyx_v_b = <span class='pyx_c_api'>__Pyx_PyObject_to_MemoryviewSlice_ds_double</span>(values[1], PyBUF_WRITABLE);<span class='error_goto'> if (unlikely(!__pyx_v_b.memview)) __PYX_ERR(0, 14, __pyx_L3_error)</span>\n }\n goto __pyx_L4_argument_unpacking_done;\n __pyx_L5_argtuple_error:;\n <span class='pyx_c_api'>__Pyx_RaiseArgtupleInvalid</span>(\"jacobian_c\", 1, 2, 2, <span class='py_macro_api'>PyTuple_GET_SIZE</span>(__pyx_args)); <span class='error_goto'>__PYX_ERR(0, 14, __pyx_L3_error)</span>\n __pyx_L3_error:;\n <span class='pyx_c_api'>__Pyx_AddTraceback</span>(\"_cython_magic_89adf65e112d455b90bc5490de5d61a1.jacobian_c\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n <span class='refnanny'>__Pyx_RefNannyFinishContext</span>();\n return NULL;\n __pyx_L4_argument_unpacking_done:;\n __pyx_r = __pyx_pf_46_cython_magic_89adf65e112d455b90bc5490de5d61a1_jacobian_c(__pyx_self, __pyx_v_A, __pyx_v_b);\n\n /* function exit code */\n <span class='refnanny'>__Pyx_RefNannyFinishContext</span>();\n return __pyx_r;\n}\n\nstatic PyObject *__pyx_pf_46_cython_magic_89adf65e112d455b90bc5490de5d61a1_jacobian_c(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_A, __Pyx_memviewslice __pyx_v_b) {\n int __pyx_v_MAX_ITER;\n int __pyx_v_N;\n double __pyx_v_epsilon;\n PyArrayObject *__pyx_v_x = 0;\n __Pyx_memviewslice __pyx_v_x0 = { 0, 0, { 0 }, { 0 }, { 0 } };\n double __pyx_v_sigma;\n double __pyx_v_delta;\n double __pyx_v_max_delta;\n CYTHON_UNUSED int __pyx_v_k;\n int __pyx_v_i;\n int __pyx_v_j;\n __Pyx_LocalBuf_ND __pyx_pybuffernd_x;\n __Pyx_Buffer __pyx_pybuffer_x;\n PyObject *__pyx_r = NULL;\n <span class='refnanny'>__Pyx_RefNannyDeclarations</span>\n <span class='refnanny'>__Pyx_RefNannySetupContext</span>(\"jacobian_c\", 0);\n __pyx_pybuffer_x.pybuffer.buf = NULL;\n __pyx_pybuffer_x.refcount = 0;\n __pyx_pybuffernd_x.data = NULL;\n __pyx_pybuffernd_x.rcbuffer = &amp;__pyx_pybuffer_x;\n/* … */\n /* function exit code */\n __pyx_L1_error:;\n <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_t_1);\n <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_t_2);\n <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_t_3);\n <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_t_4);\n __PYX_XDEC_MEMVIEW(&amp;__pyx_t_6, 1);\n { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;\n __Pyx_PyThreadState_declare\n __Pyx_PyThreadState_assign\n <span class='pyx_c_api'>__Pyx_ErrFetch</span>(&amp;__pyx_type, &amp;__pyx_value, &amp;__pyx_tb);\n <span class='pyx_c_api'>__Pyx_SafeReleaseBuffer</span>(&amp;__pyx_pybuffernd_x.rcbuffer-&gt;pybuffer);\n <span class='pyx_c_api'>__Pyx_ErrRestore</span>(__pyx_type, __pyx_value, __pyx_tb);}\n <span class='pyx_c_api'>__Pyx_AddTraceback</span>(\"_cython_magic_89adf65e112d455b90bc5490de5d61a1.jacobian_c\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n __pyx_r = NULL;\n goto __pyx_L2;\n __pyx_L0:;\n <span class='pyx_c_api'>__Pyx_SafeReleaseBuffer</span>(&amp;__pyx_pybuffernd_x.rcbuffer-&gt;pybuffer);\n __pyx_L2:;\n <span class='pyx_macro_api'>__Pyx_XDECREF</span>((PyObject *)__pyx_v_x);\n __PYX_XDEC_MEMVIEW(&amp;__pyx_v_x0, 1);\n __PYX_XDEC_MEMVIEW(&amp;__pyx_v_A, 1);\n __PYX_XDEC_MEMVIEW(&amp;__pyx_v_b, 1);\n <span class='refnanny'>__Pyx_XGIVEREF</span>(__pyx_r);\n <span class='refnanny'>__Pyx_RefNannyFinishContext</span>();\n return __pyx_r;\n}\n/* … */\n __pyx_tuple__26 = <span class='py_c_api'>PyTuple_Pack</span>(13, __pyx_n_s_A, __pyx_n_s_b, __pyx_n_s_MAX_ITER, __pyx_n_s_N, __pyx_n_s_epsilon, __pyx_n_s_x, __pyx_n_s_x0, __pyx_n_s_sigma, __pyx_n_s_delta, __pyx_n_s_max_delta, __pyx_n_s_k, __pyx_n_s_i, __pyx_n_s_j);<span class='error_goto'> if (unlikely(!__pyx_tuple__26)) __PYX_ERR(0, 14, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_tuple__26);\n <span class='refnanny'>__Pyx_GIVEREF</span>(__pyx_tuple__26);\n/* … */\n __pyx_t_1 = PyCFunction_NewEx(&amp;__pyx_mdef_46_cython_magic_89adf65e112d455b90bc5490de5d61a1_1jacobian_c, NULL, __pyx_n_s_cython_magic_89adf65e112d455b90);<span class='error_goto'> if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 14, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);\n if (<span class='py_c_api'>PyDict_SetItem</span>(__pyx_d, __pyx_n_s_jacobian_c, __pyx_t_1) &lt; 0) <span class='error_goto'>__PYX_ERR(0, 14, __pyx_L1_error)</span>\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_1); __pyx_t_1 = 0;\n __pyx_codeobj__27 = (PyObject*)<span class='pyx_c_api'>__Pyx_PyCode_New</span>(2, 0, 13, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__26, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_kvedala_ipython_cython__c, __pyx_n_s_jacobian_c, 14, __pyx_empty_bytes);<span class='error_goto'> if (unlikely(!__pyx_codeobj__27)) __PYX_ERR(0, 14, __pyx_L1_error)</span>\n</pre><pre class=\"cython line score-0\">&#xA0;<span class=\"\">15</span>: <span class=\"k\">cdef</span><span class=\"p\">:</span></pre>\n<pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">16</span>: <span class=\"nb\">int</span> <span class=\"n\">MAX_ITER</span> <span class=\"o\">=</span> <span class=\"mf\">100</span></pre>\n<pre class='cython code score-0 '> __pyx_v_MAX_ITER = 0x64;\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">17</span>: <span class=\"nb\">int</span> <span class=\"n\">N</span> <span class=\"o\">=</span> <span class=\"n\">A</span><span class=\"o\">.</span><span class=\"n\">shape</span><span class=\"p\">[</span><span class=\"mf\">0</span><span class=\"p\">]</span></pre>\n<pre class='cython code score-0 '> __pyx_v_N = (__pyx_v_A.shape[0]);\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">18</span>: <span class=\"n\">double</span> <span class=\"n\">epsilon</span> <span class=\"o\">=</span> <span class=\"mf\">1e-10</span></pre>\n<pre class='cython code score-0 '> __pyx_v_epsilon = 1e-10;\n</pre><pre class=\"cython line score-36\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">19</span>: <span class=\"n\">np</span><span class=\"o\">.</span><span class=\"n\">ndarray</span><span class=\"p\">[</span><span class=\"n\">ndim</span><span class=\"o\">=</span><span class=\"mf\">1</span><span class=\"p\">,</span><span class=\"n\">dtype</span><span class=\"o\">=</span><span class=\"n\">double</span><span class=\"p\">]</span> <span class=\"n\">x</span> <span class=\"o\">=</span> <span class=\"n\">np</span><span class=\"o\">.</span><span class=\"n\">random</span><span class=\"o\">.</span><span class=\"n\">random</span><span class=\"p\">(</span><span class=\"n\">A</span><span class=\"o\">.</span><span class=\"n\">shape</span><span class=\"p\">[</span><span class=\"mf\">1</span><span class=\"p\">])</span> <span class=\"o\">-</span> <span class=\"mf\">0.5</span></pre>\n<pre class='cython code score-36 '> <span class='pyx_c_api'>__Pyx_GetModuleGlobalName</span>(__pyx_t_2, __pyx_n_s_np);<span class='error_goto'> if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 19, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_2);\n __pyx_t_3 = <span class='pyx_c_api'>__Pyx_PyObject_GetAttrStr</span>(__pyx_t_2, __pyx_n_s_random);<span class='error_goto'> if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 19, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_3);\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_2); __pyx_t_2 = 0;\n __pyx_t_2 = <span class='pyx_c_api'>__Pyx_PyObject_GetAttrStr</span>(__pyx_t_3, __pyx_n_s_random);<span class='error_goto'> if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 19, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_2);\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_3); __pyx_t_3 = 0;\n __pyx_t_3 = <span class='py_c_api'>PyInt_FromSsize_t</span>((__pyx_v_A.shape[1]));<span class='error_goto'> if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 19, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_3);\n __pyx_t_4 = NULL;\n if (CYTHON_UNPACK_METHODS &amp;&amp; likely(<span class='py_c_api'>PyMethod_Check</span>(__pyx_t_2))) {\n __pyx_t_4 = <span class='py_macro_api'>PyMethod_GET_SELF</span>(__pyx_t_2);\n if (likely(__pyx_t_4)) {\n PyObject* function = <span class='py_macro_api'>PyMethod_GET_FUNCTION</span>(__pyx_t_2);\n <span class='pyx_macro_api'>__Pyx_INCREF</span>(__pyx_t_4);\n <span class='pyx_macro_api'>__Pyx_INCREF</span>(function);\n <span class='pyx_macro_api'>__Pyx_DECREF_SET</span>(__pyx_t_2, function);\n }\n }\n __pyx_t_1 = (__pyx_t_4) ? __Pyx_PyObject_Call2Args(__pyx_t_2, __pyx_t_4, __pyx_t_3) : <span class='pyx_c_api'>__Pyx_PyObject_CallOneArg</span>(__pyx_t_2, __pyx_t_3);\n <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_t_4); __pyx_t_4 = 0;\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_3); __pyx_t_3 = 0;\n if (unlikely(!__pyx_t_1)) <span class='error_goto'>__PYX_ERR(0, 19, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_2); __pyx_t_2 = 0;\n __pyx_t_2 = <span class='pyx_c_api'>__Pyx_PyFloat_SubtractObjC</span>(__pyx_t_1, __pyx_float_0_5, 0.5, 0, 0);<span class='error_goto'> if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 19, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_2);\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_1); __pyx_t_1 = 0;\n if (!(likely(((__pyx_t_2) == Py_None) || likely(<span class='pyx_c_api'>__Pyx_TypeTest</span>(__pyx_t_2, __pyx_ptype_5numpy_ndarray))))) <span class='error_goto'>__PYX_ERR(0, 19, __pyx_L1_error)</span>\n __pyx_t_5 = ((PyArrayObject *)__pyx_t_2);\n {\n __Pyx_BufFmt_StackElem __pyx_stack[1];\n if (unlikely(<span class='pyx_c_api'>__Pyx_GetBufferAndValidate</span>(&amp;__pyx_pybuffernd_x.rcbuffer-&gt;pybuffer, (PyObject*)__pyx_t_5, &amp;__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 1, 0, __pyx_stack) == -1)) {\n __pyx_v_x = ((PyArrayObject *)Py_None); <span class='pyx_macro_api'>__Pyx_INCREF</span>(Py_None); __pyx_pybuffernd_x.rcbuffer-&gt;pybuffer.buf = NULL;\n <span class='error_goto'>__PYX_ERR(0, 19, __pyx_L1_error)</span>\n } else {__pyx_pybuffernd_x.diminfo[0].strides = __pyx_pybuffernd_x.rcbuffer-&gt;pybuffer.strides[0]; __pyx_pybuffernd_x.diminfo[0].shape = __pyx_pybuffernd_x.rcbuffer-&gt;pybuffer.shape[0];\n }\n }\n __pyx_t_5 = 0;\n __pyx_v_x = ((PyArrayObject *)__pyx_t_2);\n __pyx_t_2 = 0;\n</pre><pre class=\"cython line score-21\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">20</span>: <span class=\"n\">double</span><span class=\"p\">[:]</span> <span class=\"n\">x0</span> <span class=\"o\">=</span> <span class=\"n\">x</span><span class=\"o\">.</span><span class=\"n\">copy</span><span class=\"p\">()</span></pre>\n<pre class='cython code score-21 '> __pyx_t_1 = <span class='pyx_c_api'>__Pyx_PyObject_GetAttrStr</span>(((PyObject *)__pyx_v_x), __pyx_n_s_copy);<span class='error_goto'> if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 20, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);\n __pyx_t_3 = NULL;\n if (CYTHON_UNPACK_METHODS &amp;&amp; likely(<span class='py_c_api'>PyMethod_Check</span>(__pyx_t_1))) {\n __pyx_t_3 = <span class='py_macro_api'>PyMethod_GET_SELF</span>(__pyx_t_1);\n if (likely(__pyx_t_3)) {\n PyObject* function = <span class='py_macro_api'>PyMethod_GET_FUNCTION</span>(__pyx_t_1);\n <span class='pyx_macro_api'>__Pyx_INCREF</span>(__pyx_t_3);\n <span class='pyx_macro_api'>__Pyx_INCREF</span>(function);\n <span class='pyx_macro_api'>__Pyx_DECREF_SET</span>(__pyx_t_1, function);\n }\n }\n __pyx_t_2 = (__pyx_t_3) ? <span class='pyx_c_api'>__Pyx_PyObject_CallOneArg</span>(__pyx_t_1, __pyx_t_3) : <span class='pyx_c_api'>__Pyx_PyObject_CallNoArg</span>(__pyx_t_1);\n <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_t_3); __pyx_t_3 = 0;\n if (unlikely(!__pyx_t_2)) <span class='error_goto'>__PYX_ERR(0, 20, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_2);\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_1); __pyx_t_1 = 0;\n __pyx_t_6 = <span class='pyx_c_api'>__Pyx_PyObject_to_MemoryviewSlice_ds_double</span>(__pyx_t_2, PyBUF_WRITABLE);<span class='error_goto'> if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 20, __pyx_L1_error)</span>\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_2); __pyx_t_2 = 0;\n __pyx_v_x0 = __pyx_t_6;\n __pyx_t_6.memview = NULL;\n __pyx_t_6.data = NULL;\n</pre><pre class=\"cython line score-0\">&#xA0;<span class=\"\">21</span>: <span class=\"c\"># double[:] x = _x</span></pre>\n<pre class=\"cython line score-0\">&#xA0;<span class=\"\">22</span>: <span class=\"c\"># double[:,:] _A = A</span></pre>\n<pre class=\"cython line score-0\">&#xA0;<span class=\"\">23</span>: <span class=\"c\"># double[:] _b = b</span></pre>\n<pre class=\"cython line score-0\">&#xA0;<span class=\"\">24</span>: <span class=\"n\">double</span> <span class=\"n\">sigma</span><span class=\"p\">,</span> <span class=\"n\">delta</span><span class=\"p\">,</span> <span class=\"n\">max_delta</span></pre>\n<pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">25</span>: <span class=\"nb\">int</span> <span class=\"n\">k</span> <span class=\"o\">=</span> <span class=\"mf\">1</span><span class=\"p\">,</span> <span class=\"n\">i</span><span class=\"p\">,</span> <span class=\"n\">j</span></pre>\n<pre class='cython code score-0 '> __pyx_v_k = 1;\n</pre><pre class=\"cython line score-0\">&#xA0;<span class=\"\">26</span>: </pre>\n<pre class=\"cython line score-4\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">27</span>: <span class=\"k\">with</span> <span class=\"k\">nogil</span><span class=\"p\">:</span></pre>\n<pre class='cython code score-4 '> {\n #ifdef WITH_THREAD\n PyThreadState *_save;\n Py_UNBLOCK_THREADS\n <span class='pyx_c_api'>__Pyx_FastGIL_Remember</span>();\n #endif\n /*try:*/ {\n/* … */\n /*finally:*/ {\n /*normal exit:*/{\n #ifdef WITH_THREAD\n <span class='pyx_c_api'>__Pyx_FastGIL_Forget</span>();\n Py_BLOCK_THREADS\n #endif\n goto __pyx_L5;\n }\n __pyx_L5:;\n }\n }\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">28</span>: <span class=\"k\">for</span> <span class=\"n\">k</span> <span class=\"ow\">in</span> <span class=\"nb\">range</span><span class=\"p\">(</span><span class=\"n\">MAX_ITER</span><span class=\"p\">):</span></pre>\n<pre class='cython code score-0 '> __pyx_t_7 = __pyx_v_MAX_ITER;\n __pyx_t_8 = __pyx_t_7;\n for (__pyx_t_9 = 0; __pyx_t_9 &lt; __pyx_t_8; __pyx_t_9+=1) {\n __pyx_v_k = __pyx_t_9;\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">29</span>: <span class=\"n\">max_delta</span> <span class=\"o\">=</span> <span class=\"mf\">0</span></pre>\n<pre class='cython code score-0 '> __pyx_v_max_delta = 0.0;\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">30</span>: <span class=\"k\">for</span> <span class=\"n\">i</span> <span class=\"ow\">in</span> <span class=\"nb\">range</span><span class=\"p\">(</span><span class=\"n\">N</span><span class=\"p\">):</span></pre>\n<pre class='cython code score-0 '> __pyx_t_10 = __pyx_v_N;\n __pyx_t_11 = __pyx_t_10;\n for (__pyx_t_12 = 0; __pyx_t_12 &lt; __pyx_t_11; __pyx_t_12+=1) {\n __pyx_v_i = __pyx_t_12;\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">31</span>: <span class=\"n\">sigma</span> <span class=\"o\">=</span> <span class=\"mf\">0.0</span></pre>\n<pre class='cython code score-0 '> __pyx_v_sigma = 0.0;\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">32</span>: <span class=\"k\">for</span> <span class=\"n\">j</span> <span class=\"ow\">in</span> <span class=\"nb\">range</span><span class=\"p\">(</span><span class=\"n\">N</span><span class=\"p\">):</span></pre>\n<pre class='cython code score-0 '> __pyx_t_13 = __pyx_v_N;\n __pyx_t_14 = __pyx_t_13;\n for (__pyx_t_15 = 0; __pyx_t_15 &lt; __pyx_t_14; __pyx_t_15+=1) {\n __pyx_v_j = __pyx_t_15;\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">33</span>: <span class=\"k\">if</span> <span class=\"n\">i</span> <span class=\"o\">!=</span> <span class=\"n\">j</span><span class=\"p\">:</span></pre>\n<pre class='cython code score-0 '> __pyx_t_16 = ((__pyx_v_i != __pyx_v_j) != 0);\n if (__pyx_t_16) {\n/* … */\n }\n }\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">34</span>: <span class=\"n\">sigma</span> <span class=\"o\">+=</span> <span class=\"n\">A</span><span class=\"p\">[</span><span class=\"n\">i</span><span class=\"p\">,</span><span class=\"n\">j</span><span class=\"p\">]</span> <span class=\"o\">*</span> <span class=\"n\">x0</span><span class=\"p\">[</span><span class=\"n\">j</span><span class=\"p\">]</span></pre>\n<pre class='cython code score-0 '> __pyx_t_17 = __pyx_v_i;\n __pyx_t_18 = __pyx_v_j;\n __pyx_t_19 = __pyx_v_j;\n __pyx_v_sigma = (__pyx_v_sigma + ((*((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_A.data + __pyx_t_17 * __pyx_v_A.strides[0]) ) + __pyx_t_18 * __pyx_v_A.strides[1]) ))) * (*((double *) ( /* dim=0 */ (__pyx_v_x0.data + __pyx_t_19 * __pyx_v_x0.strides[0]) )))));\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">35</span>: <span class=\"n\">x</span><span class=\"p\">[</span><span class=\"n\">i</span><span class=\"p\">]</span> <span class=\"o\">=</span> <span class=\"p\">(</span><span class=\"n\">b</span><span class=\"p\">[</span><span class=\"n\">i</span><span class=\"p\">]</span> <span class=\"o\">-</span> <span class=\"n\">sigma</span><span class=\"p\">)</span> <span class=\"o\">/</span> <span class=\"n\">A</span><span class=\"p\">[</span><span class=\"n\">i</span><span class=\"p\">,</span><span class=\"n\">i</span><span class=\"p\">]</span></pre>\n<pre class='cython code score-0 '> __pyx_t_20 = __pyx_v_i;\n __pyx_t_21 = __pyx_v_i;\n __pyx_t_22 = __pyx_v_i;\n __pyx_t_23 = __pyx_v_i;\n *__Pyx_BufPtrStrided1d(double *, __pyx_pybuffernd_x.rcbuffer-&gt;pybuffer.buf, __pyx_t_23, __pyx_pybuffernd_x.diminfo[0].strides) = (((*((double *) ( /* dim=0 */ (__pyx_v_b.data + __pyx_t_20 * __pyx_v_b.strides[0]) ))) - __pyx_v_sigma) / (*((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_A.data + __pyx_t_21 * __pyx_v_A.strides[0]) ) + __pyx_t_22 * __pyx_v_A.strides[1]) ))));\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">36</span>: <span class=\"n\">delta</span> <span class=\"o\">=</span> <span class=\"n\">fabs</span><span class=\"p\">(</span><span class=\"n\">x</span><span class=\"p\">[</span><span class=\"n\">i</span><span class=\"p\">]</span> <span class=\"o\">-</span> <span class=\"n\">x0</span><span class=\"p\">[</span><span class=\"n\">i</span><span class=\"p\">])</span></pre>\n<pre class='cython code score-0 '> __pyx_t_24 = __pyx_v_i;\n __pyx_t_25 = __pyx_v_i;\n __pyx_v_delta = fabs(((*__Pyx_BufPtrStrided1d(double *, __pyx_pybuffernd_x.rcbuffer-&gt;pybuffer.buf, __pyx_t_24, __pyx_pybuffernd_x.diminfo[0].strides)) - (*((double *) ( /* dim=0 */ (__pyx_v_x0.data + __pyx_t_25 * __pyx_v_x0.strides[0]) )))));\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">37</span>: <span class=\"n\">max_delta</span> <span class=\"o\">=</span> <span class=\"n\">fmax</span><span class=\"p\">(</span><span class=\"n\">delta</span><span class=\"p\">,</span> <span class=\"n\">max_delta</span><span class=\"p\">)</span></pre>\n<pre class='cython code score-0 '> __pyx_v_max_delta = fmax(__pyx_v_delta, __pyx_v_max_delta);\n }\n</pre><pre class=\"cython line score-0\">&#xA0;<span class=\"\">38</span>: <span class=\"c\"># with gil:</span></pre>\n<pre class=\"cython line score-0\">&#xA0;<span class=\"\">39</span>: <span class=\"c\"># print(k, sigma, max_delta)</span></pre>\n<pre class=\"cython line score-0\">&#xA0;<span class=\"\">40</span>: <span class=\"c\"># print(&#39;-----&#39;)</span></pre>\n<pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">41</span>: <span class=\"k\">if</span> <span class=\"n\">max_delta</span> <span class=\"o\">&lt;</span> <span class=\"n\">epsilon</span><span class=\"p\">:</span></pre>\n<pre class='cython code score-0 '> __pyx_t_16 = ((__pyx_v_max_delta &lt; __pyx_v_epsilon) != 0);\n if (__pyx_t_16) {\n/* … */\n }\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">42</span>: <span class=\"k\">break</span></pre>\n<pre class='cython code score-0 '> goto __pyx_L7_break;\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">43</span>: <span class=\"k\">for</span> <span class=\"n\">i</span> <span class=\"ow\">in</span> <span class=\"nb\">range</span><span class=\"p\">(</span><span class=\"n\">N</span><span class=\"p\">):</span></pre>\n<pre class='cython code score-0 '> __pyx_t_10 = __pyx_v_N;\n __pyx_t_11 = __pyx_t_10;\n for (__pyx_t_12 = 0; __pyx_t_12 &lt; __pyx_t_11; __pyx_t_12+=1) {\n __pyx_v_i = __pyx_t_12;\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">44</span>: <span class=\"n\">x0</span><span class=\"p\">[</span><span class=\"n\">i</span><span class=\"p\">]</span> <span class=\"o\">=</span> <span class=\"n\">x</span><span class=\"p\">[</span><span class=\"n\">i</span><span class=\"p\">]</span></pre>\n<pre class='cython code score-0 '> __pyx_t_26 = __pyx_v_i;\n __pyx_t_27 = __pyx_v_i;\n *((double *) ( /* dim=0 */ (__pyx_v_x0.data + __pyx_t_27 * __pyx_v_x0.strides[0]) )) = (*__Pyx_BufPtrStrided1d(double *, __pyx_pybuffernd_x.rcbuffer-&gt;pybuffer.buf, __pyx_t_26, __pyx_pybuffernd_x.diminfo[0].strides));\n }\n }\n __pyx_L7_break:;\n }\n</pre><pre class=\"cython line score-2\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">45</span>: <span class=\"k\">return</span> <span class=\"n\">x</span></pre>\n<pre class='cython code score-2 '> <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_r);\n <span class='pyx_macro_api'>__Pyx_INCREF</span>(((PyObject *)__pyx_v_x));\n __pyx_r = ((PyObject *)__pyx_v_x);\n goto __pyx_L0;\n</pre></div></body></html>"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "A = np.array([[2,1],[5,7]],dtype=float)\nb = np.array([11,13],dtype=float)\n\n# A = np.array([[16,3],[7,-11]],dtype=float)\n# b = np.array([11,13],dtype=float)",
"execution_count": 5,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "jacobian_c(A,b)",
"execution_count": 6,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 6,
"data": {
"text/plain": "array([ 7.11111111, -3.22222222])"
},
"metadata": {}
}
]
},
{
"metadata": {
"scrolled": true,
"trusted": true
},
"cell_type": "code",
"source": "%timeit jacobian(A, b)",
"execution_count": 7,
"outputs": [
{
"output_type": "stream",
"text": "280 µs ± 1.67 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
"name": "stdout"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "%timeit jacobian_n(A, b)",
"execution_count": 8,
"outputs": [
{
"output_type": "stream",
"text": "10.8 µs ± 6.24 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
"name": "stdout"
}
]
},
{
"metadata": {
"scrolled": true,
"trusted": true
},
"cell_type": "code",
"source": "%timeit jacobian_c(A, b)",
"execution_count": 9,
"outputs": [
{
"output_type": "stream",
"text": "7.79 µs ± 94 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
"name": "stdout"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "raw",
"source": "%lprun -f jacobian_c jacobian_c(A, b)"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "",
"execution_count": null,
"outputs": []
}
],
"metadata": {
"hide_input": false,
"kernelspec": {
"name": "python3",
"display_name": "Python 3",
"language": "python"
},
"nbTranslate": {
"hotkey": "alt-t",
"sourceLang": "en",
"targetLang": "fr",
"displayLangs": [
"*"
],
"langInMainMenu": true,
"useGoogleTranslate": true
},
"toc": {
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"base_numbering": 1,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
},
"varInspector": {
"window_display": true,
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"library": "var_list.py",
"delete_cmd_prefix": "del ",
"delete_cmd_postfix": "",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"library": "var_list.r",
"delete_cmd_prefix": "rm(",
"delete_cmd_postfix": ") ",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
]
},
"language_info": {
"name": "python",
"version": "3.7.6",
"mimetype": "text/x-python",
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"pygments_lexer": "ipython3",
"nbconvert_exporter": "python",
"file_extension": ".py"
},
"gist": {
"id": "",
"data": {
"description": "GitHub/Jacobian to solve linear system.ipynb",
"public": true
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment