text:
using JuMP
using Ipopt

include("../OptMaps.jl")
using .OptMaps

function solve_model(solver, model_log, solver_options=Dict())
    model = Model(solver)

    for (k, v) in solver_options
        set_optimizer_attribute(model, k, v)
    end

    n_variables = 10
    n_constraints = 1

    @variable(model, -3.0 <= x[1:n_variables] <= 5.0)
    @NLconstraint(model, sum((x[i]^2 for i in 1:n_variables)) >= 5.0)

    constraints_f = []
    push!(constraints_f, (v) -> sum(v.^2) >= 5.0)

    @NLobjective(model, Min, sum(x[i]^2 for i in 1:n_variables))
    obj_f(v) = sum(v.^2)

    starting_point = zeros(Float64, n_variables) .+ 3.0
    solution_point = zeros(Float64, n_variables)

    ###############################################################
    #save model textual description
    open(model_log, "w") do f
        print(f, model)

        ###############################################################
        #solve model
        elapsed_time = @elapsed optimize!(model)

        println(f, "-----------------------------------------------")
        println(f, "Elapsed time: $elapsed_time")
        println(f, "-----------------------------------------------")
        println(f, "Status:    $(termination_status(model))")
        println(f, "Objective: $(objective_value(model))")
        println(f, "-----------------------------------------------")
        println(f, "Variables:")

        for i in 1:n_variables
            println(f, "$(x[i]): $(value(x[i]))")
            solution_point[i] = value(x[i])
        end
    end

    variable_names = vcat(["x[$(string(i))]" for i in 1:n_variables])
    constraint_names = ["c1"]

    obj_val = objective_value(model)

    return obj_f, constraints_f, variable_names, constraint_names, starting_point, solution_point, obj_val
end

optimizers = [Ipopt.Optimizer]
opt_names = ["ipopt"]

problems = Dict()

for (opt, opt_name) in zip(optimizers, opt_names)
    objective_f, constraints_f, variable_names, constraint_names, starting_point, solution_point, obj_val = solve_model(opt, "sphere_$(opt_name).log")
    problems[opt_name] = (objective_f, constraints_f, variable_names, constraint_names, starting_point, solution_point, obj_val)
end

#All functions are equal, get first to initialize OptMap
objective_f, constraints_f, variable_names, constraint_names, _, _, obj_val = problems[opt_names[1]]
variable_domains = convert(Array{Tuple{Float64,Float64},1}, [])

for i in 1:length(variable_names)
    push!(variable_domains, (-4.0,5.0))
end

O = OptMap(variable_names, variable_domains, objective_f; constraint_names=constraint_names, constraints_f=constraints_f, verbose=true)
create_map!(O)

for opt_name in opt_names
    _, _, _, _, starting_point, solution_point, _ = problems[opt_name]

    path_to_solution = zeros(Float64, length(variable_domains), 2)
    path_to_solution[:,1] = starting_point
    path_to_solution[:,2] = solution_point

    draw_solution!(O, path_to_solution; line_width=3.0)
    save_all(O, "sphere_$(opt_name)", ".")
end

max_stars_repo_name: nmaterise/optmap

max_issues_repo_name: nmaterise/optmap

max_forks_repo_name: nmaterise/optmap


text:
using StanSample
using RegressionAndOtherStories
using Test

df = DataFrame(x = 1:10, y = 1:10)

@testset "test_hibbs.jl" begin

    mod_sum1 = model_summary(df, [:x, :y])
    @test mod_sum1[:x, :median] ≈ 5.5 atol=0.1

    mod_sum2 = model_summary(df, ["x", "y"])
    @test mod_sum2[:x, :median] ≈ 5.5 atol=0.1

    mod_sum3 = model_summary(df, ["x", "y", "z"])
    @test size(mod_sum3) == (2, 4)
    @test mod_sum3[:x, :median] ≈ 5.5 atol=0.1

    mod_sum4 = model_summary(df, ["z"])
    @test mod_sum4 === nothing

end

ext: jl


text:
# call_fast_LP.jl
# Wrapper for fast_LP.jl, for Sloan computing clusters
# Authored by Arthur J Delarue on 10/28/15

include("intro.jl")
include("fast_LP.jl")

function load_parameters()
	if length(ARGS) == 0
		println("Error: shell script failed to pass in arguments")
		return nothing
	else
		JSONfileName = ARGS[1]
		params = JSON.parsefile(JSONfileName)
		return JSONfileName, params
	end
end

t=time()
# Load prameters
JSONfile, params = load_parameters()
# Load Manhattan
manhattan = loadCityGraph()
ws = params["WARMSTART"]
if ws
	startTimes = load("Outputs/fast_r200_minr1_i100_wd_1214_m112_y2013_relaxed_smooth_ppc3_dataaverage_errorboth_dynConstr_st1000_add1000_every1_tc0.0/manhattan-times-100.jld", "times")
else
	startTimes = manhattan.roadTime
end
# Load data
if params["LOCALIZED_TESTING"]
	travel_times, num_rides, trainTestDf, testing_travel_times, testing_num_rides, testDf = loadInputTravelTimes(manhattan.positions, params["METHOD"], year=params["YEAR"], startTime=params["START_TIME"], endTime=params["END_TIME"], startMonth=params["START_MONTH"], endMonth=params["END_MONTH"], radius = params["RADIUS"], split_training_type = "localized")
else
	travel_times, num_rides, trainTestDf, testing_travel_times, testing_num_rides, testDf = loadInputTravelTimes(manhattan.positions, params["METHOD"], year=params["YEAR"], startTime=params["START_TIME"], endTime=params["END_TIME"], startMonth=params["START_MONTH"], endMonth=params["END_MONTH"], radius = params["RADIUS"], split_training_type = "random")
end
# Choose constraints if random
if params["RANDOM_CONSTRAINTS"]
	travel_times, num_rides = chooseConstraints(travel_times, num_rides, sample_size=params["SAMPLE_SIZE"]);
end
println("**** Launching computation ****")
# Launch ! (three cases for different error computations)
if params["ERROR_COMPUTATION"] == "single-ride"
	@time status, new_times, directory = fast_LP(manhattan, travel_times, num_rides, trainTestDf, startTimes, errorComputation=params["ERROR_COMPUTATION"], model_type=params["MODEL"], last_smooth=params["LAST_SMOOTH"], max_rounds=params["MAX_ROUNDS"], min_rides=params["MIN_RIDES"], radius=params["RADIUS"], year=params["YEAR"], startTime=params["START_TIME"], endTime=params["END_TIME"], startMonth=params["START_MONTH"], endMonth=params["END_MONTH"], method=params["METHOD"], sample_size=params["SAMPLE_SIZE"], turnCost=params["TURN_COST"], turnCostAsVariable=params["TURN_COST_AS_VARIABLE"], delta_bound=params["DELTA_BOUND"], maxNumPathsPerOD=params["MAX_NUM_PATHS_PER_OD"], computeFinalSP = params["COMPUTE_FINAL_SHORTEST_PATHS"], startWithSimpleLP=params["START_SIMPLE"], randomConstraints=params["RANDOM_CONSTRAINTS"], dynamicConstraints=params["DYNAMIC_CONSTRAINTS"], numPairsToAdd = params["NUM_OD_ADDED"], numPairsToRemove = params["NUM_OD_REMOVED"], addOnlyIfAbovePercentile=params["ADD_IF_ABOVE_PERCENTILE"], removeOnlyIfBelowPercentile=params["REMOVE_IF_BELOW_PERCENTILE"], globalConstraintUpdate=params["GLOBAL_CONSTRAINT_UPDATE"],iterationMultiple = params["UPDATE_EVERY_N_ITERATIONS"], warmstart=params["WARMSTART"], metropolis=params["METROPOLIS"], localized_testing=params["LOCALIZED_TESTING"])
elseif params["ERROR_COMPUTATION"] == "average"
	@time status, new_times, directory = fast_LP(manhattan, travel_times, num_rides, testing_travel_times, startTimes, errorComputation=params["ERROR_COMPUTATION"], model_type=params["MODEL"], last_smooth=params["LAST_SMOOTH"], max_rounds=params["MAX_ROUNDS"], min_rides=params["MIN_RIDES"], radius=params["RADIUS"], year=params["YEAR"], startTime=params["START_TIME"], endTime=params["END_TIME"], startMonth=params["START_MONTH"], endMonth=params["END_MONTH"], method=params["METHOD"], sample_size=params["SAMPLE_SIZE"], turnCost=params["TURN_COST"], turnCostAsVariable=params["TURN_COST_AS_VARIABLE"], delta_bound=params["DELTA_BOUND"], maxNumPathsPerOD=params["MAX_NUM_PATHS_PER_OD"], computeFinalSP = params["COMPUTE_FINAL_SHORTEST_PATHS"], startWithSimpleLP=params["START_SIMPLE"], randomConstraints=params["RANDOM_CONSTRAINTS"], dynamicConstraints=params["DYNAMIC_CONSTRAINTS"], numPairsToAdd = params["NUM_OD_ADDED"], numPairsToRemove = params["NUM_OD_REMOVED"], addOnlyIfAbovePercentile=params["ADD_IF_ABOVE_PERCENTILE"], removeOnlyIfBelowPercentile=params["REMOVE_IF_BELOW_PERCENTILE"], globalConstraintUpdate=params["GLOBAL_CONSTRAINT_UPDATE"],iterationMultiple = params["UPDATE_EVERY_N_ITERATIONS"], warmstart=params["WARMSTART"], metropolis=params["METROPOLIS"], localized_testing=params["LOCALIZED_TESTING"])
elseif params["ERROR_COMPUTATION"] == "both"
	@time status, new_times, directory = fast_LP(manhattan, travel_times, num_rides, trainTestDf, startTimes, testingData2 = testing_travel_times, errorComputation=params["ERROR_COMPUTATION"], model_type=params["MODEL"], last_smooth=params["LAST_SMOOTH"], max_rounds=params["MAX_ROUNDS"], min_rides=params["MIN_RIDES"], radius=params["RADIUS"], year=params["YEAR"], startTime=params["START_TIME"], endTime=params["END_TIME"], startMonth=params["START_MONTH"], endMonth=params["END_MONTH"], method=params["METHOD"], sample_size=params["SAMPLE_SIZE"], turnCost=params["TURN_COST"], turnCostAsVariable=params["TURN_COST_AS_VARIABLE"], delta_bound=params["DELTA_BOUND"], maxNumPathsPerOD=params["MAX_NUM_PATHS_PER_OD"], computeFinalSP = params["COMPUTE_FINAL_SHORTEST_PATHS"], startWithSimpleLP=params["START_SIMPLE"], randomConstraints=params["RANDOM_CONSTRAINTS"], dynamicConstraints=params["DYNAMIC_CONSTRAINTS"], numPairsToAdd = params["NUM_OD_ADDED"], numPairsToRemove = params["NUM_OD_REMOVED"], addOnlyIfAbovePercentile=params["ADD_IF_ABOVE_PERCENTILE"], removeOnlyIfBelowPercentile=params["REMOVE_IF_BELOW_PERCENTILE"], globalConstraintUpdate=params["GLOBAL_CONSTRAINT_UPDATE"],iterationMultiple = params["UPDATE_EVERY_N_ITERATIONS"], warmstart=params["WARMSTART"], metropolis=params["METROPOLIS"], localized_testing=params["LOCALIZED_TESTING"])
end
# Move log files
f=open(string(directory, "timestats.csv"), "a")
write(f, string(time()-t, '\n'))
close(f)
mv(JSONfile, string(directory, "parameters.json"), remove_destination=true)

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
# Copyright (C) 2020 Argonne National Laboratory
# Written by Alinson Santos Xavier <axavier@anl.gov>

using RELOG

@testset "Resolve" begin
    # Shoud not crash
    filename = "$(pwd())/../instances/s1.json"
    solution_old, model_old = RELOG.solve(filename, return_model = true)
    solution_new = RELOG.resolve(model_old, filename)
end

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
#   Translation of Matlab version by John Myles White
#   Translation of minpack subroutine cvsrch
#   Dianne O'Leary   July 1991
#
#     **********
#
#     Subroutine cvsrch
#
#     The purpose of cvsrch is to find a step which satisfies
#     a sufficient decrease condition and a curvature condition.
#     The user must provide a subroutine which calculates the
#     function and the gradient.
#
#     At each stage the subroutine updates an interval of
#     uncertainty with endpoints stx and sty. The interval of
#     uncertainty is initially chosen so that it contains a
#     minimizer of the modified function
#
#          f(x + stp * s) - f(x) - f_tol * stp * (gradf(x)' * s).
#
#     If a step is obtained for which the modified function
#     has a nonpositive function value and nonnegative derivative,
#     then the interval of uncertainty is chosen so that it
#     contains a minimizer of f(x + stp * s).
#
#     The algorithm is designed to find a step which satisfies
#     the sufficient decrease condition
#
#           f(x + stp * s) <= f(x) + f_tol * stp * (gradf(x)' * s),
#
#     and the curvature condition
#
#           abs(gradf(x + stp * s)' * s)) <= gtol * abs(gradf(x)' * s).
#
#     If f_tol is less than gtol and if, for example, the function
#     is bounded below, then there is always a step which satisfies
#     both conditions. If no step can be found which satisfies both
#     conditions, then the algorithm usually stops when rounding
#     errors prevent further progress. In this case stp only
#     satisfies the sufficient decrease condition.
#
#     The subroutine statement is
#
#        subroutine cvsrch(df,n,x,f,g,s,stp,f_tol,gtol,x_tol,
#                          stpmin,stpmax,maxfev,info,nfev,wa)
#
#     where
#
#	df is the name of the user-supplied subroutine which
#    calculates the function and the gradient.  df must
#    be declared in an external statement in the user
#    calling program, and should be written as follows.
#
#    function [f,g] = df(n,x) (Matlab)
#                     (10/2010 change in documentation)
#                     (derived from Fortran subroutine df(n,x,f,g))
#    integer n
#    f
#    x(n),g(n)
#
#    Calculate the function at x and
#    return this value in the variable f.
#    Calculate the gradient at x and
#    return this vector in g.
#
#  n is a positive integer input variable set to the number
#	  of variables.
#
#	x is an array of length n. On input it must contain the
#	  base point for the line search. On output it contains
#    x + stp * s.
#
#	f is a variable. On input it must contain the value of f
#    at x. On output it contains the value of f at x + stp * s.
#
#	g is an array of length n. On input it must contain the
#    gradient of f at x. On output it contains the gradient
#    of f at x + stp * s.
#
#	s is an input array of length n which specifies the
#    search direction.
#
#	stp is a nonnegative variable. On input stp contains an
#    initial estimate of a satisfactory step. On output
#    stp contains the final estimate.
#
#  f_tol and gtol are nonnegative input variables. Termination
#    occurs when the sufficient decrease condition and the
#    directional derivative condition are satisfied.
#
#	x_tol is a nonnegative input variable. Termination occurs
#    when the relative width of the interval of uncertainty
#	  is at most x_tol.
#
#	stpmin and stpmax are nonnegative input variables which
#	  specify lower and upper bounds for the step.
#
#	maxfev is a positive integer input variable. Termination
#    occurs when the number of calls to df is at least
#    maxfev by the end of an iteration.
#
#	info is an integer output variable set as follows:
#
#	  info = 0  Improper input parameters.
#
#	  info = 1  The sufficient decrease condition and the
#              directional derivative condition hold.
#
#	  info = 2  Relative width of the interval of uncertainty
#		         is at most x_tol.
#
#	  info = 3  Number of calls to df has reached maxfev.
#
#	  info = 4  The step is at the lower bound stpmin.
#
#	  info = 5  The step is at the upper bound stpmax.
#
#	  info = 6  Rounding errors prevent further progress.
#              There may not be a step which satisfies the
#              sufficient decrease and curvature conditions.
#              Tolerances may be too small.
#
#    nfev is an integer output variable set to the number of
#         calls to df.
#
#     Argonne National Laboratory. MINPACK Project. June 1983
#     Jorge J. More', David J. Thuente
#
#     **********

