Scaling example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import MultivariatePolynomials\n",
    "const MP = MultivariatePolynomials\n",
    "using DynamicPolynomials\n",
    "using HomotopyContinuation\n",
    "using JLD2\n",
    "using DataFrames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1-element Array{Symbol,1}:\n",
       " :f"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@load \"HC_result_dim3.jld2\" f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "munorm_grad (generic function with 1 method)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function weylnorm(f::MP.AbstractPolynomialLike)\n",
    "    sqrt(sum(MP.terms(f)) do t\n",
    "        exps = MP.exponents(t)\n",
    "        c = MP.coefficient(t)\n",
    "        c * conj(c) * prod(factorial, exps) / factorial(sum(exps))\n",
    "    end)\n",
    "end\n",
    "\n",
    "function weylnorm(f::Vector{<:MP.AbstractPolynomialLike})\n",
    "    sqrt(sum(fi -> weylnorm(fi)^2, f))\n",
    "end\n",
    "\n",
    "\n",
    "function jac(F::Vector, x)\n",
    "    vars = MP.variables(F)\n",
    "    [MP.differentiate(fi, vars[i])(vars=>x) for fi in F, i in 1:length(vars)]\n",
    "end\n",
    "\n",
    "function jac(F::Vector)\n",
    "    vars = MP.variables(F)\n",
    "    [MP.differentiate(fi, vars[i]) for fi in F, i in 1:length(vars)]\n",
    "end\n",
    "\n",
    "function munorm(F::Vector, x)\n",
    "    d = diagm(map(F) do fi\n",
    "        di = MP.maxdegree(fi)\n",
    "        sqrt(di) * norm(x)^(di-1)\n",
    "    end)\n",
    "    sqrt(weylnorm(F) * norm(pinv(jac(F, x)) * d))\n",
    "end\n",
    "\n",
    "function scale(F, w, λ)\n",
    "    vars = variables(F)\n",
    "    map(F) do p\n",
    "        MP.subs(p, vars => inv.(λ) .* vars)\n",
    "    end, λ .* w\n",
    "end\n",
    "\n",
    "function scale_munorm(F, w, λ)\n",
    "    vars = variables(F)\n",
    "    munorm(map(F) do p\n",
    "        MP.subs(p, vars => inv.(λ) .* vars)\n",
    "    end, λ .* w)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8-element Array{DynamicPolynomials.Polynomial{true,Float64},1}:\n",
       " 1.00663l1p1^3 + 0.926589l1p1^2p2 + 0.752524l1p1^2p3 + 0.756662l1p1^2p4 + 0.376208l1p1p2^2 - 0.140268l1p1p2p3 + 0.577388l1p1p2p4 - 0.473155l1p1p3^2 - 0.939223l1p1p3p4 - 0.215605l1p1p4^2 - 1.81154l3p1^3 - 0.194997l3p1^2p2 - 1.81281l3p1^2p3 + 1.06788l3p1^2p4 + 0.104512l3p1p2^2 + 0.83405l3p1p2p3 + 0.820324l3p1p2p4 + 0.45035l3p1p3^2 - 0.00748252l3p1p3p4 + 0.372711l3p1p4^2 - 0.118499l2p1^2 + 0.757563l2p1p2 - 0.51333l2p1p3 + 0.283454l2p1p4 - l0p1 + 51.0  \n",
       " 0.463295l1p1^2p2 + 0.752416l1p1p2^2 - 0.140268l1p1p2p3 + 0.577388l1p1p2p4 - 1.93474l1p2^3 + 0.336938l1p2^2p3 - 1.42015l1p2^2p4 + 0.722479l1p2p3^2 + 0.648254l1p2p3p4 - 0.163331l1p2p4^2 - 0.0974987l3p1^2p2 + 0.209024l3p1p2^2 + 0.83405l3p1p2p3 + 0.820324l3p1p2p4 + 2.83552l3p2^3 + 0.676585l3p2^2p3 + 1.89576l3p2^2p4 - 0.0697744l3p2p3^2 + 0.333876l3p2p3p4 + 0.879106l3p2p4^2 + 0.757563l2p1p2 + 1.05544l2p2^2 + 0.759711l2p2p3 - 0.0636985l2p2p4 - l0p2 + 59.0\n",
       " 0.376262l1p1^2p3 - 0.140268l1p1p2p3 - 0.94631l1p1p3^2 - 0.939223l1p1p3p4 + 0.168469l1p2^2p3 + 1.44496l1p2p3^2 + 0.648254l1p2p3p4 + 0.534104l1p3^3 - 0.904022l1p3^2p4 + 0.496171l1p3p4^2 - 0.906407l3p1^2p3 + 0.83405l3p1p2p3 + 0.9007l3p1p3^2 - 0.00748252l3p1p3p4 + 0.338293l3p2^2p3 - 0.139549l3p2p3^2 + 0.333876l3p2p3p4 - 2.19678l3p3^3 + 0.236926l3p3^2p4 + 0.40213l3p3p4^2 - 0.51333l2p1p3 + 0.759711l2p2p3 - 0.660021l2p3^2 - 0.567906l2p3p4 - l0p3 + 36.0   \n",
       " 0.378331l1p1^2p4 + 0.577388l1p1p2p4 - 0.939223l1p1p3p4 - 0.43121l1p1p4^2 - 0.710074l1p2^2p4 + 0.648254l1p2p3p4 - 0.326662l1p2p4^2 - 0.452011l1p3^2p4 + 0.992342l1p3p4^2 + 1.01696l1p4^3 + 0.53394l3p1^2p4 + 0.820324l3p1p2p4 - 0.00748252l3p1p3p4 + 0.745423l3p1p4^2 + 0.947882l3p2^2p4 + 0.333876l3p2p3p4 + 1.75821l3p2p4^2 + 0.118463l3p3^2p4 + 0.804259l3p3p4^2 - 1.08092l3p4^3 + 0.283454l2p1p4 - 0.0636985l2p2p4 - 0.567906l2p3p4 + 1.46964l2p4^2 - l0p4 + 41.0\n",
       " -0.335542p1^3 - 0.463295p1^2p2 - 0.376262p1^2p3 - 0.378331p1^2p4 - 0.376208p1p2^2 + 0.140268p1p2p3 - 0.577388p1p2p4 + 0.473155p1p3^2 + 0.939223p1p3p4 + 0.215605p1p4^2 + 0.644914p2^3 - 0.168469p2^2p3 + 0.710074p2^2p4 - 0.722479p2p3^2 - 0.648254p2p3p4 + 0.163331p2p4^2 - 0.178035p3^3 + 0.452011p3^2p4 - 0.496171p3p4^2 - 0.338985p4^3                                                                                                                          \n",
       " 0.0592493p1^2 - 0.757563p1p2 + 0.51333p1p3 - 0.283454p1p4 - 0.527722p2^2 - 0.759711p2p3 + 0.0636985p2p4 + 0.33001p3^2 + 0.567906p3p4 - 0.734821p4^2                                                                                                                                                                                                                                                                                                                 \n",
       " 0.603846p1^3 + 0.0974987p1^2p2 + 0.906407p1^2p3 - 0.53394p1^2p4 - 0.104512p1p2^2 - 0.83405p1p2p3 - 0.820324p1p2p4 - 0.45035p1p3^2 + 0.00748252p1p3p4 - 0.372711p1p4^2 - 0.945175p2^3 - 0.338293p2^2p3 - 0.947882p2^2p4 + 0.0697744p2p3^2 - 0.333876p2p3p4 - 0.879106p2p4^2 + 0.732261p3^3 - 0.118463p3^2p4 - 0.40213p3p4^2 + 0.360306p4^3                                                                                                                           \n",
       " -p1 - p2 - p3 - p4 + 1.0                                                                                                                                                                                                                                                                                                                                                                                                                                            "
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = f[:input_system][500]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "μ_s: 107239.0765186863\n",
      "s: Complex{Float64}[1.0+0.0im, 187.0+7.46203e-14im, 143.381+43.9299im, 130.967+16.0205im, 119.802+96.0335im, 0.126401+1.26771im, 0.27151-0.963521im, 0.176132-0.468405im, 0.425957+0.164214im]\n"
     ]
    }
   ],
   "source": [
    "result = solve(g, report_progress=false)\n",
    "r1 = finite(result)[1]\n",
    "g_hom = homogenize(g)\n",
    "s = [1; r1.solution]\n",
    "\n",
    "μ_s = munorm(g_hom, s)\n",
    "\n",
    "println(\"μ_s: $μ_s\")\n",
    "println(\"s: $s\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "munorm_grad (generic function with 1 method)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compute the gradient by finite differences since working out the derivative is really ugly\n",
    "function munorm_grad(F, w; epsilon=1e-12)\n",
    "    map(1:length(w)) do i\n",
    "        s = ones(length(w))\n",
    "        s[i] += epsilon\n",
    "        m1 = scale_munorm(F, w, s)\n",
    "        s = ones(length(w))\n",
    "        s[i] -= epsilon\n",
    "        m2 = scale_munorm(F, w, s)\n",
    "\n",
    "        (m2 .- m1) ./ (2*epsilon)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8552.50535422097"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "∇s = munorm_grad(g_hom, s)\n",
    "# ∇s is quite large, let this scale by the current norm\n",
    "# Not sure yet that this is \n",
    "g_hom2, s2 = scale(g_hom, s, 1 .+ ∇s ./ μ_s)\n",
    "\n",
    "μ_s2 = munorm(g_hom2, s2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "μ_s3: 3299.5198271093363\n",
      "s3: Complex{Float64}[3.20885+0.0im, 54.1821+2.16208e-14im, 59.3841+18.1944im, 34.0198+4.16146im, 54.6657+43.8201im, 0.159938+1.60406im, 0.42225-1.49846im, 0.226657-0.60277im, 0.512818+0.1977im]\n"
     ]
    }
   ],
   "source": [
    "# lets do a second step\n",
    "∇s2 = munorm_grad(g_hom2, s2)\n",
    "# ∇s2 is quite large, let this scale by the current norm\n",
    "g_hom3, s3 = scale(g_hom2, s2, 1 .+ ∇s2 ./ μ_s2)\n",
    "\n",
    "μ_s3 = munorm(g_hom3, s3)\n",
    "println(\"μ_s3: $μ_s3\")\n",
    "println(\"s3: $s3\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8-element Array{Complex{Float64},1}:\n",
       "  5.15143e-14+3.29958e-13im\n",
       "  2.45137e-13+2.47358e-13im\n",
       " -3.28626e-14+8.28226e-14im\n",
       " -6.26166e-14+3.73035e-14im\n",
       " -2.84495e-16+3.19189e-16im\n",
       " -2.08167e-16+0.0im        \n",
       "  4.28477e-16-9.12465e-16im\n",
       "  1.66533e-16+0.0im        "
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check that we have indeed a correct solution\n",
    "map(gi -> gi(variables(g_hom3)=>s3), g_hom3)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.1",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}