{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"../../L-BFGS/OracleFunction.jl\")\n",
    "include(\"../../L-BFGS/BFGS.jl\")\n",
    "include(\"../../utilities/genFunc.jl\")\n",
    "using .BFGS\n",
    "using .OracleFunction\n",
    "using LinearAlgebra, BenchmarkTools, CSV, DataFrames\n",
    "\n",
    "baseDir = joinpath(\"../\", \"results/LBFGS/comparison_BFGS/\")\n",
    "mkpath(baseDir);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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.6269851049934334 -0.9953855131718308 … 0.26819190323331954 -0.6780465661809252; -0.9332327939739362 0.8084805224537295 … 0.28452273460460176 0.9029306010736042; … ; 0.0 0.0 … 10.0 0.0; 0.0 0.0 … 0.0 10.0], [-10.186364190714537, -1.6231132115458327, 15.375526057866429, 1.5006156140417093, -14.546392908739245, -4.939853634188801, -0.2053418568399169, 20.47017872706546, 1.0276498922046493, -12.900436004701787  …  3.3479087325325874, 4.529977532247315, 9.395681468606902, -9.098074604739619, -8.545976276298948, 5.1912832757085265, -8.133605566489734, -7.81791180224994, 8.53647733342897, 0.2170916132884937]), Main.OracleFunction.var\"#df#2\"{Matrix{Float64}, Vector{Float64}}([0.6269851049934334 -0.9953855131718308 … 0.26819190323331954 -0.6780465661809252; -0.9332327939739362 0.8084805224537295 … 0.28452273460460176 0.9029306010736042; … ; 0.0 0.0 … 10.0 0.0; 0.0 0.0 … 0.0 10.0], [-10.186364190714537, -1.6231132115458327, 15.375526057866429, 1.5006156140417093, -14.546392908739245, -4.939853634188801, -0.2053418568399169, 20.47017872706546, 1.0276498922046493, -12.900436004701787  …  3.3479087325325874, 4.529977532247315, 9.395681468606902, -9.098074604739619, -8.545976276298948, 5.1912832757085265, -8.133605566489734, -7.81791180224994, 8.53647733342897, 0.2170916132884937])), [0.6269851049934334 -0.9953855131718308 … 0.26819190323331954 -0.6780465661809252; -0.9332327939739362 0.8084805224537295 … 0.28452273460460176 0.9029306010736042; … ; 0.0 0.0 … 10.0 0.0; 0.0 0.0 … 0.0 10.0], [-10.186364190714537, -1.6231132115458327, 15.375526057866429, 1.5006156140417093, -14.546392908739245, -4.939853634188801, -0.2053418568399169, 20.47017872706546, 1.0276498922046493, -12.900436004701787  …  3.3479087325325874, 4.529977532247315, 9.395681468606902, -9.098074604739619, -8.545976276298948, 5.1912832757085265, -8.133605566489734, -7.81791180224994, 8.53647733342897, 0.2170916132884937], [106.26650786511657 -1.0664727828949865 … 2.4445684049045426 0.0322174082238121; -1.0664727828949865 106.1318558809136 … -0.6525752258209305 1.2988114776625301; … ; 2.4445684049045426 -0.6525752258209305 … 106.71985607542891 0.16719323022752428; 0.0322174082238121 1.2988114776625301 … 0.16719323022752428 105.42561554395674], [2.0243350631400165 57.59287736844317 … 68.90017274988415 -7.460256118037204])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# parameters for exact LS vs AWLS vs quadratic function \n",
    "λ        = 10^1\n",
    "ϵ        = 10^-14\n",
    "maxIters = 1:200\n",
    "m        = 1000\n",
    "n        = 20\n",
    "num_trials = 10\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": 7,
   "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"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\"../results/LBFGS/comparison_BFGS/statisticsBFGS-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 = BroydenFletcherGoldfarbShanno(ls, ϵ=ϵ, MaxEvaluations=maxIter)\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, \"statisticsBFGS-iterations-m\" * string(m) * \"n\" * string(n) * \"--error-norm.csv\");\n",
    "\n",
    "accData = Dict(\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[: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": 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
}