# Returns x, f, g, stp, info, nfev
# TODO: Decide whether to update x, f, g and info
#       or just return step and nfev and let existing code do its job

function mt_linesearch!{T}(df,
                         x::Vector,
                         s::Vector,
                         new_x::Vector,
                         g::Vector,
                         lsr::LineSearchResults{T},
                         c::Real,
                         mayterminate::Bool;
                         n::Integer = length(x),
                         stp::Real = 1.0,
                         f_tol::Real = 1e-4,
                         gtol::Real = 0.9,
                         x_tol::Real = 1e-8,
                         stpmin::Real = 1e-16,
                         stpmax::Real = 65536.0,
                         maxfev::Integer = 100)

   info = 0
   info_cstep = 1 # Info from step

   # Count function and gradient calls
   f_calls = 0
   g_calls = 0

   f = df.fg!(x, g)
   f_calls += 1
   g_calls += 1

   #
   # Check the input parameters for errors.
   #

   if n <= 0 || stp <= 0.0 || f_tol < 0.0 || gtol < 0.0 ||
      x_tol < 0.0 || stpmin < 0.0 || stpmax < stpmin || maxfev <= 0
      throw(ArgumentError("Invalid parameters to mure_thuente_line_search"))
   end

   #
   # Compute the initial gradient in the search direction
   # and check that s is a descent direction.
   #

   dginit = vecdot(g, s)
   if dginit >= 0.0
      throw(ArgumentError("Search direction is not a direction of descent"))
   end

   #
   # Initialize local variables.
   #

   bracketed = false
   stage1 = true
   nfev = 0
   finit = f
   dgtest = f_tol * dginit
   width = stpmax - stpmin
   width1 = 2 * width
   copy!(new_x, x)
   # Keep this across calls
   # Replace with new_x?

   #
   # The variables stx, fx, dgx contain the values of the step,
   # function, and directional derivative at the best step.
   # The variables sty, fy, dgy contain the value of the step,
   # function, and derivative at the other endpoint of
   # the interval of uncertainty.
   # The variables stp, f, dg contain the values of the step,
   # function, and derivative at the current step.
   #

   stx = 0.0
   fx = finit
   dgx = dginit
   sty = 0.0
   fy = finit
   dgy = dginit

   while true
      #
      # Set the minimum and maximum steps to correspond
      # to the present interval of uncertainty.
      #

      if bracketed
         stmin = min(stx, sty)
         stmax = max(stx, sty)
      else
         stmin = stx
         stmax = stp + 4 * (stp - stx) # Why 4?
      end

      #
      # Force the step to be within the bounds stpmax and stpmin
      #

      stp = max(stp, stpmin)
      stp = min(stp, stpmax)

      #
      # If an unusual termination is to occur then let
      # stp be the lowest point obtained so far.
      #

      if (bracketed && (stp <= stmin || stp >= stmax)) ||
           nfev >= maxfev-1 || info_cstep == 0 ||
           (bracketed && stmax - stmin <= x_tol * stmax)
         stp = stx
      end

      #
      # Evaluate the function and gradient at stp
      # and compute the directional derivative.
      #

      for i in 1:n
         new_x[i] = x[i] + stp * s[i] # TODO: Use x_new here
      end
      f = df.fg!(new_x, g)
      f_calls += 1
      g_calls += 1
      nfev += 1 # This includes calls to f() and g!()
      dg = vecdot(g, s)
      ftest1 = finit + stp * dgtest

      #
      # Test for convergence.
      #

      # What's does info_cstep stand for?

      if (bracketed && (stp <= stmin || stp >= stmax)) || info_cstep == 0
         info = 6
      end
      if stp == stpmax && f <= ftest1 && dg <= dgtest
         info = 5
      end
      if stp == stpmin && (f > ftest1 || dg >= dgtest)
         info = 4
      end
      if nfev >= maxfev
         info = 3
      end
      if bracketed && stmax - stmin <= x_tol * stmax
         info = 2
      end
      if f <= ftest1 && abs(dg) <= -gtol * dginit
         info = 1
      end

      #
      # Check for termination.
      #

      if info != 0
         return stp, f_calls, g_calls
      end

      #
      # In the first stage we seek a step for which the modified
      # function has a nonpositive value and nonnegative derivative.
      #

      if stage1 && f <= ftest1 && dg >= min(f_tol, gtol) * dginit
         stage1 = false
      end

      #
      # A modified function is used to predict the step only if
      # we have not obtained a step for which the modified
      # function has a nonpositive function value and nonnegative
      # derivative, and if a lower function value has been
      # obtained but the decrease is not sufficient.
      #

      if stage1 && f <= fx && f > ftest1
         #
         # Define the modified function and derivative values.
         #
         fm = f - stp * dgtest
         fxm = fx - stx * dgtest
         fym = fy - sty * dgtest
         dgm = dg - dgtest
         dgxm = dgx - dgtest
         dgym = dgy - dgtest
         #
         # Call cstep to update the interval of uncertainty
         # and to compute the new step.
         #
         stx, fxm, dgxm,
         sty, fym, dgym,
         stp, fm, dgm,
         bracketed, info_cstep =
           cstep(stx, fxm, dgxm, sty, fym, dgym,
                 stp, fm, dgm, bracketed, stmin, stmax)
         #
         # Reset the function and gradient values for f.
         #
         fx = fxm + stx * dgtest
         fy = fym + sty * dgtest
         dgx = dgxm + dgtest
         dgy = dgym + dgtest
      else
         #
         # Call cstep to update the interval of uncertainty
         # and to compute the new step.
         #
         stx, fx, dgx,
         sty, fy, dgy,
         stp, f, dg,
         bracketed, info_cstep =
           cstep(stx, fx, dgx, sty, fy, dgy,
                 stp, f, dg, bracketed, stmin, stmax)
      end

      #
      # Force a sufficient decrease in the size of the
      # interval of uncertainty.
      #

      if bracketed
         if abs(sty - stx) >= 0.66 * width1
            stp = stx + 0.5 * (sty - stx)
         end
         width1 = width
         width = abs(sty - stx)
      end
   end # while
end # function


text:
using GameZero
using Test

@testset "basic" begin
    global g
    @test_nowarn begin
        g = GameZero.initgame(joinpath("..","example","BasicGame","basic.jl"))
        GameZero.quitSDL(g)
    end

end

@testset "basic2" begin
    @test_nowarn begin
        g = GameZero.initgame(joinpath("..","example","BasicGame","basic2.jl"))
        GameZero.quitSDL(g)
    end
end


text:
using OptimNewton
using Test

@testset "newton" begin include("newtont.jl")
    # Write your own tests here.
end


text:
@testset "Utils" begin
  @testset "sumar_columna" begin
    df = cargar_csv("./src/cities.csv")
    @test Reporstat.sumar_columna(df, 1) == 4969
    @test Reporstat.sumar_columna(df, "LatD") == 4969
  end;

  @testset "sumar_fila" begin
    df = DataFrame(A = 1:4, B = 4.0:-1.0:1.0)
    @test Reporstat.sumar_fila(df, 1) == 5.0
  end;

  @testset "cargar_csv url" begin
    datos_url = "https://raw.githubusercontent.com/mucinoab/Reporstat/master/src/cities.csv"
    dlocal = cargar_csv("./src/cities.csv")
    dremoto = cargar_csv(datos_url, "URL")
    @test dlocal == dremoto
  end;

  @testset "cargar_csv local" begin
    file = "./src/cities.csv"
    df_custom = cargar_csv(file)
    df = DataFrame(CSV.File(file))
    @test df_custom == df
  end;

  @testset "API INEGI" begin
    @test Reporstat.poblacion_mexico().lugar == ["México"]

    @test Reporstat.poblacion_entidad("32").lugar == ["Zacatecas"]
    @test Reporstat.poblacion_entidad("06").densidad_poblacion == [129.981519801061]

    @test Reporstat.poblacion_municipio("01", "001").lugar == ["Aguascalientes, Aguascalientes"]

    cdmx = Reporstat.poblacion_municipio("09", "016")
    @test cdmx.lugar == ["Ciudad de México, Miguel Hidalgo"]
    @test cdmx.hombres == [195467]
    @test cdmx.mujeres == [219003]
    @test cdmx.porcentaje_hombres == [47.1607112698145]
    @test cdmx.porcentaje_mujeres == [52.8392887301855]
    @test cdmx.porcentaje_indigena == [5.0142822]
    @test cdmx.densidad_poblacion == [8927.81062496375]
  end;
end;

max_stars_repo_name: yollotltamayo/Covid

max_issues_repo_name: yollotltamayo/Covid

max_forks_repo_name: yollotltamayo/Covid


text:
module Test
export generate_dot, generate_sum

using LinearAlgebra, Random
using ..EFT: two_prod

"""
    (x, y, d, c) = generate_dot(n, c)

Generate two Float64 vectors whose dot product is ill-conditioned.

Inputs:
  n -- vectors size
  c -- target condition number

Results:
  x, y -- vectors of size n
  d    -- accurate dot product, rounded to nearest
  c    -- actual condition number of the dot product
"""
function generate_dot(n, c::T) where {T}
    R = Rational{BigInt}

    # Initialization
    x = zeros(T, n)
    y = zeros(T, n)

    # First half of the vectors:
    #   random numbers within a large exponent range
    n2 = div(n, 2)
    b = log2(c)

    e = rand(n2) .* b/2
    e[1]  = b/2 + 1           # Make sure exponents b/2
    e[n2] = 0                 # and 0 actually occur
    for i in 1:n2
        x[i] = (2*rand(T)-1) * 2^(e[i])
        y[i] = (2*rand(T)-1) * 2^(e[i])
    end

    # Second half of the vectors such that
    #   (*) log2( dot (x[1:i], y[1:i]) ) decreases from b/2 to 0
    δe = -b/(2*(n-n2-1))
    e = b/2:δe:0
    for i in eachindex(e)
        # Random x[i]
        cx = (2*rand()-1) * 2^(e[i])
        x[i+n2] = cx

        # y[i] chosen according to (*)
        cy = (2*rand()-1) * 2^(e[i])
        y[i+n2] = (cy - T(dot(R.(x), R.(y)))) / cx
    end

    # Random permutation of x and y
    perm = randperm(n)
    X = x[perm]
    Y = y[perm]

    # Dot product, rounded to nearest
    d = T(dot(R.(X), R.(Y)))

    # Actual condition number
    c = 2 * dot(abs.(X), abs.(Y)) / abs(d)

    (X,Y,d,c)
end

"""
    (x, s, c) = generate_sum(n, c)

Generate a Float64 vectors whose sum is ill-conditioned.

Inputs:
  n -- vectors size
  c -- target condition number

Results:
  x -- vectors of size n
  s -- accurate sum, rounded to nearest
  c -- actual condition number of the sum
"""
function generate_sum(n, c::T) where {T}
    R = Rational{BigInt}

    (x, y, _, _) = generate_dot(n÷2, c)

    z = (two_prod.(x, y)
         |> Iterators.flatten
         |> collect)

    # Complete if necessary
    if length(z) < n
        push!(z, rand())
    end

    z = shuffle(z)

    # Sum, rounded to nearest
    s = T(sum(R.(z)))

    # Actual condition number
    c = sum(abs.(z)) / abs(s)

    (z, s, c)
end

end


text:
module PoseComposition

import Base: @kwdef
import LinearAlgebra: dot, norm, cross
import Rotations: AngleAxis, Rotation, QuatRotation, RotZYX, RotMatrix
import StaticArrays: StaticVector, SVector, @SVector

include("docstring_extensions.jl")

"""
Struct representing the pose (position and orientation) of an object.

Can represent either a relative pose (relative to some parent coordinate frame
which must be supplied by context) or an absolute pose.  An absolute pose is,
by definition, a pose relative to the world coordinate frame.

See [Operations on Poses](@ref) for documentation of this data structure and
the operations defined on it.
"""
@kwdef struct Pose
  """Origin of the object's coordinate frame."""
  pos::StaticVector{3, <:Real}
  """Orientation of the object's coordinate frame."""
  orientation::Rotation{3}
end

### Constructors ###

function Pose(x::Real, y::Real, z::Real, orientation::Rotation{3})::Pose
  Pose(@SVector([x, y, z]), orientation)
end

function Pose(pos::AbstractVector{<:Real}, orientation::Rotation{3})::Pose
  (x, y, z) = pos
  Pose(@SVector([x, y, z]), orientation)
end

function Pose(pos::AbstractVector{<:Real}, ypr::NamedTuple{(:yaw, :pitch, :roll)})::Pose
  Pose(pos, RotZYX(ypr.yaw, ypr.pitch, ypr.roll))
end

function Pose(x::Real, y::Real, z::Real, ypr::NamedTuple{(:yaw, :pitch, :roll)})::Pose
  Pose(x, y, z, RotZYX(ypr.yaw, ypr.pitch, ypr.roll))
end

function Base.isapprox(a::Pose, b::Pose; kwargs...)::Bool
  (isapprox(a.pos, b.pos; kwargs...) &&
   isapprox(a.orientation, b.orientation; kwargs...))
end

### Constants ###

"""
The identity quaternion, representing the identity orientation.
"""
IDENTITY_ORN = one(QuatRotation)

"""
Identity pose, a.k.a. the relative pose of any coordinate frame relative to itself.

This is the identity element for the pose group (the identity element is the same in both
group structures).
"""
const IDENTITY_POSE = Pose([0, 0, 0], IDENTITY_ORN)

### Pretty-printing ###

function Base.show(io::IO, pose::Pose)
  # Not sure whether most people will want (yaw, pitch, roll) or quaternion
  # components here.
  print(io, strWithQuat(pose))
end

