tarea14.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bibliotecas\n",
    "import numpy as np\n",
    "import math\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,\n",
    "     0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5]\n",
    "v = [4.74, 5.12, 5.39, 5.49, 5.43, 5.21, 4.85,\n",
    "     4.42, 3.97, 3.56, 3.26, 3.11, 3.13, 3.31, 3.63]\n",
    "\n",
    "\n",
    "def f(x):\n",
    "    \"\"\"Funición a evaluar\"\"\"\n",
    "    a = 0\n",
    "    for i in range(0, 15):\n",
    "        aux = x[0]*math.sin(x[1]*u[i])+x[2]-v[i]\n",
    "        a += aux**2\n",
    "    return a\n",
    "\n",
    "\n",
    "def jHess(x):\n",
    "    jH = 0\n",
    "    for i in range(0, 15):\n",
    "        j = jacobian(x, i)\n",
    "        jH += np.outer(j, j)\n",
    "    return jH\n",
    "\n",
    "\n",
    "def jGrad(x):\n",
    "    jG = 0\n",
    "    for i in range(0, 15):\n",
    "        j = jacobian(x, i)\n",
    "        jG += j*fxerror(x, i)\n",
    "    return jG\n",
    "\n",
    "def fxerror(x, i):\n",
    "    return x[0]*(math.sin(x[1]*u[i]))+x[2]-v[i]\n",
    "\n",
    "\n",
    "def jacobian(x, i):\n",
    "    return np.array([math.sin(x[1]*u[i]),\n",
    "                     u[i]*x[0]*math.cos(x[1]*u[i]),\n",
    "                     1\n",
    "                     ])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "# gauss newton\n",
    "\n",
    "def gaussNewton(x0, tol=1e-4):\n",
    "    k = 0\n",
    "    print(k, x0)\n",
    "    while np.linalg.norm(jGrad(x0)) > tol:\n",
    "        a = jHess(x0)\n",
    "        s =  np.dot(np.linalg.inv(a), -jGrad(x0))\n",
    "        x0 = x0 + s\n",
    "        k = k+1\n",
    "        print(k,  x0, f(x0))\n",
    "    return x0\n",
    "\n",
    "     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 [1.  2.  1.3]\n",
      "1 [-0.93568029  3.49735647  5.39809436] 53.37613144510951\n",
      "2 [1.14538958 3.08262573 4.32987585] 3.2740392361556117\n",
      "3 [0.8072063  3.95417301 4.54756344] 2.3393236847142225\n",
      "4 [1.18849459 3.72984109 4.30197552] 0.03666322402337587\n",
      "5 [1.19729029 3.80013066 4.29694345] 0.0002045916233588216\n",
      "6 [1.19974186 3.79842054 4.29544624] 0.00010264046529553048\n"
     ]
    }
   ],
   "source": [
    "#Evaluar Gauss Newton\n",
    "x0 = np.array([1.0, 2.0, 1.3])\n",
    "# grad = jHess(x0)\n",
    "gn = gaussNewton(x0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "def jHessFull(x):\n",
    "    jH = 0\n",
    "    for i in range(0, 15):\n",
    "        j = jacobian(x, i)\n",
    "        t = tk(x, i)\n",
    "        fxe = fxerror(x, i)\n",
    "        jH += np.outer(j, j) + np.dot(t, fxe)\n",
    "    return jH\n",
    "\n",
    "def tk(x, i):\n",
    "    return np.array([\n",
    "            [0, u[i]*math.cos(x[1]*u[i]), 0],\n",
    "            [u[i]*math.cos(x[1]*u[i]), -1*x[0]*(u[i]**2)*math.sin(x[1]*u[i]), 0],\n",
    "            [0, 0, 0]\n",
    "        ])\n",
    "\n",
    "\n",
    "def levenbergMarquardtMC(x0, tol=1e-6):\n",
    "    k = 0\n",
    "    mu = 1e-4\n",
    "    v = 10\n",
    "    I = np.identity(x0.size)\n",
    "    print(k, x0, f(x0))\n",
    "    while np.linalg.norm(jGrad(x0)) >= tol:\n",
    "        mu = mu/v\n",
    "        while True:\n",
    "            delta = np.dot(np.linalg.inv(jHess(x0) + mu*I), jGrad(x0))\n",
    "            x1 = x0 - delta\n",
    "            mu = mu * v\n",
    "            # k = k+1\n",
    "            print(k,  x0, round(f(x0), 6), x1, round(f(x1), 6), mu)\n",
    "            if (f(x1) < f(x0)):\n",
    "                x0 = x1\n",
    "                break\n",
    "        k = k+1\n",
    "    return x0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 [1.  2.  1.3] 94.15754849549766\n",
      "0 [1.  2.  1.3] 94.157548 [-0.93562588  3.49733997  5.39805025] 53.372779 0.0001\n",
      "1 [-0.93562588  3.49733997  5.39805025] 53.372779 [1.14538051 3.08256099 4.32988084] 3.274627 0.0001\n",
      "2 [1.14538051 3.08256099 4.32988084] 3.274627 [0.80712493 3.95419439 4.54761399] 2.340253 0.0001\n",
      "3 [0.80712493 3.95419439 4.54761399] 2.340253 [1.18849118 3.72981044 4.30197739] 0.036695 0.0001\n",
      "4 [1.18849118 3.72981044 4.30197739] 0.036695 [1.19728804 3.80013201 4.29694481] 0.000205 0.0001\n",
      "5 [1.19728804 3.80013201 4.29694481] 0.000205 [1.19974185 3.79842054 4.29544624] 0.000103 0.0001\n",
      "6 [1.19974185 3.79842054 4.29544624] 0.000103 [1.19974328 3.79842307 4.29544531] 0.000103 0.0001\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([1.19974328, 3.79842307, 4.29544531])"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([1.0, 2.0, 1.3])\n",
    "# jHessFull(x0)\n",
    "levenbergMarquardtMC(x0)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "# newtonNDim\n",
    "def newtonNDim(x0, tol=1e-4):\n",
    "    k = 0\n",
    "    while np.linalg.norm(jGrad(x0)) > tol:\n",
    "        x0 = x0 - np.dot(np.linalg.inv(jHessFull(x0)), jGrad(x0))\n",
    "        k = k+1\n",
    "        print(k, f(x0), x0)\n",
    "    return x0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 54.017252177768384 [3.99580259 1.31260472 1.43307402]\n",
      "2 8.564631058236923 [-1.30516037  1.55750348  5.12558769]\n",
      "3 1274.6327100242845 [ 9.86327327  4.93423432 -3.2325974 ]\n",
      "4 942.0638629707175 [-0.4929472   0.52657729 12.38857608]\n",
      "5 9.84720393540005 [-0.46540121  0.50124599  4.48503307]\n",
      "6 26.760342877078983 [ 1.01592982 -0.53718744  3.59041409]\n",
      "7 5.887812686621256 [ 1.31759996 -0.69634203  4.95010312]\n",
      "8 4.385025031366485 [ 3.10649552 -0.99218937  6.10389088]\n",
      "9 3.4176903027688565 [ 2.52439716 -1.02930125  5.99144875]\n",
      "10 3.125172208826077 [ 2.70551128 -0.81162292  5.90479549]\n",
      "11 3.002247617753201 [ 2.91097729 -0.739984    5.86287679]\n",
      "12 2.9569746664736765 [ 3.22339738 -0.62415701  5.82879228]\n",
      "13 2.9220905985036865 [ 3.42279675 -0.59167454  5.82282501]\n",
      "14 2.9199622970592403 [ 3.8478207  -0.49635976  5.79664164]\n",
      "15 2.893917399085221 [ 3.99653985 -0.48904989  5.79901873]\n",
      "16 3.007176060307242 [ 4.80749608 -0.36879838  5.76958427]\n",
      "17 2.878508171965549 [ 5.01466233 -0.38004975  5.78751216]\n",
      "18 2.89095334989235 [ 4.44940979 -0.42401811  5.78356055]\n",
      "19 2.8820331258129706 [ 4.58714969 -0.41571777  5.782439  ]\n",
      "20 2.9359638232636582 [ 5.40265506 -0.33129887  5.7651284 ]\n",
      "21 2.874081627219626 [ 5.55813272 -0.33735404  5.7735395 ]\n",
      "22 2.9750451695978204 [ 4.39890199 -0.41258712  5.78043699]\n",
      "23 2.882855972247002 [ 4.56797402 -0.42075467  5.7909828 ]\n",
      "24 7.871469087362141 [ 2.0501762  -0.68411953  5.83659281]\n",
      "25 165.57752102235185 [ 8.76660878 -2.12629312 10.41463946]\n",
      "26 53.00610706563754 [ 7.10607753 -1.82540823  9.70457072]\n",
      "27 16.35852647229084 [-1.9891641  -2.02192043  2.58870305]\n",
      "28 7.73981815931943 [-0.68358595 -2.52864241  4.18269156]\n",
      "29 12.85349182645464 [ 1.03857795 -0.60923646  4.16396918]\n",
      "30 5.006562368378349 [ 1.45216754 -0.93095138  5.17541062]\n",
      "31 3.144402337055223 [ 2.70979089 -0.82994914  5.93665706]\n",
      "32 3.0095416695567536 [ 2.88454367 -0.74812459  5.86294245]\n",
      "33 2.9593906221786517 [ 3.18965388 -0.63432712  5.83187072]\n",
      "34 2.9243780842217206 [ 3.39549798 -0.59735384  5.82404089]\n",
      "35 2.9171917842938777 [ 3.79298651 -0.50686607  5.79925353]\n",
      "36 2.895355068842953 [ 3.9481618  -0.49592502  5.80010505]\n",
      "37 2.956806088184937 [ 4.65066822 -0.38871999  5.77385737]\n",
      "38 2.879746080382678 [ 4.81382189 -0.39603611  5.78508671]\n",
      "39 4.817245241081642 [ 2.71652298 -0.58860147  5.81524574]\n",
      "40 3.5850452706799802 [ 3.69571186 -0.73718632  6.18976055]\n",
      "41 2.9767203359538437 [ 3.11418619 -0.70378832  5.88463368]\n",
      "42 1121.1043027154333 [9.67138831 1.56092858 5.16206171]\n",
      "43 11.256945372394021 [-0.22382973  1.67040036  4.51888775]\n",
      "44 11.790776856652124 [0.09499013 1.96529342 4.23010561]\n",
      "45 11.843277474960784 [5.45411333e-04 1.91407185e+00 4.30619798e+00]\n",
      "46 11.84324000001304 [-8.50610872e-07  1.91498460e+00  4.30800072e+00]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-8.50610872e-07,  1.91498460e+00,  4.30800072e+00])"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([1.0, 2.0, 1.3])\n",
    "newtonNDim(x0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "def broyden(x0, tol=1e-4):\n",
    "    k = 0\n",
    "    # a = np.identity(x0.size)\n",
    "    a = jHessFull(x0)\n",
    "    while np.linalg.norm(jGrad(x0)) > tol:\n",
    "        s = np.dot(np.linalg.inv(a), -jGrad(x0))\n",
    "        x1 = x0 + s\n",
    "        y = jGrad(x1)-jGrad(x0)\n",
    "        x0 = x1\n",
    "        a = a + np.outer(y-np.dot(a, s), s)/np.dot(s, s)\n",
    "        print(k, x0, f(x0))\n",
    "        k = k+1\n",
    "        # print(a)\n",
    "    return x0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 [3.99580259 1.31260472 1.43307402] 54.017252177768384\n",
      "1 [0.12661468 2.39969539 4.3641327 ] 11.188289591093318\n",
      "2 [0.97517877 2.35773684 3.67572858] 8.12335348532393\n",
      "3 [0.82291731 2.55533109 3.92182163] 6.52009318244472\n",
      "4 [0.96472794 2.84546574 3.9601331 ] 4.006487163169713\n",
      "5 [1.59635521 4.91802645 4.70673353] 20.325730117359218\n",
      "6 [1.11476687 3.13913286 4.00338175] 2.211347463249108\n",
      "7 [1.06339624 3.39078854 4.18029776] 0.9834566334938528\n",
      "8 [1.18632157 3.64078034 4.20677046] 0.21865227184736577\n",
      "9 [ 1.14158132 11.98852512  8.3535281 ] 265.6630076453519\n",
      "10 [1.20022345 3.92630501 4.2779131 ] 0.14112998749486014\n",
      "11 [1.16197706 3.79283684 4.30672858] 0.013701155287514604\n",
      "12 [1.2047975  3.73992301 4.27164062] 0.02845506539357117\n",
      "13 [1.19454165 3.7924672  4.28816816] 0.0011665406386229464\n",
      "14 [1.19657858 3.80300768 4.29385073] 0.0004397136158992153\n",
      "15 [1.19801838 3.80066011 4.29458852] 0.00019127713685014125\n",
      "16 [1.19967131 3.79846254 4.29538823] 0.00010275926514522527\n",
      "17 [1.1997263  3.79844572 4.29543703] 0.00010264918015628212\n",
      "18 [1.19975064 3.79841332 4.29544892] 0.00010264203169308395\n",
      "19 [1.19974329 3.79842306 4.29544531] 0.00010264038551143695\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([1.19974329, 3.79842306, 4.29544531])"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([1.0, 2.0, 1.3])\n",
    "broyden(x0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "#descenso de gradiente\n",
    "\n",
    "def dirgrad(x):\n",
    "    vgrad = jGrad(x)\n",
    "    magGrad = np.sqrt(vgrad.dot(vgrad))\n",
    "    p = -vgrad/magGrad\n",
    "    return p\n",
    "\n",
    "\n",
    "def phiAlpha(x0, alpha, p):\n",
    "    paX = x0 + p * alpha\n",
    "    return f(paX)\n",
    "\n",
    "\n",
    "def phipAlpha(x0, alpha, p):\n",
    "    x = x0 + alpha * p\n",
    "    vgrad = jGrad(x)\n",
    "    return (np.dot(vgrad, p))\n",
    "\n",
    "\n",
    "def phipp(x0, alpha, p):\n",
    "    x = x0 + alpha * p\n",
    "    ahess = jHessFull(x)\n",
    "    return np.dot(np.dot(ahess, p), p)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def exhaustivoRefinado(p, xini, alpha=0, h=0.1, tol=1e-9):\n",
    "    \"\"\"Busqueda de minimo con metodo exhaustivo refinado. puedes cambiar el paso\n",
    "    Retorna f(a) y alpha\n",
    "    \"\"\"\n",
    "    k = 0\n",
    "    while h > tol:\n",
    "        while phiAlpha(xini, alpha+h, p) < phiAlpha(xini, alpha, p):\n",
    "            alpha = alpha + h\n",
    "            fnow = phiAlpha(xini, alpha, p)\n",
    "            # print(k, h, fnow)\n",
    "            k += 1\n",
    "        alpha = alpha-h\n",
    "        h = h / 10\n",
    "    return alpha\n",
    "\n",
    "\n",
    "def gradDescent(x0, k=0, tol=1e-6):\n",
    "    print(\"k, x^(k), f(x^(k), pk, θ,\")\n",
    "    op = dirgrad(x0)\n",
    "    while np.linalg.norm(jGrad(x0)) >= tol:\n",
    "        p = dirgrad(x0)\n",
    "        alpha = exhaustivoRefinado(p, x0)\n",
    "        # print(f\"a: {alpha}\")\n",
    "        x0 = x0 + alpha*p\n",
    "        if k >= 1:\n",
    "            angulo = np.arccos(np.dot(op, p))\n",
    "            op = p\n",
    "            print(f\"{k}, {x0}, {f(x0)}, {p},  {round(np.degrees(angulo),6)} \")\n",
    "        else:\n",
    "            print(f\"{k}, {x0}, {f(x0)}, {p},  - \")\n",
    "\n",
    "        k = k+1\n",
    "    return x0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k, x^(k), f(x^(k), pk, θ,\n",
      "0, [1.94914572 1.79160024 2.76408284], 18.74744377395421, [ 0.54013837 -0.11859581  0.83317799],  - \n",
      "1, [1.88687047 2.38495533 2.88891415], 11.5212870332878, [-0.10216867  0.97345741  0.20479804],  89.999998 \n",
      "2, [1.95980422 2.33139573 3.1798815 ], 9.384534187552484, [ 0.23935232 -0.17577069  0.95489011],  89.999998 \n",
      "3, [1.88242504 2.5826336  3.24552378], 8.294732937054158, [-0.28557885  0.92722908  0.24226217],  89.999994 \n",
      "4, [1.83521421 2.51673855 3.44207696], 7.397429706178855, [-0.22205081 -0.30992996  0.92446572],  89.999999 \n",
      "5, [1.59940207 2.84449563 3.49531798], 6.350047855403266, [-0.57901071  0.80477135  0.13072746],  89.999994 \n",
      "6, [1.54715709 2.764177   3.75836703], 4.889099039147942, [-0.18661846 -0.28689721  0.93960819],  89.999998 \n",
      "7, [1.04683873 3.42038571 3.85936219], 2.8392608494774683, [-0.60182065  0.78933732  0.12148459],  89.999995 \n",
      "8, [1.13580647 3.43235639 4.22231967], 0.7577182034806647, [0.23794913 0.03201624 0.9707498 ],  89.999997 \n",
      "9, [1.15947276 3.76005127 4.20571096], 0.12904806079908762, [ 0.07194104  0.99613027 -0.05048733],  89.999995 \n",
      "10, [1.18108238 3.76292369 4.29317682], 0.010771012348283978, [0.23972977 0.03186565 0.97031655],  89.999997 \n",
      "11, [1.19726094 3.79523746 4.2881185 ], 0.0009329899329476622, [ 0.44337091  0.88555408 -0.13862251],  89.999986 \n",
      "12, [1.19853878 3.79571744 4.29527175], 0.00016221706478745037, [0.1754723  0.06590977 0.98227561],  89.999973 \n",
      "13, [1.19951171 3.79821804 4.29493016], 0.00010690901018059429, [ 0.35969337  0.92448488 -0.12628691],  89.999985 \n",
      "14, [1.19963518 3.79823929 4.29543737], 0.00010294917843339105, [0.23633963 0.04066337 0.97081928],  89.99996 \n",
      "15, [1.19972942 3.79840642 4.29540743], 0.00010266281886255482, [ 0.48529545  0.86064868 -0.1541823 ],  89.999526 \n",
      "16, [1.19973661 3.79840902 4.29544452], 0.00010264203143953922, [0.18966241 0.06853676 0.97945438],  89.999207 \n",
      "17, [1.19974208 3.79842193 4.29544256], 0.00010264050649110953, [ 0.38647842  0.91185476 -0.1384028 ],  89.986434 \n",
      "18, [1.19974272 3.79842207 4.29544527], 0.00010264039447382076, [0.22878328 0.05100538 0.97214025],  89.978092 \n",
      "19, [1.1997432  3.79842298 4.29544511], 0.00010264038617197361, [ 0.46397918  0.87237992 -0.15387202],  89.939177 \n",
      "20, [1.19974324 3.79842299 4.29544531], 0.0001026403855596459, [0.20111919 0.06944013 0.97710242],  89.796911 \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([1.19974324, 3.79842299, 4.29544531])"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([1.0, 2.0, 1.3])\n",
    "gradDescent(x0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gD(x0):\n",
    "    p = dirgrad(x0)\n",
    "    alpha = exhaustivoRefinado(p, x0)\n",
    "    # TODO: buscar alpha con newton para mayor precisión ?\n",
    "    x0 = x0 + alpha*p\n",
    "    return x0\n",
    "\n",
    "\n",
    "def forsyte(x0, k=0, m=0, tol=1e-2):\n",
    "    \"\"\"Algoritmo de forsyte.\"\"\"\n",
    "    print(\"k, x ^ (k), f(x ^ (k), pk\")\n",
    "    while np.linalg.norm(jGrad(x0)) >= tol:\n",
    "        x1 = gD(x0)\n",
    "        x2 = gD(x1)\n",
    "        x3 = gD(x2)\n",
    "        y = x3\n",
    "        d = (y - x0)/np.linalg.norm(y - x0)\n",
    "        alpha = exhaustivoRefinado(d, x0)\n",
    "        # TODO: buscar alpha con newton para mayor precisión ?\n",
    "        # print(f\"alpha: {alpha}\")\n",
    "        x0 = x0 + alpha*d\n",
    "        # itTime = timer()\n",
    "        print(f\"{k}, {x0}, {f(x0)}, {d} \")\n",
    "        k = k + 1\n",
    "    return x0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k, x ^ (k), f(x ^ (k), pk\n",
      "0, [1.98768275 2.34102147 3.23448466], 9.323111447940448, [0.4492232  0.15510523 0.87985276] \n",
      "1, [0.92034716 3.72185805 4.0734375 ], 1.3464740630673415, [-0.5511876   0.71308407  0.43324743] \n",
      "2, [1.19160467 3.7793176  4.29915715], 0.0035638437979520383, [0.75868704 0.16071008 0.63132103] \n",
      "3, [1.19929015 3.79788801 4.29452828], 0.00011683237100980668, [ 0.3726457   0.90042314 -0.22444008] \n",
      "4, [1.19972604 3.79838614 4.29545396], 0.00010265425538361158, [0.38303801 0.4377257  0.81343598] \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([1.19972604, 3.79838614, 4.29545396])"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([1.0, 2.0, 1.3])\n",
    "forsyte(x0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gradientes Conjugados NL\n",
    "def newton(xo, p, ao, itmax=100, tol=1e-12):\n",
    "    k = 0\n",
    "    # ak = 0\n",
    "    while abs(phipAlpha(xo, ao, p)) > tol:\n",
    "        phiap = phipAlpha(xo, ao, p)\n",
    "        phiapp = phipp(xo, ao, p)\n",
    "        ak = ao - phiap/phiapp\n",
    "        # print(f\"\\t {k}, {ak}\")\n",
    "        k = k+1\n",
    "        ao = ak\n",
    "        if k >= itmax:\n",
    "            print(\"Iteraciones exedidas\")\n",
    "            break\n",
    "    return ak\n",
    "\n",
    "\n",
    "def gradientesConjugados(x0, flavor=\"FR\", k=0, tol=1e-6):\n",
    "    beta = 0\n",
    "    p = -jGrad(x0)\n",
    "    print(\"k, fx, x0, alpha, beta, \")\n",
    "    while np.linalg.norm(jGrad(x0)) >= tol:\n",
    "        alpha = newton(x0, p, 0)\n",
    "        xi = x0 + alpha*p\n",
    "        gradxi = jGrad(xi)\n",
    "        gradx0 = jGrad(x0)\n",
    "        if flavor == \"FR\":\n",
    "            beta = np.dot(gradxi, gradxi)/np.dot(gradx0, gradx0)\n",
    "        elif flavor == \"PR\":\n",
    "            beta = np.dot(gradxi, (gradxi - gradx0)) / \\\n",
    "                (np.linalg.norm(gradx0))**2\n",
    "        elif flavor == \"PR+\":\n",
    "            beta = np.dot(gradxi, (gradxi - gradx0)) / \\\n",
    "                (np.linalg.norm(gradx0))**2\n",
    "            if beta < 0:\n",
    "                beta = 0\n",
    "        elif flavor == \"HS\":\n",
    "            beta = np.dot(gradxi, (gradxi - gradx0)) / \\\n",
    "                np.dot((gradxi - gradx0), p)\n",
    "        p = -gradxi + beta*p\n",
    "        print(k, f(x0), xi, alpha, beta)\n",
    "        x0 = xi\n",
    "        k += 1\n",
    "    return x0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k, fx, x0, alpha, beta, \n",
      "0 94.15754849549766 [1.94914573 1.79160024 2.76408285] 0.04167297878840474 0.060878374343717706\n",
      "1 18.747443773954203 [2.02662322 4.14362042 3.78633343] 0.23942760378580574 1.0382856625602845\n",
      "2 10.399532329801419 [1.8379291  4.36389756 4.17489594] 0.03155052291676524 0.7463486246789234\n",
      "3 6.822440548056421 [1.51220548 4.29250593 4.53774657] 0.033573948870974674 0.6738290215457408\n",
      "4 3.994497845539156 [1.22713232 4.06238271 4.65127131] 0.03087289892251813 0.5548701086255411\n",
      "5 2.3052148287807377 [1.00300198 3.7866028  4.52150219] 0.04259656196489717 0.4461148244815673\n",
      "6 1.0774794719304255 [0.97069997 3.67944351 4.32426414] 0.04161969424279546 0.2816530262639268\n",
      "7 0.5032591487096914 [1.05202087 3.674519   4.22127262] 0.05233542691105991 0.5734709026011945\n",
      "8 0.29711358788433234 [1.18352891 3.70290328 4.20418379] 0.06517048905782614 0.660791964289599\n",
      "9 0.14986331399089345 [1.25610216 3.73937625 4.25497767] 0.052075024906906485 0.4277116631468989\n",
      "10 0.07205523071623406 [1.26009922 3.77610085 4.30150843] 0.05297491628932087 0.5560401310266453\n",
      "11 0.037876732598443104 [1.23360447 3.81090862 4.32023823] 0.055166842084049345 0.5554642331586579\n",
      "12 0.018171391017569573 [1.20852238 3.82655823 4.31257721] 0.04674610698721564 0.38768140276028357\n",
      "13 0.008977753933909033 [1.19619364 3.82421676 4.29965077] 0.048041691178374125 0.5269570847945885\n",
      "14 0.005303246666887423 [1.19093581 3.81041192 4.28981356] 0.0630045373919861 0.7679846841845113\n",
      "15 0.002776392705487479 [1.19200353 3.79567579 4.28922515] 0.05302866841463807 0.3810710087843801\n",
      "16 0.0011399980393954907 [1.19492317 3.79128393 4.29277626] 0.04184569396054978 0.31440575168082363\n",
      "17 0.0006460974065228107 [1.19802515 3.79187769 4.29576639] 0.056194446967369666 0.7270188332421548\n",
      "18 0.0004377262364834571 [1.20130201 3.79570051 4.29720569] 0.06839693281192386 0.7155622784723503\n",
      "19 0.00025313305738510253 [1.20216118 3.79873613 4.29641787] 0.046309074688708846 0.31074210866906177\n",
      "20 0.00016376792448722827 [1.20154197 3.79972062 4.29554856] 0.044274651872627864 0.41258577287199716\n",
      "21 0.00013723184972340812 [1.20035617 3.79979894 4.29509697] 0.061277541539271056 0.7495903535064209\n",
      "22 0.00012207661674596204 [1.19931412 3.79925197 4.29522984] 0.0572994408634437 0.5051689536604583\n",
      "23 0.00011145413944798035 [1.19908698 3.7986919  4.2955168 ] 0.0469941428692301 0.40989235790812584\n",
      "24 0.00010705302681774116 [1.19928729 3.79828745 4.29562965] 0.05465964929541966 0.5887306726592832\n",
      "25 0.00010495467763943222 [1.1995997  3.79810265 4.2955348 ] 0.05431577083553281 0.4694423303609594\n",
      "26 0.00010372712257374668 [1.19976699 3.79814316 4.29540455] 0.04696421354751994 0.4181660962935094\n",
      "27 0.00010322885432030028 [1.19982681 3.79828298 4.29534056] 0.0578641329450995 0.7109304282399471\n",
      "28 0.00010297213566833621 [1.19982257 3.79844412 4.29537067] 0.06084878179208638 0.557459514447771\n",
      "29 0.00010278020480990309 [1.19979357 3.79849887 4.2954333 ] 0.043964577910472785 0.2956363962966918\n",
      "30 0.00010270290233145304 [1.19976854 3.79849044 4.29546809] 0.046809149679536116 0.5005846457282779\n",
      "31 0.0001026785705529929 [1.19974023 3.79845429 4.29547898] 0.06822993109305887 0.8768415285308672\n",
      "32 0.00010266081666228658 [1.1997229  3.79841917 4.29546227] 0.05514755865997308 0.41620984924073323\n",
      "33 0.00010264823429999872 [1.19972418 3.79840992 4.29544619] 0.04176768717164093 0.30605399116223003\n",
      "34 0.0001026442679053047 [1.19973231 3.79841072 4.29543805] 0.05286614495799619 0.6354657434996448\n",
      "35 0.000102642731407581 [1.19974386 3.7984158  4.29543729] 0.06512304806790695 0.6975417478665766\n",
      "36 0.00010264152864834488 [1.19974942 3.79842042 4.29544209] 0.0490819501688343 0.38216418116837564\n",
      "37 0.0001026408963307581 [1.19974897 3.79842309 4.29544592] 0.04823471451558744 0.4764916136812723\n",
      "38 0.00010264065885317713 [1.19974601 3.79842499 4.29544753] 0.05763787887704425 0.5931794827925861\n",
      "39 0.00010264052363759499 [1.19974342 3.7984256  4.29544685] 0.05037926437206052 0.419719186811674\n",
      "40 0.00010264045353121046 [1.19974238 3.79842505 4.29544565] 0.050421972217695205 0.5390809326702103\n",
      "41 0.00010264042408125207 [1.19974226 3.79842374 4.29544478] 0.06249888695631179 0.7000647869630585\n",
      "42 0.00010264040440279846 [1.19974266 3.79842264 4.29544476] 0.05063448156601354 0.3665280995273554\n",
      "43 0.00010264039324178274 [1.19974298 3.79842237 4.29544508] 0.042952017568219134 0.35681697738983414\n",
      "44 0.00010264038977163102 [1.19974325 3.7984225  4.29544536] 0.058971874096517524 0.7804510073643888\n",
      "45 0.00010264038807160657 [1.19974346 3.79842291 4.29544548] 0.06553675946212835 0.6365031856984068\n",
      "46 0.00010264038659712402 [1.19974348 3.79842316 4.2954454 ] 0.04432395101571281 0.2880331009866186\n",
      "47 0.00010264038596238502 [1.19974341 3.79842322 4.29544533] 0.04490722050529163 0.4438667939425934\n",
      "48 0.00010264038577715696 [1.19974331 3.7984232  4.29544528] 0.06470843878205525 0.8208184197411821\n",
      "49 0.00010264038565868185 [1.19974323 3.79842312 4.29544529] 0.057006736997202656 0.47700127455299574\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([1.19974323, 3.79842312, 4.29544529])"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([1.0, 2.0, 1.3])\n",
    "gradientesConjugados(x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [],
   "source": [
    "#bfgs\n",
    "def bfgs(x0, tol=1e-6):\n",
    "    k = 0\n",
    "    # a = np.identity(x0.size)\n",
    "    a = jHessFull(x0)\n",
    "    while np.linalg.norm(jGrad(x0)) > tol:\n",
    "        s = np.dot(np.linalg.inv(a), -jGrad(x0))\n",
    "        x1 = x0 + s\n",
    "        y = jGrad(x1)-jGrad(x0)\n",
    "        a = a + (np.dot(y, y)/np.dot(y, s)) - \\\n",
    "            np.dot(np.dot(a, s), np.dot(a, s))/np.dot(np.dot(a, s), s)\n",
    "        x0 = x1\n",
    "        print(k, x0, f(x0))\n",
    "        k = k+1\n",
    "        # print(a)\n",
    "    return x0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 [3.99580259 1.31260472 1.43307402] 54.017252177768384\n",
      "1 [18.79081704 -2.81902499 -9.74084996] 8253.61395184246\n",
      "2 [ 82.22103485 -57.15376288 -19.07819159] 51405.274630641215\n",
      "3 [ 2137.96603235  -576.67781623 -1555.20948726] 76975123.40641178\n",
      "4 [-1101057.36698927   245827.62368476   855381.31390444] 20358593481069.32\n",
      "5 [-6.03398909e+11  1.36580658e+11  4.66818202e+11] 5.94602835341919e+24\n",
      "6 [-1.27215841e+23  2.87952961e+22  9.84205450e+22] 2.089188635863388e+47\n"
     ]
    },
    {
     "ename": "LinAlgError",
     "evalue": "Singular matrix",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mLinAlgError\u001b[0m                               Traceback (most recent call last)",
      "Cell \u001b[1;32mIn [109], line 2\u001b[0m\n\u001b[0;32m      1\u001b[0m x0 \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39marray([\u001b[39m1.0\u001b[39m, \u001b[39m2.0\u001b[39m, \u001b[39m1.3\u001b[39m])\n\u001b[1;32m----> 2\u001b[0m bfgs(x0)\n",
      "Cell \u001b[1;32mIn [108], line 7\u001b[0m, in \u001b[0;36mbfgs\u001b[1;34m(x0, tol)\u001b[0m\n\u001b[0;32m      5\u001b[0m a \u001b[39m=\u001b[39m jHessFull(x0)\n\u001b[0;32m      6\u001b[0m \u001b[39mwhile\u001b[39;00m np\u001b[39m.\u001b[39mlinalg\u001b[39m.\u001b[39mnorm(jGrad(x0)) \u001b[39m>\u001b[39m tol:\n\u001b[1;32m----> 7\u001b[0m     s \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mdot(np\u001b[39m.\u001b[39;49mlinalg\u001b[39m.\u001b[39;49minv(a), \u001b[39m-\u001b[39mjGrad(x0))\n\u001b[0;32m      8\u001b[0m     x1 \u001b[39m=\u001b[39m x0 \u001b[39m+\u001b[39m s\n\u001b[0;32m      9\u001b[0m     y \u001b[39m=\u001b[39m jGrad(x1)\u001b[39m-\u001b[39mjGrad(x0)\n",
      "File \u001b[1;32m<__array_function__ internals>:180\u001b[0m, in \u001b[0;36minv\u001b[1;34m(*args, **kwargs)\u001b[0m\n",
      "File \u001b[1;32mc:\\Users\\Daniel\\AppData\\Local\\Programs\\Python\\Python38-32\\lib\\site-packages\\numpy\\linalg\\linalg.py:552\u001b[0m, in \u001b[0;36minv\u001b[1;34m(a)\u001b[0m\n\u001b[0;32m    550\u001b[0m signature \u001b[39m=\u001b[39m \u001b[39m'\u001b[39m\u001b[39mD->D\u001b[39m\u001b[39m'\u001b[39m \u001b[39mif\u001b[39;00m isComplexType(t) \u001b[39melse\u001b[39;00m \u001b[39m'\u001b[39m\u001b[39md->d\u001b[39m\u001b[39m'\u001b[39m\n\u001b[0;32m    551\u001b[0m extobj \u001b[39m=\u001b[39m get_linalg_error_extobj(_raise_linalgerror_singular)\n\u001b[1;32m--> 552\u001b[0m ainv \u001b[39m=\u001b[39m _umath_linalg\u001b[39m.\u001b[39;49minv(a, signature\u001b[39m=\u001b[39;49msignature, extobj\u001b[39m=\u001b[39;49mextobj)\n\u001b[0;32m    553\u001b[0m \u001b[39mreturn\u001b[39;00m wrap(ainv\u001b[39m.\u001b[39mastype(result_t, copy\u001b[39m=\u001b[39m\u001b[39mFalse\u001b[39;00m))\n",
      "File \u001b[1;32mc:\\Users\\Daniel\\AppData\\Local\\Programs\\Python\\Python38-32\\lib\\site-packages\\numpy\\linalg\\linalg.py:89\u001b[0m, in \u001b[0;36m_raise_linalgerror_singular\u001b[1;34m(err, flag)\u001b[0m\n\u001b[0;32m     88\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39m_raise_linalgerror_singular\u001b[39m(err, flag):\n\u001b[1;32m---> 89\u001b[0m     \u001b[39mraise\u001b[39;00m LinAlgError(\u001b[39m\"\u001b[39m\u001b[39mSingular matrix\u001b[39m\u001b[39m\"\u001b[39m)\n",
      "\u001b[1;31mLinAlgError\u001b[0m: Singular matrix"
     ]
    }
   ],
   "source": [
    "x0 = np.array([1.0, 2.0, 1.3])\n",
    "bfgs(x0)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.1 32-bit",
   "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.8.1"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "dbd870ce3ecbb50051e36c03c034f91bd0ffd37935fd5f5a0057175bdb6e15f6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}