Christine Model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"We define functions to convert Duals, the variable types used by ForwardDiff, \\nto Floats. This is useful if we want to print out the value of a variable \\n(since print doesn't know how to Duals). Note that after being converted to a Float, no\\ndifferentiation by ForwardDiff can happen!  e.g. after\\n    x = convert(Float64, y)\\nForwardDiff can still differentiate y, but it can't differentiate x\\n\""
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using PyCall\n",
    "using PyPlot\n",
    "using ForwardDiff\n",
    "using DiffBase\n",
    "using MAT\n",
    "\n",
    "pygui(true)\n",
    "\n",
    "import Base.convert\n",
    "convert(::Type{Float64}, x::ForwardDiff.Dual) = Float64(x.value)\n",
    "function convert(::Array{Float64}, x::Array{ForwardDiff.Dual}) \n",
    "    y = zeros(size(x)); \n",
    "    for i in 1:prod(size(x)) \n",
    "        y[i] = convert(Float64, x[i]) \n",
    "    end\n",
    "    return y\n",
    "end\n",
    "\n",
    "include(\"general_utils.jl\")\n",
    "include(\"hessian_utils.jl\")\n",
    "\n",
    "\"\"\"\n",
    "We define functions to convert Duals, the variable types used by ForwardDiff, \n",
    "to Floats. This is useful if we want to print out the value of a variable \n",
    "(since print doesn't know how to Duals). Note that after being converted to a Float, no\n",
    "differentiation by ForwardDiff can happen!  e.g. after\n",
    "    x = convert(Float64, y)\n",
    "ForwardDiff can still differentiate y, but it can't differentiate x\n",
    "\"\"\"\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Method definition prob_weight(Any) in module Main at In[21]:2 overwritten at In[60]:2.\n",
      "WARNING: Method definition #prob_weight(Array{Any, 1}, Main.#prob_weight, Any) in module Main overwritten.\n",
      "WARNING: Method definition prob_right(Any, Any, Any, Any) in module Main at In[21]:12 overwritten at In[60]:12.\n",
      "WARNING: Method definition #prob_right(Array{Any, 1}, Main.#prob_right, Any, Any, Any, Any) in module Main overwritten.\n",
      "WARNING: Method definition one_trial_cost(Any, Any, Any, Any, Any) in module Main at In[21]:21 overwritten at In[60]:21.\n",
      "WARNING: Method definition #one_trial_cost(Array{Any, 1}, Main.#one_trial_cost, Any, Any, Any, Any, Any) in module Main overwritten.\n",
      "WARNING: Method definition J(Any, Any, Any, Any, Any) in module Main at In[21]:34 overwritten at In[60]:36.\n",
      "WARNING: Method definition #J(Array{Any, 1}, Main.#J, Any, Any, Any, Any, Any) in module Main overwritten.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "J (generic function with 1 method)"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function prob_weight(nflashes; delta=0.5, gamma=3, other_unused_pars...)\n",
    "    nflashes = nflashes/10\n",
    "    return delta.*nflashes.^gamma./(delta.*nflashes.^gamma + (1-nflashes).^gamma)\n",
    "end\n",
    "# From Christine: delta = .5, gamma = 3, \n",
    "# From Christine: lapse = .008, softmax/lambda = 1.8, alpha = .93\n",
    "\n",
    "\n",
    "function prob_right(clickrateL, clickrateR, nflashesL, nflashesR; \n",
    "    lapse=0.008, lambda=1.8, alpha=0.93, bias=0.5, pars...)\n",
    "\n",
    "    answer = exp(-lambda*(clickrateR.^alpha*prob_weight(nflashesR; pars...) - \n",
    "        clickrateL.^alpha*prob_weight(nflashesL; pars...)) + bias)\n",
    "    \n",
    "    return lapse + (1 - 2*lapse)./(1 + answer)    \n",
    "end\n",
    "\n",
    "\n",
    "function one_trial_cost(rat_went_right, clickrateL, clickrateR, nflashesL, nflashesR; pars...)\n",
    "\n",
    "    PR = prob_right(clickrateL, clickrateR, nflashesL, nflashesR; pars...)\n",
    "    \n",
    "    if abs(1-PR) < 1e-16; PR = 1-1e-16; end\n",
    "    if PR < 1e-16;        PR =16      ; end\n",
    "    if rat_went_right==1\n",
    "        return log(PR)\n",
    "    else\n",
    "        return log(1-PR)\n",
    "    end\n",
    "end\n",
    "\n",
    "\n",
    "function J(rat_went_right, clickrateL, clickrateR, nflashesL, nflashesR; \n",
    "            nderivs=0, difforder=0, pars...)\n",
    "     \n",
    "    answer = ForwardDiffZeros(length(rat_went_right), 1, nderivs=nderivs, difforder=difforder)\n",
    "\n",
    "    for i=1:length(rat_went_right)\n",
    "        answer[i] = \n",
    "            one_trial_cost(rat_went_right[i], clickrateL[i], clickrateR[i], \n",
    "                nflashesL[i], nflashesR[i]; pars...)\n",
    "    end    \n",
    "    return sum(answer), answer\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "F = matread(\"J266.mat\")[\"J266_data\"]\n",
    "nflashesL = F[\"lflashes\"]\n",
    "nflashesR = F[\"rflashes\"]\n",
    "clickrateL = F[\"lclickHz\"]\n",
    "clickrateR = F[\"rclickHz\"]\n",
    "rat_went_right = F[\"went_right\"];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "J(rat_went_right, clickrateL, clickrateR, nflashesL, nflashesR)[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ntrials = 1000\n",
    "rat_went_right[ntrials]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0: eta=0.001 ps=[0.400, 0.500]\n",
      "1: eta=0.0011 cost=-37388.2 jtype=constrained costheta=-1.000 ps=[0.400951, 0.500027]\n",
      "2: eta=0.00121 cost=-37397.6 jtype=constrained costheta=-1.000 ps=[0.401998, 0.500056]\n",
      "3: eta=0.001331 cost=-37408.2 jtype=constrained costheta=-1.000 ps=[0.403151, 0.500088]\n",
      "4: eta=0.0014641 cost=-37420.1 jtype=constrained costheta=-1.000 ps=[0.40442, 0.500122]\n",
      "5: eta=0.00161051 cost=-37433.7 jtype=constrained costheta=-1.000 ps=[0.405818, 0.50016]\n",
      "6: eta=0.00177156 cost=-37449 jtype=constrained costheta=-1.000 ps=[0.407358, 0.500201]\n",
      "7: eta=0.00194872 cost=-37466.5 jtype=constrained costheta=-1.000 ps=[0.409055, 0.500247]\n",
      "8: eta=0.00214359 cost=-37486.4 jtype=constrained costheta=-1.000 ps=[0.410923, 0.500296]\n",
      "9: eta=0.00235795 cost=-37509.2 jtype=constrained costheta=-1.000 ps=[0.412982, 0.500349]\n",
      "10: eta=0.00259374 cost=-37535.3 jtype=constrained costheta=-1.000 ps=[0.415251, 0.500408]\n",
      "11: eta=0.00285312 cost=-37565.2 jtype=constrained costheta=-1.000 ps=[0.417752, 0.500471]\n",
      "12: eta=0.00313843 cost=-37599.8 jtype=constrained costheta=-1.000 ps=[0.420507, 0.50054]\n",
      "13: eta=0.00345227 cost=-37639.8 jtype=constrained costheta=-1.000 ps=[0.423545, 0.500615]\n",
      "14: eta=0.0037975 cost=-37686.1 jtype=constrained costheta=-1.000 ps=[0.426894, 0.500697]\n",
      "15: eta=0.00417725 cost=-37740.1 jtype=constrained costheta=-1.000 ps=[0.430585, 0.500786]\n",
      "16: eta=0.00459497 cost=-37803.1 jtype=constrained costheta=-1.000 ps=[0.434656, 0.500883]\n",
      "17: eta=0.00505447 cost=-37877 jtype=constrained costheta=-1.000 ps=[0.439145, 0.500988]\n",
      "18: eta=0.00555992 cost=-37963.8 jtype=constrained costheta=-1.000 ps=[0.444095, 0.501102]\n",
      "19: eta=0.00611591 cost=-38066.2 jtype=constrained costheta=-1.000 ps=[0.449554, 0.501226]\n",
      "20: eta=0.0067275 cost=-38187.5 jtype=constrained costheta=-1.000 ps=[0.455574, 0.501361]\n",
      "21: eta=0.00740025 cost=-38331.4 jtype=constrained costheta=-1.000 ps=[0.462212, 0.501508]\n",
      "22: eta=0.00814027 cost=-38502.9 jtype=constrained costheta=-1.000 ps=[0.46953, 0.501667]\n",
      "23: eta=0.0089543 cost=-38707.7 jtype=constrained costheta=-1.000 ps=[0.477596, 0.501841]\n",
      "24: eta=0.00984973 cost=-38952.9 jtype=constrained costheta=-1.000 ps=[0.486483, 0.502031]\n",
      "25: eta=0.0108347 cost=-39246.9 jtype=constrained costheta=-1.000 ps=[0.49627, 0.502238]\n",
      "26: eta=0.0119182 cost=-39599.9 jtype=constrained costheta=-1.000 ps=[0.507039, 0.502464]\n",
      "27: eta=0.01311 cost=-40023.7 jtype=constrained costheta=-1.000 ps=[0.518878, 0.502711]\n",
      "28: eta=0.014421 cost=-40531.9 jtype=constrained costheta=-1.000 ps=[0.531877, 0.502982]\n",
      "29: eta=0.0158631 cost=-41139.5 jtype=constrained costheta=-1.000 ps=[0.546125, 0.50328]\n",
      "30: eta=0.0174494 cost=-41862.2 jtype=constrained costheta=-1.000 ps=[0.56171, 0.50361]\n",
      "31: eta=0.0191943 cost=-42715.1 jtype=constrained costheta=-1.000 ps=[0.578712, 0.503975]\n",
      "32: eta=0.0211138 cost=-43710.8 jtype=constrained costheta=-1.000 ps=[0.597196, 0.504382]\n",
      "33: eta=0.0232252 cost=-44856.7 jtype=constrained costheta=-1.000 ps=[0.617207, 0.504839]\n",
      "34: eta=0.0255477 cost=-46152.1 jtype=constrained costheta=-1.000 ps=[0.63876, 0.505356]\n",
      "35: eta=0.0281024 cost=-47587.2 jtype=constrained costheta=-1.000 ps=[0.661825, 0.505948]\n",
      "36: eta=0.0309127 cost=-49142 jtype=constrained costheta=-1.000 ps=[0.686314, 0.506633]\n",
      "37: eta=0.0340039 cost=-50788.3 jtype=constrained costheta=-1.000 ps=[0.712072, 0.507436]\n",
      "38: eta=0.0374043 cost=-52491.2 jtype=constrained costheta=-1.000 ps=[0.738857, 0.508393]\n",
      "39: eta=0.0411448 cost=-54209.4 jtype=constrained costheta=-1.000 ps=[0.766335, 0.509552]\n",
      "40: eta=0.0452593 cost=-55896.2 jtype=constrained costheta=-1.000 ps=[0.794072, 0.510982]\n",
      "41: eta=0.0497852 cost=-57502.4 jtype=constrained costheta=-1.000 ps=[0.821543, 0.512788]\n",
      "42: eta=0.0547637 cost=-58986.2 jtype=constrained costheta=-0.999 ps=[0.848149, 0.51512]\n",
      "43: eta=0.0602401 cost=-60324 jtype=constrained costheta=-0.999 ps=[0.873257, 0.518192]\n",
      "44: eta=0.0662641 cost=-61511.9 jtype=constrained costheta=-0.998 ps=[0.896256, 0.522285]\n",
      "45: eta=0.0728905 cost=-62559.6 jtype=constrained costheta=-0.996 ps=[0.916623, 0.527773]\n",
      "46: eta=0.0801795 cost=-63484.3 jtype=constrained costheta=-0.994 ps=[0.933992, 0.535121]\n",
      "47: eta=0.0881975 cost=-64308.6 jtype=constrained costheta=-0.992 ps=[0.948229, 0.544846]\n",
      "48: eta=0.0970172 cost=-65058.7 jtype=constrained costheta=-0.990 ps=[0.959479, 0.557442]\n",
      "49: eta=0.106719 cost=-65761.7 jtype=constrained costheta=-0.990 ps=[0.968124, 0.57332]\n",
      "50: eta=0.117391 cost=-66441.9 jtype=constrained costheta=-0.992 ps=[0.974667, 0.59283]\n",
      "51: eta=0.12913 cost=-67118.9 jtype=constrained costheta=-0.994 ps=[0.979607, 0.616329]\n",
      "52: eta=0.142043 cost=-67806.9 jtype=constrained costheta=-0.996 ps=[0.983361, 0.644262]\n",
      "53: eta=0.156247 cost=-68515.5 jtype=constrained costheta=-0.997 ps=[0.986249, 0.677214]\n",
      "54: eta=0.171872 cost=-69249.8 jtype=constrained costheta=-0.998 ps=[0.988505, 0.715959]\n",
      "55: eta=0.189059 cost=-70010 jtype=constrained costheta=-0.999 ps=[0.990295, 0.761495]\n",
      "56: eta=0.207965 cost=-70791.2 jtype=constrained costheta=-0.999 ps=[0.991739, 0.815101]\n",
      "57: eta=0.228762 cost=-71581.3 jtype=constrained costheta=-0.999 ps=[0.992922, 0.878403]\n",
      "58: eta=0.251638 cost=-72361 jtype=constrained costheta=-1.000 ps=[0.993906, 0.953454]\n",
      "59: eta=0.276801 cost=-73105.2 jtype=constrained costheta=-1.000 ps=[0.994738, 1.04286]\n",
      "60: eta=0.304482 cost=-73788.8 jtype=constrained costheta=-1.000 ps=[0.995454, 1.1499]\n",
      "61: eta=0.33493 cost=-74395.1 jtype=constrained costheta=-1.000 ps=[0.99608, 1.27875]\n",
      "62: eta=0.368423 cost=-74921.3 jtype=constrained costheta=-1.000 ps=[0.996634, 1.43466]\n",
      "63: eta=0.405265 cost=-75376.3 jtype=constrained costheta=-1.000 ps=[0.997122, 1.62428]\n",
      "64: eta=0.445792 cost=-75775.8 jtype=constrained costheta=-1.000 ps=[0.997547, 1.85584]\n",
      "65: eta=0.490371 cost=-76135.6 jtype=constrained costheta=-1.000 ps=[0.997913, 2.13931]\n",
      "66: eta=0.539408 cost=-76463.6 jtype=constrained costheta=-1.000 ps=[0.998225, 2.48633]\n",
      "67: eta=0.593349 cost=-76755.7 jtype=constrained costheta=-1.000 ps=[0.998486, 2.90955]\n",
      "68: eta=0.652683 cost=-77000.7 jtype=constrained costheta=-1.000 ps=[0.998705, 3.42105]\n",
      "69: eta=0.717952 cost=-77194.4 jtype=constrained costheta=-1.000 ps=[0.998892, 4.02918]\n",
      "70: eta=0.789747 cost=-77347.1 jtype=constrained costheta=-1.000 ps=[0.999071, 4.73353]\n",
      "71: eta=0.868722 cost=-77473.5 jtype=constrained costheta=-1.000 ps=[0.99926, 5.51911]\n",
      "72: eta=0.955594 cost=-77581 jtype=constrained costheta=-1.000 ps=[0.99945, 6.35204]\n",
      "73: eta=1.05115 cost=-77669.9 jtype=constrained costheta=-1.000 ps=[0.999623, 7.18086]\n",
      "74: eta=1.15627 cost=-77739.9 jtype=constrained costheta=-1.000 ps=[0.999759, 7.94628]\n",
      "75: eta=1.2719 cost=-77792.6 jtype=constrained costheta=-1.000 ps=[0.999854, 8.59702]\n",
      "76: eta=1.39908 cost=-77830.5 jtype=constrained costheta=-1.000 ps=[0.999916, 9.1041]\n",
      "77: eta=1.53899 cost=-77856.2 jtype=constrained costheta=-1.000 ps=[0.999953, 9.46608]\n",
      "78: eta=1.69289 cost=-77872.6 jtype=constrained costheta=-1.000 ps=[0.999975, 9.70332]\n",
      "79: eta=1.86218 cost=-77882.3 jtype=constrained costheta=-1.000 ps=[0.999987, 9.84651]\n",
      "80: eta=2.0484 cost=-77887.6 jtype=constrained costheta=-1.000 ps=[0.999994, 9.92626]\n",
      "81: eta=2.25324 cost=-77890.4 jtype=constrained costheta=-1.000 ps=[0.999997, 9.96723]\n",
      "82: eta=2.47856 cost=-77891.6 jtype=constrained costheta=-1.000 ps=[0.999999, 9.98661]\n",
      "83: eta=2.72642 cost=-77892.2 jtype=constrained costheta=-1.000 ps=[1, 9.995]\n",
      "84: eta=2.99906 cost=-77892.4 jtype=Newton costheta=-1.000 ps=[1, 9.99816]\n",
      "85: eta=3.29897 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 9.99932]\n",
      "86: eta=3.62887 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 9.99975]\n",
      "87: eta=3.99175 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 9.99991]\n",
      "88: eta=4.39093 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 9.99997]\n",
      "89: eta=4.83002 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 9.99999]\n",
      "90: eta=5.31302 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "91: eta=5.84432 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "92: eta=6.42876 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "93: eta=7.07163 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "94: eta=7.7788 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "95: eta=8.55668 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "96: eta=9.41234 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "97: eta=10.3536 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "98: eta=11.3889 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "99: eta=12.5278 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "100: eta=13.7806 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "101: eta=15.1587 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "102: eta=16.6745 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "103: eta=18.342 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "104: eta=20.1762 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "105: eta=22.1938 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "106: eta=24.4132 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "107: eta=26.8545 cost=-77892.5 jtype=Newton costheta=-1.000 ps=[1, 10]\n",
      "About to break -- tol=1e-12, new_cost-cost=0, eta=26.8545\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "([1.0,10.0],\n",
       "[0.001 0.0011 … 24.4132 26.8545; -37379.8 -37388.2 … -77892.5 -77892.5; 0.4 0.400951 … 1.0 1.0; 0.5 0.500027 … 10.0 10.0],\n",
       "\n",
       "-77892.53968849515,\n",
       "[148.0 148.0 … 5.0 55.0; 0.0 0.0 … -2.84343e-25 1.13899e-25])"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# From Christine: alpha = .93, delta = .5, gamma = 3, softmax/lambda = 1.8, lapse = .008\n",
    "args = [\"alpha\", \"gamma\", \"lambda\", \"lapse\"]\n",
    "seed = [0.93, 0.5]\n",
    "bbox = Dict(:alpha=>[0.0001 1], :lambda=>[0.0001 10])\n",
    "\n",
    "# ntrials = 1000\n",
    "ntrials = length(rat_went_right)\n",
    "\n",
    "params, traj = bbox_Hessian_keyword_minimization(seed, args, bbox, (;params...) -> J(rat_went_right[1:ntrials], \n",
    "clickrateL[1:ntrials], clickrateR[1:ntrials], nflashesL[1:ntrials], nflashesR[1:ntrials]; params...)[1];\n",
    "verbose=true, start_eta=0.001, tol=1e-12)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1: eta=0.001 cost=-8.0895 costheta=-1.000 ps=[0.410, 0.507]\n",
      "2: eta=0.001 cost=-8.2630 costheta=-1.000 ps=[0.421, 0.515]\n",
      "3: eta=0.001 cost=-8.4924 costheta=-1.000 ps=[0.435, 0.525]\n",
      "4: eta=0.001 cost=-8.8029 costheta=-1.000 ps=[0.453, 0.536]\n",
      "5: eta=0.002 cost=-9.2312 costheta=-1.000 ps=[0.474, 0.550]\n",
      "6: eta=0.002 cost=-9.8227 costheta=-1.000 ps=[0.501, 0.567]\n",
      "7: eta=0.002 cost=-10.5967 costheta=-1.000 ps=[0.534, 0.587]\n",
      "8: eta=0.002 cost=-11.4519 costheta=-1.000 ps=[0.569, 0.609]\n",
      "9: eta=0.002 cost=-12.1965 costheta=-1.000 ps=[0.601, 0.631]\n",
      "10: eta=0.003 cost=-12.8181 costheta=-0.999 ps=[0.631, 0.653]\n",
      "11: eta=0.003 cost=-13.4380 costheta=-0.999 ps=[0.661, 0.678]\n",
      "12: eta=0.003 cost=-14.2051 costheta=-1.000 ps=[0.697, 0.707]\n",
      "13: eta=0.003 cost=-15.3089 costheta=-1.000 ps=[0.745, 0.744]\n",
      "14: eta=0.004 cost=-16.9658 costheta=-1.000 ps=[0.810, 0.790]\n",
      "15: eta=0.004 cost=-19.2418 costheta=-1.000 ps=[0.890, 0.843]\n",
      "16: eta=0.005 cost=-21.4741 costheta=-1.000 ps=[0.973, 0.895]\n",
      "17: eta=0.005 cost=-22.7657 costheta=-1.000 ps=[1.036, 0.932]\n",
      "18: eta=0.006 cost=-23.3753 costheta=-1.000 ps=[1.079, 0.957]\n",
      "19: eta=0.006 cost=-23.6824 costheta=-1.000 ps=[1.111, 0.975]\n",
      "20: eta=0.007 cost=-23.8514 costheta=-1.000 ps=[1.137, 0.989]\n",
      "21: eta=0.007 cost=-23.9515 costheta=-1.000 ps=[1.157, 1.001]\n",
      "22: eta=0.008 cost=-24.0146 costheta=-1.000 ps=[1.175, 1.010]\n",
      "23: eta=0.009 cost=-24.0565 costheta=-1.000 ps=[1.189, 1.018]\n",
      "24: eta=0.010 cost=-24.0854 costheta=-1.000 ps=[1.203, 1.026]\n",
      "25: eta=0.011 cost=-24.1060 costheta=-1.000 ps=[1.214, 1.032]\n",
      "26: eta=0.012 cost=-24.1213 costheta=-1.000 ps=[1.225, 1.038]\n",
      "27: eta=0.013 cost=-24.1327 costheta=-1.000 ps=[1.235, 1.043]\n",
      "28: eta=0.014 cost=-24.1415 costheta=-1.000 ps=[1.244, 1.048]\n",
      "29: eta=0.016 cost=-24.1484 costheta=-1.000 ps=[1.252, 1.052]\n",
      "30: eta=0.017 cost=-24.1539 costheta=-1.000 ps=[1.260, 1.056]\n",
      "31: eta=0.019 cost=-24.1583 costheta=-1.000 ps=[1.267, 1.060]\n",
      "32: eta=0.021 cost=-24.1619 costheta=-1.000 ps=[1.274, 1.064]\n",
      "33: eta=0.023 cost=-24.1648 costheta=-1.000 ps=[1.281, 1.067]\n",
      "34: eta=0.026 cost=-24.1673 costheta=-1.000 ps=[1.288, 1.071]\n",
      "35: eta=0.028 cost=-24.1693 costheta=-1.000 ps=[1.294, 1.074]\n",
      "36: eta=0.031 cost=-24.1710 costheta=-1.000 ps=[1.300, 1.077]\n",
      "37: eta=0.034 cost=-24.1724 costheta=-1.000 ps=[1.305, 1.080]\n",
      "38: eta=0.037 cost=-24.1737 costheta=-1.000 ps=[1.311, 1.083]\n",
      "39: eta=0.041 cost=-24.1747 costheta=-1.000 ps=[1.316, 1.085]\n",
      "40: eta=0.045 cost=-24.1756 costheta=-1.000 ps=[1.321, 1.088]\n",
      "41: eta=0.050 cost=-24.1764 costheta=-1.000 ps=[1.327, 1.091]\n",
      "42: eta=0.055 cost=-24.1770 costheta=-1.000 ps=[1.331, 1.093]\n",
      "43: eta=0.060 cost=-24.1776 costheta=-1.000 ps=[1.336, 1.096]\n",
      "44: eta=0.066 cost=-24.1781 costheta=-1.000 ps=[1.341, 1.098]\n",
      "45: eta=0.073 cost=-24.1785 costheta=-1.000 ps=[1.346, 1.100]\n",
      "46: eta=0.080 cost=-24.1789 costheta=-1.000 ps=[1.350, 1.103]\n",
      "47: eta=0.088 cost=-24.1792 costheta=-1.000 ps=[1.355, 1.105]\n",
      "48: eta=0.097 cost=-24.1795 costheta=-1.000 ps=[1.359, 1.107]\n",
      "49: eta=0.107 cost=-24.1798 costheta=-1.000 ps=[1.363, 1.109]\n",
      "50: eta=0.117 cost=-24.1800 costheta=-1.000 ps=[1.368, 1.111]\n",
      "51: eta=0.129 cost=-24.1802 costheta=-1.000 ps=[1.372, 1.114]\n",
      "52: eta=0.142 cost=-24.1804 costheta=-1.000 ps=[1.376, 1.116]\n",
      "53: eta=0.156 cost=-24.1805 costheta=-1.000 ps=[1.380, 1.118]\n",
      "54: eta=0.172 cost=-24.1807 costheta=-1.000 ps=[1.384, 1.120]\n",
      "55: eta=0.189 cost=-24.1808 costheta=-1.000 ps=[1.388, 1.122]\n",
      "56: eta=0.208 cost=-24.1809 costheta=-1.000 ps=[1.392, 1.123]\n",
      "57: eta=0.229 cost=-24.1810 costheta=-1.000 ps=[1.395, 1.125]\n",
      "58: eta=0.252 cost=-24.1811 costheta=-1.000 ps=[1.399, 1.127]\n",
      "59: eta=0.277 cost=-24.1811 costheta=-1.000 ps=[1.403, 1.129]\n",
      "60: eta=0.304 cost=-24.1812 costheta=-1.000 ps=[1.407, 1.131]\n",
      "61: eta=0.335 cost=-24.1813 costheta=-1.000 ps=[1.410, 1.133]\n",
      "62: eta=0.368 cost=-24.1813 costheta=-1.000 ps=[1.414, 1.134]\n",
      "63: eta=0.405 cost=-24.1814 costheta=-1.000 ps=[1.417, 1.136]\n",
      "64: eta=0.446 cost=-24.1814 costheta=-1.000 ps=[1.421, 1.138]\n",
      "65: eta=0.490 cost=-24.1814 costheta=-1.000 ps=[1.424, 1.140]\n",
      "66: eta=0.539 cost=-24.1815 costheta=-1.000 ps=[1.428, 1.141]\n",
      "67: eta=0.593 cost=-24.1815 costheta=-1.000 ps=[1.431, 1.143]\n",
      "68: eta=0.653 cost=-24.1815 costheta=-1.000 ps=[1.435, 1.145]\n",
      "69: eta=0.718 cost=-24.1815 costheta=-1.000 ps=[1.438, 1.146]\n",
      "70: eta=0.790 cost=-24.1816 costheta=-1.000 ps=[1.441, 1.148]\n",
      "71: eta=0.869 cost=-24.1816 costheta=-1.000 ps=[1.445, 1.150]\n",
      "72: eta=0.956 cost=-24.1816 costheta=-1.000 ps=[1.448, 1.151]\n",
      "73: eta=1.051 cost=-24.1816 costheta=-1.000 ps=[1.451, 1.153]\n",
      "74: eta=1.156 cost=-24.1816 costheta=-1.000 ps=[1.455, 1.154]\n",
      "75: eta=1.272 cost=-24.1816 costheta=-1.000 ps=[1.458, 1.156]\n",
      "76: eta=1.399 cost=-24.1816 costheta=-1.000 ps=[1.461, 1.157]\n",
      "77: eta=1.539 cost=-24.1817 costheta=-1.000 ps=[1.464, 1.159]\n",
      "78: eta=1.693 cost=-24.1817 costheta=-1.000 ps=[1.467, 1.160]\n",
      "79: eta=1.862 cost=-24.1817 costheta=-1.000 ps=[1.470, 1.162]\n",
      "80: eta=2.048 cost=-24.1817 costheta=-1.000 ps=[1.473, 1.163]\n",
      "81: eta=2.253 cost=-24.1817 costheta=-1.000 ps=[1.476, 1.165]\n",
      "82: eta=2.479 cost=-24.1817 costheta=-1.000 ps=[1.479, 1.166]\n",
      "83: eta=2.726 cost=-24.1817 costheta=-1.000 ps=[1.482, 1.168]\n",
      "84: eta=2.999 cost=-24.1817 costheta=-1.000 ps=[1.485, 1.169]\n",
      "85: eta=3.299 cost=-24.1817 costheta=-1.000 ps=[1.488, 1.170]\n",
      "86: eta=3.629 cost=-24.1817 costheta=-1.000 ps=[1.491, 1.172]\n",
      "87: eta=3.992 cost=-24.1817 costheta=-1.000 ps=[1.494, 1.173]\n",
      "88: eta=4.391 cost=-24.1817 costheta=-1.000 ps=[1.497, 1.175]\n",
      "89: eta=4.830 cost=-24.1817 costheta=-1.000 ps=[1.500, 1.176]\n",
      "90: eta=5.313 cost=-24.1817 costheta=-1.000 ps=[1.503, 1.177]\n",
      "91: eta=5.844 cost=-24.1817 costheta=-1.000 ps=[1.506, 1.179]\n",
      "92: eta=6.429 cost=-24.1817 costheta=-1.000 ps=[1.508, 1.180]\n",
      "93: eta=7.072 cost=-24.1817 costheta=-1.000 ps=[1.511, 1.181]\n",
      "94: eta=7.779 cost=-24.1817 costheta=-1.000 ps=[1.514, 1.183]\n",
      "95: eta=8.557 cost=-24.1817 costheta=-1.000 ps=[1.517, 1.184]\n",
      "96: eta=9.412 cost=-24.1817 costheta=-1.000 ps=[1.519, 1.185]\n",
      "97: eta=10.354 cost=-24.1817 costheta=-1.000 ps=[1.522, 1.187]\n",
      "98: eta=11.389 cost=-24.1817 costheta=-1.000 ps=[1.525, 1.188]\n",
      "99: eta=12.528 cost=-24.1817 costheta=-1.000 ps=[1.528, 1.189]\n",
      "100: eta=13.781 cost=-24.1817 costheta=-1.000 ps=[1.530, 1.190]\n",
      "101: eta=15.159 cost=-24.1817 costheta=-1.000 ps=[1.533, 1.192]\n",
      "102: eta=16.675 cost=-24.1817 costheta=-1.000 ps=[1.536, 1.193]\n",
      "103: eta=18.342 cost=-24.1817 costheta=-1.000 ps=[1.538, 1.194]\n",
      "104: eta=20.176 cost=-24.1817 costheta=-1.000 ps=[1.541, 1.195]\n",
      "105: eta=22.194 cost=-24.1817 costheta=-1.000 ps=[1.543, 1.196]\n",
      "106: eta=24.413 cost=-24.1817 costheta=-1.000 ps=[1.546, 1.198]\n",
      "107: eta=26.855 cost=-24.1817 costheta=-1.000 ps=[1.548, 1.199]\n",
      "108: eta=29.540 cost=-24.1817 costheta=-1.000 ps=[1.551, 1.200]\n",
      "109: eta=32.494 cost=-24.1817 costheta=-1.000 ps=[1.554, 1.201]\n",
      "110: eta=35.743 cost=-24.1817 costheta=-1.000 ps=[1.556, 1.202]\n",
      "111: eta=39.318 cost=-24.1817 costheta=-1.000 ps=[1.559, 1.203]\n",
      "112: eta=43.249 cost=-24.1817 costheta=-1.000 ps=[1.561, 1.205]\n",
      "113: eta=47.574 cost=-24.1817 costheta=-1.000 ps=[1.564, 1.206]\n",
      "114: eta=52.332 cost=-24.1817 costheta=-1.000 ps=[1.566, 1.207]\n",
      "115: eta=57.565 cost=-24.1817 costheta=-1.000 ps=[1.568, 1.208]\n",
      "116: eta=63.322 cost=-24.1817 costheta=-1.000 ps=[1.571, 1.209]\n",
      "117: eta=69.654 cost=-24.1817 costheta=-1.000 ps=[1.573, 1.210]\n",
      "118: eta=76.619 cost=-24.1817 costheta=-1.000 ps=[1.576, 1.211]\n",
      "119: eta=84.281 cost=-24.1817 costheta=-1.000 ps=[1.578, 1.212]\n",
      "120: eta=92.709 cost=-24.1817 costheta=-1.000 ps=[1.580, 1.214]\n",
      "121: eta=101.980 cost=-24.1817 costheta=-1.000 ps=[1.583, 1.215]\n",
      "122: eta=112.178 cost=-24.1817 costheta=-1.000 ps=[1.585, 1.216]\n",
      "123: eta=123.396 cost=-24.1817 costheta=-1.000 ps=[1.587, 1.217]\n",
      "124: eta=135.735 cost=-24.1817 costheta=-1.000 ps=[1.590, 1.218]\n",
      "125: eta=149.309 cost=-24.1817 costheta=-1.000 ps=[1.592, 1.219]\n",
      "126: eta=164.240 cost=-24.1817 costheta=-1.000 ps=[1.594, 1.220]\n",
      "127: eta=180.664 cost=-24.1817 costheta=-1.000 ps=[1.597, 1.221]\n",
      "128: eta=198.730 cost=-24.1817 costheta=-1.000 ps=[1.599, 1.222]\n",
      "129: eta=218.603 cost=-24.1817 costheta=-1.000 ps=[1.601, 1.223]\n",
      "130: eta=240.463 cost=-24.1817 costheta=-1.000 ps=[1.604, 1.224]\n",
      "131: eta=264.510 cost=-24.1817 costheta=-1.000 ps=[1.606, 1.225]\n",
      "132: eta=290.961 cost=-24.1817 costheta=-1.000 ps=[1.608, 1.226]\n",
      "133: eta=320.057 cost=-24.1817 costheta=-1.000 ps=[1.610, 1.227]\n",
      "134: eta=352.063 cost=-24.1817 costheta=-1.000 ps=[1.612, 1.228]\n",
      "135: eta=387.269 cost=-24.1817 costheta=-1.000 ps=[1.615, 1.229]\n",
      "136: eta=425.996 cost=-24.1817 costheta=-1.000 ps=[1.617, 1.230]\n",
      "137: eta=468.595 cost=-24.1817 costheta=-1.000 ps=[1.619, 1.231]\n",
      "138: eta=515.455 cost=-24.1817 costheta=-1.000 ps=[1.621, 1.232]\n",
      "139: eta=567.000 cost=-24.1817 costheta=-1.000 ps=[1.623, 1.233]\n",
      "140: eta=623.700 cost=-24.1817 costheta=-1.000 ps=[1.625, 1.234]\n",
      "141: eta=686.070 cost=-24.1817 costheta=-1.000 ps=[1.628, 1.235]\n",
      "142: eta=754.677 cost=-24.1817 costheta=-1.000 ps=[1.630, 1.236]\n",
      "143: eta=830.145 cost=-24.1817 costheta=-1.000 ps=[1.632, 1.237]\n",
      "144: eta=913.160 cost=-24.1817 costheta=-1.000 ps=[1.634, 1.238]\n",
      "145: eta=1004.475 cost=-24.1817 costheta=-1.000 ps=[1.636, 1.239]\n",
      "146: eta=1104.923 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "147: eta=1215.415 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "148: eta=1336.957 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "149: eta=1470.653 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "150: eta=1617.718 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "151: eta=1779.490 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "152: eta=1957.439 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "153: eta=2153.182 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "154: eta=2368.501 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "155: eta=2605.351 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "156: eta=2865.886 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "157: eta=3152.474 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "158: eta=3467.722 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "159: eta=3814.494 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "160: eta=4195.943 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "161: eta=4615.538 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "162: eta=5077.092 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "163: eta=5584.801 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "164: eta=6143.281 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "165: eta=6757.609 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "166: eta=7433.370 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "167: eta=8176.707 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "168: eta=8994.377 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "169: eta=9893.815 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "170: eta=10883.197 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "171: eta=11971.516 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "172: eta=13168.668 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "173: eta=14485.535 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "174: eta=15934.088 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "175: eta=17527.497 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "176: eta=19280.247 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "177: eta=21208.271 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "178: eta=23329.099 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "179: eta=25662.008 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "180: eta=28228.209 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "181: eta=31051.030 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "182: eta=34156.133 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "183: eta=37571.747 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "184: eta=41328.921 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "185: eta=45461.813 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "186: eta=50007.995 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "187: eta=55008.794 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "188: eta=60509.674 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "189: eta=66560.641 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "190: eta=73216.705 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "191: eta=80538.375 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "192: eta=88592.213 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "193: eta=97451.434 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "194: eta=107196.578 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "195: eta=117916.236 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "196: eta=129707.859 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "197: eta=142678.645 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "198: eta=156946.509 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "199: eta=172641.160 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "200: eta=189905.276 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "201: eta=208895.804 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "202: eta=229785.385 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "203: eta=252763.923 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "204: eta=278040.315 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "205: eta=305844.347 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "206: eta=336428.781 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "207: eta=370071.660 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "208: eta=407078.826 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "209: eta=447786.708 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "210: eta=492565.379 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "211: eta=541821.917 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "212: eta=596004.109 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "213: eta=655604.519 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "214: eta=721164.971 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "215: eta=793281.468 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "216: eta=872609.615 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "217: eta=959870.577 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "218: eta=1055857.635 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "219: eta=1161443.398 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "220: eta=1277587.738 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "221: eta=1405346.512 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "222: eta=1545881.163 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "223: eta=1700469.279 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "224: eta=1870516.207 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "225: eta=2057567.828 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "226: eta=2263324.610 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "227: eta=2489657.071 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "228: eta=2738622.778 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "229: eta=3012485.056 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "230: eta=3313733.562 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "231: eta=3645106.918 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "232: eta=4009617.610 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "233: eta=4410579.371 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "234: eta=4851637.308 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "235: eta=5336801.039 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "236: eta=5870481.143 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "237: eta=6457529.257 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "238: eta=7103282.183 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "239: eta=7813610.401 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "240: eta=8594971.441 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "241: eta=9454468.585 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "242: eta=10399915.444 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "243: eta=11439906.988 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "244: eta=12583897.687 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "245: eta=13842287.456 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "246: eta=15226516.201 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "247: eta=16749167.821 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "248: eta=18424084.603 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "249: eta=20266493.064 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "250: eta=22293142.370 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "251: eta=24522456.607 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "252: eta=26974702.268 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "253: eta=29672172.495 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "254: eta=32639389.744 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "255: eta=35903328.718 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "256: eta=39493661.590 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "257: eta=43443027.749 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "258: eta=47787330.524 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "259: eta=52566063.577 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "260: eta=57822669.934 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "261: eta=63604936.928 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "262: eta=69965430.620 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "263: eta=76961973.682 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "264: eta=84658171.051 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "265: eta=93123988.156 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "266: eta=102436386.971 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "267: eta=112680025.669 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "268: eta=123948028.235 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "269: eta=136342831.059 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "270: eta=149977114.165 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "271: eta=164974825.581 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "272: eta=181472308.139 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "273: eta=199619538.953 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "274: eta=219581492.849 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "275: eta=241539642.134 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "276: eta=265693606.347 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "277: eta=292262966.982 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "278: eta=321489263.680 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "279: eta=353638190.048 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "280: eta=389002009.053 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "281: eta=427902209.958 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "282: eta=470692430.954 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "283: eta=517761674.049 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "284: eta=569537841.454 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "285: eta=626491625.599 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "286: eta=689140788.159 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "287: eta=758054866.975 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "288: eta=833860353.673 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "289: eta=917246389.040 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "290: eta=1008971027.944 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "291: eta=1109868130.738 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "292: eta=1220854943.812 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "293: eta=1342940438.193 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "294: eta=1477234482.012 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "295: eta=1624957930.214 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "296: eta=1787453723.235 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "297: eta=1966199095.559 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "298: eta=2162819005.114 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "299: eta=2379100905.626 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "300: eta=2617010996.188 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "301: eta=2878712095.807 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "302: eta=3166583305.388 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "303: eta=3483241635.927 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "304: eta=3831565799.520 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "305: eta=4214722379.471 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "306: eta=4636194617.419 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "307: eta=5099814079.160 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "308: eta=5609795487.077 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "309: eta=6170775035.784 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "310: eta=6787852539.363 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "311: eta=7466637793.299 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "312: eta=8213301572.629 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "313: eta=9034631729.892 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "314: eta=9938094902.881 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "315: eta=10931904393.169 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "316: eta=12025094832.486 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "317: eta=13227604315.734 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "318: eta=14550364747.308 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "319: eta=16005401222.039 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "320: eta=17605941344.242 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "321: eta=19366535478.667 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "322: eta=21303189026.533 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "323: eta=23433507929.187 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "324: eta=25776858722.105 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "325: eta=28354544594.316 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "326: eta=31189999053.748 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "327: eta=34308998959.122 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "328: eta=37739898855.035 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "329: eta=41513888740.538 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "330: eta=45665277614.592 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "331: eta=50231805376.051 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "332: eta=55254985913.656 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "333: eta=60780484505.022 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "334: eta=66858532955.524 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "335: eta=73544386251.076 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "336: eta=80898824876.184 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "337: eta=88988707363.802 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "338: eta=97887578100.182 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "339: eta=107676335910.201 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "340: eta=118443969501.221 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "341: eta=130288366451.343 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "342: eta=143317203096.477 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "343: eta=157648923406.125 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "344: eta=173413815746.737 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "345: eta=190755197321.411 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "346: eta=209830717053.552 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "347: eta=230813788758.908 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "348: eta=253895167634.798 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "349: eta=279284684398.278 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "350: eta=307213152838.106 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "351: eta=337934468121.917 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "352: eta=371727914934.109 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "353: eta=408900706427.519 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "354: eta=449790777070.271 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "355: eta=494769854777.299 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "356: eta=544246840255.028 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "357: eta=598671524280.531 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "358: eta=658538676708.585 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "359: eta=724392544379.443 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "360: eta=796831798817.387 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "361: eta=876514978699.126 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "362: eta=964166476569.039 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "363: eta=1060583124225.943 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "364: eta=1166641436648.537 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "365: eta=1283305580313.391 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "366: eta=1411636138344.730 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "367: eta=1552799752179.203 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "368: eta=1708079727397.124 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "369: eta=1878887700136.836 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "370: eta=2066776470150.520 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "371: eta=2273454117165.572 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "372: eta=2500799528882.130 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "373: eta=2750879481770.343 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "374: eta=3025967429947.378 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "375: eta=3328564172942.116 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "376: eta=3661420590236.328 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "377: eta=4027562649259.961 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "378: eta=4430318914185.958 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "379: eta=4873350805604.555 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "380: eta=5360685886165.011 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "381: eta=5896754474781.513 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "382: eta=6486429922259.664 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "383: eta=7135072914485.631 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "384: eta=7848580205934.194 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "385: eta=8633438226527.614 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "386: eta=9496782049180.377 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "387: eta=10446460254098.416 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "388: eta=11491106279508.258 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "389: eta=12640216907459.084 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "390: eta=13904238598204.994 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "391: eta=15294662458025.494 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "392: eta=16824128703828.045 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "393: eta=18506541574210.852 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "394: eta=20357195731631.938 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "395: eta=22392915304795.133 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "396: eta=24632206835274.648 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "397: eta=27095427518802.117 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "398: eta=29804970270682.332 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "399: eta=32785467297750.566 cost=NaN costheta=NaN ps=[NaN, NaN]\n",
      "400: eta=36064014027525.625 cost=NaN costheta=NaN ps=[NaN, NaN]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "([NaN,NaN],NaN,\n",
       "[0.0011 0.00121 … 3.27855e13 3.6064e13; -8.08949 -8.263 … NaN NaN; 0.409587 0.421136 … NaN NaN; 0.50686 0.514938 … NaN NaN])"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# From Christine: alpha = .93, delta = .5, gamma = 3, softmax/lambda = 1.8, lapse = .008\n",
    "args = [\"alpha\", \"lambda\"]\n",
    "seed = [0.4, 0.5]\n",
    "\n",
    "ntrials = 10\n",
    "# ntrials = length(rat_went_right)\n",
    "\n",
    "params, cost, traj = trust_region_keyword_Hessian_minimization(seed, args, (;params...) -> J(rat_went_right[1:ntrials], \n",
    "clickrateL[1:ntrials], clickrateR[1:ntrials], nflashesL[1:ntrials], nflashesR[1:ntrials]; params...)[1];\n",
    "verbose=true, start_eta=0.001, tol=1e-12)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Array{Float64,2}:\n",
       " 830.145    913.16     1004.48     1104.92\n",
       " -24.1817   -24.1817    -24.1817    NaN   \n",
       "   1.63179    1.63389     1.63599   NaN   \n",
       "   1.23695    1.2379      1.23884   NaN   "
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ps = traj[3:end,145]; traj[:,143:146]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-24.181729489660793,[-2.30333e-6,-1.03767e-6],\n",
       "[NaN NaN; NaN NaN])"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "func = (x) -> J(rat_went_right[1:ntrials], \n",
    "clickrateL[1:ntrials], clickrateR[1:ntrials], nflashesL[1:ntrials], nflashesR[1:ntrials]; \n",
    "nderivs=length(x), difforder=2, make_dict(args, x)...)[1]\n",
    "\n",
    "cost, grad, hess = vgh(func, ps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dict{Any,Any} with 2 entries:\n",
       "  :alpha  => 1.63599\n",
       "  :lambda => 1.23884"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "make_dict(args, traj[3:end,145])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "InterruptException:",
     "output_type": "error",
     "traceback": [
      "InterruptException:",
      "",
      " in (::#kw##bbox_Hessian_keyword_minimization)(::Array{Any,1}, ::#bbox_Hessian_keyword_minimization, ::Array{Float64,1}, ::Array{String,1}, ::Dict{Symbol,Array{Float64,2}}, ::Function) at ./<missing>:0"
     ]
    }
   ],
   "source": [
    "# From Christine: alpha = .93, delta = .5, gamma = 3, softmax/lambda = 1.8, lapse = .008\n",
    "args = [\"alpha\", \"delta\"]\n",
    "seed = [0.4, 0.5]\n",
    "bbox = Dict(:alpha=>[0.0001 2], :delta=>[0.0001 10])\n",
    "\n",
    "ntrials = 10\n",
    "# ntrials = length(rat_went_right)\n",
    "\n",
    "params, cost, traj = bbox_Hessian_keyword_minimization(seed, args, bbox, (;params...) -> J(rat_went_right[1:ntrials], \n",
    "clickrateL[1:ntrials], clickrateR[1:ntrials], nflashesL[1:ntrials], nflashesR[1:ntrials]; params...)[1];\n",
    "verbose=true, start_eta=0.001)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Method definition trust_region_keyword_Hessian_minimization(Any, Any, Any) in module Main at In[45]:67 overwritten at In[51]:67.\n",
      "WARNING: Method definition #trust_region_keyword_Hessian_minimization(Array{Any, 1}, Main.#trust_region_keyword_Hessian_minimization, Any, Any, Any) in module Main overwritten.\n",
      "\u001b[1m\u001b[31mWARNING: replacing docs for 'trust_region_keyword_Hessian_minimization :: Tuple{Any,Any,Any}' in module 'Main'.\u001b[0m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "trust_region_keyword_Hessian_minimization"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "#############################################################################\n",
    "#                                                                           #\n",
    "#                   TRUST_REGION_KEYWORD_HESSIAN_MINIMIZATION               #\n",
    "#                                                                           #\n",
    "#############################################################################\n",
    "\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "function trust_region_keyword_Hessian_minimization(seed, args, func; start_eta=10, tol=1e-6, maxiter=400,\n",
    "    verbose=false)\n",
    "\n",
    "(below, x stands for delta_x, the step from the current x=x0 position at which the cost = const)\n",
    "\n",
    "cost = 0.5*x'*H*x + grad*x + const\n",
    "\n",
    "dcost/dx = H*x + grad  ;   dcost/dx = 0  ==> x =  - inv(H)*grad\n",
    "\n",
    "Trust-region says have a parameter lambda, and replace H with hat{H} = H +  I/eta.  \n",
    "When eta is very large, this is equivalent to a straight Newton method jump, \n",
    "because hat{H} ~= H.  But when eta is small, this is more like a small gradient\n",
    "descent step, because for small eta inv(hat{H}) ~= eta and therefore the delta x is like \n",
    "-eta*grad.  So, if the cost function is going down, make eta larger, and if it is going\n",
    "up, make eta a lot smaller. Just like we do in other adaptive methods\n",
    "\n",
    "PARAMETERS:\n",
    "===========\n",
    "\n",
    "seed        column vector, representing the starting value of the parameters.\n",
    "\n",
    "args        List of strings, representing args to be minimized.\n",
    "\n",
    "func        Function that takes a vector and returns a scalar.  If you want to\n",
    "            work with a function that tales more parameterrs and returns more than one \n",
    "            output, you can use something like\n",
    "\n",
    "                    x -> orig_func(x, other_params)[1]\n",
    "\n",
    "            You only need the \"[1]\" part if the orig_func returns more outputs than a scalar. \n",
    "\n",
    "OPTIONAL PARAMETERS:\n",
    "====================\n",
    "\n",
    "start_eta=10    Starting value of eta.  It's good to start with somethibg biggish, if it is\n",
    "                too much, it'll quickly get cut down.\n",
    "\n",
    "tol=1e-15       Numerical tolerance. If a proposed jump produces a change in func that is less than\n",
    "               this, the minimization stops.\n",
    "\n",
    "maxiter=400    Maximum number of iterations to do before stopping\n",
    "\n",
    "verbose=false   If true, print out a report on each iteration of iteration number, radius size (eta),\n",
    "                what type jump was proposed (\"Newton\" means going straight to global min, \"constrained\" means jump has \n",
    "                norm eta, failed means that finding the minimum at a given radius somehow didn't work). Will also\n",
    "                print out the cosine of the angle between the proposed jump and the gradient.\n",
    "\n",
    "RETURNS:\n",
    "========\n",
    "\n",
    "params       A vector the size of seed that has the last values of the minimizing parameters for func\n",
    "\n",
    "\"\"\"\n",
    "function trust_region_keyword_Hessian_minimization(seed, args, func; start_eta=10, tol=1e-15, maxiter=400,\n",
    "    verbose=false, verbose_level=1)\n",
    "\n",
    "    params = seed\n",
    "    eta = start_eta\n",
    "\n",
    "    cost, grad, hess = keyword_vgh(func, args, params)\n",
    "    if verbose && verbose_level >= 2\n",
    "        @printf(\"Initial cost=%g, grad=\", cost); print_vector_g(grad); print(\"\\n\")\n",
    "        @printf(\"hess = \"); print_vector_g(hess[:]); print(\"\\n\")\n",
    "    end\n",
    "    \n",
    "    traj = zeros(2+length(params),maxiter)\n",
    "\n",
    "    for i in [1:maxiter;]\n",
    "        hathess    = hess + eye(length(grad), length(grad))/eta        \n",
    "        new_params = params - inv(hathess)*grad\n",
    "        new_cost, new_grad, new_hess = keyword_vgh(func, args, new_params)\n",
    "            \n",
    "        if abs(new_cost - cost) < tol\n",
    "            break\n",
    "        end\n",
    "\n",
    "        if new_cost >= cost\n",
    "            eta = eta/2\n",
    "            costheta = NaN\n",
    "        else\n",
    "            eta = eta*1.1\n",
    "            costheta = dot(new_params-params, grad)/(norm(new_params-params)*norm(grad))\n",
    "\n",
    "            params = new_params\n",
    "            cost = new_cost\n",
    "            grad = new_grad\n",
    "            hess = new_hess\n",
    "        end\n",
    "            \n",
    "        traj[1,i] = eta\n",
    "        traj[2,i] = cost\n",
    "        traj[3:end,i] = params\n",
    "\n",
    "        if verbose\n",
    "            @printf \"%d: eta=%.3f cost=%.4f costheta=%.3f ps=\" i eta cost  costheta\n",
    "            print_vector(params)\n",
    "            @printf \"\\n\"\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    return params, cost, traj\n",
    "end\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "func = (x;params...) -> J(rat_went_right[1:ntrials], \n",
    "clickrateL[1:ntrials], clickrateR[1:ntrials], nflashesL[1:ntrials], nflashesR[1:ntrials]; make_dict(args, x)...)[1]\n",
    "\n",
    "e = 0.000001; d=[0,e]; (func(seed+d) - func(seed))/e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 219,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.4736842105263158,[0.277008],\n",
       "[-0.291588])"
      ]
     },
     "execution_count": 219,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "args = [\"alpha\"]\n",
    "seed = [0.9]\n",
    "func = (x;pars...) -> J(rat_went_right, clickrateL, \n",
    "    clickrateR, nflashesL, nflashesR; make_dict(args, x)...)[1]\n",
    "\n",
    "\n",
    "args = [\"delta\"]\n",
    "seed = [0.9]\n",
    "ntrials = 2783 # 2783\n",
    "func = (x;pars...) -> J(rat_went_right[ntrials], clickrateL[ntrials], clickrateR[ntrials], nflashesL[ntrials], nflashesR[ntrials]; nderivs=length(args), difforder=2, make_dict(args, x)...)[1]\n",
    "func = (x;pars...) -> one_trial_cost(rat_went_right[ntrials], clickrateL[ntrials], clickrateR[ntrials], nflashesL[ntrials], nflashesR[ntrials]; make_dict(args, x)...)\n",
    "func = (x;pars...) -> prob_right(clickrateL[ntrials], clickrateR[ntrials], nflashesL[ntrials], nflashesR[ntrials]; make_dict(args, x)...)\n",
    "func = (x;pars...) -> prob_weight(nflashesR[ntrials]; make_dict(args, x)...)\n",
    "\n",
    "\n",
    "cost, grad, hess = vgh(func, seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-4.8283137373023015,\n",
       "[-4.82831])"
      ]
     },
     "execution_count": 194,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "strial = 27\n",
    "gu = J(rat_went_right[strial], clickrateL[strial], clickrateR[strial], nflashesL[strial], nflashesR[strial])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "56265×1 Array{Float64,2}:\n",
       " 10.0\n",
       " 10.0\n",
       " 10.0\n",
       " 10.0\n",
       "  0.0\n",
       " 10.0\n",
       "  0.0\n",
       " 10.0\n",
       " 10.0\n",
       "  0.0\n",
       " 10.0\n",
       "  0.0\n",
       "  0.0\n",
       "  ⋮  \n",
       "  3.0\n",
       " 10.0\n",
       " 10.0\n",
       "  1.0\n",
       "  1.0\n",
       "  2.0\n",
       " 10.0\n",
       " 10.0\n",
       "  2.0\n",
       "  0.0\n",
       "  0.0\n",
       "  0.0"
      ]
     },
     "execution_count": 196,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nflashesR\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0: eta=0.001 ps=[0.930]\n",
      "Constrained parabolic minimization failed with error ArgumentError(\"collection must be non-empty\")\n",
      "\n",
      "eta was 0.001\n",
      "grad was\n",
      "[NaN]\n",
      "\n",
      "hess was\n",
      "[NaN]\n",
      "\n",
      "1: eta=0.0005 cost=-176622.7834 jtype=failed costheta=NaN ps=[0.930]\n",
      "Constrained parabolic minimization failed with error ArgumentError(\"collection must be non-empty\")\n",
      "\n",
      "eta was 0.0005\n",
      "grad was\n",
      "[NaN]\n",
      "\n",
      "hess was\n",
      "[NaN]\n",
      "\n",
      "2: eta=0.00025 cost=-176622.7834 jtype=failed costheta=NaN ps=[0.930]\n",
      "Constrained parabolic minimization failed with error ArgumentError(\"collection must be non-empty\")\n",
      "\n",
      "eta was 0.00025\n",
      "grad was\n",
      "[NaN]\n",
      "\n",
      "hess was\n",
      "[NaN]\n",
      "\n",
      "3: eta=0.000125 cost=-176622.7834 jtype=failed costheta=NaN ps=[0.930]\n",
      "Constrained parabolic minimization failed with error ArgumentError(\"collection must be non-empty\")\n",
      "\n",
      "eta was 0.000125\n",
      "grad was\n",
      "[NaN]\n",
      "\n",
      "hess was\n",
      "[NaN]\n",
      "\n",
      "4: eta=6.25e-05 cost=-176622.7834 jtype=failed costheta=NaN ps=[0.930]\n",
      "Constrained parabolic minimization failed with error ArgumentError(\"collection must be non-empty\")\n",
      "\n",
      "eta was 6.25e-05\n",
      "grad was\n",
      "[NaN]\n",
      "\n",
      "hess was\n",
      "[NaN]\n",
      "\n",
      "5: eta=3.125e-05 cost=-176622.7834 jtype=failed costheta=NaN ps=[0.930]\n",
      "Constrained parabolic minimization failed with error ArgumentError(\"collection must be non-empty\")\n",
      "\n",
      "eta was 3.125e-05\n",
      "grad was\n",
      "[NaN]\n",
      "\n",
      "hess was\n",
      "[NaN]\n",
      "\n",
      "6: eta=1.5625e-05 cost=-176622.7834 jtype=failed costheta=NaN ps=[0.930]\n",
      "Constrained parabolic minimization failed with error ArgumentError(\"collection must be non-empty\")\n",
      "\n",
      "eta was 1.5625e-05\n",
      "grad was\n",
      "[NaN]\n",
      "\n",
      "hess was\n",
      "[NaN]\n",
      "\n",
      "7: eta=7.8125e-06 cost=-176622.7834 jtype=failed costheta=NaN ps=[0.930]\n",
      "Constrained parabolic minimization failed with error ArgumentError(\"collection must be non-empty\")\n",
      "\n",
      "eta was 7.8125e-06\n",
      "grad was\n",
      "[NaN]\n",
      "\n",
      "hess was\n",
      "[NaN]\n",
      "\n",
      "8: eta=3.90625e-06 cost=-176622.7834 jtype=failed costheta=NaN ps=[0.930]\n",
      "Constrained parabolic minimization failed with error ArgumentError(\"collection must be non-empty\")\n",
      "\n",
      "eta was 3.90625e-06\n",
      "grad was\n",
      "[NaN]\n",
      "\n",
      "hess was\n",
      "[NaN]\n",
      "\n",
      "9: eta=1.95313e-06 cost=-176622.7834 jtype=failed costheta=NaN ps=[0.930]\n",
      "Constrained parabolic minimization failed with error ArgumentError(\"collection must be non-empty\")\n",
      "\n",
      "eta was 1.95313e-06\n",
      "grad was\n",
      "[NaN]\n",
      "\n",
      "hess was\n",
      "[NaN]\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1-element Array{Float64,1}:\n",
       " 0.93"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seed = [0.93]\n",
    "args = [\"alpha\"]\n",
    "bbox = [0 3]\n",
    "func = (;pars...) -> J(rat_went_right, clickrateL, clickrateR, nflashesL, nflashesR; pars...)\n",
    "\n",
    "params, traj = bbox_Hessian_keyword_minimization(seed, args, bbox, func, wallwidth_factor=0.001, verbose=true, start_eta=0.001)\n",
    "params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Base.KeyIterator for a Dict{String,Any} with 12 entries. Keys:\n",
       "  \"rightvolume\"\n",
       "  \"went_right\"\n",
       "  \"rclickHz\"\n",
       "  \"lclickHz\"\n",
       "  \"lefthit\"\n",
       "  \"rflashes\"\n",
       "  \"lflashes\"\n",
       "  \"leftvolume\"\n",
       "  \"righthit\"\n",
       "  \"leftprob\"\n",
       "  \"rightprob\"\n",
       "  \"cues_lowtohigh\""
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "keys(matread(\"J266.mat\")[\"J266_data\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 228,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([4018.0,2215.0,2045.0,1973.0,1892.0,2749.0,1980.0,1997.0,1951.0,35445.0],[0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0],Any[PyObject <matplotlib.patches.Rectangle object at 0x32690b510>,PyObject <matplotlib.patches.Rectangle object at 0x327b005d0>,PyObject <matplotlib.patches.Rectangle object at 0x327b1ead0>,PyObject <matplotlib.patches.Rectangle object at 0x33943a510>,PyObject <matplotlib.patches.Rectangle object at 0x33943ab90>,PyObject <matplotlib.patches.Rectangle object at 0x33944d250>,PyObject <matplotlib.patches.Rectangle object at 0x33944d8d0>,PyObject <matplotlib.patches.Rectangle object at 0x33944df50>,PyObject <matplotlib.patches.Rectangle object at 0x339459610>,PyObject <matplotlib.patches.Rectangle object at 0x339459c90>])"
      ]
     },
     "execution_count": 228,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "figure(1); clf()\n",
    "h = plt[:hist](nflashesL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Method definition trust_region_Hessian_minimization(Any, Any) in module Main at In[42]:65 overwritten at In[50]:65.\n",
      "WARNING: Method definition #trust_region_Hessian_minimization(Array{Any, 1}, Main.#trust_region_Hessian_minimization, Any, Any) in module Main overwritten.\n",
      "\u001b[1m\u001b[31mWARNING: replacing docs for 'trust_region_Hessian_minimization :: Tuple{Any,Any}' in module 'Main'.\u001b[0m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "trust_region_Hessian_minimization"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "#############################################################################\n",
    "#                                                                           #\n",
    "#                   TRUST_REGION_HESSIAN_MINIMIZATION                       #\n",
    "#                                                                           #\n",
    "#############################################################################\n",
    "\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "function trust_region_Hessian_minimization(seed, func; start_eta=10, tol=1e-6, maxiter=400,\n",
    "    verbose=false)\n",
    "\n",
    "(below, x stands for delta_x, the step from the current x=x0 position at which the cost = const)\n",
    "\n",
    "cost = 0.5*x'*H*x + grad*x + const\n",
    "\n",
    "dcost/dx = H*x + grad  ;   dcost/dx = 0  ==> x =  - inv(H)*grad\n",
    "\n",
    "Trust-region says have a parameter lambda, and replace H with hat{H} = H +  I/eta.  \n",
    "When eta is very large, this is equivalent to a straight Newton method jump, \n",
    "because hat{H} ~= H.  But when eta is small, this is more like a small gradient\n",
    "descent step, because for small eta inv(hat{H}) ~= eta and therefore the delta x is like \n",
    "-eta*grad.  So, if the cost function is going down, make eta larger, and if it is going\n",
    "up, make eta a lot smaller. Just like we do in other adaptive methods\n",
    "\n",
    "PARAMETERS:\n",
    "===========\n",
    "\n",
    "seed        column vector, representing the starting value of the parameters.\n",
    "\n",
    "func        Function that takes a vector and returns a scalar.  If you want to\n",
    "            work with a function that tales more parameterrs and returns more than one \n",
    "            output, you can use something like\n",
    "\n",
    "                    x -> orig_func(x, other_params)[1]\n",
    "\n",
    "            You only need the \"[1]\" part if the orig_func returns more outputs than a scalar. \n",
    "\n",
    "OPTIONAL PARAMETERS:\n",
    "====================\n",
    "\n",
    "start_eta=10    Starting value of eta.  It's good to start with somethibg biggish, if it is\n",
    "                too much, it'll quickly get cut down.\n",
    "\n",
    "tol=1e-15       Numerical tolerance. If a proposed jump produces a change in func that is less than\n",
    "               this, the minimization stops.\n",
    "\n",
    "maxiter=400    Maximum number of iterations to do before stopping\n",
    "\n",
    "verbose=false   If true, print out a report on each iteration of iteration number, radius size (eta),\n",
    "                what type jump was proposed (\"Newton\" means going straight to global min, \"constrained\" means jump has \n",
    "                norm eta, failed means that finding the minimum at a given radius somehow didn't work). Will also\n",
    "                print out the cosine of the angle between the proposed jump and the gradient.\n",
    "\n",
    "RETURNS:\n",
    "========\n",
    "\n",
    "params       A vector the size of seed that has the last values of the minimizing parameters for func\n",
    "\n",
    "\"\"\"\n",
    "function trust_region_Hessian_minimization(seed, func; start_eta=10, tol=1e-15, maxiter=400,\n",
    "    verbose=false, verbose_level=1)\n",
    "\n",
    "    params = seed\n",
    "    eta = start_eta\n",
    "\n",
    "    cost, grad, hess = vgh(func, params)\n",
    "    if verbose && verbose_level >= 2\n",
    "        @printf(\"Initial cost, grad, hess:\\n\")\n",
    "        print_vector_g(:cost)\n",
    "        print_vector_g(:grad)\n",
    "        print_vector_g(:hess)\n",
    "    end\n",
    "    \n",
    "    traj = zeros(2+length(params),maxiter)\n",
    "\n",
    "    for i in [1:maxiter;]\n",
    "        hathess    = hess + eye(length(grad), length(grad))/eta        \n",
    "        new_params = params - inv(hathess)*grad\n",
    "        new_cost, new_grad, new_hess = vgh(func, new_params)\n",
    "            \n",
    "        if abs(new_cost - cost) < tol\n",
    "            break\n",
    "        end\n",
    "\n",
    "        if new_cost >= cost\n",
    "            eta = eta/2\n",
    "            costheta = NaN\n",
    "        else\n",
    "            eta = eta*1.1\n",
    "            costheta = dot(new_params-params, grad)/(norm(new_params-params)*norm(grad))\n",
    "\n",
    "            params = new_params\n",
    "            cost = new_cost\n",
    "            grad = new_grad\n",
    "            hess = new_hess\n",
    "        end\n",
    "\n",
    "        traj[1,i] = eta\n",
    "        traj[2,i] = cost\n",
    "        traj[3:end,i] = params\n",
    "        \n",
    "        if verbose\n",
    "            @printf \"%d: eta=%.3f cost=%.4f costheta=%.3f ps=\" i eta cost  costheta\n",
    "            print_vector(params)\n",
    "            @printf \"\\n\"\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    return params, cost, traj\n",
    "end\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.2",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}