function strWithYPR(pose::Pose)::String
  (yaw, pitch, roll) = _yawPitchRoll(pose.orientation)
  "Pose⟨pos=$(pose.pos), orientation=(yaw=$yaw, pitch=$pitch, roll=$roll)⟩"
end

function strWithQuat(pose::Pose)::String
  q = QuatRotation(pose.orientation)
  (w, x, y, z) = componentsWXYZ(q)
  "Pose⟨pos=$(pose.pos), orientation=(w=$w, x=$x, y=$y, z=$z)⟩"
end

function _yawPitchRoll(orn::Rotation{3})
  ypr = RotZYX(orn)
  (yaw=ypr.theta1, pitch=ypr.theta2, roll=ypr.theta3)
end

componentsWXYZ(q::QuatRotation) = @SVector([q.q.s, q.q.v1, q.q.v2, q.q.v3])

"""
Like `isapprox`, but does not consider a quaternion to be equivalent to its
negative (even though they correspond to the same rotation matrix).  Note that
this is stricter than `Base.isapprox`, since for a `Rotations.QuatRotation`
`q`, we have `-q ≈ q` and in fact `-q == q`.
"""
function isapproxIncludingQuaternionSign(a::Pose, b::Pose; kwargs)::Bool
  (isapprox(a.pos, b.pos; kwargs...) &&
   isapprox(componentsWXYZ(QuatRotation(a.orientation)),
            componentsWXYZ(QuatRotation(b.orientation));
            kwargs...))
end

# Operations for combining poses.
Base.:(*)(a::Pose, b::Pose)::Pose = Pose(
    a.pos + a.orientation * b.pos,
    a.orientation * b.orientation)

Base.:(/)(a::Pose, b::Pose)::Pose = Pose(
    a.pos - (a.orientation / b.orientation) * b.pos,
    a.orientation / b.orientation)

Base.:(\)(a::Pose, b::Pose)::Pose = Pose(
    a.orientation \ (-a.pos + b.pos),
    a.orientation \ b.orientation)

Base.:(^)(pose::Pose, t::Real) = Pose(t * pose.pos,
                                      # quaternion exponentiation = SLERP
                                      quatPow(QuatRotation(pose.orientation),
                                              t))

function Base.inv(a::Pose)::Pose
  Pose(a.orientation \ -a.pos,
       inv(a.orientation))
end

# Action of a pose on a point.
function Base.:(*)(a::Pose, bpos::StaticVector{3, <:Real})
  a.pos + a.orientation * bpos
end

function Base.:(\)(a::Pose, bpos::StaticVector{3, <:Real})
  a.orientation \ (-a.pos + bpos)
end

# Convenience wrappers when the user supplies a Vector instead of a StaticVector
Base.:(*)(a::Pose, bpos::AbstractVector{<:Real}) = a * SVector{3}(bpos)
Base.:(\)(a::Pose, bpos::AbstractVector{<:Real}) = a \ SVector{3}(bpos)

"""
Vectorized pose–point multiplication.  Returns the matrix whose `i`th column is
`a * bpoints[:, i]`.

The matrix `bpoints` must have 3 rows, as each column represents a point in 3D
space.
"""
function Base.:(*)(a::Pose, bpoints::AbstractMatrix{<:Real})
  size(bpoints, 1) == 3 || error(
      "Must pass a 3×N matrix (one column per point)")
  a.pos .+ a.orientation * bpoints
end

"""
Vectorized version of pose–point left division.  Returns the matrix whose `i`th
column is `a \\ bpoints[:, i]`.

The matrix `bpoints` must have 3 rows, as each column represents a point in 3D
space.
"""
function Base.:(\)(a::Pose, bpoints::AbstractMatrix{<:Real})
  size(bpoints, 1) == 3 || error(
      "Must pass a 3×N matrix (one column per point)")
  a.orientation \ (-a.pos .+ bpoints)
end

"""
!!! note "TODO"
    This code mostly duplicates `GenDirectionalStats.hopf`.  The two should
    probably be consolidated into one.

Returns a rotation that carries the z-axis to `newZ`, with the remaining degree
of freedom determined by `planarAngle` as described below.

Start with the case `planarAngle = 0`.  In that case, the returned rotation is
the unique (except at singularities) rotation that carries `[0, 0, 1]` to
`newZ` "along a great circle" (more precisely: the unique rotation that carries
`[0, 0, 1]` to `newZ` and whose equator contains `[0, 0, 1]` and `newZ`;
"equator" means the unique great circle that is fixed setwise by the rotation).

Next consider the general case.  This works the same as the above special case,
except that we precede that rotation with a rotation by angle `planarAngle`
around `[0, 0, 1]` (or equivalently, we follow that rotation with a rotation by
`planarAngle` around `newZ`).

The name of this function comes from the fact that we are using geodesics
(great circles) to define a coordinate chart on the fiber over `newZ`  in the
Hopf fibration.

See also: [`invGeodesicHopf`](@ref)
"""
function geodesicHopf(newZ::StaticVector{3, <:Real}, planarAngle::Real)
  @assert norm(newZ) ≈ 1
  zUnit = @SVector([0, 0, 1])
  if newZ ≈ -zUnit
    @warn "Singularity: anti-parallel z-axis, rotation has an undetermined degree of freedom"
    axis = @SVector([1, 0, 0])
    geodesicAngle = π
  elseif newZ ≈ zUnit
    # Choice of axis doesn't matter here as long as it's nonzero
    axis = @SVector([1, 0, 0])
    geodesicAngle = 0
  else
    axis = cross(zUnit, newZ)
    @assert !(axis ≈ zero(axis)) || newZ ≈ zUnit
    geodesicAngle = let θ = asin(clamp(norm(axis), -1, 1))
      dot(zUnit, newZ) > 0 ? θ : π - θ
    end
  end
  return (AngleAxis(geodesicAngle, axis...) *
          AngleAxis(planarAngle, zUnit...))
end

geodesicHopf(newZ::AbstractVector{<:Real}, planarAngle::Real) = geodesicHopf(
    SVector{3}(newZ), planarAngle)

"""
Inverse function of [`geodesicHopf`](@ref).

Satisfies the round-trip conditions

    geodesicHopf(invGeodesicHopf(r)...) == r

and

    invGeodesicHopf(geodesicHopf(newZ, planarAngle))
    == (newZ=newZ, planarAngle=planarAngle)
"""
function invGeodesicHopf(r::Rotation{3})::NamedTuple{(:newZ, :planarAngle)}
  zUnit = @SVector([0, 0, 1])
  newZ = r * zUnit
  if newZ ≈ -zUnit
    @warn "Singularity: anti-parallel z-axis, planarAngle is undetermined"
    planarAngle = 0
  else
    # Solve `planarRot == AngleAxis(planarAngle, zUnit...)` for `planarAngle`
    planarRot = AngleAxis(geodesicHopf(newZ, 0) \ r)
    axis = @SVector([planarRot.axis_x, planarRot.axis_y, planarRot.axis_z])
    # `axis` is either `zUnit` or `-zUnit`, and we need to ensure that it's
    # `zUnit`.  (Exception: the degenerate case `planarAngle == 0`)
    if axis[3] < 0
      axis = -axis
      planarAngle = -planarRot.theta
    else
      planarAngle = planarRot.theta
    end
    atol = 1e-14
    @assert isapprox(axis, zUnit; atol=atol) ||
            abs(rem2pi(planarAngle, RoundNearest)) < atol
  end
  return (newZ=newZ, planarAngle=planarAngle)
end

# Second group structure on poses: Direct product ``ℝ^3 × SO_3(ℝ)``.
⊗(a::Pose, b::Pose) = Pose(a.pos + b.pos, a.orientation * b.orientation)
⊘(a::Pose, b::Pose) = Pose(a.pos - b.pos, a.orientation / b.orientation)
⦸(a::Pose, b::Pose) = Pose(b.pos - a.pos, a.orientation \ b.orientation)

"""
Interpolates between the identity pose and `b`.

Namely, `interp(b, 0) == IDENTITY_POSE` and `interp(b, 1) == b`.  The position
is interpolated linearly and the orientation is interpolated by quaternion
SLERP.  That is, this interpolation treats position and orientation
independently, as in the [`⊗`](@ref) operation (not the [`*`](@ref
Base.:*(::Pose, ::Pose)) operation).
"""
interp(b::Pose, t::Real) = Pose(t * b.pos,
                                quatPow(QuatRotation(b.orientation), t))

"""
Like [`interp`](@ref interp(::Pose, ::Real)), but interpolates between two
given poses rather than always starting at the identity.  That is,

    interp(a, b, 0) == a
    interp(a, b, 1) == b

and as a special case, we have

    interp(b, t) == interp(IDENTITY_POSE, b, t)
"""
interp(a::Pose, b::Pose, t::Real) = a * interp(a ⦸ b, t)

function quatPow(q::QuatRotation, t::Real)
  # TODO: Once https://github.com/JuliaGeometry/Rotations.jl/issues/126 is
  # fixed, this special case won't be necessary
  if t == 0 || q == one(QuatRotation) || q == -one(QuatRotation)
    return one(QuatRotation)
  end
  return RotMatrix{3}(exp(t * log(q)))
end

end  # module PoseComposition


text:
@testset "Test Scalar Function Spaces" begin

using HDiscontinuousGalerkin
using Tensors

root_file = "mesh/figure2.1"
mesh = parse_mesh_triangle(root_file)

finiteElement = ContinuousLagrange{2,RefTetrahedron,1}()
Wh = ScalarFunctionSpace(mesh, finiteElement)

# Basic Test
sq2 = sqrt(2)
sq3 = sqrt(3)
invs = (Tensor{2,2}([1.0 1.0;-2.0 0.0]),
        Tensor{2,2}([-1.0 -1.0;2.0 0.0]),
        Tensor{2,2}([1.0 1.0;-1.0 1.0]),
        Tensor{2,2}([1.0 -1.0;0.0 2.0]))
detsJf = ([sq2/2,sq2/2,1],
          [sq2/2,sq2/2,1],
          [1,sq2/2,sq2/2],
          [sq2/2,sq2/2,1])
@test getnlocaldofs(Wh) == 3
@inbounds for (cellcount, cell) in enumerate(CellIterator(mesh))
    reinit!(Wh, cell)
    @test Wh.detJ[] ≈ 0.5
    @test getdetJdV(Wh,1)/Wh.qr_weights[1] ≈ 0.5
    @test Wh.Jinv[] ≈ invs[cellcount]
    # Face Data
    @test Wh.detJf ≈ detsJf[cellcount]
    if cellcount == 1
        @test Wh.normals ≈ [Vec{2}([-sq2/2,sq2/2]), Vec{2}([-sq2/2,-sq2/2]), Vec{2}([1.0,0.0])]
    end
end
end


text:
module Inspect
using ..Salsa
using ..Salsa: AbstractSalsaStorage

"""
    Salsa.Inspect.all_inputs_and_derived_vals(st::SalsaStorage)

Should be implemented by salsa storage backends, to support generating a dependency graph
from the current state of a Salsa Runtime.

Returns two dicts (inputs: K=>V and derived_vals: K=>(V, K[deps...]))
"""
function all_inputs_and_derived_vals end

# Crude code to dump a .dot (graphviz) version of the Salsa Node tree.
# For debugging / illustration purposes - not intended to be reachable from production
# code.
function dump_graph(rt::Runtime; module_boxes=false)
    println(build_graph(rt; module_boxes=module_boxes))
end

# There is no IdSet in julia, so we build our own from IdDict.
# We need an IdSet to prevent comparing all the contents of the various maps, which is
# expensive. We just want to compare via `===` which is what an IdDict does.
# (Used the name _IdSet to indicate that this is internal only, and not a full-featured set)
struct _IdSet{T}
    d::IdDict{T,Nothing}  # Only using the Keys
    _IdSet{T}() where T = new(IdDict{T,Nothing}())
end
Base.in(k, s::_IdSet) = haskey(s.d, k)
Base.push!(s::_IdSet, k) = s.d[k] = nothing

function build_graph(rt::Salsa.Runtime; module_boxes=false)
    io = IOBuffer()
    println(io, """digraph G {""")
    println(io, """edge [dir="back"];""")
    seen = _IdSet{Any}()
    modules_map = Dict{Module,Set}()
    edges = Dict{Pair, Int}()
    inputs = Dict{Any, String}()  # Note, there might be duplicate strings, Any must be the key

    _build_graph(io, rt.storage, seen, modules_map, edges, inputs)

    if module_boxes
        for (m,keys) in modules_map
            mname = nameof(m)
            println(io, """subgraph cluster_$mname {""")
            #println(io, """node [style=filled];""")
            # First print non-input keys
            println(io, """ $(join((repr(vertex_name(k)) for k in keys
                                   if !haskey(inputs, k)), " ")) """);
            begin  # Then print input keys all on the bottom
                println(io, "{")
                println(io, "rank=sink;")
                println(io, """ $(join((repr(vertex_name(k)) for k in keys
                                    if haskey(inputs, k)), " ")) """);
                println(io, "}")
            end
            println(io, """ label = "Module `$mname`"; """)
            println(io, """ fontsize = 25; """)
            println(io, """ color=blue; """)
            println(io, """}""")
        end
    end

    if !module_boxes
        println(io, "{")
        println(io, "rank=sink;")
    end
    for (input_key, name) in inputs
        println(io, """$(vertex_name(input_key)) [label=$(repr(name))]""")
    end
    if !module_boxes
        println(io, "}")
    end

    max_count = length(edges) == 0 ? 0 : maximum(values(edges))
    maxwidth = 10
    for ((a,b), count) in edges
        normwidth = 1 + maxwidth * (count / max_count)
        println(io, """ "$(vertex_name(a))" -> "$(vertex_name(b))"  [penwidth=$normwidth, weight=$normwidth]""")
    end
    @show max_count
    println(io, "}")
    return String(take!(io))
end

function _build_graph(io::IO, st::AbstractSalsaStorage, seen::_IdSet, modules_map::Dict, out_edges::Dict, out_inputs::Dict)
    inputs, derived_vals = all_inputs_and_derived_vals(st)
    for (k,v) in inputs
        F = key_func_type(k)
        m = F.name.module
        out_inputs[k] = "@input $(nameof(F.instance))"
        push!(get!(modules_map, m, Set([])), k)
    end

    seen_key_types = Set{Any}()
    for (k,(v, deps)) in derived_vals
        K = key_identifier(k)
        if !(K ∈ seen_key_types)
            push!(seen_key_types, K)
            key_str = _derived_key_as_call_str(k)
            println(io, "$(vertex_name(k)) [shape=rect,label=$(repr(key_str))]")
            F = key_func_type(k)
            m = F.name.module
            push!(get!(modules_map, m, Set([])), k)
        end
        _build_graph(io, st, k, v, deps, seen, modules_map, out_edges)
    end
end

