{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"../../L-BFGS/OracleFunction.jl\")\n",
    "include(\"../../L-BFGS/LBFGS.jl\")\n",
    "include(\"../../utilities/genFunc.jl\")\n",
    "using .LBFGS\n",
    "using .OracleFunction\n",
    "using LinearAlgebra, BenchmarkTools, CSV, DataFrames\n",
    "\n",
    "baseDir = joinpath(\"../\", \"results/LBFGS/ill_conditioned/\")\n",
    "mkpath(baseDir);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LeastSquaresF{Float64, Main.OracleFunction.var\"#f#1\"{Matrix{Float64}, Vector{Float64}}, Main.OracleFunction.var\"#df#2\"{Matrix{Float64}, Vector{Float64}}}(OracleF{Float64, Main.OracleFunction.var\"#f#1\"{Matrix{Float64}, Vector{Float64}}, Main.OracleFunction.var\"#df#2\"{Matrix{Float64}, Vector{Float64}}}([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], Main.OracleFunction.var\"#f#1\"{Matrix{Float64}, Vector{Float64}}([0.2592586022640595 0.6963535943217778 … 0.2686596970120587 -0.36678046257922436; 0.514702105566236 -0.7473274556601393 … 0.19887814981786733 0.8663661527822621; … ; 0.0 0.0 … 1.0e-12 0.0; 0.0 0.0 … 0.0 1.0e-12], [7.689547521112595, -23.188865118086323, 13.727391619018281, -4.242513852908478, 23.90743053982707, 1.8343747758964255, -1.159884201591351, 2.216694337124031, 2.730975021441503, -19.180979473592664  …  -6.088266825619777e-13, 1.1943252953272254e-13, 4.048371886383304e-13, -6.926115820555791e-13, 3.8752674617747763e-13, 6.919446457660259e-13, 5.147174879596801e-13, -5.887351028568868e-13, 6.537120155767028e-13, 3.946146858948962e-13]), Main.OracleFunction.var\"#df#2\"{Matrix{Float64}, Vector{Float64}}([0.2592586022640595 0.6963535943217778 … 0.2686596970120587 -0.36678046257922436; 0.514702105566236 -0.7473274556601393 … 0.19887814981786733 0.8663661527822621; … ; 0.0 0.0 … 1.0e-12 0.0; 0.0 0.0 … 0.0 1.0e-12], [7.689547521112595, -23.188865118086323, 13.727391619018281, -4.242513852908478, 23.90743053982707, 1.8343747758964255, -1.159884201591351, 2.216694337124031, 2.730975021441503, -19.180979473592664  …  -6.088266825619777e-13, 1.1943252953272254e-13, 4.048371886383304e-13, -6.926115820555791e-13, 3.8752674617747763e-13, 6.919446457660259e-13, 5.147174879596801e-13, -5.887351028568868e-13, 6.537120155767028e-13, 3.946146858948962e-13])), [0.2592586022640595 0.6963535943217778 … 0.2686596970120587 -0.36678046257922436; 0.514702105566236 -0.7473274556601393 … 0.19887814981786733 0.8663661527822621; … ; 0.0 0.0 … 1.0e-12 0.0; 0.0 0.0 … 0.0 1.0e-12], [7.689547521112595, -23.188865118086323, 13.727391619018281, -4.242513852908478, 23.90743053982707, 1.8343747758964255, -1.159884201591351, 2.216694337124031, 2.730975021441503, -19.180979473592664  …  -6.088266825619777e-13, 1.1943252953272254e-13, 4.048371886383304e-13, -6.926115820555791e-13, 3.8752674617747763e-13, 6.919446457660259e-13, 5.147174879596801e-13, -5.887351028568868e-13, 6.537120155767028e-13, 3.946146858948962e-13], [6.711447540940819 -2.307272578298896 … 1.2541719054073133 -1.8639946218858157; -2.307272578298896 8.480411930046433 … -1.627313075365743 0.5078496705429104; … ; 1.2541719054073133 -1.627313075365743 … 7.4470442813402595 -4.4848817373554315; -1.8639946218858157 0.5078496705429104 … -4.4848817373554315 6.7106397848856], [-3.33043813332588 57.71972704197855 … 26.336102875181652 -62.14554952332052])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# parameters for exact LS vs AWLS vs quadratic function \n",
    "λ        = 10^-12\n",
    "ϵ        = 10^-14\n",
    "memsize  = 7\n",
    "maxIters = 1:200\n",
    "m        = 1000\n",
    "n        = 20\n",
    "num_trials = 20\n",
    "\n",
    "gf = genFunc(:exactRandDataset, λ=λ, m=m, n=n)\n",
    "non_quadratic = OracleF(ones(size(gf[:X_hat], 2)),\n",
    "        (x) -> norm(gf[:X_hat] * x - gf[:y_hat]),\n",
    "        (x) -> inv(norm(gf[:X_hat] * x - gf[:y_hat])) * gf[:X_hat]' * (gf[:X_hat] * x - gf[:y_hat])\n",
    "    )\n",
    "quadratic = OracleF(ones(size(gf[:X_hat], 2)),\n",
    "        (x) -> norm(gf[:X_hat] * x - gf[:y_hat])^2,\n",
    "        (x) -> 2 * gf[:X_hat]' * (gf[:X_hat] * x - gf[:y_hat])\n",
    "    )\n",
    "ls = LeastSquaresF(gf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done trial 1\n",
      "Done trial 2\n",
      "Done trial 3\n",
      "Done trial 4\n",
      "Done trial 5\n",
      "Done trial 6\n",
      "Done trial 7\n",
      "Done trial 8\n",
      "Done trial 9\n",
      "Done trial 10\n",
      "Done trial 11\n",
      "Done trial 12\n",
      "Done trial 13\n",
      "Done trial 14\n",
      "Done trial 15\n",
      "Done trial 16\n",
      "Done trial 17\n",
      "Done trial 18\n",
      "Done trial 19\n",
      "Done trial 20\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\"../results/LBFGS/ill_conditioned/statisticsLBFGS-iterations-m1000n20--error-norm.csv\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "### residual, relative error and norm of gradient with respect to iterations with exact line search\n",
    "using Statistics, CSV\n",
    "\n",
    "# Preallocate arrays\n",
    "num_iterations = length(maxIters)\n",
    "\n",
    "gradients = zeros(num_trials, num_iterations)\n",
    "residuals = zeros(num_trials, num_iterations)\n",
    "relative_errors = zeros(num_trials, num_iterations)\n",
    "\n",
    "for trial in 1:num_trials\n",
    "\n",
    "    gf = genFunc(:exactRandDataset, λ=λ, m=m, n=n)\n",
    "    ls = LeastSquaresF(gf)\n",
    "\n",
    "    for (i, maxIter) in enumerate(maxIters)\n",
    "        t = LimitedMemoryBFGS(ls, ϵ=ϵ, MaxEvaluations=maxIter, m=memsize)\n",
    "\n",
    "        relative_errors[trial, i] = norm(t[:x] - gf[:w_star]) / norm(gf[:w_star])\n",
    "        residuals[trial, i] = norm(gf[:X_hat] * t[:x] - gf[:y_hat]) / norm(gf[:y_hat])\n",
    "        gradients[trial, i] = norm(t[:grad])\n",
    "    end\n",
    "\n",
    "    println(\"Done trial \", trial)\n",
    "    \n",
    "end\n",
    "\n",
    "# Calculate mean and standard deviation\n",
    "mean_relative = mean(relative_errors, dims=1)'\n",
    "std_relative = std(relative_errors, dims=1)'\n",
    "mean_residual = mean(residuals, dims=1)'\n",
    "std_residual = std(residuals, dims=1)'\n",
    "mean_gradient = mean(gradients, dims=1)'\n",
    "std_gradient = std(gradients, dims=1)'\n",
    "\n",
    "\n",
    "# Write results to CSV\n",
    "outputvsc = joinpath(baseDir, \"statisticsLBFGS-iterations-m\" * string(m) * \"n\" * string(n) * \"--error-norm.csv\");\n",
    "\n",
    "accData = Dict(\n",
    "    :lambda        => Array{Float64}(undef, 0),\n",
    "    :epsilon       => Array{Float64}(undef, 0),\n",
    "    :memsize       => Array{Int64}(undef, 0),\n",
    "    :maxiterations => Array{Int64}(undef, 0),\n",
    "    :mean_relative => Array{Float64}(undef, 0),\n",
    "    :std_relative  => Array{Float64}(undef, 0),\n",
    "    :mean_residual => Array{Float64}(undef, 0),\n",
    "    :std_residual  => Array{Float64}(undef, 0),\n",
    "    :mean_gradient => Array{Float64}(undef, 0),\n",
    "    :std_gradient  => Array{Float64}(undef, 0)\n",
    ")\n",
    "\n",
    "# create dataframe with columns from arrays\n",
    "for maxIter ∈ maxIters\n",
    "    push!(accData[:lambda], λ)\n",
    "    push!(accData[:epsilon], ϵ)\n",
    "    push!(accData[:memsize], memsize)\n",
    "    push!(accData[:maxiterations], maxIter)\n",
    "    push!(accData[:mean_relative], mean_relative[maxIter])\n",
    "    push!(accData[:std_relative], std_relative[maxIter])\n",
    "    push!(accData[:mean_residual], mean_residual[maxIter])\n",
    "    push!(accData[:std_residual], std_residual[maxIter])\n",
    "    push!(accData[:mean_gradient], mean_gradient[maxIter])\n",
    "    push!(accData[:std_gradient], std_gradient[maxIter])\n",
    "\n",
    "end\n",
    "\n",
    "\n",
    "CSV.write(outputvsc, accData)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "### relative error with respect to iterations and memeory size\n",
    "# output csv\n",
    "outputvsc = joinpath(baseDir, \"statisticsLBFGS-iterations-m\" * string(m) * \"n\" * string(n) * \"--memsize.csv\");\n",
    "accData = Dict(\n",
    "\n",
    "\n",
    "    :memsize       => Array{Int64}(undef, 0),\n",
    "    :iterations    => Array{Int64}(undef, 0),\n",
    "    :maxiterations => Array{Int64}(undef, 0),\n",
    "    :relative      => Array{Float64}(undef, 0),\n",
    "    :residual      => Array{Float64}(undef, 0),\n",
    "    :gradient      => Array{Float64}(undef, 0)\n",
    "    )\n",
    "\n",
    "memsizes = [1 3 5 7 9 11]\n",
    "\n",
    "for memsize ∈ memsizes\n",
    "    gf = genFunc(:exactRandDataset, λ=λ, m=m, n=n)\n",
    "    ls = LeastSquaresF(gf)\n",
    "    for maxIter ∈ maxIters\n",
    "\n",
    "        t = LimitedMemoryBFGS(ls, ϵ=ϵ, MaxEvaluations=maxIter, m=memsize)\n",
    "\n",
    "        relative_error = norm(t[:x] - gf[:w_star]) / norm(gf[:w_star])\n",
    "        residual = norm(gf[:X_hat] * t[:x] - gf[:y_hat]) / norm(gf[:y_hat])\n",
    "\n",
    "        iterations = maxIter - t[:RemainingEvaluations]\n",
    "\n",
    "        push!(accData[:memsize], memsize)\n",
    "        push!(accData[:iterations], iterations)\n",
    "        push!(accData[:maxiterations], maxIter)\n",
    "        push!(accData[:relative], relative_error)\n",
    "        push!(accData[:residual], residual)\n",
    "        push!(accData[:gradient], norm(t[:grad]))\n",
    "        flush(stdout)\n",
    "    end\n",
    "end\n",
    "\n",
    "CSV.write(outputvsc, DataFrame(accData));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.9.3",
   "language": "julia",
   "name": "julia-1.9"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.9.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}