key_func_type(::Union{Salsa.InputKey{F},Salsa.DerivedKey{F}}) where F = F
# TODO: Change this to use key values, not just types, for unrolling the graph.
key_identifier(k::Union{Salsa.InputKey,Salsa.DerivedKey}) = typeof(k)

# TODO: Change these to use key values, not just types, for unrolling the graph.
function vertex_name(k::Salsa.InputKey)
    return vertex_name(typeof(k))
end
function vertex_name(k::Salsa.DerivedKey)
    return vertex_name(typeof(k))
end
function vertex_name(x::Any)::String
    return "v$(objectid(x))"
end

function _build_graph(io, st::AbstractSalsaStorage,
        derived_key, v, deps::Vector,  # Vector of salsa keys
        seen::_IdSet{Any}, modules_map::Dict{Module,Set}, edges::Dict{Pair, Int}) where {F,TT}
    in(derived_key, seen) && return
    push!(seen, derived_key)
    for d in deps
        edge = key_identifier(derived_key) => key_identifier(d)
        count = get!(edges, edge, 0) + 1
        edges[edge] = count
    end
    #_build_graph(io, s.leaf_node, seen)
end

function _derived_key_as_call_str(key::Salsa.DerivedKey{F,TT})::String where {F,TT}
    f = isdefined(F, :instance) ? nameof(F.instance) : nameof(F)
    argsexprs = [
        :rt,
        (Expr(:(::), fieldtype(TT, i)) for i = 1:fieldcount(TT))...,
    ]
    # Display f in a copy/paste executable way, by exploiting Expr printing.
    f_str = string(:(($f,)))[2:end-2]  # (wraps f name in var"" if needed)
    return "$f_str($(join(argsexprs, ", ")))"
end

end # module


text:
# ---
# title: 1041. Robot Bounded In Circle
# id: problem1041
# author: Indigo
# date: 2021-07-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/robot-bounded-in-circle/description/>
# hidden: true
# ---
#
# On an infinite plane, a robot initially stands at `(0, 0)` and faces north.
# The robot can receive one of three instructions:
#
#   * `"G"`: go straight 1 unit;
#   * `"L"`: turn 90 degrees to the left;
#   * `"R"`: turn 90 degress to the right.
#
# The robot performs the `instructions` given in order, and repeats them
# forever.
#
# Return `true` if and only if there exists a circle in the plane such that the
# robot never leaves the circle.
#
#
#
# **Example 1:**
#
#
#
#     Input: "GGLLGG"
#     Output: true
#     Explanation:
#     The robot moves from (0,0) to (0,2), turns 180 degrees, and then returns to (0,0).
#     When repeating these instructions, the robot remains in the circle of radius 2 centered at the origin.
#
#
# **Example 2:**
#
#
#
#     Input: "GG"
#     Output: false
#     Explanation:
#     The robot moves north indefinitely.
#
#
# **Example 3:**
#
#
#
#     Input: "GL"
#     Output: true
#     Explanation:
#     The robot moves from (0, 0) -> (0, 1) -> (-1, 1) -> (-1, 0) -> (0, 0) -> ...
#
#
#
#
# **Note:**
#
#   1. `1 <= instructions.length <= 100`
#   2. `instructions[i]` is in `{'G', 'L', 'R'}`
#
#
## @lc code=start
using LeetCode

function is_robot_bounded(instructions::String)::Bool
    x, y = 0, 0
    dirs = ((0, 1), (1, 0), (0, -1), (-1, 0))
    st = 1
    for ch in instructions
        if ch == 'G'
            x += dirs[st][1]
            y += dirs[st][2]
        elseif ch == 'L'
            st = mod1(st + 1, 4)
        else
            st = mod1(st - 1, 4)
        end
    end
    return (x, y) == (0, 0) || st != 1
end
## @lc code=end


text:
include("./dream6_m1/dream6_m1_@macro.jl")

# using JSON

using NLopt, LikelihoodProfiler

# Scale types
mixed_scale = fill(true, length(p0))
mixed_scale[[15,17,19,21,23,25,27,29]] = false

# Bounds
bounds_params = fill([1e-12, 1e12], length(p0))
bounds_params[[15,17,19,21,23,25,27,29]] = fill([0.5, 9.], 8)

# fit optimum
opt_0 = Opt(
    :LN_NELDERMEAD, # :LN_SBPLX
    length(p0)
)

min_objective!(opt_0, obj)
ftol_abs!(opt_0, 1e-3)
# bound constraints
lb = minimum.(bounds_params)
ub = maximum.(bounds_params)
lower_bounds!(opt_0, lb)
upper_bounds!(opt_0, ub)
(obj_0, p_0, ret_0) = optimize(opt_0, p0)

# calculate critical level
obj_crit = obj_0 + chisqinvcdf(1, 0.95)

# scan bounds
scan_bounds = fill([1e-6, 1e6], length(p0))
scan_bounds[[15,17,19,21,23,25,27,29]] = fill([0.6, 8.], 8)

# calculate intervals for all
max_iter = 100000
result = []
for i = 1:length(p0)
    println()
    println("START # $i")
    out_inq = params_intervals(
        p_0,
        i,
        obj_crit,
        obj,
        logscale = mixed_scale,
        bounds = bounds_params,
        scan_bound = scan_bounds[i],
        max_iter = max_iter
    )
    println(
        "FINISH #$i, iteration count: ",
        round(out_inq.count_evals[1]/max_iter*100, 1), "%   ",
        round(out_inq.count_evals[2]/max_iter*100, 1), "%"
    )
    push!(result, out_inq)
end

result

open("result.json", "w") do f
    JSON.print(f, result, 4)
end

open("result.csv", "w") do f
    CSV.print(f, result)
end

[println(r.loss_final) for r in result]


text:
# Autogenerated wrapper script for libaio_jll for armv7l-linux-musleabihf
export libaio

JLLWrappers.@generate_wrapper_header("libaio")
JLLWrappers.@declare_library_product(libaio, "libaio.so.1")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libaio,
        "lib/libaio.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

__precompile__(true)

module Unicode

export graphemes

"""
    Unicode.normalize(s::AbstractString, normalform::Symbol)

Normalize the string `s` according to one of the four "normal forms" of the Unicode
standard: `normalform` can be `:NFC`, `:NFD`, `:NFKC`, or `:NFKD`.  Normal forms C
(canonical composition) and D (canonical decomposition) convert different visually identical
representations of the same abstract string into a single canonical form, with form C being
more compact.  Normal forms KC and KD additionally canonicalize "compatibility equivalents":
they convert characters that are abstractly similar but visually distinct into a single
canonical choice (e.g. they expand ligatures into the individual characters), with form KC
being more compact.

Alternatively, finer control and additional transformations may be be obtained by calling
`Unicode.normalize(s; keywords...)`, where any number of the following boolean keywords
options (which all default to `false` except for `compose`) are specified:

* `compose=false`: do not perform canonical composition
* `decompose=true`: do canonical decomposition instead of canonical composition
  (`compose=true` is ignored if present)
* `compat=true`: compatibility equivalents are canonicalized
* `casefold=true`: perform Unicode case folding, e.g. for case-insensitive string comparison
* `newline2lf=true`, `newline2ls=true`, or `newline2ps=true`: convert various newline
  sequences (LF, CRLF, CR, NEL) into a linefeed (LF), line-separation (LS), or
  paragraph-separation (PS) character, respectively
* `stripmark=true`: strip diacritical marks (e.g. accents)
* `stripignore=true`: strip Unicode's "default ignorable" characters (e.g. the soft hyphen
  or the left-to-right marker)
* `stripcc=true`: strip control characters; horizontal tabs and form feeds are converted to
  spaces; newlines are also converted to spaces unless a newline-conversion flag was
  specified
* `rejectna=true`: throw an error if unassigned code points are found
* `stable=true`: enforce Unicode Versioning Stability

For example, NFKC corresponds to the options `compose=true, compat=true, stable=true`.

# Examples
```jldoctest
julia> "μ" == Unicode.normalize("µ", compat=true) #LHS: Unicode U+03bc, RHS: Unicode U+00b5
true

julia> Unicode.normalize("JuLiA", casefold=true)
"julia"

julia> Unicode.normalize("JúLiA", stripmark=true)
"JuLiA"
```
"""
function normalize end
normalize(s::AbstractString, nf::Symbol) = Base.Unicode.normalize(s, nf)
normalize(s::AbstractString; kwargs...) = Base.Unicode.normalize(s; kwargs...)

"""
    Unicode.isassigned(c) -> Bool

Returns `true` if the given char or integer is an assigned Unicode code point.

# Examples
```jldoctest
julia> Unicode.isassigned(101)
true

julia> Unicode.isassigned('\\x01')
true
```
"""
isassigned(c) = Base.Unicode.isassigned(c)

"""
    graphemes(s::AbstractString) -> GraphemeIterator

Returns an iterator over substrings of `s` that correspond to the extended graphemes in the
string, as defined by Unicode UAX #29. (Roughly, these are what users would perceive as
single characters, even though they may contain more than one codepoint; for example a
letter combined with an accent mark is a single grapheme.)
"""
graphemes(s::AbstractString) = Base.Unicode.GraphemeIterator{typeof(s)}(s)

# BEGIN 0.7 deprecations

@deprecate is_assigned_char(c::Char) Unicode.isassigned(c)
@deprecate normalize_string(s::AbstractString, nf::Symbol; kwargs...) Unicode.normalize(s, nf; kwargs...)
@deprecate normalize_string(s::AbstractString; kwargs...) Unicode.normalize(s; kwargs...)

# END 0.7 deprecations

end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
function vermeer_simple_shear_update_du(u,p,t)
    # unpacking
    r, du_prev, params = p
    σ̇₁₁_prev, σ̇₁₂_prev, σ̇ₒₒₚ_prev, ϵ̇₂₂_prev, _ = du_prev
    σ₁₁, σ₁₂, σₒₒₚ, ϵ₂₂, D = u
    # get stress tensor
    σᵢⱼ = SymmetricTensor{2,3}(SA[σ₁₁ σ₁₂ 0 ; σ₁₂ params.σ₃ 0 ; 0 0 σₒₒₚ])
    # compute damage growth rate
    KI = compute_KI(r,σᵢⱼ,D)
    Ḋ = compute_subcrit_damage_rate(r,KI,D ; vmax = params.dldt_max)

    # solve for non linear derivatives
    nl_u = SA[σ̇₁₁_prev, σ̇₁₂_prev, σ̇ₒₒₚ_prev, ϵ̇₂₂_prev] # first guess is previous iter
    nl_u = nl_solve(residual_simple_shear_old,nl_u,u,p,σᵢⱼ,Ḋ)
    # (du,u,σᵢⱼ,Ḋ,p)->residual_simple_shear(du,u,σᵢⱼ,Ḋ,p ; R1_allowed=params.regime1)
    # build back unknown vector
    du = SA[nl_u[1], nl_u[2], nl_u[3], nl_u[4], Ḋ]
    #println("du after update : ", du)
    p[2]=du

    return du
end

function nl_solve(res_func::F,nl_u,u,p,σᵢⱼ,Ḋ) where F<:Function
    # unpack
    r, du_prev, params = p
    # scale strain term
    nl_u_scaled = setindex(nl_u,nl_u[4]*r.G, 4)

    #result = DiffResults.JacobianResult(nl_u) no need for DiffResults with static vectors
    for i in 1:params.newton_maxiter
        # get residual and its gradient with respect to u
        ∇res = ForwardDiff.jacobian(nl_u -> res_func(nl_u,u,σᵢⱼ,Ḋ,p), nl_u_scaled)
        res  = res_func(nl_u_scaled,u,σᵢⱼ,Ḋ,p)

        # update u with Newton algo
        nl_u_scaled = nl_u_scaled - ∇res\res

        if (norm(res) <= params.newton_abstol) && (abs(res[3]) > 1e-10) #test plain strain condition
            @show(norm(res),res[3])
        end
        (norm(res) <= params.newton_abstol) && (break) # @debug("Newton iter $i ending norm res = $(norm(res))") ;

        (i == params.newton_maxiter) && @debug("Newton maxiter ending norm res = $(norm(res))")#("max newton iteration reached ($i), residual still higher than abstol with $(norm(res))")
    end
    # rescale strain term :
    nl_u = setindex(nl_u_scaled,nl_u_scaled[4]/r.G, 4)
    return nl_u
end

function residual_simple_shear(du,u,σᵢⱼ,Ḋ,p ; R1_allowed=false)
    # unpacking
    r, du_prev, params = p
    σ̇₁₁, σ̇₁₂, σ̇ₒₒₚ, Gϵ̇₂₂ = du
    _, _, _, _, D = u
    ϵ̇₂₂ = Gϵ̇₂₂/r.G # rescale ϵ̇₂₂

    # build σ̇ and ϵ̇ from du
    #zeroTdu = zero(eltype(du))
    σ̇ᵢⱼ = SymmetricTensor{2,3}([ σ̇₁₁    σ̇₁₂   0.0  ;
                                 σ̇₁₂    0.0   0.0  ;
                                 0.0    0.0   σ̇ₒₒₚ ])
    σ = tr(σᵢⱼ)/3
    sᵢⱼ = dev(σᵢⱼ)
    τ = get_τ(sᵢⱼ)

    σ̇ = tr(σ̇ᵢⱼ)/3
    τ̇ = sᵢⱼ ⊡ σ̇ᵢⱼ / (2*τ)

    # if R1_allowed && (Ḋ == 0)
    #     KI = DSB.compute_KI(r,σ,τ,D)
    #     if KI <= 0
    #         res = SA[σ̇ᵢⱼ[1,1] - (3r.ν/(1+r.ν))*σ̇,
    #                  -2*r.G*ϵ̇₂₂ + σ̇ᵢⱼ[2,2] - (3r.ν/(1+r.ν))*σ̇,
    #                  σ̇ᵢⱼ[3,3] - (3r.ν/(1+r.ν))*σ̇,
    #                  -2*r.G*params.ϵ̇₁₂ + σ̇ᵢⱼ[1,2]]
    #         return res
    #     end
    # end

    c1, c2, c3 = compute_c1c2c3(r,D)
    A, B = compute_AB(r,c1,c2,c3)
    A1, B1 = compute_A1B1(r,A,B)

    dc1dD = compute_dc1dD(r,D)
    dc2dD = compute_dc2dD(r,D)
    dc3dD = compute_dc3dD(r,D)
    dA1dD = compute_dA1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)
    dB1dD = compute_dB1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)

    la1 = λ₁(A1,B1,σ,τ)
    la2 = λ₂(r,A1,B1,σ,τ)
    la3 = λ₃(A1,B1)
    dla1dD  = dλ₁dD(A1,B1,dA1dD,dB1dD,σ,τ)
    dla2dD  = dλ₂dD(A1,B1,dA1dD,dB1dD,σ,τ)
    dla3dD  = dλ₃dD(A1,B1,dA1dD,dB1dD)
    dla1dσ  = dλ₁dσ(A1,B1,τ)
    dla1dτ  = dλ₁dτ(A1,B1,σ,τ)
    # build residual
    res = SA[ la1*σ̇₁₁ - la2*σ̇ + la3*τ̇ + Ḋ * (dla1dD*σᵢⱼ[1,1] - dla2dD*σ + dla3dD*τ) + (dla1dσ*σ̇ + dla1dτ*τ̇)*sᵢⱼ[1,1],
              -2*r.G*ϵ̇₂₂ - la2*σ̇ + la3*τ̇ + Ḋ * (dla1dD*params.σ₃ - dla2dD*σ + dla3dD*τ) + (dla1dσ*σ̇ + dla1dτ*τ̇)*sᵢⱼ[2,2],
              la1*σ̇ₒₒₚ - la2*σ̇ + la3*τ̇ + Ḋ * (dla1dD*σᵢⱼ[3,3] - dla2dD*σ + dla3dD*τ) + (dla1dσ*σ̇ + dla1dτ*τ̇)*sᵢⱼ[3,3],
              -2*r.G*params.ϵ̇₁₂ + la1*σ̇₁₂ + Ḋ*dla1dD*σᵢⱼ[1,2] + (dla1dσ*σ̇ + dla1dτ*τ̇)*sᵢⱼ[1,2] ]
    return res
end

function residual_simple_shear_old(du,u,σᵢⱼ,Ḋ,p ; keep_2G=false)
    # unpacking
    r, du_prev, params = p
    σ̇₁₁, σ̇₁₂, σ̇ₒₒₚ, Gϵ̇₂₂ = du
    _, _, _, _, D = u
    ϵ̇₂₂ = Gϵ̇₂₂/r.G # rescale ϵ̇₂₂

    # build σ̇ and ϵ̇ from du
    #zeroTdu = zero(eltype(du))
    σ̇ᵢⱼ = SymmetricTensor{2,3}([ σ̇₁₁    σ̇₁₂   0.0  ;
                                 σ̇₁₂    0.0   0.0  ;
                                 0.0    0.0   σ̇ₒₒₚ ])
    ϵ̇ᵢⱼ = compute_ϵ̇ij_2(r,D,Ḋ,σᵢⱼ,σ̇ᵢⱼ; damaged_allowed=true, keep_2G)

    # build residual
    if keep_2G
        res = SA[ϵ̇ᵢⱼ[1,1],
                -2*r.G*ϵ̇₂₂ + ϵ̇ᵢⱼ[2,2],
                ϵ̇ᵢⱼ[3,3],
                -2*r.G*params.ϵ̇₁₂ + ϵ̇ᵢⱼ[1,2] ]
    else
        res = SA[ϵ̇ᵢⱼ[1,1],
                -ϵ̇₂₂ + ϵ̇ᵢⱼ[2,2],
                ϵ̇ᵢⱼ[3,3],
                -params.ϵ̇₁₂ + ϵ̇ᵢⱼ[1,2] ]*2*r.G
    end
    return res
end
# function vermeer_simple_shear_update_du(u,p,t)
#     # unpacking
#     r, mp, σ₃, ϵ̇₁₂, du_prev, vmax = @view(p[1:6])
#     σ̇₁₁_prev, σ̇₁₂_prev, σ̇ₒₒₚ_prev, ϵ̇₂₂_prev, _ = du_prev
#     σ₁₁, σ₁₂, σₒₒₚ, ϵ₂₂, D = u
#     # get stress tensor
#     σᵢⱼ = SymmetricTensor{2,3}(SA[σ₁₁ σ₁₂ 0 ; σ₁₂ σ₃ 0 ; 0 0 σₒₒₚ])

#     # correct flaws orientation according to stress rotation
#     if length(p) == 8
#         ψ_cor = compute_flaws_angle_wrt_σ1(σᵢⱼ,p)
#         r = Rheology(r,(ψ=ψ_cor,))
#     end

#     # compute damage growth rate
#     KI = compute_KI(r,σᵢⱼ,D)
#     #@show r.ψ KI
#     Ḋ = compute_subcrit_damage_rate(r,KI,D ; vmax)

#     # solve for non linear derivatives
#     nl_u = SA[σ̇₁₁_prev, σ̇₁₂_prev, σ̇ₒₒₚ_prev, ϵ̇₂₂_prev] # first guess is previous iter
#     nl_u = nl_solve(residual_simple_shear,nl_u,u,p,σᵢⱼ,Ḋ)

#     # build back unknown vector
#     du = SA[nl_u[1], nl_u[2], nl_u[3], nl_u[4], Ḋ]
#     #println("du after update : ", du)
#     p[5]=du

#     return du
# end

# function nl_solve(res_func,nl_u,u,p,σᵢⱼ,Ḋ)
#     # unpack
#     r, mp, σ₃, ϵ̇₁₂, du_prev = @view(p[1:5])
#     # scale strain term
#     nl_u_scaled = setindex(nl_u,nl_u[4]*r.G, 4)
#     if all(nl_u .== 0) # in case du_prev is initialized to zero (should be able to remove that)
#         nl_u_scaled = SA[r.G*ϵ̇₁₂, r.G*ϵ̇₁₂, r.G*ϵ̇₁₂, r.G*ϵ̇₁₂]
#     end
#     #result = DiffResults.JacobianResult(nl_u) no need for DiffResults with static vectors
#     for i in 1:mp.solver.newton_maxiter
#         # get residual and its gradient with respect to u
#         ∇res = ForwardDiff.jacobian(nl_u -> res_func(nl_u,u,σᵢⱼ,Ḋ,p), nl_u_scaled)
#         res  = res_func(nl_u_scaled,u,σᵢⱼ,Ḋ,p)
#         #(i == 1) && (norm_res0 = norm(res))
#         #println(∇res)
#         # update u with Newton algo
#         δnl_u = - ∇res\res
#         nl_u_scaled = nl_u_scaled + δnl_u
#         #@debug "δu = $δu"
#         #@debug "typeof(u) = $(typeof(u))"
#         if (norm(res) <= mp.solver.newton_abstol) && (abs(res[3]) > 1e-10) #test plain strain condition
#             @show(norm(res),res[3])
#         end
#         (norm(res) <= mp.solver.newton_abstol) && (break) # @debug("Newton iter $i ending norm res = $(norm(res))") ;

#         (i == mp.solver.newton_maxiter) && @debug("Newton maxiter ending norm res = $(norm(res))")#("max newton iteration reached ($i), residual still higher than abstol with $(norm(res))")
#     end
#     # rescale strain term :
#     nl_u = setindex(nl_u_scaled,nl_u_scaled[4]/r.G, 4)
#     return nl_u
# end

# function residual_simple_shear(du,u,σᵢⱼ,Ḋ,p)
#     # unpacking
#     r, mp, σ₃, ϵ̇₁₂ = @view(p[1:4])
#     σ̇₁₁, σ̇₁₂, σ̇ₒₒₚ, ϵ̇₂₂ = du
#     _, _, _, _, D = u
#     ϵ̇₂₂ = ϵ̇₂₂/r.G # rescale ϵ̇₂₂

#     # build σ̇ and ϵ̇ from du
#     #zeroTdu = zero(eltype(du))
#     σ̇ᵢⱼ = SymmetricTensor{2,3}([ σ̇₁₁    σ̇₁₂   0.0  ;
#                                  σ̇₁₂    0.0   0.0  ;
#                                  0.0    0.0   σ̇ₒₒₚ ])
#     ϵ̇ᵢⱼ = compute_ϵ̇ij_2(r,D,Ḋ,σᵢⱼ,σ̇ᵢⱼ; damaged_allowed=true)

#     # build residual
#     res = SA[ϵ̇ᵢⱼ[1,1],
#              ϵ̇₂₂ - ϵ̇ᵢⱼ[2,2],
#              ϵ̇ᵢⱼ[3,3],
#              ϵ̇₁₂ - ϵ̇ᵢⱼ[1,2]] / ϵ̇₁₂
#     return res
# end

# # SHEAR STRESS CONTROL

# function vermeer_simple_shear_sigma_update_du(u,p,t ; log_eps=false)
#     # unpacking
#     r, mp, σ₃, σ̇₁₂, du_prev = @view(p[1:5])
#     σ̇₁₁_p, σ̇ₒₒₚ_p, ϵ̇₂₂_p, ϵ̇₁₂_p, _ = du_prev
#     σ₁₁, σₒₒₚ, ϵ₂₂, ϵ₁₂, D = u
#     # get stress tensor
#     σᵢⱼ = SymmetricTensor{2,3}(SA[σ₁₁ σ̇₁₂*t 0 ; σ̇₁₂*t σ₃ 0 ; 0 0 σₒₒₚ])

#     # correct flaws orientation according to stress rotation
#     if length(p) == 7
#         ψ_cor = compute_flaws_angle_wrt_σ1(σᵢⱼ,p)
#         r = Rheology(r,(ψ=ψ_cor,))
#     end

#     # compute damage growth rate
#     KI = compute_KI(r,σᵢⱼ,D)
#     #@show r.ψ KI
#     Ḋ = compute_subcrit_damage_rate(r,KI,D)

#     # solve for non linear derivatives
#     nl_u = SA[σ̇₁₁_p, σ̇ₒₒₚ_p, ϵ̇₂₂_p, ϵ̇₁₂_p] # first guess is previous iter
#     nl_u = nl_solve_sigma(residual_simple_shear_sigma,nl_u,u,p,σᵢⱼ,Ḋ ; log_eps)

#     # build back unknown vector
#     du = SA[nl_u[1], nl_u[2], nl_u[3], nl_u[4], Ḋ]
#     #println("du after update : ", du)
#     p[5]=du

#     return du
# end

# function nl_solve_sigma(res_func,nl_u,u,p,σᵢⱼ,Ḋ ; log_eps=true)
#     # unpack
#     r, mp, σ₃, σ̇₁₂, du_prev = @view(p[1:5])
#     # scale strain term
#     nl_u_scaled = nl_u
#     if log_eps
#         #@show nl_u[4]
#         log_ϵ̇₁₂ = log10(sign(σ̇₁₂)*nl_u[4]) # take the log to force sign(ϵ̇₁₂) = sign(σ̇₁₂)
#         #@show logged_ϵ̇₁₂
#         # scale everything on the order of 1
#         nl_u_scaled = SA[nl_u[1]/σ̇₁₂, nl_u[2]/σ̇₁₂, nl_u[3]*r.G/σ̇₁₂, log_ϵ̇₁₂]
#     else
#         nl_u_scaled = SA[nl_u[1], nl_u[2], nl_u[3]*r.G, nl_u[4]*r.G]
#     end

#     if all(nl_u .== 0) # in case du_prev is initialized to zero (should be able to remove that)
#         nl_u_scaled = SA[σ̇₁₂, σ̇₁₂, σ̇₁₂, σ̇₁₂]
#     end
#     #result = DiffResults.JacobianResult(nl_u) no need for DiffResults with static vectors
#     for i in 1:mp.solver.newton_maxiter
#         # get residual and its gradient with respect to u
#         ∇res = ForwardDiff.jacobian(nl_u -> res_func(nl_u,u,σᵢⱼ,Ḋ,p ; log_eps), nl_u_scaled)
#         res  = res_func(nl_u_scaled,u,σᵢⱼ,Ḋ,p; log_eps)
#         #(i == 1) && (norm_res0 = norm(res))
#         #println(∇res)
#         # update u with Newton algo
#         δnl_u = - ∇res\res
#         nl_u_scaled = nl_u_scaled + δnl_u
#         #@debug "δu = $δu"
#         #@debug "typeof(u) = $(typeof(u))"
#         #println(norm(res))
#         (norm(res) <= mp.solver.newton_abstol) && (break) # @debug("Newton iter $i ending norm res = $(norm(res))") ;

#         (i == mp.solver.newton_maxiter) && @debug("Newton maxiter ending norm res = $(norm(res))")#("max newton iteration reached ($i), residual still higher than abstol with $(norm(res))")
#     end
#     # rescale everything :
#     if log_eps
#         nl_u = SA[nl_u_scaled[1]*σ̇₁₂, nl_u_scaled[2]*σ̇₁₂, nl_u_scaled[3]*σ̇₁₂/r.G, sign(σ̇₁₂)*exp10(nl_u_scaled[4])]
#     else
#         nl_u = SA[nl_u_scaled[1], nl_u_scaled[2], nl_u_scaled[3]/r.G, nl_u_scaled[4]/r.G]
#     end

#     return nl_u
# end

# function residual_simple_shear_sigma(du,u,σᵢⱼ,Ḋ,p ; log_eps=false)
#     # unpacking
#     r, mp, σ₃, σ̇₁₂ = @view(p[1:4])
#     scaled_σ̇₁₁, scaled_σ̇ₒₒₚ, scaled_ϵ̇₂₂, scaled_ϵ̇₁₂ = du
#     σ₁₁, σₒₒₚ, ϵ₂₂, ϵ₁₂, D = u

#     # rescale
#     if log_eps
#         du = SA[scaled_σ̇₁₁*σ̇₁₂, scaled_σ̇ₒₒₚ*σ̇₁₂, scaled_ϵ̇₂₂*σ̇₁₂/r.G, sign(σ̇₁₂)*exp10(scaled_ϵ̇₁₂)]
#     else
#         du = SA[scaled_σ̇₁₁, scaled_σ̇ₒₒₚ, scaled_ϵ̇₂₂/r.G, scaled_ϵ̇₁₂/r.G]
#     end
#     σ̇₁₁, σ̇ₒₒₚ, ϵ̇₂₂, ϵ̇₁₂ = du

#     # build σ̇ and ϵ̇ from du
#     #zeroTdu = zero(eltype(du))
#     σ̇ᵢⱼ = Tensor{2,3}([σ̇₁₁      σ̇₁₂     0  ;
#                         σ̇₁₂    0   0  ;
#                         0  0     σ̇ₒₒₚ  ])
#     ϵ̇ᵢⱼ = compute_ϵ̇ij_2(r,D,Ḋ,σᵢⱼ,σ̇ᵢⱼ; damaged_allowed=true)

#     # build residual
#     res = SA[ϵ̇ᵢⱼ[1,1],
#              ϵ̇₂₂ - ϵ̇ᵢⱼ[2,2],
#              ϵ̇ᵢⱼ[3,3],
#              ϵ̇₁₂ - ϵ̇ᵢⱼ[1,2]] / (σ̇₁₂/r.G)
#     return res
# end

# function compute_flaws_angle_wrt_σ1(σᵢⱼ,p)
#     _, _, _, control_rate, _, ψᵢ, σ1_initial_orientation = p
#     shear_strain_sign = (control_rate > 0) ? :positive : :negative
#     σ1_orientation = get_stress_deviation_from_y(σᵢⱼ ; shear_mode=:simple, shear_strain_sign)
#     Δσ1_orientation = σ1_orientation - σ1_initial_orientation
#     ψ_cor = ψᵢ - Δσ1_orientation
#     # ψ_cor has to be inside [-pi/2,pi/2] therefore:
#     (ψ_cor > 90) && (ψ_cor -= 180)
#     (ψ_cor < -90) && (ψ_cor += 180)
#     return ψ_cor
# end
# function vermeer_simple_shear_update_du_2(du,u,p,t)
#     # unpacking
#     r, mp, σ₁₁, σ₃, ϵ̇₁₂ = p
#     σ̇₁₂, σ̇ₒₒₚ, ϵ̇₁₁, ϵ̇₂₂, Ḋ = du
#     σ₁₂, σₒₒₚ, ϵ̇₁₁, ϵ₂₂, D = u

#     # build σ and D from u
#     σᵢⱼ = Tensor{2,3}([σ₁₁ σ₁₂ 0 ; σ₁₂ σ₃ 0 ; 0 0 σₒₒₚ])
#     KI = compute_KI(r,σᵢⱼ,D)
#     Ḋ = compute_subcrit_damage_rate(r,KI,D)

#     nl_u = SA[σ̇₁₂, σ̇ₒₒₚ, ϵ̇₁₁, ϵ̇₂₂]
#     if all(nl_u .== 0)
#         nl_u = SA[ϵ̇₁₂*2r.G, ϵ̇₁₂*2r.G, ϵ̇₁₂, ϵ̇₁₂]
#     end
#     #result = DiffResults.JacobianResult(nl_u) no need for DiffResults with static vectors
#     for i in 1:mp.solver.newton_maxiter
#         # get residual and its gradient with respect to u
#         ∇res = ForwardDiff.jacobian(nl_u -> residual_simple_shear_2(nl_u,u,σᵢⱼ,Ḋ,p), nl_u)
#         res  = residual_simple_shear_2(nl_u,u,σᵢⱼ,Ḋ,p)

#         #println(∇res)
#         # update u with Newton algo
#         δnl_u = - ∇res\res
#         dnl_u = nl_u + δnl_u
#         #@debug "δu = $δu"
#         #@debug "typeof(u) = $(typeof(u))"
#         #println(norm(res))
#         (norm(res) <= mp.solver.newton_abstol) && break

#         (i == mp.solver.newton_maxiter) && @debug("ending norm res = $(norm(res))")#("max newton iteration reached ($i), residual still higher than abstol with $(norm(res))")
#     end
#     du[:] .= [nl_u[1], nl_u[2], nl_u[3], nl_u[4], Ḋ]
#     println("du after update : ", du)
#     return du
# end

# function residual_simple_shear_2(du,u,σᵢⱼ,Ḋ,p)
#     # unpacking
#     r, mp, σ₁₁, σ₃, ϵ̇₁₂ = p
#     σ̇₁₂, σ̇ₒₒₚ, ϵ̇₁₁, ϵ̇₂₂ = du
#     σ₁₁, σ₁₂, σₒₒₚ, ϵ₂₂, D = u

#     # build σ̇ and ϵ̇ from du
#     #zeroTdu = zero(eltype(du))
#     σ̇ᵢⱼ = Tensor{2,3}([  0    σ̇₁₂   0  ;
#                         σ̇₁₂    0    0  ;
#                          0     0   σ̇ₒₒₚ ])
#     ϵ̇ᵢⱼ = compute_ϵ̇ij_2(r,D,Ḋ,σᵢⱼ,σ̇ᵢⱼ; damaged_allowed=true)

#     # build residual
#     res = SA[ϵ̇₁₁ - ϵ̇ᵢⱼ[1,1],
#              ϵ̇₂₂ - ϵ̇ᵢⱼ[2,2],
#              ϵ̇ᵢⱼ[3,3],
#              ϵ̇₁₂ - ϵ̇ᵢⱼ[1,2]] * 2*r.G
#     return res
# end

# function vermeer_simple_shear_update_du!(du,u,p,t)
#     # unpacking
#     r, mp, σ₃, ϵ̇₁₂ = p
#     σ̇₁₁, σ̇₁₂, σ̇ₒₒₚ, ϵ̇₂₂, Ḋ = du
#     σ₁₁, σ₁₂, σₒₒₚ, ϵ₂₂, D = u

#     # build σ and D from u
#     σᵢⱼ = Tensor{2,3}([σ₁₁ σ₁₂ 0 ; σ₁₂ σ₃ 0 ; 0 0 σₒₒₚ])
#     KI = compute_KI(r,σᵢⱼ,D)
#     Ḋ = compute_subcrit_damage_rate(r,KI,D)

#     nl_u = SA[σ̇₁₁, σ̇₁₂, σ̇ₒₒₚ, ϵ̇₂₂]
#     if all(nl_u .== 0)
#         nl_u = SA[1e-6 .* rand(4)...]
#     end
#     #result = DiffResults.JacobianResult(nl_u) no need for DiffResults with static vectors
#     res0 = 0.0
#     for i in 1:mp.solver.newton_maxiter
#         # get residual and its gradient with respect to u
#         ∇res = ForwardDiff.jacobian(nl_u -> residual_simple_shear(nl_u,u,σᵢⱼ,Ḋ,p), nl_u)
#         res  = residual_simple_shear(nl_u,u,σᵢⱼ,Ḋ,p)
#         (i == 1) && (norm_res0 = norm(res))
#         #println(∇res)
#         # update u with Newton algo
#         δnl_u = - ∇res\res
#         nl_u = nl_u + δnl_u
#         #@debug "δu = $δu"
#         #@debug "typeof(u) = $(typeof(u))"
#         #println(norm(res))
#         (norm(res) <= mp.solver.newton_abstol) && (@debug("Newton iter $i ending norm res = $(norm(res))") ; break)

#         (i == mp.solver.newton_maxiter) && @debug("Newton maxiter ending norm res = $(norm(res))")#("max newton iteration reached ($i), residual still higher than abstol with $(norm(res))")
#     end
#     du[:] .= [nl_u[1], nl_u[2], nl_u[3], nl_u[4], Ḋ]
#     #println("du after update : ", du)
#     return du
# end

### old pure shear related time integration ###
function time_integration(r,p,σᵢⱼ_i,ϵᵢⱼ_i,D_i,ϵ̇11,Δt,tspan)
    flags = p.flags
    t_vec = tspan[1]:Δt:tspan[2]
    σᵢⱼ_vec = Vector{SymmetricTensor{2,3}}(undef,length(t_vec))
    ϵᵢⱼ_vec = similar(σᵢⱼ_vec)
    D_vec = Vector{Float64}(undef,length(t_vec))
    σᵢⱼ_vec[1] = σᵢⱼ_i
    ϵᵢⱼ_vec[1] = ϵᵢⱼ_i
    D_vec[1] = D_i
    last_tsim_printed = 0.0
    for i in 2:length(t_vec)
      set_print_flag!(p,i,t_vec[i-1])
      flags.print && print_time_iteration(i,t_vec[i-1])

      σᵢⱼ_vec[i], ϵᵢⱼ_vec[i], D_vec[i], u = solve(r,p,σᵢⱼ_vec[i-1],ϵᵢⱼ_vec[i-1],D_vec[i-1],ϵ̇11,Δt)
    end
    return t_vec, σᵢⱼ_vec, ϵᵢⱼ_vec, D_vec
  end


text:
#
# Algorithm of:
#
# NONMONOTONE SPECTRAL PROJECTED GRADIENT METHODS ON CONVEX SETS
# ERNESTO G. BIRGIN, JOSÉ MARIO MARTÍNEZ, AND MARCOS RAYDAN
# SIAM J. O. PTIM. Vol. 10, No. 4, pp. 1196-1211
#
# Implemented by J. M. Martínez (IMECC - UNICAMP)
# Iinitally translated and adapted to Julia by L. Martínez (IQ-UNICAMP)
#
"""
```
spgbox!(f, g!, x::AbstractVecOrMat; lower=..., upper=..., options...)
```

```
spgbox!(f, g!, lower::AbstractVeOrMat, upper::AbstractVecOrMat, x::AbstractVecOrMat; options...)
```

Minimizes function `f` starting from initial point `x`, given the function to compute the gradient, `g!`. `f` must be of the form `f(x)`, and `g!` of the form `g!(g,x)`, where `g` is the gradient vector to be modified. It modifies the `x` vector, which will contain the best solution found (see `spgbox` for a non-mutating alternative).

Optional lower and upper box bounds can be provided using optional arguments `lower` and `upper`, which can be provided as the fourth and fifth arguments or with keyword parameters.

Returns a structure of type `SPGBoxResult`, containing the best solution found in `x` and the final objective function in `f`.

Alternativelly, a single function that computes the function value and the gradient can be provided, using:

```
spgbox(fg!, x; lower=..., upper=..., options...)
```
or
```
spgbox(fg!, lower upper, x; options...)
```

The `fg!` must be of the form `fg!(g,x)` where `x` is the current point and `g` the array that stores the gradient. And it must return
the function value.

# Examples
```jldocstest
julia> f(x) = x[1]^2 + x[2]^2

julia> function g!(g,x)
           g[1] = 2*x[1]
           g[2] = 2*x[2]
       end
```
## Without bounds

```jldocstest
julia> x = rand(2)

julia> spgbox!(f,g!,x)

SPGBOX RESULT:

Convergence achieved.

Final objective function value = 0.0
Sample of best point = Vector{Float64}[ 0.0, 0.0]
Projected gradient norm = 0.0

Number of iterations = 3
Number of function evaluations = 3

```

## With bounds
```jldocstest
julia> x = 2 .+ rand(2)

julia> spgbox!(f,g!,x,lower=[2.,-Inf])

SPGBOX RESULT:

Convergence achieved.

Final objective function value = 4.0
Sample of best point = Vector{Float64}[ 2.0, 0.0]
Projected gradient norm = 0.0

Number of iterations = 3
Number of function evaluations = 3

```

## With a single function to compute the function and the gradient

```julia-repl
julia> function fg!(g,x)
           g[1] = 2*x[1]
           g[2] = 2*x[2]
           fx = x[1]^2 + x[2]^2
           return fx
       end
fg! (generic function with 1 method)

julia> x = rand(2);

julia> spgbox(fg!,x)

 SPGBOX RESULT:

 Convergence achieved.

 Final objective function value = 0.0
 Sample of best point = Vector{Float64}[ 0.0, 0.0]
 Projected gradient norm = 0.0

 Number of iterations = 3
 Number of function evaluations = 3
```

"""
function spgbox!(f::F, g!::G, x::AbstractVecOrMat{T}; kargs...) where {F<:Function,G<:Function,T}
    spgbox!((g, x) -> begin
        g!(g, x)
        return f(x)
    end, x, f; kargs...)
end

#
# Call with a single function to compute the function and the gradient
#
function spgbox!(
    fg!::Function,
    x::AbstractVecOrMat{T},
    func_only = nothing;
    lower::Union{Nothing,AbstractVecOrMat{T}} = nothing,
    upper::Union{Nothing,AbstractVecOrMat{T}} = nothing,
    eps = oneunit(T) / 100_000,
    nitmax::Int = 100,
    nfevalmax::Int = 1000,
    m::Int = 10,
    vaux::VAux = VAux(x, (isnothing(func_only) ? fg!(similar(x), x) : func_only(x)), m = m),
    iprint::Int = 0,
    project_x0::Bool = true,
) where {T}
    # Adimentional variation of T (base Number type)
    adT = typeof(one(T))

    # Number of variables
    n = length(x)

    # Auxiliary arrays (associate names and check dimensions)
    g = vaux.g
    xn = vaux.xn
    gn = vaux.gn
    fprev = vaux.fprev
    @assert length(g) == n "Auxiliar gradient vector `g` must be of the same length as `x`"
    @assert length(xn) == n "Auxiliar vector `xn` must be of the same length as `x`"
    @assert length(gn) == n "Auxiliar vector `gn` must be of the same length as `x`"
    @assert length(fprev) == m "Auxiliar vector `fprev` must be of length `m`"

    # Check if bounds are defined, project or not the initial point on them
    if !isnothing(lower)
        @assert length(lower) == n "Lower bound vector `lower` must be of the same length than x, got: $(length(lower))"
        if project_x0
            @. x = max(x, lower)
        else
            for i in eachindex(x)
                if x[i] < lower[i]
                    error(
                        " Initial value of variable $i smaller than lower bound, and `project_x0` is set to `false`. ",
                    )
                end
            end
        end
    end
    if !isnothing(upper)
        @assert length(upper) == n "Upper bound vector `upper` must be of the same length than `x`, got: $(length(upper))"
        if project_x0
            @. x = min(x, upper)
        else
            for i in eachindex(x)
                if x[i] > lower[i]
                    error(
                        " Initial value of variable $i greater than upper bound, and `project_x0` is set to `false`. ",
                    )
                end
            end
        end
    end

    # Iteration counter
    nit = 1

    # Objective function and gradient at initial point
    nfeval = 1
    fcurrent = fg!(g, x)
    gnorm = pr_gradnorm(g, x, lower, upper)
    gnorm <= eps && return SPGBoxResult(x, fcurrent, gnorm, nit, nfeval, 0)

    # Do a consertive initial step
    small = T(sqrt(Base.eps(T)))
    tspg = small / max(T(1), gnorm)

    # Initialize array of previous function values
    # Allow slight nonmonotonicity
    for i in eachindex(fprev)
        fprev[i] = fcurrent + abs(fcurrent) / 10
    end

    while nit < nitmax

        if iprint > 0
            println("----------------------------------------------------------- ")
            println(" Iteration: ", nit)
            println(" x = ", x[begin], " ... ", x[end])
            println(" Objective function value = ", fcurrent)
            println(" ")
            println(" Norm of the projected gradient = ", gnorm)
            println(" Number of function evaluations = ", nfeval)
        end

        fref = maximum(fprev)
        if iprint > 2
            println(" fref = ", fref)
            println(" fprev = ", fprev)
            println(" t = ", t)
        end

        compute_xn!(xn, x, tspg, g, lower, upper)

        lsfeval, fn =
            safequad_ls(xn, gn, x, g, fcurrent, tspg, fref, nfevalmax - nfeval, iprint, func_only, fg!, lower, upper)

        if lsfeval < 0
            return SPGBoxResult(x, fcurrent, gnorm, nit, nfeval - lsfeval, 2)
        else
            nfeval += lsfeval
        end

        # Trial point accepted
        num = zero(T)
        den = zero(T)
        for i in eachindex(x)
            num = num + (xn[i] - x[i])^2 / oneunit(T)
            den = den + (xn[i] - x[i]) * (gn[i] - g[i]) / oneunit(T)
        end
        if den <= zero(T)
            tspg = adT(100)
        else
            tspg = max(min(adT(1.0e30), num / den), adT(1.0e-30))
        end
        fcurrent = fn
        for i in eachindex(x)
            x[i] = xn[i]
            g[i] = gn[i]
        end
        for i = firstindex(fprev):lastindex(fprev)-1
            fprev[i] = fprev[i+1]
        end
        fprev[end] = fcurrent
        nit = nit + 1

        # Compute projected gradient norm
        gnorm = pr_gradnorm(g, x, lower, upper)
        gnorm <= eps && return SPGBoxResult(x, fcurrent, gnorm, nit, nfeval, 0)
    end

    # Maximum number of iterations achieved
    return SPGBoxResult(x, fcurrent, gnorm, nit, nfeval, 1)
end

"Perform a safeguarded quadratic line search"
function safequad_ls(
    xn::AbstractVecOrMat{T},
    gn::AbstractVecOrMat{T},
    x::AbstractVecOrMat{T},
    g::AbstractVecOrMat{T},
    fcurrent::Number,
    tspg::Number,
    fref::Number,
    nfevalmax::Int,
    iprint::Int,
    func_only::Union{Nothing,Function},
    fg!::Function,
    lower::Union{Nothing,AbstractVecOrMat{T}},
    upper::Union{Nothing,AbstractVecOrMat{T}},
) where {T}
    one_T = one(T)
    # Armijo parameter
    gamma = one_T / 10_000

    gtd = zero(fref)
    for i in eachindex(x)
        gtd += (xn[i] - x[i])*g[i]
    end

    # Compute function values at initial point
    fn = fg!(gn, xn)
    nfeval = 1
    nfeval > nfevalmax && return -nfeval, fn

    alpha, trials = one_T, 0
    while true
        trials += 1
        if trials > 1
            # Compute a new trial point and its function values
            compute_xn!(xn, x, alpha * tspg, g, lower, upper)
            gtd = zero(T)
            for i in eachindex(x)
                gtd += (xn[i] - x[i]) * g[i]
            end
            if iprint > 2
                println(" xn = ", xn[begin], " ... ", xn[end])
                println(" f[end] = ", fn, " fref = ", fref)
            end
            if !isnothing(func_only)
                fn = func_only(xn)
            else
                fn = fg!(gn, xn)
            end
            nfeval += 1
            nfeval > nfevalmax && return -nfeval, fn
        end

        # If the point is not acceptable
        if fn >= fref + gamma * gtd
            # Perform a safeguarded quadratic interpolation step
            if alpha <= one_T / 10
                alpha /= 2
            else
                atemp = -gtd * alpha^2 / (2 * (fn - fcurrent - alpha * gtd))
                if atemp <= one_T / 10 || atemp >= 9 * one_T / 10 * alpha
                    atemp = alpha / 2
                end
                alpha = atemp
            end
            # If the point is acceptable
        else
            # Update the gradient at the accepted point, if necessary
            if trials > 1 && !isnothing(func_only)
                fn = fg!(gn, xn)
                nfeval += 1
                nfeval > nfevalmax && return -nfeval, fn
            end
            break
        end
    end

    return nfeval, fn
end

#
# Call with lower and upper as positional arguments
#
spgbox!(f::F, g!::G, lower, upper, x; kargs...) where {F,G} = spgbox!(f, g!, x, lower = lower, upper = upper, kargs...)
spgbox!(fg!::F, lower, upper, x; kargs...) where {F,G} = spgbox!(fg!, x, lower = lower, upper = upper, kargs...)

"""

See `spgbox!` for additional help.

```
spgbox(f, g!, x::AbstractVecOrMat; lower=..., upper=..., options...)`
spgbox(f, g!, lower::AbstractVecOrMat, upper::AbstractVecOrMat, x::AbstractVecOrMat; options...)`
```
Minimizes function `f` starting from initial point `x`, given the function to compute the gradient, `g!`.
`f` must be of the form `f(x)`, and `g!` of the form `g!(g,x)`, where `g` is the gradient vector to be modified.

Optional lower and upper box bounds can be provided using optional arguments `lower` and `upper`, which can be provided as the fourth and fifth arguments or with keyword parameters.

```
spgbox(fg!, x::AbstractVecOrMat; lower=..., upper=..., options...)`
spgbox(fg!, lower::AbstractVecOrMat, upper::AbstractVecOrMat, x::AbstractVecOrMat; options...)`
```

Given a single function `fg!(g,x)` that updates a gradient vector `g` and returns the function value, minimizes the function.

These functions return a structure of type `SPGBoxResult`, containing the best solution found in `x` and the final objective function in `f`.

These functions *do not* mutate the `x` vector, instead it will create a (deep)copy of it (see `spgbox!` for the in-place alternative).

"""
function spgbox(f::F, g!::G, x::AbstractVecOrMat{T}; kargs...) where {F<:Function,G<:Function,T}
    x0 = deepcopy(x)
    return spgbox!(f, g!, x0; kargs...)
end
# With a single function to compute the function and the gradient
function spgbox(fg::FG, x::AbstractVecOrMat{T}; kargs...) where {FG<:Function,T}
    x0 = deepcopy(x)
    return spgbox!(fg, x0; kargs...)
end
#
# call with lower and upper as positional arguments
#
function spgbox(f::F, g!::G, lower, upper, x::AbstractVecOrMat{T}; kargs...) where {F<:Function,G<:Function,T}
    x0 = deepcopy(x)
    return spgbox!(f, g!, lower, upper, x0; kargs...)
end
# With a single function to compute the function and the gradient
function spgbox(fg::FG, lower, upper, x::AbstractVecOrMat{T}; kargs...) where {FG<:Function,T}
    x0 = deepcopy(x)
    return spgbox!(fg, lower, upper, x0; kargs...)
end


text:
## Plot recipe
## brute force this to avoid vertical asymptotes

@recipe function f(pq::AbstractRationalFunction{T}, a=nothing, b=nothing) where {T <: Real}

    p,q = lowest_terms(convert(RationalFunction,pq), method=:numerical)

    if a==nothing && b==nothing
        u=Polynomials.poly_interval(p)
        v=Polynomials.poly_interval(q)
        a,b = min(first(u), first(v)), max(last(u), last(v))
    end

    n = 251
    Δ = (b-a)/n

    # identify vertical asumptotes
    zs = filter(z -> a+Δ < z < b-Δ, sort(real.(filter(isreal, roots(q)))))
    xs = Any[]
    ys = Any[]
    if length(zs) > 0
        z₁ = popfirst!(zs)
        xs′ = range(a, stop=z₁ - Δ, length=50)
        ys′ = pq.(xs′)
        push!(xs, xs′)
        push!(ys, ys′)
        push!(xs, z₁)
        push!(ys, NaN)
        for z₂ ∈ zs
            xs′ = range(z₁ + Δ, stop=z₂ - Δ, length=50)
            ys′ = pq.(xs′)
            push!(xs, xs′)
            push!(ys, ys′)
            push!(xs, z₂)
            push!(ys, NaN)
            z₁ = z₂
        end
        xs′ = range(z₁+Δ, stop=b, length=50)
        ys′ = pq.(xs′)
        push!(xs, xs′)
        push!(ys, ys′)

        xs = vcat(xs...)
        ys = vcat(ys...)
    else
        xs = range(a, stop=b, length=251)
        ys = pq.(xs)
    end

    xs, ys
end


text:
"""
#Column descriptions:

* 1st: 5-digit Nutrient Databank number that uniquely identifies a food item. If this field is defined as numeric, the leading zero will be lost
* 2nd: 4-digit code indicating food group to which a food item belongs. (for details see FD_GROUP.txt).
* 3rd: long description. 200-character description of food item
* 4th: short description
* 5th: Other names commonly used to describe a food, including local or regional names for various foods, for example, “soda” or “pop” for “carbonated beverages.

For more info and description of other fields, see p. 30 of sr28_doc.pdf
"""
function getfooddes(data_folder; filename="FOOD_DES.txt")
  food_des_file = joinpath(data_folder, filename)

  food_des = CSV.read(food_des_file, DataFrame, header=0, delim='^', quotechar='~');
  food_des = food_des[:, 1:5];  # only the columns that I need
  colnames = [:food, :group, :longdes, :shortdes, :othernames]
  for cn in 1:size(food_des)[2]
    rename!(food_des, names(food_des)[cn] => colnames[cn])
  end
  return food_des
end

"""
The Nutrient Data file contains mean nutrient values per 100 g of the edible portion of food, along with fields to further describe the mean value.

#  Column descriptions:

* 1st: 5-digit Nutrient Databank number that uniquely identifies a food item.
* 2nd: Unique 3-digit identifier code for a nutrient
* 3rd: Amount in 100 grams, edible portion
* 4th:
* 5th: Standard error of the mean.

For more info and description of other fields, see p. 32 of sr28_doc.pdf
"""
function getnutdata(data_folder; filename="NUT_DATA.txt")
  nut_data_file = joinpath(data_folder, filename)
  nut_data = CSV.read(nut_data_file, DataFrame, header=0, delim='^', quotechar='~');
  nut_data = nut_data[:, [1,2,3,5]];
  colnames = [:food, :nut, :amount, :stderr]
  for cn in 1:size(nut_data)[2]
    rename!(nut_data, names(nut_data)[cn] => colnames[cn])
  end
  return nut_data
end

"""
This file (Table 9) is a support file to the Nutrient Data file. It provides the 3-digit nutrient code, unit of measure, INFOODS tagname, and description.

# Column descriptions:

* 1st: Unique 3-digit identifier code for a nutrient.
* 2nd: Units of measure (mg, g, μg, and so on).

For more info and description of other fields, see p. 34 of sr28_doc.pdf
"""
function getnutrdef(data_folder; filename="NUTR_DEF.txt")
  nutr_def_file = joinpath(data_folder, filename)

  nutr_def = CSV.read(nutr_def_file, DataFrame, header=0, delim='^', quotechar='~');
  nutr_def = nutr_def[:, 1:4];
  # convert values same as  nutr_def[31, 2] to ug (i.e. micrograms)
  udefvar = nutr_def[31, 2]
  udefrows = findall(a -> a == udefvar, nutr_def[:, 2])
  nutr_def[udefrows, 2] .= "ug"

  colnames = [:nut, :unit, :nuttag, :nutname]
  for cn in 1:size(nutr_def)[2]
    rename!(nutr_def, names(nutr_def)[cn] => colnames[cn])
  end
  return nutr_def
end

function get_food_groups(data_folder)
  group_file = joinpath(data_folder, "FD_GROUP.txt")
  groups = CSV.read(group_file, DataFrame, header=0, delim='^', quotechar='~');
  return groups
end

"""
applies the updates to the data files.
"""
function apply_updates(food_des, nut_data, updatedir)
  # add new data
  newfood = getfooddes(updatedir, filename="ADD_FOOD.txt")
  newnut = getnutdata(updatedir, filename="ADD_NUTR.txt")
  food_des = vcat(food_des, newfood)
  nut_data = vcat(nut_data, newnut)

  # apply changes to the data
  foodc = getfooddes(updatedir, filename="CHG_FOOD.txt")
  nutc = getnutdata(updatedir, filename="CHG_NUTR.txt")

  for rr in 1:size(foodc)[1]
    food_des[.&(food_des[1] .== foodc[rr, 1], food_des[2] .== foodc[rr, 2]), 3:end] = DataFrame(foodc[rr, 3:end])
  end

  for rr in 1:size(nutc)[1]
    nut_data[.&(nut_data[1] .== nutc[rr, 1], nut_data[2] .== nutc[rr, 2]), 3:end] = DataFrame(nutc[rr, 3:end])
  end

  # remove deleted entries in the update
  # foodd = getfooddes(updatedir, filename="DEL_FOOD.txt")
  nutd_file = joinpath(updatedir, "DEL_NUTR.txt")
  nutd = CSV.read(nutd_file, DataFrame, header=0, delim='^', quotechar='~');
  for rr in 1:size(nutd)[1]
    nut_data = nut_data[.~.&(nut_data[1] .== nutd[rr, 1], nut_data[2] .== nutd[rr, 2]), :]
  end

  # TODO: update weights too
  "ADD_WGT.txt"
  "CHG_WGT.txt"
  "DEL_WGT.txt"

  return food_des, nut_data
end

function fix_DRI_table(;dritable::String="dri.csv", data_folder::String="sr28asc")
  counter = 0
  newfields = []
  for line in eachline(open(dritable))
    counter += 1
    if counter > 1
      fields = split(line, ",")
      dri = fields[5]
      if endswith(dri, "*") || ismatch(r"\d+\D", dri)
        dri = dri[1:end-1]
      end
      entry = vcat(fields[1:4], [dri], fields[6:end])
      push!(newfields, entry)
    end
  end

  outfile = "dri_new.csv"
  j = open(outfile, "w")
  println(j, readline(open(dritable)))
  close(j)
  open(outfile, "a") do ff
    for item in newfields
      println(ff, join(item, ","))
    end
  end
  dd = CSV.read("dri_new.csv", DataFrame, header=1, delim=',')

  nutr_def = getnutrdef(data_folder);
  dd = CSV.read("dri_new.csv", DataFrame, header=1, delim=',');
  availables = Dict()
  notavailables = Set()
  nuts = levels(dd[:Nutrient])
  for rr in 1:size(nuts)[1]
    item = nuts[rr]
    found = false
    for r2 in 1:size(nutr_def)[1]
      if contains(lowercase(nutr_def[r2, :nutname]), lowercase(item))
        if haskey(availables, item)
          push!(availables[item], nutr_def[r2, :nut])
        else
          availables[item] = [nutr_def[r2, :nut]]
          found = true
        end
      end
    end
    if ~found
      push!(notavailables, item)
    end
  end

end

"""
Load the cleaned and read dri table
"""
function load_dri(;filepath="dri_new.csv")
  dri = CSV.read(filepath, DataFrame, header=1, delim=',', missingstring="NA")
  return dri
end

"""
getnutdata(data_folder)

Return `nutrs` (a dataframe of the data), `dri` (a dataframe of recommended dosages), and `groups` (food groups).

Extracts data from SR28 database (https://www.ars.usda.gov/northeast-area/beltsville-md/beltsville-human-nutrition-research-center/nutrient-data-laboratory/docs/sr28-download-files/)

The four principal database files are the Food Description file, Nutrient Data file, Gram Weight file, and Footnote file.
"""
function allnutdata(;data_folder::String="./data/sr28asc", dri_file::String="./data/dri_new.csv", updatedir::String="")
  if ~isdir(data_folder)
    throw("Data directory does not exist.")
  end

  data_folder = abspath(data_folder)

  food_des = getfooddes(data_folder);
  nut_data = getnutdata(data_folder);
  nutr_def = getnutrdef(data_folder);
  # Option: get data from Weight File (file name = WEIGHT). This file (Table 12) contains the weight in grams of a number of common measures for each food item (p. 36)

  # apply updates
  if length(updatedir) > 0
    food_des, nut_data = apply_updates(food_des, nut_data, updatedir)
  end

  df1 = join(food_des, nut_data, on=:food, kind=:inner);
  nutrs = join(df1, nutr_def, on=:nut, kind=:inner);

  dri = load_dri(filepath=dri_file);
  groups = get_food_groups(data_folder);

  # convert all units to grams
  dri = convert_units_to_g(dri, Symbol("Unit(/d)"), [:DRI, :Upper_intake]);
  nutrs = convert_units_to_g(nutrs, :unit, [:amount]);

  return nutrs, dri, groups
end

"""
Returns lists of min and max required nutrients and lists of foods and their nutrients.

# Parameters:

* weight: in kg.
"""
function get_min_max(dri; age=22, sex="male", weight=70)
  if sex == "male" || sex == "m" || sex == "M"
    sex = "Male"
  elseif sex == "female" || sex == "f" || sex == "F"
    sex = "Female"
  end

  select_dri = @from i in dri begin
    @where getproperty(i, Symbol("AgeMin(year)")) <= age
    @where getproperty(i, Symbol("AgeMax(year)")) > age
    # @where typeof(i.id) <: Number
    @where i.Sex == sex || i.Sex == "Neutral"  # TODO: add pregnancy and lactation
    @select i #{i.DRI, i.Unit_d_, i.Upper_intake, i.id}
    @collect DataFrame
  end

  select_dri = select_dri[.~ismissing.(select_dri[:id]), :];

  minlist = select_dri[:DRI]
  maxlist = select_dri[:Upper_intake]
  dri_ids = select_dri[:id]

  # Remove the upper levels where the refernce intakes (minlist) are more than the upper levels (max list). This is from the data. In these cases, the upper intake is for intaking from supplements not from food. So it is safe to remove the upper intake in those cases.
  for rr in 1:length(minlist)
    if !ismissing(minlist[rr]) && !ismissing(maxlist[rr]) && minlist[rr] > maxlist[rr]
      maxlist[rr] = missing
    end
  end

  # fill the missing values
  mins_refilled = Float64[]
	for mm in minlist
		if ismissing(mm)
			push!(mins_refilled, 0.0)
		else
			push!(mins_refilled, mm)
		end
	end
	maxs_refilled = Float64[]
	for mm in maxlist
		if ismissing(mm)
			push!(maxs_refilled, 1e8)
		else
			push!(maxs_refilled, mm)
		end
  end

  # adjust those nutrients whose unit is g/kg to the body weight.
  for index in 1:length(mins_refilled)
    if select_dri[index, Symbol("Unit(/d)")] == "g/kg"
      if mins_refilled[index] != 0
        mins_refilled[index] = weight * mins_refilled[index]
      end
      if maxs_refilled[index] != 1e8
        maxs_refilled[index] = weight * maxs_refilled[index]
      end
    end
  end

  return dri_ids, mins_refilled, maxs_refilled
end

function pushtolist!(list::T, topush::Y, index::Z) where {T<:AbstractArray, Y<:AbstractArray, Z<:Integer}
  if length(topush) >0
    push!(list[index], topush[1])
  else
    push!(list[index], 0.0)
  end
end

"Returns a list of lists, where each inner list is the amount of nutrients in a food"
function getnutamounts(nfoods, foodids, dri_ids, nutrs)
  nutamounts = [Float64[] for i in 1:nfoods]
  prog = ProgressMeter.Progress(nfoods)
  for index in 1:nfoods
    foodid = foodids[index]
    for nutid in dri_ids
      nutamount = nutrs[.&(nutrs[:food] .== foodid, nutrs[:nut] .== nutid), :amount]
      pushtolist!(nutamounts, nutamount, index)
      # if length(nutamount) > 1
      #   warn("More than one nutrient ($nutid) for food ($foodid)")
      # end
    end
    next!(prog)
  end
  return nutamounts
end

"""
    loaddata(your_age, your_sex, your_weight; only_groups=[])

Return `nutrs`, `dri`, `groups`, `foodids`, `foodnames`, `nfoods`, `nutamounts`, `calories`, `dri_ids`, `mins`, `maxs`. They containt the following information:

* `nutrs:` list of foods. A table with the following columns: food (name of food), group (food group), longdes (long food description), shortdes (short food description), othernames, nut (id of a nutrition in the food), amount (amount of the nutrition in 100 g of edible food), stderr (standard error of the amount of the nutrition), unit (unit of the amount of the nutrition), nuttag (nutrition id), nutname (nutrition name).
* `dri:` recommended range of nutritions for various groups. A tale with the following columns: sex, AgeMin(year), AgeMax(year), Nutrient (nutrient name), DRI (recommended amound), Unit(/d)(unit of the recommended amount per day), DRI available (whether there is strong evidence for DRI), Upper intake (maximum amount before it's dangerous), id (nutrition id)
* `groups:` A list of food groups (second column) and their ids (first column).
* `foodids:` A list of all food IDs.
* `foodnames:` A list of all food names, with the same order as `foodids`.
* `nutamounts:` A list of lists, where each inner list has all the nutrient amounts per food. It has the same order as `foodids`.
* `calories:` calories in each food.
* `dri_ids:` A list of nutrient IDs.
* `mins:` The minimum recommended amount for each nutrient in grams per day. It has the same order as `dri_ids`.
* `maxs:` The maximum recommended amount for each nutrient in grams per day. It has the same order as `dri_ids`.
"""
function loaddata(your_age, your_sex, your_weight; only_groups=[])
  # mkpath("../data")
  if !isfile("all_opt_inputs_$(your_age)_$(your_sex).serialized")
    nutrs, dri, groups = allnutdata(data_folder="./data/sr28asc", dri_file="./data/dri_new.csv", updatedir="./data/SR28upd0516");
    foodids = unique(nutrs[:food]);
    foodnames = unique(nutrs[:shortdes]);
    nfoods = length(foodids)
    dri_ids, mins, maxs = get_min_max(dri, age=your_age, sex=your_sex, weight=your_weight);

    if length(only_groups) > 0
      select_rows = Array{Int64}(undef, 0)
      for rr in 1:size(nutrs)[1]
          if in(nutrs[rr, :group], only_groups) && (in(nutrs[rr, :nut], dri_ids) || nutrs[rr, :nut] .== 208)
          push!(select_rows, rr)
        end
      end

      foodids_org = deepcopy(nutrs[:food]);
      nutrs = nutrs[select_rows, :];
      foodids = unique(nutrs[:food]);
      foodnames = unique(nutrs[:shortdes]);
      nfoods = length(foodids);

      calories = Float64[] # calories in each food
      for index in 1:nfoods
        foodid = foodids[index]
        calorie = nutrs[.&(nutrs[:food] .== foodid, nutrs[:nut] .== 208), :amount]
        push!(calories, calorie[1])
      end
      nutamounts = getnutamounts(nfoods, foodids, dri_ids, nutrs);
    else
      nutamounts = getnutamounts(nfoods, foodids, dri_ids, nutrs)
      calories = Float64[] # calories in each food
      for index in 1:nfoods
        foodid = foodids[index]
        calorie = nutrs[.&(nutrs[:food] .== foodid, nutrs[:nut] .== 208), :amount]
        push!(calories, calorie[1])
      end
    end

    serialize("all_opt_inputs_$(your_age)_$(your_sex).serialized", [nutrs, dri, groups, foodids, foodnames, nfoods, nutamounts, calories, dri_ids, mins, maxs])
  else
    nutrs, dri, groups, foodids, foodnames, nfoods, nutamounts, calories, dri_ids, mins, maxs = deserialize("all_opt_inputs_$(your_age)_$(your_sex).serialized");
  end

  return nutrs, dri, groups, foodids, foodnames, nfoods, nutamounts, calories, dri_ids, mins, maxs
end

"Converts mg to g"
function mg_to_g(value::T) where T<:Real
  return value/T(1000)
end

"Converts ug to g"
function ug_to_g(value::T) where T<:Real
  return value/T(1000_000)
end

"Converts kg to g"
function kg_to_g(value::T) where T<:Real
  return value*T(1000)
end

"Converts mg to g"
function mg_to_g!(dd::AbstractDataFrame, unitcolumn::Symbol, columns::AbstractArray{Symbol, 1}, unit::AbstractString)
  for col in columns
    dd[dd[unitcolumn] .== unit, col] = mg_to_g.(dd[dd[unitcolumn] .== unit, col])
  end
end

"Converts ug to g"
function ug_to_g!(dd::AbstractDataFrame, unitcolumn::Symbol, columns::AbstractArray{Symbol, 1}, unit::AbstractString)
  for col in columns
    dd[dd[unitcolumn] .== unit, col] = ug_to_g.(dd[dd[unitcolumn] .== unit, col])
  end
end

"Converts kg to g"
function kg_to_g!(dd::AbstractDataFrame, unitcolumn::Symbol, columns::AbstractArray{Symbol, 1}, unit::AbstractString)
  for col in columns
    dd[dd[unitcolumn] .== unit, col] = kg_to_g.(dd[dd[unitcolumn] .== unit, col])
  end
end

"Converts mg to g"
function mg_to_g(value::Missing)
  return value
end

"Converts ug to g"
function ug_to_g(value::Missing)
  return value
end

"Converts kg to g"
function kg_to_g(value::Missing)
  return value
end

"""
    convert_units_to_g!(d::AbstractDataFrame, unitcolumn::Symbol, columns::AbstractArray{Symbol, 1})

Converts the values in columns `columns` to g or ml.
"""
function convert_units_to_g(d::AbstractDataFrame, unitcolumn::Symbol, columns::AbstractArray{Symbol, 1})
  current_values = unique(d[unitcolumn])
  d = DataFrame(d);
  for vvs in current_values
    if vvs == "g"
      continue
    end
    if vvs == "\xb5g"
      vv = "ug"
    else
      vv = lowercase(vvs)
    end

    if vv == "kg" || vv == "l"
      kg_to_g!(d, unitcolumn, columns, vvs)
      d[d[unitcolumn] .== vvs, unitcolumn] = "g"
    elseif vv == "ug"
      ug_to_g!(d, unitcolumn, columns, vvs)
      d[d[unitcolumn] .== vvs, unitcolumn] = "g"
    elseif vv == "mg" || vv == "mg/kg"
      mg_to_g!(d, unitcolumn, columns, vvs)
      if vv == "mg/kg"
        d[d[unitcolumn] .== vvs, unitcolumn] = "g/kg"
      else
        d[d[unitcolumn] .== vvs, unitcolumn] = "g"
      end
    end
  end
  return d
end

"""
    return_excluded_food_indices(nutrs, exclude_list, foodids)

Return a list of `indices` and `ids` of foods in which the words in `exclude_list` exist.
Use this list to set their amount to zero in the optimization process.
"""
function return_excluded_food_indices(nutrs, exclude_list, foodids)
  ids = Set{Int64}()
  for ex in exclude_list
    for row in 1:size(nutrs, 1)
      if occursin(ex, lowercase(nutrs[row, :longdes]))
        push!(ids, nutrs[row, :food])
      end
    end
  end
  ids = sort(collect(ids))
  indices = Int64[]
  for (index, id) in enumerate(foodids)
    if in(id, ids)
      push!(indices, index)
    end
  end
  return indices, ids
end

"""
Returns foods' `longdes` given their ids
"""
function find_longdes(nutrs, foodids)
  foodlongnames = AbstractString[]
  for id in foodids
    longdes = nutrs[nutrs[:food] .== id, :longdes][1]
    push!(foodlongnames, longdes)
  end
  return foodlongnames
end

"""
Returns foods groups given their ids
"""
function find_groups(nutrs, foodids)
  foodgroup = Int64[]
  for id in foodids
    groupid = nutrs[nutrs[:food] .== id, :group][1]
    push!(foodgroup, groupid)
  end
  return foodgroup
end

function constraints(dri, dri_ids, mins, maxs)
  nutrient_name_ids = [findfirst(x -> .&(!ismissing(dri.id[x]), dri.id[x] == id), 1:size(dri, 1)) for id in dri_ids]
  nutrient_names = dri[nutrient_name_ids, :Nutrient]
  return DataFrame(:nutrient => nutrient_names, :nutrient_id => dri_ids, Symbol("min(g)") => mins, Symbol("max (g)") => maxs)
end

function exclude_nutrients(dri_ids, mins, maxs, nutamounts, exclude_ids)
  new_rows = findall(x-> !in(x, exclude_ids), dri_ids)
  dri_ids = dri_ids[new_rows]
  mins = mins[new_rows]
  maxs = maxs[new_rows]
  for food in 1:length(nutamounts)
    nutamounts[food] = nutamounts[food][new_rows]
  end
  return dri_ids, mins, maxs, nutamounts
end


text:
# Message handling

export handle

handlers(o) = Dict()

handle_message(o, m) = get(handlers(o), m["type"], identity)(m)

handle(f, o, t) = (handlers(o)[t] = f)

# Callback Tasks

const callbacks = Dict{Int,Condition}()

const counter = [0]

function callback!()
  id = (counter[1] += 1)
  cb = Condition()
  callbacks[id] = cb
  return id, cb
end

function callback!(id, value = nothing)
  haskey(callbacks, id) || return
  notify(callbacks[id], value)
  delete!(callbacks, id)
  return
end

function enable_callbacks!(o)
  handle(o, "callback") do m
    callback!(m["callback"], get(m, "result", nothing))
  end
end


text:
# IoU implementation with Julia
# Ref: https://github.com/MartinThoma/algorithms/blob/master/CV/IoU/IoU.py

"""
    IoU(bb1, bb2)

Calculate the Intersection over Union (IoU) of two axis-aligned bounding boxes `bb1` and `bb2`.

Here, `bb1` and `bb2` are provided as `Dict` with keys = {"x1", "x2", "y1", "y2"}, where `x1`, `y1` are coordinates of top-left corner, and
`x2` and `y2` are coordinates of bottom-right corner.
"""
function IoU(bb1, bb2)
    @assert bb1["x1"] < bb1["x2"]
    @assert bb1["y2"] < bb1["y1"]
    @assert bb2["x1"] < bb2["x2"]
    @assert bb2["y2"] < bb2["y1"]

    x_left = max(bb1["x1"], bb2["x1"])
    y_top = min(bb1["y1"], bb2["y1"])
    x_right = min(bb1["x2"], bb2["x2"])
    y_bottom = max(bb1["y2"], bb2["y2"])

    if x_right < x_left || y_bottom > y_top
        return 0.0
    end

    # compute the intersection area
    intersection_area = (x_right - x_left) * (y_top - y_bottom)

    # compute the area of both boxes
    bb1_area = (bb1["x2"] - bb1["x1"]) * (bb1["y1"] - bb1["y2"])
    bb2_area = (bb2["x2"] - bb2["x1"]) * (bb2["y1"] - bb2["y2"])

    # compute the intersection over union by taking the intersection
    # area and dividing it by the sum of prediction + ground-truth
    # areas - the interesection area
    iou = intersection_area / Float64(bb1_area + bb2_area - intersection_area)

    @assert iou >= 0.0
    @assert iou <= 1.0
    return iou
end