text:
#=
function eqn_deepcopy(eqn::AbstractSolutionData)

  error("Generic fallback eqn_deepcopy() called erroneously. eqn_deepcopy() should be extended for each physics module's corresponding eqn object.")

  return nothing

end
=#

"""
  ODLCommonTools.eqn_deepcopy

  This function performs a proper deepcopy (unlike julia's builtin deepcopy)
    on an equation object.
  It preserves reference topology (i.e. q & q_vec pointing to same array in DG schemes).

    Inputs:
      eqn
      mesh
      sbp
      opts

    Outputs:
      eqn_copy

    One reason for doing this is this case:
      a = rand(2,2)
      b = a
      a[3] = 8
      b[3] == 8
      this is because 'a[3] =' is actually setindex!

"""
# function eqn_deepcopy{Tmsh, Tsol, Tres}(mesh::AbstractMesh{Tmsh}, sbp, eqn::AbstractSolutionData{Tsol, Tres}, opts::Dict)
function eqn_deepcopy(mesh::AbstractMesh, sbp, eqn::AbstractSolutionData, opts::Dict)

  # The eqn object has over 100 fields, so it is necessary to write a better approach for
  #   copying than explicitly copying every named field
  # This is the second write of eqn_deepcopy; the first version explicitly copied each field.
  #   Was done for SimpleODE and Advection before Euler made it obvious that that was intractable.

  error("eqn_deepcopy fallback reached. A physics-specific method should be defined.")

  return nothing

end

"""
  This function copies the data from one array of arrays to another array
  of arrays with the same structure.  It recurses down to the deepest array
  (whose element type is not a subtype of AbstractArray) and calls copy! on it.

  **Inputs**

   * src: an AbstractArray (possibly an array of arrays)

  **Inputs/Outputs**

   * dest: an AbstractArray (must be similar to src)

  Note that this function will not work correctly if any array (including
  the nested ones) has an element type of Any.

  Aliasing: aliasing is allowed, however cycles in the recursion tree are not
            allowed.  This function does not attempt to handle this case, and
            it will likely result in a StackOverflow error.
"""
function copy_array_recursive!(dest::AbstractArray{Td}, src::AbstractArray{Ts}) where {Td <: AbstractArray, Ts <: AbstractArray}

  @assert length(dest) == length(src)

  for i=1:length(src)
    # if the element type of src[i] is an AbstractArray, call this method
    # again, otherwise call the other method
    copy_array_recursive!(dest[i], src[i])
  end

  return nothing
end

function copy_array_recursive!(dest::AbstractArray{Td},
                                src::AbstractArray{Ts}) where {Td, Ts}

  copy!(dest, src)

  return nothing
end

"""
  ODLCommonTools.eqn_deepcopy_fields

  This function is physics-agnostic and copies over all fields.
  It assumes the fields of the destination object are already initialized
  and (in the case of arrays) the right size.

    **Inputs**:

     * eqn

    **Outputs:**

     * eqn_copy

"""
function copy!(eqn_copy::AbstractSolutionData, eqn::AbstractSolutionData)
  # 2: copy over fields

  for fdnm in fieldnames(eqn)    # loop over first level fieldnames in eqn

    if fdnm == :file_dict
      setfield!(eqn_copy, fdnm, getfield(eqn, fdnm))
    end

    fdnm_type = typeof(getfield(eqn, fdnm))    # get the super type of the current field

    # ------- handle params
    # if this first level fieldname is of type ParamType; ex: eqn.params
    if issubtype(fdnm_type, AbstractParamType)
      params_copy = getfield(eqn_copy, fdnm)
      params = getfield(eqn, fdnm)

      # loop over eqn.params, eqn.params_conservative, or eqn.params_entropy
      # println(" is a subtype of AbstractParamType, fdnm: ", fdnm)

      for fdnm_lvl2 in fieldnames(params)      # loop over 2nd level fieldnames

        # get the super type of the current 2nd level field
        fdnm_lvl2_type = typeof(getfield(params, fdnm_lvl2))

        # if the 2nd level fieldname is of type Array; ex: eqn.params.q_vals
        if issubtype(fdnm_lvl2_type, AbstractArray)

          # this does not work: setfield!(getfield(eqn_copy, a), b , getfield(getfield(eqn, a),b))
          #   must use copy, or else changing eqn's value changes eqn_copy
          copy!(getfield(params_copy, fdnm_lvl2), getfield(params, fdnm_lvl2))

          # Note: this is assuming that there are no Arrays of Arrays inside
          #       an eqn.params (or params_entropy, etc)

        else            # if the 2nd level fieldname is not of type Array; ex: eqn.params.gamma

          # because copy is not defined for all non-array types, such as functions
          setfield!(params_copy, fdnm_lvl2, getfield(params, fdnm_lvl2))

        end

      end

    # -------- handle arrays
    # if this first level fieldname is of type Array; ex: eqn.q or eqn.q_face_send
    elseif issubtype(fdnm_type, AbstractArray)

      arr_copy = getfield(eqn_copy, fdnm)
      arr = getfield(eqn, fdnm)
      copy_array_recursive!(arr_copy, arr)
    else        # handle non-arrays

      # copy is not defined for many of these non-array types: use assignment
      setfield!(eqn_copy, fdnm, getfield(eqn, fdnm))

    end

  end     # end of loop over first level fieldnames

  return eqn_copy

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -19.58980225336104 <= q <= 16.332536463331582)

                     add_NL_constraint(m, :(1/(1 + exp(-(0.1097699546178803 + 0.5720816873107069*1/(1 + exp(-(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])))) + 0.11808627755560508*1/(1 + exp(-(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])))) + -0.12730937872071246*1/(1 + exp(-(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])))) + 0.6166642606025916*1/(1 + exp(-(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])))) + 0.45996993866072433*1/(1 + exp(-(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))))))) + 1/(1 + exp(-(-0.8288330741680832 + -0.5580861219922526*1/(1 + exp(-(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])))) + -0.8786358511296881*1/(1 + exp(-(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])))) + -0.6915444540383735*1/(1 + exp(-(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])))) + 0.5554450106282389*1/(1 + exp(-(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])))) + -0.20937059223327026*1/(1 + exp(-(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))))))) + 1/(1 + exp(-(0.13615042538554878 + 0.08502572267265318*1/(1 + exp(-(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])))) + -0.9980047609363689*1/(1 + exp(-(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])))) + -0.5678523898140919*1/(1 + exp(-(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])))) + 0.29618367003382495*1/(1 + exp(-(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])))) + -0.7256587322261439*1/(1 + exp(-(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))))))) + 1/(1 + exp(-(-0.16139551949919406 + -0.5886735107938796*1/(1 + exp(-(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])))) + -0.02630736114123211*1/(1 + exp(-(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])))) + -0.03441572289730743*1/(1 + exp(-(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])))) + 0.3935924617025486*1/(1 + exp(-(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])))) + 0.3425056157482871*1/(1 + exp(-(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))))))) + 1/(1 + exp(-(-0.5329713534130147 + -0.3183971019485696*1/(1 + exp(-(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])))) + 0.5847814320797045*1/(1 + exp(-(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])))) + -0.7330877632484647*1/(1 + exp(-(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])))) + -0.37731199331975374*1/(1 + exp(-(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])))) + -0.08817423127621415*1/(1 + exp(-(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

ext: jl

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# ======= ALL ROTATIONS ALLOWED START =======
#lx + ly + lz = 1
@constraints(model, begin
    [i=1:length(boxes)], (lpx_length_of_box_parallel_to_x[i] +
    lpy_length_of_box_parallel_to_y[i] +
    lpz_length_of_box_parallel_to_z[i] == 1)
end)

#wx + wy + wz = 1
@constraints(model, begin
    [i=1:length(boxes)], (wpx_width_of_box_parallel_to_x[i] +
    wpy_width_of_box_parallel_to_y[i] +
    wpz_width_of_box_parallel_to_z[i] == 1)
end)

#hx + hy + hz = 1
@constraints(model, begin
    [i=1:length(boxes)], (hpx_height_of_box_parallel_to_x[i] +
    hpy_height_of_box_parallel_to_y[i] +
    hpz_height_of_box_parallel_to_z[i] == 1)
end)

#lx + wx + hx = 1
@constraints(model, begin
    [i=1:length(boxes)], (lpx_length_of_box_parallel_to_x[i] +
    wpx_width_of_box_parallel_to_x[i] +
    hpx_height_of_box_parallel_to_x[i] == 1)
end)

#ly + wy + hy = 1
@constraints(model, begin
    [i=1:length(boxes)], (lpy_length_of_box_parallel_to_y[i] +
    wpy_width_of_box_parallel_to_y[i] +
    hpy_height_of_box_parallel_to_y[i] == 1)
end)

#lz + wz + hz = 1
@constraints(model, begin
    [i=1:length(boxes)], (lpz_length_of_box_parallel_to_z[i] +
    wpz_width_of_box_parallel_to_z[i] +
    hpz_height_of_box_parallel_to_z[i] == 1)
end)
# ======= ALL ROTATIONS ALLOWED END =======


text:
macro primitive(ex)
    @capture(shortdef(ex), (logpdf(argname_::name_, args__) = body_) | (logpdf(argname_::name_, args__) where {Ts__} = body_)) || error("PrimitiveError: defining new contextual primitive calls requires a logpdf definition for the call.")
    argname = gensym(argname)
    expr = quote
        $ex

        @inline function (ctx::Jaynes.SimulateContext)(call::typeof(trace), addr::Jaynes.Address, $argname::$name, args...)
            Jaynes.visit!(ctx.visited, addr)
            s = $argname(args...)
            Jaynes.add_choice!(ctx, addr, logpdf($argname, args..., s), s)
            return s
        end

        function (ctx::Jaynes.GenerateContext)(call::typeof(trace), addr::T, $argname::$name, args...) where {T <: Jaynes.Address, K}
            Jaynes.visit!(ctx.visited, addr)
            if Jaynes.haskey(ctx.target, addr)
                s = Jaynes.getindex(ctx.target, addr)
                score = logpdf($argname, args..., s)
                Jaynes.add_choice!(ctx, addr, score, s)
                Jaynes.increment!(ctx, score)
            else
                s = $argname(args...)
                score = logpdf($argname, args..., s)
                Jaynes.add_choice!(ctx, addr, score, s)
            end
            return s
        end

        @inline function (ctx::Jaynes.ProposeContext)(call::typeof(trace), addr::T, $argname::$name, args...) where {T <: Jaynes.Address, K}
            s = $argname(args...)
            score = logpdf($argname, args..., s)
            Jaynes.add_choice!(ctx, addr, score, s)
            Jaynes.increment!(ctx, score)
            return s
        end

        @inline function (ctx::Jaynes.RegenerateContext)(call::typeof(trace), addr::T, $argname::$name, args...) where {T <: Jaynes.Address, K}
            Jaynes.visit!(ctx, addr)
            in_prev_chm = Jaynes.has_value(Jaynes.get_trace(ctx.prev), addr)
            in_sel = Jaynes.haskey(ctx.target, addr)

            if in_prev_chm
                prev = Jaynes.get_sub(Jaynes.get_trace(ctx.prev), addr)
            end

            if in_sel && in_prev_chm
                ret = rand(d)
                Jaynes.set_sub!(ctx.discard, addr, prev)
            elseif in_prev_chm
                ret = prev.val
            else
                ret = rand(d)
            end

            score = logpdf(d, ret)
            in_prev_chm && Jaynes.increment!(ctx, score - Jaynes.get_score(prev))
            Jaynes.add_choice!(ctx, addr, score, ret)
            return ret
        end

        @inline function (ctx::Jaynes.UpdateContext)(call::typeof(trace),
                                                     addr::T,
                                                     $argname::$name,
                                                     args...) where {T <: Jaynes.Address, K}
            in_prev_chm = Jaynes.haskey(ctx.prev, addr)
            in_prev_chm && begin
                prev = Jaynes.getindex(ctx.prev, addr)
                prev_ret = prev.val
                prev_score = prev.score
            end
            in_target = Jaynes.haskey(ctx.target, addr)
            if in_target
                ret = Jaynes.getindex(ctx.target, addr)
                in_prev_chm && begin
                    Jaynes.add_choice!(ctx.discard, addr, prev)
                end
                Jaynes.visit!(ctx.visited, addr)
            elseif in_prev_chm
                ret = prev_ret
            else
                ret = $argname(args...)
            end
            score = logpdf($argname, args..., ret)
            if in_prev_chm
                Jaynes.increment!(ctx, score - prev_score)
            elseif in_target
                Jaynes.increment!(ctx, score)
            end
            Jaynes.add_choice!(ctx, addr, score, ret)
            return ret
        end

        @inline function (ctx::Jaynes.AssessContext)(call::typeof(trace),
                                                     addr::T,
                                                     $argname::$name,
                                                     args...) where {T <: Jaynes.Address, K}
            Jaynes.haskey(ctx.target, addr) || error("AssessError: constrained target must provide constraints for all possible addresses in trace. Missing at address $addr.")
            val = Jaynes.getindex(ctx.target, addr)
            Jaynes.increment!(ctx, logpdf($argname, args..., val))
            return val

        end

        @inline function (ctx::Jaynes.ParameterBackpropagateContext)(call::typeof(trace),
                                                                     addr::T,
                                                                     $argname::$name,
                                                                     args...) where {T <: Jaynes.Address, K}
            if Jaynes.haskey(ctx.fixed, addr)
                s = Jaynes.getindex(ctx.fixed, addr)
            else
                s = Jaynes.get_value(Jaynes.get_sub(ctx.call, addr))
            end
            Jaynes.increment!(ctx, logpdf($argname, args..., s))
            return s
        end

        @inline function (ctx::Jaynes.ChoiceBackpropagateContext)(call::typeof(trace),
                                                                  addr::T,
                                                                  $argname::$name,
                                                                  args...) where {T <: Jaynes.Address, K}
            Jaynes.haskey(ctx.target, addr) || return Jaynes.get_value(Jaynes.get_sub(ctx.call, addr))
            s = Jaynes.getindex(ctx.call, addr)
            Jaynes.increment!(ctx, logpdf($argname, args..., s))
            return s
        end
    end
    expr = MacroTools.prewalk(unblock ∘ rmlines, expr)
    esc(expr)
end


text:
export iterate, length, haskey, keys, values, getindex

## Overloads for transparent iteration on `Exponent`, `Polyomial`, `Point`

## Exponent iterator
iterate(expo::Exponent) = iterate(expo.expo)
iterate(expo::Exponent, state) = iterate(expo.expo, state)
length(expo::Exponent) = length(expo.expo)
haskey(expo::Exponent, key) = haskey(expo.expo, key)
keys(expo::Exponent) = keys(expo.expo)
values(expo::Exponent) = values(expo.expo)
getindex(expo::Exponent, var::Variable) = expo.expo[var]

## Polynomial
iterate(poly::Polynomial) = iterate(poly.poly)
iterate(poly::Polynomial, state) = iterate(poly.poly, state)
length(poly::Polynomial) = length(poly.poly)
haskey(poly::Polynomial, key) = haskey(poly.poly, key)
keys(poly::Polynomial) = keys(poly.poly)
values(poly::Polynomial) = values(poly.poly)
getindex(poly::Polynomial, expo::Exponent) = poly.poly[expo]

## Point
iterate(pt::Point) = iterate(pt.coords)
iterate(pt::Point, state) = iterate(pt.coords, state)
length(pt::Point) = length(pt.coords)
haskey(pt::Point, key) = haskey(pt.coords, key)
keys(pt::Point) = keys(pt.coords)
values(pt::Point) = values(pt.coords)
getindex(pt::Point, var::Variable) = pt.coords[var]


text:
######################################################################################################################
# This file does a two-stage optimization of NICE (with and without revenue recycling) and saves the results. Users
# can change some of the core model parameters and optimization settings.
######################################################################################################################

# Load required Julia packages.
using NLopt

# Load RICE+AIR source code.
include("MimiNICE_recycle_time_varying.jl")

# ------------------------------------------------------------------------------------------------
# NICE + REVENUE RECYCLING PARAMETERS TO CHANGE
# ------------------------------------------------------------------------------------------------

# Pure rate of time preference.
ρ = 0.015

# Elasticity of marginal utility of consumption.
η = 1.5

# Income elasticity of climate damages (1 = proportional to income, -1 = inversely proportional to income).
damage_elasticity = 1.0

# Share of carbon tax revenue that is lost and cannot be recycled (1 = 100% of revenue lost, 0 = nothing lost)
lost_revenue_share = 0.0

# Shares of regional revenues that are recycled globally as international transfers (1 = 100% of revenue recycled globally).
global_recycle_share = zeros(12)

# Share of recycled carbon tax revenue that each region-quintile pair receives (row = region, column = quintile)
quintile_recycle_share = ones(12,5) .* 0.2

# Should the time-varying elasticity values only change across the range of GDP values from the studies?
# true = limit calculations to study gdp range, false = allow calculations for 0 to +Inf GDP.
bound_gdp_elasticity = false

# Quintile income distribution scenario (options = "constant", "SSP1", "SSP2", "SSP3", "SSP4", or "SSP5")
quintile_income_scenario = "constant"

# Type of slope for regression analysis (options are :central, :steeper, :flatter, :percentile)
regression_slope_type = :central

# Value if setting the regression slope type to a specific elasticity percentile value.
elasticity_percentile = 0.9

# ------------------------------------------------------------------------------------------------
# CHOICES ABOUT YOUR ANALYSIS & OPTIMZATION
# ------------------------------------------------------------------------------------------------

# Name of folder to store your results in (a folder will be created with this name).
results_folder = "My Results"

# Do you also want to perform a reference case optimization run (with no revenue recycling)?
run_reference_case = true

# Number of 10-year timesteps to find optimal carbon tax for (after which model assumes full decarbonization).
n_objectives = 45

# Global optimization algorithm (:symbol) used for initial result. See options at http://ab-initio.mit.edu/wiki/index.php/NLopt_Algorithms
global_opt_algorithm = :GN_DIRECT_L

# Local optimization algorithm (:symbol) used to "polish" global optimum solution.
local_opt_algorithm = :LN_SBPLX

# Maximum time in seconds to run global optimization (in case optimization does not converge).
global_stop_time = 600

# Maximum time in seconds to run local optimization (in case optimization does not converge).
local_stop_time = 300

# Relative tolerance criteria for global optimization convergence (will stop if |Δf| / |f| < tolerance from one iteration to the next.)
global_tolerance = 1e-8

# Relative tolerance criteria for global optimization convergence (will stop if |Δf| / |f| < tolerance from one iteration to the next.)
local_tolerance = 1e-12

# ------------------------------------------------------------------------------------------------
# ------------------------------------------------------------------------------------------------
# RUN EVERYTHING & SAVE KEY RESULTS
# ------------------------------------------------------------------------------------------------
# ------------------------------------------------------------------------------------------------

# Create a folder to save results.
output_directory = joinpath((@__DIR__), "..", "results", results_folder)
mkpath(output_directory)

# Load quintile income distribution scenario data.
income_distribution_raw = DataFrame(load(joinpath(@__DIR__, "..", "data", quintile_income_scenario*"_quintile_distributions_consumption.csv")))

# Clean up and organize time-varying income distribution data into required NICE format (time × regions × quintiles).
income_distributions = get_quintile_income_shares(income_distribution_raw)

#---------------------------------------------------------------------------------------------------
# Create a baseline version of the model without CO₂ mitigation.
#---------------------------------------------------------------------------------------------------

# This includes the user-specifications but has no CO₂ mitigation policy (will be used to calculte global CO₂ policy).
bau_model = create_nice_recycle(slope_type=regression_slope_type, percentile=elasticity_percentile)
n_steps   = length(dim_keys(bau_model, :time))

# Run model and extract some generic values needed to set up optimizations below.
run(bau_model)
rice_backstop_prices  = bau_model[:emissions, :pbacktime] .* 1000

# Now set user-specified parameter settings.
update_param!(bau_model, :MIU, zeros(n_steps, 12))
update_param!(bau_model, :damage_elasticity, damage_elasticity)
update_param!(bau_model, :quintile_income_shares, income_distributions)
update_param!(bau_model, :recycle_share, quintile_recycle_share)
update_param!(bau_model, :lost_revenue_share, lost_revenue_share)
update_param!(bau_model, :global_recycle_share, global_recycle_share)
update_param!(bau_model, :global_carbon_tax, zeros(n_steps))
update_param!(bau_model, :rho, ρ)
update_param!(bau_model, :eta, η)

# If selected, allow elasticities to be calculated for all GDP values (not just those observed in studies).
if bound_gdp_elasticity == false
    update_param!(bau_model, :min_study_gdp, 1e-10)
    update_param!(bau_model, :max_study_gdp, +Inf)
end

run(bau_model)

##############################################################################################################
# Optimization Run for NICE + Revenue Recycling
##############################################################################################################

# Create an instance of NICE with revenue recycling.
nice_rev_recycle = create_nice_recycle(slope_type=regression_slope_type, percentile=elasticity_percentile)

# Set user-specified parameter settings.
update_param!(nice_rev_recycle, :damage_elasticity, damage_elasticity)
update_param!(nice_rev_recycle, :quintile_income_shares, income_distributions)
update_param!(nice_rev_recycle, :recycle_share, quintile_recycle_share)
update_param!(nice_rev_recycle, :lost_revenue_share, lost_revenue_share)
update_param!(nice_rev_recycle, :global_recycle_share, global_recycle_share)
update_param!(nice_rev_recycle, :rho, ρ)
update_param!(nice_rev_recycle, :eta, η)

# If selected, allow elasticities to be calculated for all GDP values (not just those observed in studies).
if bound_gdp_elasticity == false
    update_param!(nice_rev_recycle, :min_study_gdp, 1e-10)
    update_param!(nice_rev_recycle, :max_study_gdp, +Inf)
end

#---------------------------------------------------------------------
# Set up first-stage global optimization for NICE + Revenue Recycling.
#---------------------------------------------------------------------

println("Beginning global optimization for NICE + Revenue Recycling.")

# Create NICE + Revenue Recycling objective function with user-specified settings.
recycle_objective = construct_nice_recycle_objective(nice_rev_recycle, revenue_recycling = true)

# Create an NLopt optimization object.
recycle_opt_global = Opt(global_opt_algorithm, n_objectives)

# Set upper and lower bounds for global CO₂ tax (upper bound = maximum RICE backstop prices).
upper_bound = maximum(rice_backstop_prices, dims=2)[2:(n_objectives+1)]

lower_bounds!(recycle_opt_global, zeros(n_objectives))
upper_bounds!(recycle_opt_global, upper_bound)

# Set maximum run time.
maxtime!(recycle_opt_global, global_stop_time)

# Set convergence tolerance.
ftol_rel!(recycle_opt_global, global_tolerance)

# Set objective function.
max_objective!(recycle_opt_global, (x, grad) -> recycle_objective(x))

# Optimize model using global optimization algorithm, initialized at a fraction of the upper bound.
recycle_max_welfare_global, recycle_opt_tax_global, recycle_convergence_global = optimize(recycle_opt_global, upper_bound ./ 100)

println("Convergence result = ", recycle_convergence_global)

# Calculate full tax time series (just to have to cross-check with local solution).
recycle_full_global_opt_tax, recycle_global_opt_regional_mitigation = mu_from_tax(recycle_opt_tax_global, rice_backstop_prices)

#----------------------------------------------------------------------
# Set up second-stage local optimization for NICE + Revenue Recycling.
#----------------------------------------------------------------------

println("Beginning local optimization for NICE + Revenue Recycling.")

# Create an NLopt optimization object for second-stage local optimization.
recycle_opt = Opt(local_opt_algorithm, n_objectives)

# Set upper and lower bounds for global CO₂ tax (upper bound = maximum RICE backstop prices).
lower_bounds!(recycle_opt, zeros(n_objectives))
upper_bounds!(recycle_opt, upper_bound)

# Set maximum run time.
maxtime!(recycle_opt, local_stop_time)

# Set convergence tolerance.
ftol_rel!(recycle_opt, local_tolerance)

# Set objective function.
max_objective!(recycle_opt, (x, grad) -> recycle_objective(x))

# Optimize model using local optimization algorithm, initialized at global optimization solution.
recycle_max_welfare, recycle_opt_tax, recycle_convergence = optimize(recycle_opt, recycle_opt_tax_global)

println("Convergence result = ", recycle_convergence)

#---------------------------------------------------
# Run Model Under Optimal Policy and Save Results
#---------------------------------------------------

# Get optimal regional CO₂ mitigation rates from final optimal tax.
recycle_full_opt_tax, recycle_opt_regional_mitigation = mu_from_tax(recycle_opt_tax, rice_backstop_prices)

# Run model under optimal policy.
update_param!(nice_rev_recycle, :MIU, recycle_opt_regional_mitigation)
update_param!(nice_rev_recycle, :global_carbon_tax, recycle_full_opt_tax)
run(nice_rev_recycle)

# Save results.
save_nice_recycle_results(nice_rev_recycle, bau_model, recycle_full_opt_tax, recycle_full_global_opt_tax, output_directory, revenue_recycling=true)

##############################################################################################################
# Optimization Run for NICE Without Revenue Recycling (Reference Case)
##############################################################################################################

if run_reference_case == true

    # Create an instance of NICE to run without revenue recycling.
    nice_reference = create_nice_recycle(slope_type=regression_slope_type, percentile=elasticity_percentile)

    # Set user-specified parameter settings.
    update_param!(nice_reference, :damage_elasticity, damage_elasticity)
    update_param!(nice_reference, :quintile_income_shares, income_distributions)
    update_param!(nice_reference, :recycle_share, quintile_recycle_share)
    update_param!(nice_reference, :lost_revenue_share, lost_revenue_share)
    update_param!(nice_reference, :global_recycle_share, global_recycle_share)
    update_param!(nice_reference, :rho, ρ)
    update_param!(nice_reference, :eta, η)

    # If selected, allow elasticities to be calculated for all GDP values (not just those observed in studies).
    if bound_gdp_elasticity == false
        update_param!(nice_reference, :min_study_gdp, 1e-10)
        update_param!(nice_reference, :max_study_gdp, +Inf)
    end

    #-------------------------------------------------------------------------------------
    # Set up first-stage global optimization for reference case without revenue recycling.
    #-------------------------------------------------------------------------------------

    println("Beginning global optimization for reference case without revenue recycling.")

    # Create NICE objective function with user parameter settings and revenue recycling switched off.
    reference_objective = construct_nice_recycle_objective(nice_reference, revenue_recycling = false)

    # Create an NLopt optimization object.
    reference_opt_global = Opt(global_opt_algorithm, n_objectives)

    # Set upper and lower bounds for global CO₂ tax (upper bound = maximum RICE backstop prices).
    lower_bounds!(reference_opt_global, zeros(n_objectives))
    upper_bounds!(reference_opt_global, upper_bound)

    # Set max run time.
    maxtime!(reference_opt_global, global_stop_time)

    # Set convergence tolerance.
    ftol_rel!(reference_opt_global, global_tolerance)

    # Set objective function.
    max_objective!(reference_opt_global, (x, grad) -> reference_objective(x))

    # Optimize model using global optimization algorithm, initialized at a fraction of the upper bound.
    reference_max_welfare_global, reference_opt_tax_global, reference_convergence_global = optimize(reference_opt_global, upper_bound ./ 100)

    # Calculate full tax time series (just to have to cross-check with local solution).
    reference_full_global_opt_tax, reference_global_opt_regional_mitigation = mu_from_tax(reference_opt_tax_global, rice_backstop_prices)

    println("Convergence result = ", reference_convergence_global)

    #--------------------------------------------------------------------------------------
    # Set up second-stage local optimization for reference case without revenue recycling.
    #--------------------------------------------------------------------------------------

    println("Beginning local optimization for reference case without revenue recycling.")

    # Create an NLopt optimization object for second-stage local optimization.
    reference_opt = Opt(local_opt_algorithm, n_objectives)

    # Set upper and lower bounds for global CO₂ tax (upper bound = maximum RICE backstop prices).
    lower_bounds!(reference_opt, zeros(n_objectives))
    upper_bounds!(reference_opt, upper_bound)

    # Set max run time.
    maxtime!(reference_opt, local_stop_time)

    # Set convergence tolerance.
    ftol_rel!(reference_opt, local_tolerance)

    # Set objective function.
    max_objective!(reference_opt, (x, grad) -> reference_objective(x))

    # Optimize model using local optimization algorithm, initialized at global optimization solution.
    reference_max_welfare, reference_opt_tax, reference_convergence = optimize(reference_opt, reference_opt_tax_global)

    println("Convergence result = ", reference_convergence)

    #-----------------------------------------------------------
    # Run Reference Model Under Optimal Policy and Save Results
    #-----------------------------------------------------------

    # Get optimal regional CO₂ mitigation rates from final optimal tax.
    reference_full_opt_tax, reference_opt_regional_mitigation = mu_from_tax(reference_opt_tax, rice_backstop_prices)

    # Run model under optimal policy.
    update_param!(nice_reference, :MIU, reference_opt_regional_mitigation)
    update_param!(nice_reference, :global_carbon_tax, zeros(n_steps))
    run(nice_reference)

    # Save results.
    save_nice_recycle_results(nice_reference, bau_model, reference_full_opt_tax, reference_full_global_opt_tax, output_directory, revenue_recycling=false)
end

println("Model runs complete.")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule helpme_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("helpme")
JLLWrappers.@generate_main_file("helpme", UUID("f33c8c47-f3e6-5e52-bf5f-5eacb4024144"))
end  # module helpme_jll


text:
function _random_normal(shape...)
  return u_type.(rand(Normal(0.f0,0.02f0),shape...))|> gpu
end

function extract_bboxes(mask)
  nth = last(size(mask))
  boxes = zeros(Integer, nth, 4)
  for i =1:nth
    m = mask[:,:,i]
    cluster = findall(!iszero, m)
    if length(cluster) > 0
      Is = map(x -> [x.I[1], x.I[2]], cluster) |> x -> hcat(x...)'
      x1, x2 = extrema(Is[:,1])
      y1, y2 = extrema(Is[:,2])
    else
      x1 ,x2, y1, y2 = 0, 0, 0, 0
    end
      boxes[i,:] = [y1, x1, y2, x2]
  end
  boxes
end

function extract_bboxes(masks::AbstractArray{T,4}) where T
  bs = []
  for i in 1:size(masks, 4)
    b = extract_bboxes(masks[:,:,:,i])
    push!(bs, b)
  end
  reduce(vcat, bs)
end

expand_dims(x,n::Int) = reshape(x,ones(Int32,n)...,size(x)...)
function squeeze(x)
    if size(x)[end] != 1
        return dropdims(x, dims = tuple(findall(size(x) .== 1)...))
    else
        # For the case BATCH_SIZE = 1
        int_val = dropdims(x, dims = tuple(findall(size(x) .== 1)...))
        return reshape(int_val,size(int_val)...,1)
    end
end

function bce(ŷ, y; ϵ=gpu(fill(eps(first(ŷ)), size(ŷ)...)))
  l1 = -y.*log.(ŷ .+ ϵ)
  l2 = (1 .- y).*log.(1 .- ŷ .+ ϵ)
  l1 .- l2
end

function loss_bce(x, y)
  op = clamp.(model(x), 0.001f0, 1.f0)
  mean(bce(op, y))
end

function complex_grid_to_channels!(grid)
    grid_channels = CuArray{r_type}(undef, size(grid,1), size(grid,2), 2, 1)
    grid_channels[:, :, 1, :] = real(grid)
    grid_channels[:, :, 2, :] = imag(grid)
    return grid_channels
end

function complex_helmholtz_to_channels!(helmholtz_matrix)
    helmholtz_channels = CuArray{r_type}(undef, size(grid,1), size(grid,2), 2, 2)
    helmholtz_channels[:,:,1,1] = real(helmholtz_matrix)
    helmholtz_channels[:,:,2,1] = -imag(helmholtz_matrix)
    helmholtz_channels[:,:,1,2] = imag(helmholtz_matrix)
    helmholtz_channels[:,:,2,2] = real(helmholtz_matrix)
    return helmholtz_channels
end

function check_helmholtz_channels!(helmholtz_matrix, x, n)
    helmholtz_channels = complex_helmholtz_to_channels!(helmholtz_matrix, n)|> gpu
    x_channels = complex_grid_to_channels!(x, n)|> gpu

    original_result = helmholtz_chain!(x, helmholtz_matrix; h=h)|> gpu
    channels_result = helmholtz_chain_channels!(x_channels, helmholtz_channels, n; h=h)
    channels_result = channels_result[:,:,1,:]+im*channels_result[:,:,2,:]

    @info "$(Dates.format(now(), "HH:MM:SS")) - Check Helmholtz Channels $(norm(channels_result .- original_result)), $(original_result[1,1,1,1]), $(channels_result[1,1,1,1])"
    println(norm(channels_result - original_result))
end


text:
# define initial simplex
abstract type AbstractSimplexer end
struct SSH{T} <: AbstractSimplexer
    L::T
end
function initial_simplex(is::SSH, x0)
    T = eltype(x0)
    L = T(is.L)
    n = length(x0)
    p = (n - 1 + sqrt(n + 1))/(n*sqrt(2))
    q = (sqrt(n + 1) - 1)/(n*sqrt(2))

    S = [copy(x0)]
    for i = 1:n
        vₙ = [j == i ?  x0[j] + L*p : x0[j] + L*q for j = 1:n]
        push!(S, vₙ)
    end
    S
end

struct ABA{T} <: AbstractSimplexer
    L::T
end
function initial_simplex(is::ABA{T}, x0) where T<:Real
    Tx = eltype(x0)
    iss = ABA(x0*Tx(0).+Tx(is.L))
    initial_simplex(iss, x0)
end
function initial_simplex(is::ABA, x0)
    T = eltype(x0)
    L = T.(is.L)
    n = length(x0)

    S = [copy(x0)]
    for i = 1:n
        vₙ = x0 .+ ((j == i)*L[j] for j = 1:n)
        push!(S, vₙ)
    end
    S
end

struct RB <: AbstractSimplexer
    lb
    ub
end
function simplex_rb(x0)
end
#
# struct LPS{A, R} <: AbstractSimplexer
#     abs::A
#     rel::R
# end
# function initial_simplex(is::LPS, x0)
#     T = eltype(x0)
#     small = abs.(x0) .<= sqrt(eps(T))
#
#     ...
# end

nmparameters(t) = t
struct GaoHan{Tx}
    n::Int
    T::Tx
end
function (gh::GaoHan)(t)
    n = gh.n
    T = gh.T
    nT = T(n)
    t .+ (0, 2/nT, inv(2*nT), 1/nT)
end

struct NelderMead{Ts, Tp, Tsi}
    initial_simplex::Ts
    parameters::Tp
    shrink_index::Tsi
end

"""
# NelderMead
## Constructor
```julia
NelderMead(; parameters = nmparameters,
             initial_simplex = AffineSimplex())
```

The constructor takes 2 keywords:

* `parameters`, an instance of either `AdaptiveParameters` or `FixedParameters`,
and is used to generate parameters for the Nelder-Mead Algorithm
* `initial_simplex`, an instance of `AffineSimplexer`

## Description
Our current implementation of the Nelder-Mead algorithm is based on [1] and [3].
Gradient-free methods can be a bit sensitive to starting values and tuning parameters,
so it is a good idea to be careful with the defaults provided in Optim.jl.

Instead of using gradient information, Nelder-Mead is a direct search method. It keeps
track of the function value at a number of points in the search space. Together, the
points form a simplex. Given a simplex, we can perform one of four actions: reflect,
expand, contract, or shrink. Basically, the goal is to iteratively replace the worst
point with a better point. More information can be found in [1], [2] or [3].

## References
- [1] Nelder, John A. and R. Mead (1965). "A simplex method for function minimization". Computer Journal 7: 308–313. doi:10.1093/comjnl/7.4.308
- [2] Lagarias, Jeffrey C., et al. "Convergence properties of the Nelder–Mead simplex method in low dimensions." SIAM Journal on Optimization 9.1 (1998): 112-147
- [3] Gao, Fuchang and Lixing Han (2010). "Implementing the Nelder-Mead simplex algorithm with adaptive parameters". Computational Optimization and Applications. doi:10.1007/s10589-010-9329-3
"""
function NelderMead()
    return NelderMead(0, x->(1.0, 2.0, 0.5, 0.5), nothing)
end
Base.summary(::NelderMead) = "Nelder-Mead"
function solve(prob::OptimizationProblem, x0, method::NelderMead, options::OptimizationOptions)
    solve(mstyle(prob), prob, x0, method, ABA(x0.*0 .+ 1), options)
end

# centroid except h-th vertex
function centroid!(c, simplex_vector, h=0)
    T = eltype(c)
    n = length(c)
    fill!(c, zero(T))
    @inbounds for i in 1:n+1
        if i != h
            xi = simplex_vector[i]
            c .+= xi
        end
    end
    c .= c ./ n
end

centroid(simplex_vector, h) = centroid!(similar(simplex_vector[1]), simplex_vector, h)
function nmobjective(y::Vector)
    a = sqrt(var(y) * (length(y) / (length(y) - 1)))
    return a
end
struct ValuedSimplex{TS, TV, TO}
    S::TS
    V::TV
    O::TO
end
ValuedSimplex(S, V) = ValuedSimplex(S, V, sortperm(V))

function solve(mstyle::InPlace, prob::OptimizationProblem, x0, method::NelderMead, as::AbstractSimplexer, options)
    simplex_vector = initial_simplex(as, x0)
    simplex_value = batched_value(prob, simplex_vector) # this could be batched
    order = sortperm(simplex_value)
    simplex = ValuedSimplex(simplex_vector, simplex_value, order)
    res = solve(mstyle, prob, simplex, method, options)
    x0 .= res.info.minimizer
    return res
end

function NMCaches(simplex)
    x_reflect = copy(first(simplex.S))
    x_cache = copy(first(simplex.S))
    x_centroid = copy(x_cache)
    return (x_reflect=x_reflect, x_cache=x_cache, x_centroid=x_centroid)
end

function solve(mstyle::InPlace, prob::OptimizationProblem, simplex::ValuedSimplex, method::NelderMead, options::OptimizationOptions, nmcache=NMCaches(simplex))
    t0 = time()
    simplex_vector, simplex_value, i_order = simplex.S, simplex.V, simplex.O
    f0 = minimum(simplex.V)
    n = length(first(simplex_vector))
    m = length(simplex_vector)

    # We only need two caches
    x_reflect = nmcache.x_reflect
    x_cache = nmcache.x_cache

    x_centroid = centroid!(nmcache.x_centroid, simplex_vector, i_order[end])

    α, β, γ, δ = method.parameters(n)
    step_type = "none"
    iter = 0
    nm_obj = f0
    is_converged = false
    while iter <= options.maxiter && !any(is_converged)
        iter += 1
        nm_obj, x_centroid = iterate!(prob, method, simplex_vector, simplex_value, i_order, x_cache, x_centroid, x_reflect, α, β, γ, δ)
        print_trace(method, options, iter, t0, simplex_value)
        if nm_obj ≤ options.nm_tol
            is_converged = true
        end
    end
    f_centroid_min = value(prob, x_centroid)
    f_min, i_f_min = findmin(simplex_value)
    x_min = simplex_vector[i_f_min]
    if f_centroid_min < f_min
        x_min = x_centroid
        f_min = f_centroid_min
    end
    ConvergenceInfo(method, (nm_obj=nm_obj, centroid=x_centroid, simplex=simplex, minimizer=x_min, minimum=f_min, f0=f0, iter=iter, time=time()-t0), options)
end
function print_trace(::NelderMead, options, iter, t0, simplex_value)
end

function iterate!(prob, method::NelderMead, simplex_vector, simplex_value, i_order, x_cache, x_centroid, x_reflect, α, β, γ, δ)
    # Augment the iteration counter
    shrink = false
    n = length(first(simplex_vector))
    m = length(simplex_vector)

    # Compute a reflection
    x_highest = simplex_vector[i_order[end]]
    @. x_reflect .= x_centroid + α * (x_centroid - x_highest)

    f_reflect = value(prob, x_reflect)
    if f_reflect < simplex_value[i_order[1]] # f_lowest(simplex)
        # Compute an expansion
        x_expand = x_cache
        @. x_expand .= x_centroid + β * (x_reflect - x_centroid)

        f_expand = value(prob, x_expand)

        if f_expand < f_reflect
            simplex_vector[i_order[end]] .= x_expand
            simplex_value[i_order[end]] = f_expand
            step_type = "expansion"
        else
            simplex_vector[i_order[end]] .= x_reflect
            simplex_value[i_order[end]] = f_reflect
            step_type = "reflection"
        end
    elseif f_reflect < simplex_value[i_order[end - 1]]  # f_second_highest(simplex)
        simplex_vector[i_order[end]] .= x_reflect
        simplex_value[i_order[end]] = f_reflect
        step_type = "reflection"
    else
        if f_reflect < simplex_value[i_order[end]] # f_highest(simplex)
            # Outside contraction
            x_contract = x_cache
            @. x_contract .= x_centroid + γ * (x_reflect-x_centroid)

            f_outside_contraction = value(prob, x_contract)
            if f_outside_contraction < f_reflect
                simplex_vector[i_order[end]] .= x_contract
                simplex_value[i_order[end]] = f_outside_contraction
                step_type = "outside contraction"
            else
                shrink = true
            end
        else # f_reflect > f_highest
            # Inside contraction
            x_inside_contract = x_cache
            @. x_inside_contract .= x_centroid - γ *(x_reflect - x_centroid)
            f_inside_contraction = value(prob, x_inside_contract)

            if f_inside_contraction < simplex_value[i_order[end]] #  f_highest(simplex)
                simplex_vector[i_order[end]] .= x_inside_contract
                simplex_value[i_order[end]] = f_inside_contraction
                step_type = "inside contraction"
            else
                shrink = true
            end
        end
    end

    if shrink
        if isa(method.shrink_index, Nothing)
            # shrink all vertices indexed 2 and up
            low_range_index = 2
        elseif isa(method.shrink_index, Number)
            # only shrink the vertices indexed shrink_index and up
            low_range_index = method.shrink_index
        else
            # shrink the indeces returned by shrink_index ( could be random)
            low_range_index = method.shrink_index(m)
        end
        shrink_range = low_range_index:m
        for i = shrink_range
            ord = i_order[i]
            x_lowest = simplex_vector[i_order[1]]
            @. simplex_vector[ord] .= x_lowest + δ*(simplex_vector[ord] - x_lowest)
            simplex_value[ord] = value(prob, simplex_vector[ord])
        end
        step_type = "shrink"
    end

    sortperm!(i_order, simplex_value)

    x_centroid = centroid!(x_centroid, simplex_vector, i_order[end])
    nm_obj = nmobjective(simplex_value)

    nm_obj, x_centroid
end

#####################
#    out-of-place   #
#####################
function solve(mstyle::OutOfPlace, prob::OptimizationProblem, x0, method::NelderMead, as::AbstractSimplexer, options::OptimizationOptions)
    simplex_vector = initial_simplex(as, copy(x0))
    simplex_value = batched_value(prob, simplex_vector) # this could be batched
    order = sortperm(simplex_value)
    simplex = ValuedSimplex(simplex_vector, simplex_value, order)
    res = solve(mstyle, prob, simplex, method, options)
    return res
end
function solve(mstyle::OutOfPlace, prob::OptimizationProblem, simplex::ValuedSimplex, method::NelderMead, options::OptimizationOptions)
    t0 = time()
    simplex_vector, simplex_value = simplex.S, simplex.V
    n = length(first(simplex_vector))
    m = length(simplex_vector)
    f0 = minimum(simplex.V)
    # Get the indices that correspond to the ordering of the f values
    # at the vertices. i_order[1] is the index in the simplex of the vertex
    # with the lowest function value, and i_order[end] is the index in the
    # simplex of the vertex with the highest function value
    i_order = sortperm(simplex_value)

    x_centroid = centroid(simplex_vector, i_order[end])

    α, β, γ, δ = method.parameters(n)
    step_type = "none"
    is_converged = false
    iter = 0
    nm_obj = f0

    while iter <= options.maxiter && !any(is_converged)
        iter += 1

        # Augment the iteration counter
        shrink = false

        # Compute a reflection
        x_highest = simplex_vector[i_order[end]]
        x_reflect = @. x_centroid + α * (x_centroid - x_highest)

        f_reflect = value(prob, x_reflect)
        if f_reflect < simplex_value[i_order[1]] # f_lowest(simplex)
            # Compute an expansion
            x_expand = @. x_centroid + β * (x_reflect - x_centroid)

            f_expand = value(prob, x_expand)

            if f_expand < f_reflect
                simplex_vector[i_order[end]] = x_expand
                simplex_value[i_order[end]] = f_expand
                step_type = "expansion"
            else
                simplex_vector[i_order[end]] = x_reflect
                simplex_value[i_order[end]] = f_reflect
                step_type = "reflection"
            end
        elseif f_reflect < simplex_value[i_order[end - 1]]  # f_second_highest(simplex)
            simplex_vector[i_order[end]] = x_reflect
            simplex_value[i_order[end]] = f_reflect
            step_type = "reflection"
        else
            if f_reflect < simplex_value[i_order[end]] # f_highest(simplex)
                # Outside contraction
                x_contract = @. x_centroid + γ * (x_reflect-x_centroid)

                f_outside_contraction = value(prob, x_contract)
                if f_outside_contraction < f_reflect
                    simplex_vector[i_order[end]] = x_contract
                    simplex_value[i_order[end]] = f_outside_contraction
                    step_type = "outside contraction"
                else
                    shrink = true
                end
            else # f_reflect > f_highest
                # Inside contraction
                x_inside_contract = @. x_centroid - γ *(x_reflect - x_centroid)

                f_inside_contraction = value(prob, x_inside_contract)
                if f_inside_contraction < simplex_value[i_order[end]] #  f_highest(simplex)
                    simplex_vector[i_order[end]] = x_inside_contract
                    simplex_value[i_order[end]] = f_inside_contraction
                    step_type = "inside contraction"
                else
                    shrink = true
                end
            end
        end

        if shrink
            # only shirk
            if isa(method.shrink_index, Nothing)
                low_range_index = 2
            elseif isa(method.shrink_index, Number)
                low_range_index = method.shrink_index
            else
                low_range_index = method.shrink_index(m)
            end
            shrink_range = low_range_index:m
            for i = shrink_range
                ord = i_order[i]
                x_lowest = simplex_vector[i_order[1]]
                simplex_vector[ord] = @. x_lowest + δ*(simplex_vector[ord] - x_lowest)
            end
            step_type = "shrink"
        end

        sortperm!(i_order, simplex_value)

        x_centroid = centroid(simplex_vector, i_order[end])

        nm_obj = nmobjective(simplex_value)
        # if nm_x < 1e-18
        #     break
        # end
        # check conv
    end
    f_centroid_min = value(prob, x_centroid)
    f_min, i_f_min = findmin(simplex_value)
    x_min = simplex_vector[i_f_min]
    if f_centroid_min < f_min
        x_min = x_centroid
        f_min = f_centroid_min
    end
    ConvergenceInfo(method, (nm_obj=nm_obj, centroid=x_centroid, simplex=simplex, minimizer=x_min, minimum=f_min, f0=f0, iter=iter, time=time()-t0), options)
end


text:
import NMRDataSetup
import NMRSpectraSimulator
import NMRSpecifyRegions

include("../src/NMRCalibrate.jl")
import .NMRCalibrate

using LinearAlgebra
using FFTW
import PyPlot
import BSON
import JSON

import NLopt
import Statistics
import Random

PyPlot.close("all")
fig_num = 1

Random.seed!(25)
PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])
#

# First, run prep_full.jl without optimization.

r = 1
T = Float64

y_cost = y[cost_inds_set[r]]
U_cost = U_y[cost_inds_set[r]]
P_cost = P_y[cost_inds_set[r]]
# LS_inds = 1:length(U_cost)

Es = collect( NMRSpectraSimulator.καFIDModelType(Bs[i]) for i = 1:length(Bs) )

Es[1].κs_α[1][:] = [1.0, 1.0, 1.0]
#Es[1].core.λ0 = 2.3752439850648
Es[1].core.ss_params.κs_λ[1] = 1.0
#Es[1].core.ss_params.κs_β[1] = [1.7457371427914214, 2.9647985946238675, 2.5466218556915283] # calibration answer.
#Es[1].core.ss_params.κs_β[1] = [1.5664070671305295, 2.839984010713529, 2.3915878840313702]

# 1.5432104948436969
#  2.8279246306255725
#  2.3711578754117753

#Es[1].core.ss_params.d[1] = 5.4534685586024585 # calibration answer?
#Es[1].core.ss_params.d[1] = 4.527727288274466

#tmp = (1.7457371427914214 + 2.9647985946238675)/2
#Es[1].core.ss_params.κs_β[1] = [tmp; tmp; 2.5466218556915283]
#Es[1].core.ss_params.κs_β[1][2] = 1.7457371427914214
Es[1].core.ss_params.κs_β[1][2] = Es[1].core.ss_params.κs_β[1][1]

# # originals.
# As[1].Δc_bar[1][1] = [-0.04550785321752465, -0.9527933956628394, -0.0016987511196359661]
# As[1].Δc_bar[1][2] = [-0.9542311847361354, -0.04525616042734937, -0.0005126548365153027]
# As[1].Δc_bar[1][3] = [-0.0002609620463399076, -0.0019504439098113058, -0.9977885940438487]

## modify/test.
#As[1].Δc_bar[1][1] = As[1].Δc_bar[1][2]

# http://gissmo.nmrfam.wisc.edu/entry/bmse000867/simulation_1
# http://gissmo.nmrfam.wisc.edu/entry/bmse000048/simulation_1

βs = collect( dot(Es[1].core.ss_params.κs_β[1], As[1].Δc_bar[1][k]) for k = 1:length(As[1].Δc_bar[1]))
vs = collect( Es[1].core.ss_params.κs_β[1] .* As[1].Δc_bar[1][k] for k = 1:length(As[1].Δc_bar[1]))

# for calibration, set w to all ones.
q2 = uu->NMRSpectraSimulator.evalitpproxymixture(uu, As, Es; w = ws[1])

# eval.
N_viz = 50000
U = LinRange(u_min, u_max, N_viz)
P = hz2ppmfunc.(U)
U_rad = U .* (2*π)

q_U = q2.(U_rad)

# viz.

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(P, real.(q_U), label = "fit")
PyPlot.plot(P_cost, real.(y_cost), label = "data")
PyPlot.plot(P_cost, real.(y_cost), "x")

PyPlot.gca().invert_xaxis()
PyPlot.legend()
PyPlot.xlabel("ppm")
PyPlot.ylabel("real")
PyPlot.title("data vs. bandpassed")


text:
module Utilities

using Test
using MLJBase
using Random
using ComputationalResources

@test MLJBase.finaltypes(Union{Missing,Int}) == [Union{Missing,Int64}]
@test MLJBase.finaltypes(Float64) == [Float64]

abstract type Foo end
struct Bar <: Foo end
struct Baz <: Foo end

@test MLJBase.finaltypes(Foo) == [Bar, Baz]

@testset "flat_values" begin
    t = (X = (x = 1, y = 2), Y = 3)
    @test flat_values(t) == (1, 2, 3)
end

mutable struct M
    a1
    a2
end
mutable struct A1
    a11
    a12
end
mutable struct A2
    a21
end
mutable struct A21
    a211
    a212
end

@testset "recursive getproperty, setproperty!" begin
    m = (a1 = (a11 = 10, a12 = 20), a2 = (a21 = (a211 = 30, a212 = 40),))

    @test MLJBase.recursive_getproperty(m, :(a1.a12)) == 20
    @test MLJBase.recursive_getproperty(m, :a1) == (a11 = 10, a12 = 20)
    @test MLJBase.recursive_getproperty(m, :(a2.a21.a212)) == 40

    m = M(A1(10, 20), A2(A21(30, 40)))
    MLJBase.recursive_setproperty!(m, :(a2.a21.a212), 42)
    @test MLJBase.recursive_getproperty(m, :(a1.a11)) == 10
    @test MLJBase.recursive_getproperty(m, :(a1.a12)) == 20
    @test MLJBase.recursive_getproperty(m, :(a2.a21.a211)) == 30
    @test MLJBase.recursive_getproperty(m, :(a2.a21.a212)) == 42
    @test MLJBase.recursive_getproperty(
        MLJBase.recursive_getproperty(m, :(a2.a21)), :a212) == 42
end

@testset "shuffle rows" begin
    # check dims
    x = randn(5)
    y = randn(5, 5)
    z = randn(5, 5)
    @test MLJBase.check_dimensions(x, y) === nothing
    @test MLJBase.check_dimensions(z, x) === nothing
    @test MLJBase.check_dimensions(y, z) === nothing
    @test_throws DimensionMismatch MLJBase.check_dimensions(x, randn(4))

    x = 1:5 |> collect
    y = 1:5 |> collect
    rng = 555
    Random.seed!(rng)
    perm = randperm(5)
    @test MLJBase.shuffle_rows(x, y; rng=rng) == (x[perm], y[perm])
    y = randn(5,5)
    @test MLJBase.shuffle_rows(x, y; rng=rng) == (x[perm], y[perm,:])
    @test MLJBase.shuffle_rows(z, y; rng=rng) == (z[perm,:], y[perm,:])
    @test MLJBase.shuffle_rows(x, x; rng=rng) == (x[perm], x[perm])
end

@testset "unwind" begin
    iterators = ([1, 2], ["a","b"], ["x", "y", "z"])
    @test unwind(iterators...) ==
        [1  "a"  "x";
         2  "a"  "x";
         1  "b"  "x";
         2  "b"  "x";
         1  "a"  "y";
         2  "a"  "y";
         1  "b"  "y";
         2  "b"  "y";
         1  "a"  "z";
         2  "a"  "z";
         1  "b"  "z";
         2  "b"  "z"]
end

@testset "comp resources" begin
    @test default_resource() == CPU1()
    default_resource(ComputationalResources.CPUProcesses())
    @test default_resource() == ComputationalResources.CPUProcesses()
    default_resource(CPU1())
end

end # module
true


text:
using DataStructures

function execute(src)
    pointers = Dict{Int,Int}()
    stack    = Int[]
    for (ptr, opcode) in enumerate(src)
        if opcode == '[' push!(stack, ptr) end
        if opcode == ']'
            if isempty(stack)
                src = src[1:ptr]
                break
            end
            sptr = pop!(stack)
            pointers[ptr], pointers[sptr] = sptr, ptr
        end
    end
    if ! isempty(stack) error("unclosed loops at $stack") end
    tape = DefaultDict{Int,Int}(0)
    cell, ptr = 0, 1
    while ptr ≤ length(src)
        opcode = src[ptr]
        if     opcode == '>' cell += 1
        elseif opcode == '<' cell -= 1
        elseif opcode == '+' tape[cell] += 1
        elseif opcode == '-' tape[cell] -= 1
        elseif opcode == ',' tape[cell] = Int(read(STDIN, 1))
        elseif opcode == '.' print(STDOUT, Char(tape[cell]))
        elseif (opcode == '[' && tape[cell] == 0) ||
               (opcode == ']' && tape[cell] != 0) ptr = pointers[ptr]
        end
        ptr += 1
    end
end

const src = """\
    >++++++++[<+++++++++>-]<.>>+>+>++>[-]+<[>[->+<<++++>]<<]>.+++++++..+++.>
    >+++++++.<<<[[-]<[-]>]<+++++++++++++++.>>.+++.------.--------.>>+.>++++."""
execute(src)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
f(x) = 1*x*x*x+2*x*x+3*x+4
f1(x) = 3*x*x+4*x+3
f2(x)= 6*x+4
function Horner(x,list)
	a = pop!(list)
	b = zero(x)
	b = a + x*b
	a = pop!(list) + x*a
	c = zero(x)
	for  i = 1:length(list)
		c = 2*b + x*c
		b = a + x*b
		a = pop!(list) + x*a
	end
	return (a,b,c)
end
@show test = Horner(2,[4,3,2,1])
@show f(2)
@show f1(2)
@show f2(2)

max_stars_repo_name: Magikis/Uniwersity

max_issues_repo_name: Magikis/Uniwersity

max_forks_repo_name: Magikis/Uniwersity


text:
# using Pkg
# pkg"activate /home/fredrikb/.julia/dev/BlobTracking/"
using BlobTracking
using Test, Statistics, ImageDraw, Images, VideoIO, StaticArrays
img = Gray.(zeros(100,100))
locs = [CartesianIndex(10,10), CartesianIndex(30,30)]
draw!(img,locs[1], c=Gray(1.0))
draw!(img,locs[2], c=Gray(1.0))

img2 = Gray.(zeros(100,100))
locs2 = [CartesianIndex(15,15), CartesianIndex(35,35)]
draw!(img2,locs2[1], c=Gray(1.0))
draw!(img2,locs2[2], c=Gray(1.0))

img3 = Gray.(zeros(100,100))
locs3 = [CartesianIndex(20,20), CartesianIndex(40,40)]
draw!(img3,locs3[1], c=Gray(1.0))
draw!(img3,locs3[2], c=Gray(1.0))

##
@testset "BlobTracking.jl" begin

    @testset "to_static" begin
        @info "Testing to_static"
        @test to_static(1) == 1
        @test to_static(randn(2,2)) isa SMatrix{2,2}
        @test to_static(randn(2)) isa SVector{2}

    end

    @testset "Blob" begin
        @info "Testing Blob"
        b = Blob()
        blobs = [b]
        corr = HungarianCorrespondence(dist_th = 2)
        measurement = BlobTracking.assign(corr, blobs, locs)
        BlobTracking.correct!(blobs, measurement)
        pos1 = to_static(location(b))
        @test all(to_static(location(b)) .> 0)
        @test all(to_static(location(b)) .< to_static(locs[1]))

        b = Blob()
        blobs = [b]
        corr = MCCorrespondence(HungarianCorrespondence(dist_th = 2), 20)
        measurement = BlobTracking.assign(corr, blobs, locs)
        BlobTracking.correct!(blobs, measurement)
        pos2 = to_static(location(b))
        @test pos1 == pos2

    end

    @testset "Correspondence" begin
        @info "Testing Correspondence"

        @testset "NearestNeighborCorrespondence" begin
            @info "Testing NearestNeighborCorrespondence"
            bt = BlobTracker(2:2, 10, 5.0)
            coordinates = locs
            blobs = Blob.(bt.params, coordinates)
            c = NearestNeighborCorrespondence(dist_th = 2)

            meas = BlobTracking.assign(c, blobs, coordinates)
            @test meas.assi == [1,2]

            meas = BlobTracking.assign(c, blobs, [CartesianIndex(1000,1000); coordinates])
            @test meas.assi == [2,3]

            meas = BlobTracking.assign(c, [blobs; blobs], coordinates)
            @test meas.assi == [1,2,1,2]

        end

        @testset "HungarianCorrespondence" begin
            @info "Testing HungarianCorrespondence"
            bt = BlobTracker(2:2, 10, 5.0)
            coordinates = locs
            blobs = Blob.(bt.params, coordinates)
            c = HungarianCorrespondence(dist_th = 2)

            meas = BlobTracking.assign(c, blobs, coordinates)
            @test meas.assi == [1,2]

            meas = BlobTracking.assign(c, blobs, [CartesianIndex(1000,1000); coordinates])
            @test meas.assi == [2,3]

            meas = BlobTracking.assign(c, [blobs; blobs], coordinates)
            @test meas.assi == [1,2,0,0]

        end

        @testset "MCCorrespondence" begin
            @info "Testing MCCorrespondence"
            bt = BlobTracker(2:2, 10, 5.0)
            coordinates = locs
            blobs = Blob.(bt.params, coordinates)
            c = MCCorrespondence()

            meas = BlobTracking.assign(c, blobs, coordinates)
            @test median(meas[i].assi[1] for i in eachindex(meas)) == 1
            @test median(meas[i].assi[2] for i in eachindex(meas)) == 2

            meas = BlobTracking.assign(c, blobs, [CartesianIndex(1000,1000); coordinates])
            @test median(meas[i].assi[1] for i in eachindex(meas)) == 2
            @test median(meas[i].assi[2] for i in eachindex(meas)) == 3

            meas = BlobTracking.assign(c, [blobs; blobs], coordinates)
            @test median(meas[i].assi[1] for i in eachindex(meas)) <= 1
            @test median(meas[i].assi[2] for i in eachindex(meas)) <= 2

        end

    end

    @testset "FrameBuffer" begin
        @info "Testing FrameBuffer"
        b = FrameBuffer{Float64}(2,2,2)
        @test length(b) == 0
        @test !isready(b)
        @test_throws BoundsError b[1]
        @test_throws BoundsError b[1,1,1]
        b1 = randn(2,2)
        push!(b, b1)
        @test length(b) == 1
        @test isready(b)
        @test b[1] == b1
        @test b[1,1,1] == b1[1,1]
        @test length(collect(b)) == 1
        for f in (median, mean, sum)
            @test f(b) == f(b.b[:,:,1:1], dims=3)[:,:,1]
        end

        push!(b, randn(2,2))

        for f in (median, mean, sum, std, var)
            @test f(b) == f(b.b[:,:,1:2], dims=3)[:,:,1]
        end
        @test b[1] == b1
        @test size(Matrix(b)) == (4,2)
        @test length(b) == 2
        @test length(collect(b)) == 2
        @test BlobTracking.imgsize(b) == (2,2)

        db = diff(b)
        @test length(db) == 1
        @test db.b == diff(b.b, dims=3)

    end

    @testset "Background extraction" begin
        @info "Testing Background extraction"
        @testset "MedianBackground" begin
            @info "Testing MedianBackground"
            be = MedianBackground(img,3)
            update!(be, img2)
            update!(be, img3)
            @test BlobTracking.background(be, img) == 0 .* img
            @test BlobTracking.foreground(be, img) == img
        end

        @testset "DiffBackground" begin
            @info "Testing DiffBackground"
            be = DiffBackground(img)
            update!(be, img2)
            @test imadjustintensity(BlobTracking.background(be, img)) == complement.(img)
            @test imadjustintensity(BlobTracking.foreground(be, img)) == img
        end
    end

    @testset "Tracking" begin
        @info "Testing Tracking"
        bt = BlobTracker(2:2, 10, 5.0)
        result = TrackingResult()
        ws = BlobTracking.Workspace(img, bt)
        BlobTracking.prepare_image!(ws,bt,img)
        measurement = BlobTracking.Measurement(ws, bt, img, result)
        @test measurement.assi == Int[]
        @test measurement.coordinates == locs
        @test isempty(result.blobs)
        BlobTracking.spawn_blobs!(result, bt, measurement)
        @test length(result.blobs) == 2
        @test isempty(result.dead)
        @test location.(result.blobs) == locs
        @test tracem(result.blobs[1])[1] == locs[1]
        @test tracem(result.blobs[2])[1] == locs[2]

        blobs = result.blobs
        BlobTracking.prepare_image!(ws,bt,img2)
        measurement = BlobTracking.Measurement(ws, bt, img2, result)
        @test measurement.coordinates == locs2
        @test measurement.assi == 1:2
        BlobTracking.predict!(result)
        @test location.(blobs) == locs
        filter!(result, bt, measurement)
        @test length(measurement.coordinates) == 2
        @test measurement.assi == 1:2
        BlobTracking.correct!(result, measurement)
        @test location.(blobs) == locs2

        @test result.blobs[1].tracem[2] == locs2[1]
        @test result.blobs[2].tracem[2] == locs2[2]

        @test all(>(0), blobs[1].kf.x[3:4]) # Positive velocities
        @test all(>(0), blobs[2].kf.x[3:4])

        @test blobs[1].counter == 0
        @test blobs[2].counter == 0

        measurement = BlobTracking.Measurement(ws, bt, 0*img2, result)
        @test isempty(measurement.coordinates)
        @test all(iszero, measurement.assi)
        filter!(result, bt, measurement)
        @test isempty(measurement.coordinates)
        @test all(iszero, measurement.assi)
        BlobTracking.correct!(result, measurement)
        @test location.(blobs) == locs2

        @test result.blobs[1].tracem[3] == OOB
        @test result.blobs[2].tracem[3] == OOB

        @testset "Coordinate iterator" begin
            @info "Testing Coordinate iterator"
            bt = BlobTracker(2:2, 10, 5.0)
            coords = BlobTracking.get_coordinates(bt, [img,img2,img,img2], threads=true)
            @test length(coords) == 4
            @test coords == [locs,locs2,locs,locs2]

            coords = BlobTracking.get_coordinates(bt, [img,img2,img,img2], threads=false)
            @test length(coords) == 4
        end

        @testset "track_blobs" begin
            @info "Testing track_blobs"
            bt = BlobTracker(2:2, 10, 5.0)
            result = track_blobs(bt,[img,img2])

            @test result.blobs[1].tracem[1] == locs[1]
            @test result.blobs[2].tracem[1] == locs[2]
            @test result.blobs[1].tracem[2] == locs2[1]
            @test result.blobs[2].tracem[2] == locs2[2]

            result = track_blobs(bt,[img,img2], threads=false)

            @test result.blobs[1].tracem[1] == locs[1]
            @test result.blobs[2].tracem[1] == locs[2]
            @test result.blobs[1].tracem[2] == locs2[1]
            @test result.blobs[2].tracem[2] == locs2[2]

            bt = BlobTracker(2:2, 10, 5.0, correspondence=MCCorrespondence())
            result = track_blobs(bt,[img,img3])

            @test result.blobs[1].tracem[1] == locs[1]
            @test result.blobs[2].tracem[1] == locs[2]
            @test result.blobs[1].tracem[2] == locs3[1]
            @test result.blobs[2].tracem[2] == locs3[2]
            # @test tracem.(allblobs(result))
            @test location(result) == [locs, locs3]

            bt = BlobTracker(2:2, 10, 5.0)
            coords = BlobTracking.get_coordinates(bt, [img,img2,img,img2], threads=true)
            result = track_blobs(bt,coords)

            @test result.blobs[1].tracem[1] == locs[1]
            @test result.blobs[2].tracem[1] == locs[2]
            @test result.blobs[1].tracem[2] == locs2[1]
            @test result.blobs[2].tracem[2] == locs2[2]
        end
    end
    @testset "display" begin
        @info "Testing display"
        bt = BlobTracker(2:2, 10, 5.0)
        recorder = Recorder()
        result = track_blobs(bt,[N0f8.(Gray.(img)),N0f8.(Gray.(img2)),N0f8.(Gray.(img)),N0f8.(Gray.(img2)),N0f8.(Gray.(img)),N0f8.(Gray.(img2)),N0f8.(Gray.(img)),N0f8.(Gray.(img2)),N0f8.(Gray.(img)),N0f8.(Gray.(img2))], display=img->println("displaying image"), recorder=recorder)
        traces =  trace(result, minlife=2)
        tracems =  tracem(result, minlife=2)

        drawimg = RGB.(img)
        draw!(drawimg, traces, c=RGB(0,0,0.5))
        draw!(drawimg, tracems, c=RGB(0.5,0,0))
        @test isfile(recorder.filename)
        io = VideoIO.open(recorder.filename)
        vid = VideoIO.openvideo(io)
        frame1 = first(vid)
        @test size(frame1) == size(img)

        if isfile(recorder.filename)
            rm(recorder.filename)
        end

        @testset "tune_sizes" begin
            @info "Testing tune_sizes"
            bt = BlobTracker(2:2, 10, 5.0)
            tune_sizes(bt,img)

        end
    end
end


text:
module TestUpdateMethod

using MLJBase
using Test

import LightGBM

# Setup stuff
nrows = 3_000
nfeatures = 10

x = randn(nrows, nfeatures)
y = randn(nrows)
k = Dict(Symbol(i) => Int64(i) for i in 1:nfeatures)
X = MLJBase.Tables.MatrixTable(collect(keys(k)), k, x)
w = abs.(randn(nrows))
w = w / sum(w)

# setup MLJ machinery -- will be used in several testsets
r = LightGBM.MLJInterface.LGBMRegressor(num_iterations=1, min_data_in_leaf=20)
m = MLJBase.machine(r, X, y)

# initial fit
MLJBase.fit!(m; verbosity=0)

@testset "MLJ update" begin
    # machine keeps a field `fitresult` which is the the fitresult,
    # the first element of this tuple is the LGBMEstimator
    iteration_count = LightGBM.get_iter_number(first(m.fitresult))
    cached_iteration_count = sum(m.cache.num_boostings_done)
    @test iteration_count == 1
    @test iteration_count == cached_iteration_count
    @test m.cache.num_boostings_done == [1]

    # check an update chaining on from a fit works
    m.model.num_iterations += 10

    MLJBase.fit!(m; verbosity=0)

    iteration_count = LightGBM.get_iter_number(first(m.fitresult))
    cached_iteration_count = sum(m.cache.num_boostings_done)
    @test iteration_count == 11
    @test iteration_count == cached_iteration_count
    @test m.cache.num_boostings_done == [1, 10]
end

@testset "MLJ update chain" begin
    # check an update chaining on from an update works
    m.model.num_iterations += 5

    MLJBase.fit!(m; verbosity=0)

    iteration_count = LightGBM.get_iter_number(first(m.fitresult))
    cached_iteration_count = sum(m.cache.num_boostings_done)
    @test iteration_count == 16
    @test iteration_count == cached_iteration_count
    @test m.cache.num_boostings_done == [1, 10, 5]

end

@testset "MLJ update less iterations triggers refit" begin
    m.model.num_iterations -= 1
    MLJBase.fit!(m; verbosity=0)

    iteration_count = LightGBM.get_iter_number(first(m.fitresult))
    cached_iteration_count = sum(m.cache.num_boostings_done)
    @test iteration_count == 15
    @test iteration_count == cached_iteration_count
    @test m.cache.num_boostings_done == [15]
end

@testset "MLJ update non-updatable parameters" begin
    # next, check that  changing parameters other than num_iterations results in entirely refitting
    m.model.num_iterations += 1
    m.model.min_data_in_leaf += 5

    MLJBase.fit!(m; verbosity=0)

    iteration_count = LightGBM.get_iter_number(first(m.fitresult))
    cached_iteration_count = sum(m.cache.num_boostings_done)
    @test iteration_count == 16
    @test iteration_count == cached_iteration_count
    @test m.cache.num_boostings_done == [16]
end

@testset "MLJ User provided weights" begin
# check with user provided weights
    r = LightGBM.MLJInterface.LGBMRegressor(num_iterations=1, min_data_in_leaf=20)
    weights_machine = MLJBase.machine(r, X, y, w)

    MLJBase.fit!(weights_machine; verbosity=0)

    iteration_count = LightGBM.get_iter_number(first(weights_machine.fitresult))
    cached_iteration_count = sum(weights_machine.cache.num_boostings_done)
    @test iteration_count == 1
    @test iteration_count == cached_iteration_count
    @test weights_machine.cache.num_boostings_done == [1]

    weights_machine.model.num_iterations += 1
    MLJBase.fit!(weights_machine; verbosity=0)

    iteration_count = LightGBM.get_iter_number(first(weights_machine.fitresult))
    cached_iteration_count = sum(weights_machine.cache.num_boostings_done)
    @test iteration_count == 2
    @test iteration_count == cached_iteration_count
    @test weights_machine.cache.num_boostings_done == [1, 1]

end

end # module


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

module FMICore

# check float size (32 or 64 bits)
juliaArch = Sys.WORD_SIZE
@assert (juliaArch == 64 || juliaArch == 32) "FMICore: Unknown Julia Architecture with $(juliaArch)-bit, must be 64- or 32-bit."
Creal = Cdouble
if juliaArch == 32
    Creal = Cfloat
end

"""
The mutable struct representing an abstract (version unknown) FMU.
"""
abstract type FMU end
abstract type FMUSolution end
abstract type FMUExecutionConfig end
abstract type FMUEvent end

include("FMI2_c.jl")
include("FMI2.jl")

include("FMI3_c.jl")
include("FMI3.jl")

### EXPORTING LISTS START ###

export FMU
export FMU2ExecutionConfiguration, FMU2_EXECUTION_CONFIGURATION_RESET, FMU2_EXECUTION_CONFIGURATION_NO_RESET, FMU2_EXECUTION_CONFIGURATION_NO_FREEING
export FMU3ExecutionConfiguration, FMU3_EXECUTION_CONFIGURATION_RESET, FMU3_EXECUTION_CONFIGURATION_NO_RESET, FMU3_EXECUTION_CONFIGURATION_NO_FREEING

# FMI2.jl
export FMU2, FMU2Component, FMU2ComponentEnvironment, FMU2Solution, FMU2Event
export fmi2StatusToString
export fmi2CausalityToString, fmi2StringToCausality
export fmi2VariabilityToString, fmi2StringToVariability
export fmi2InitialToString, fmi2StringToInitial
export fmi2DependencyKindToString, fmi2StringToDependencyKind

# FMI2_c.jl

# data types
export fmi2ValueReference, fmi2Component
export fmi2Char, fmi2String, fmi2Boolean, fmi2Real, fmi2Integer, fmi2Byte, fmi2FMUstate, fmi2ComponentEnvironment, fmi2Enum
export fmi2Status, fmi2StatusOK, fmi2StatusWarning, fmi2StatusDiscard, fmi2StatusError, fmi2StatusFatal, fmi2StatusPending
export fmi2CallbackFunctions
export fmi2Annotation
export fmi2Causality, fmi2CausalityParameter, fmi2CausalityCalculatedParameter, fmi2CausalityInput, fmi2CausalityOutput, fmi2CausalityLocal, fmi2CausalityIndependent
export fmi2Variability, fmi2VariabilityConstant, fmi2VariabilityFixed, fmi2VariabilityTunable, fmi2VariabilityDiscrete, fmi2VariabilityContinuous
export fmi2Initial, fmi2InitialExact, fmi2InitialApprox, fmi2InitialCalculated
export fmi2True, fmi2False
export fmi2Type, fmi2TypeModelExchange, fmi2TypeCoSimulation
export fmi2DatatypeVariable
export fmi2ScalarVariable
export fmi2StatusKind, fmi2StatusKindDoStepStatus, fmi2StatusKindPendingStatus, fmi2StatusKindLastSuccessfulTime, fmi2StatusKindTerminated
export fmi2EventInfo
export fmi2Unit, fmi2SimpleType
export fmi2DependencyKind, fmi2DependencyKindDependent, fmi2DependencyKindConstant, fmi2DependencyKindFixed, fmi2DependencyKindTunable, fmi2DependencyKindUnknown
export fmi2VariableNamingConvention, fmi2VariableNamingConventionFlat, fmi2VariableNamingConventionStructured
export fmi2VariableDependency
export fmi2DefaultExperiment, fmi2Unknown, fmi2ModelStructure
export fmi2ModelDescription
export fmi2ComponentState, fmi2ComponentStateInstantiated, fmi2ComponentStateInitializationMode, fmi2ComponentStateEventMode, fmi2ComponentStateContinuousTimeMode, fmi2ComponentStateTerminated, fmi2ComponentStateError, fmi2ComponentStateFatal

# functions
export fmi2Instantiate, fmi2FreeInstance!, fmi2GetTypesPlatform, fmi2GetVersion
export fmi2SetDebugLogging, fmi2SetupExperiment, fmi2EnterInitializationMode, fmi2ExitInitializationMode, fmi2Terminate, fmi2Reset
export fmi2GetReal!, fmi2SetReal, fmi2GetInteger!, fmi2SetInteger, fmi2GetBoolean!, fmi2SetBoolean, fmi2GetString!, fmi2SetString
export fmi2GetFMUstate!, fmi2SetFMUstate, fmi2FreeFMUstate!, fmi2SerializedFMUstateSize!, fmi2SerializeFMUstate!, fmi2DeSerializeFMUstate!
export fmi2GetDirectionalDerivative!, fmi2SetRealInputDerivatives, fmi2GetRealOutputDerivatives!
export fmi2DoStep, fmi2CancelStep, fmi2GetStatus!, fmi2GetRealStatus!, fmi2GetIntegerStatus!, fmi2GetBooleanStatus!, fmi2GetStringStatus!
export fmi2SetTime, fmi2SetContinuousStates, fmi2EnterEventMode, fmi2NewDiscreteStates!, fmi2EnterContinuousTimeMode, fmi2CompletedIntegratorStep!
export fmi2GetDerivatives!, fmi2GetEventIndicators!, fmi2GetContinuousStates!, fmi2GetNominalsOfContinuousStates!

# FMI3.jl
export FMU3, FMU3Instance, FMU3InstanceEnvironment, FMU3Solution, FMU3Event
export fmi3StatusToString
export fmi3CausalityToString, fmi3StringToCausality
export fmi3VariabilityToString, fmi3StringToVariability
export fmi3InitialToString, fmi3StringToInitial
export fmi3TypeToString, fmi3StringToType
export fmi3IntervalQualifierToString, fmi3StringToIntervalQualifier
export fmi3DependencyKindToString, fmi3StringToDependencyKind
export fmi3VariableNamingConventionToString, fmi3StringToVariableNamingConvention

# FMI3_c.jl

# data types
export fmi3ValueReference, fmi3ModelDescription, fmi3Instance, fmi3ModelVariable, fmi3DatatypeVariable
export fmi3False, fmi3True
export fmi3Status, fmi3StatusOK, fmi3StatusWarning, fmi3StatusDiscard, fmi3StatusError, fmi3StatusFatal
export fmi3Causality, fmi3CausalityParameter, fmi3CausalityCalculatedParameter, fmi3CausalityInput, fmi3CausalityOutput, fmi3CausalityLocal, fmi3CausalityIndependent, fmi3CausalityStructuralParameter
export fmi3Initial, fmi3InitialExact, fmi3InitialApprox, fmi3InitialCalculated
export fmi3Variability, fmi3VariabilityConstant, fmi3VariabilityFixed, fmi3VariabilityTunable, fmi3VariabilityDiscrete, fmi3VariabilityContinuous
export fmi3Type, fmi3TypeModelExchange, fmi3TypeCoSimulation, fmi3TypeScheduledExecution
export fmi3Clock
export fmi3IntervalQualifier, fmi3IntervalQualifierIntervalChanged, fmi3IntervalQualifierIntervalNotYetKnown, fmi3IntervalQualifierIntervalUnchanged
export fmi3Char, fmi3String, fmi3Boolean, fmi3Binary, fmi3Float32, fmi3Float64, fmi3Int8, fmi3UInt8, fmi3Int16, fmi3UInt16, fmi3Int32, fmi3UInt32, fmi3Int64, fmi3UInt64, fmi3Byte, fmi3FMUState, fmi3InstanceEnvironment, fmi3Enum
export fmi3DependencyKind, fmi3DependencyKindDependent, fmi3DependencyKindIndependent, fmi3DependencyKindConstant, fmi3DependencyKindTunable, fmi3DependencyKindDiscrete, fmi3DependencyKindFixed, fmi3DependencyKindGetVariableDependencies
export fmi3VariableNamingConvention, fmi3VariableNamingConventionFlat, fmi3VariableNamingConventionStructured
export fmi3VariableDependency
export fmi3InstanceState, fmi3InstanceStateInstantiated, fmi3InstanceStateInitializationMode, fmi3InstanceStateEventMode, fmi3InstanceStateStepMode, fmi3InstanceStateClockActivationMode, fmi3InstanceStateContinuousTimeMode
export fmi3ConfigurationMode, fmi3ReconfigurationMode, fmi3InstanceStateTerminated, fmi3InstanceStateError, fmi3InstanceStateFatal

# functions
export fmi3InstantiateCoSimulation, fmi3InstantiateModelExchange, fmi3InstantiateScheduledExecution, fmi3FreeInstance!, fmi3GetVersion
export fmi3SetDebugLogging, fmi3EnterInitializationMode, fmi3ExitInitializationMode, fmi3Terminate, fmi3Reset
export fmi3GetFloat32!, fmi3SetFloat32, fmi3GetFloat64!, fmi3SetFloat64
export fmi3GetInt8!, fmi3SetInt8, fmi3GetInt16!, fmi3SetInt16,fmi3GetInt32!, fmi3SetInt32, fmi3GetInt64!, fmi3SetInt64
export fmi3GetUInt8!, fmi3SetUInt8, fmi3GetUInt16!, fmi3SetUInt16,fmi3GetUInt32!, fmi3SetUInt32, fmi3GetUInt64!, fmi3SetUInt64
export fmi3GetBoolean!, fmi3SetBoolean, fmi3GetString!, fmi3SetString, fmi3GetBinary!, fmi3SetBinary, fmi3GetClock!, fmi3SetClock
export fmi3GetFMUState!, fmi3SetFMUState, fmi3FreeFMUState!, fmi3SerializedFMUStateSize!, fmi3SerializeFMUState!, fmi3DeSerializeFMUState!
export fmi3SetIntervalDecimal, fmi3SetIntervalFraction, fmi3GetIntervalDecimal!, fmi3GetIntervalFraction!, fmi3GetShiftDecimal!, fmi3GetShiftFraction!, fmi3ActivateModelPartition
export fmi3GetNumberOfVariableDependencies!, fmi3GetVariableDependencies!
export fmi3GetDirectionalDerivative!, fmi3GetAdjointDerivative!, fmi3GetOutputDerivatives!
export fmi3EnterConfigurationMode, fmi3ExitConfigurationMode
export fmi3GetNumberOfContinuousStates!, fmi3GetNumberOfEventIndicators!
export fmi3DoStep!, fmi3EnterStepMode
export fmi3SetTime, fmi3SetContinuousStates, fmi3EnterEventMode, fmi3UpdateDiscreteStates, fmi3EnterContinuousTimeMode, fmi3CompletedIntegratorStep!
export fmi3GetContinuousStateDerivatives, fmi3GetEventIndicators!, fmi3GetContinuousStates!, fmi3GetNominalsOfContinuousStates!, fmi3EvaluateDiscreteStates

### EXPORTING LISTS END ###

end # module


text:
using SpecialFunctions
using LsqFit

import Distributions: TDist, cdf

import Base: show, <<, push!, count, isempty, merge, merge!, zero, zeros, sum

if VERSION < v"0.7.0-beta.85"
    import Base: mean, var
else
    import Statistics: mean, var
end

export MCObservable, ScalarObservable, VectorObservable
export mean, var, stderror, confidence_interval
export p_value
export show, dump_plot
export merge, merge!

export confidence_interval

abstract type MCObservable end
abstract type ScalarObservable <: MCObservable end
abstract type VectorObservable <: MCObservable end

function p_value(X::MCObservable, y::Real; verbose::Bool=false)
    MX = mean(X)
    mid = 0.5(MX+y)
    er = max(stderror(X), abs(mid)*sqrt(eps()), sqrt(nextfloat(0.0)))
    t = (MX-y)/er
    d = TDist(count(X)-1)
    if verbose
        @show MX, y, er, t
    end
    return 2cdf(d,-abs(t))
end

function p_value(X::MCObservable, Y::MCObservable)
    NX = count(X)
    NY = count(Y)
    MX = mean(X)
    MY = mean(Y)
    mid = 0.5(MX+MY)
    N = 1.0/(1.0/NX + 1.0/NY)
    V = ( (NX-1)*var(X) + (NY-1)*var(Y) )/(NX+NY-2)
    se = max(sqrt(V/N), abs(mid)*sqrt(eps()), sqrt(nextfloat(0.0)))
    t = (MX-MY)/se
    d = TDist(NX+NY-2)
    return 2cdf(d,-abs(t))
end

isempty(obs::MCObservable) = count(obs) == 0

zero(::Type{Obs}) where (Obs<:MCObservable) = Obs()
zero(obs::MCObservable) = zero(typeof(obs))
function zeros(::Type{Obs}, dim...) where (Obs<:MCObservable)
    reshape(Obs[Obs() for i in 1:prod(dim)], dim)
end

function show(io::IO, obs::MCObservable)
    if !isempty(obs)
        print(io, mean(obs), " +/- ", stderror(obs))
    else
        print(io, "No Entries")
    end
end

function dump_plot(io::IO, obs::MCObservable; put_following_space::Bool=false, newline::Bool=false)
    print(io, mean(obs), " ", stderror(obs))
    if put_following_space
        print(io, " ")
    end
    if newline
        println(io)
    end
end

const confidence_rate_1sigma = 0.5erf(0.5sqrt(2.0))

include("util.jl")
include("observableset.jl")
include("parsesigma.jl")
include("tiny.jl")
include("tinyvector.jl")
include("simple.jl")
include("simplevector.jl")
include("binning.jl")
include("binningvector.jl")
include("jackknife.jl")
include("jackknifevector.jl")

## these three definitions are needed to resolve ambiguousness with Base.<<(Any,Integer)
<<(obs::MCObservable, x::Int32) = push!(obs,x)
<<(obs::MCObservable, x::Int64) = push!(obs,x)
<<(obs::MCObservable, x::Integer) = push!(obs,x)

<<(obs::MCObservable, x) = push!(obs,x)


text:
using RunMetrics
using Test, SafeTestsets

@time @safetestset "TCX Parser Tests" begin include("test_tcx_parser.jl") end
@time @safetestset "RunSummary Tests" begin include("test_runsummary.jl") end
@time @safetestset "Regression Tests" begin include("test_multi_regression.jl") end


text:
module MonteCarloCollisions

using StaticArrays
using PhysicalConstants.CODATA2018
using Unitful
using Interpolations
using Optim
using Random: rand, randn, randexp
using LinearAlgebra: norm

export Neutrals, Particles, Interaction, Interactions, add_interaction!
export load_interaction_lxcat, load_interactions_lxcat, svmax_find!, init_rates!, make_interactions
export init_time, init_monoenergetic, init_thermal
export advance, advance!, energy
export m_e, q_e, amu

const m_e = ustrip(u"kg", CODATA2018.m_e)
const amu = ustrip(u"kg", CODATA2018.AtomicMassConstant)
const q_e = ustrip(u"C", CODATA2018.e)
const k_B = ustrip(u"J/K", CODATA2018.k_B)

@inline function random_direction(l=1.)
    """generate a 3D vector with random direction and length l"""
    cos_theta =  1. - 2*rand()
    sin_theta = sqrt(1. - cos_theta^2)
    sin_cos_phi = sincos(2*pi*rand())
    SVector(
        l*cos_theta,
        l*sin_theta * sin_cos_phi[1],
        l*sin_theta * sin_cos_phi[2]
    )
end

@inline function random_maxwell_v(vthermal)
    SVector(randn(),randn(), randn())*vthermal # sqrt(k_B*T/m)
end

@inline function choose(probs)
    # equivalent to dist=Categorical(probs); rand(dist),
    # > 30% faster, but no error checking (sum(probs) == 1 ?)
    l = length(probs)
    i = 1
    if l > 1
        U = rand()
        @inbounds while (U -= probs[i]) > 0 && i < l
            i += 1
        end
    end
    i
end

abstract type Species end

struct Neutrals <: Species
    name::String
    T::Float64
    m::Float64
    n::Float64
    vthermal::Float64
end

Neutrals(name, T, m, n) = Neutrals(name, T, m, n, sqrt(k_B*T/m))

@inline function random_sample(s::Neutrals)
    random_maxwell_v(s.vthermal)
end

mutable struct Particle
    v::SVector{3, Float64}
    t::Float64
    tau::Float64
end

struct Particles <: Species
    name::String
    m::Float64
    q::Float64
    n::Int64
    list::Array{Particle, 1}
end

function Particles(name, m, q, n)
    Particles(
        name, m, q, n,
        [Particle(SVector(0.,0.,0.), 0., 0.) for i in 1:n]
    )
end

function init_monoenergetic(particles::Particles, Eev)
    vtot = sqrt(2*Eev*particles.q/particles.m)
    for p in particles.list
        p.v = random_direction(vtot)
    end
end

function init_thermal(particles::Particles, T)
    vthermal = sqrt(k_B*T/particles.m)
    for p in particles.list
        p.v = random_maxwell_v(vthermal)
    end
end

function init_time(particles::Particles)
    for p in particles.list
        p.t = 0.
    end
end

function energy(species::Species, particle::Particle)
    0.5*species.m*sum(particle.v.^2)
end

function energy(particles::Particles)
    map(p -> energy(particles, p), particles.list)
end

InterpolationType = typeof(interpolate(([0.,1.],), [0.,1.], Gridded(Linear())))

struct Interaction
    name::String
    DE::Float64
    species1::Particles
    species2::Neutrals
    mu::Float64
    sigmav::InterpolationType
    #Interpolations.GriddedInterpolation{Float64, 1, Float64, Gridded{Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}}}
    vthresh::Float64
end

import Interpolations: interpolate, Gridded, Linear
function Interaction(name, sigma, E, DE, species1, species2)
    mu = species1.m*species2.m / (species1.m + species2.m)

    v = sqrt.((2*q_e/mu).*E)

    sigmav_interp = interpolate((v,), sigma.*v, Gridded(Linear()))

    vthresh = sigma[2]==0. ? v[2] : v[1]

    Interaction(
        name, DE, species1, species2, mu,
        sigmav_interp, vthresh
    )
end

@inline function sigmav(interaction, v)
    # simple optimization for fast eval of below-treshold collisions
    if v < interaction.vthresh
        sv = 0.
    else
        sv = interaction.sigmav(v)
    end
    sv
end

function load_interaction_lxcat(filename, process, species)
    # load cross section from text file from lxcat database
    nl = "\r\n"
    open(filename) do f
        readuntil(f, process*nl*species*nl)
        DE = parse(Float64, readline(f))

        readuntil(f, "-----"*nl)
        data_string = readuntil(f, "-----"*nl)
        data_numbers = map(split(data_string, nl)[1:end-1]) do line
            map(x->parse(Float64, x), split(line, "\t"))
        end
        DE, hcat(data_numbers...)'
    end
end

function load_interactions_lxcat(filename, species1, species2)
    collision_types = ["ELASTIC", "EFFECTIVE", "EXCITATION", "IONIZATION", "ATTACHMENT"]
    result = Interaction[]
    # load cross section from text file from lxcat database
    nl = "\r\n"
    open(filename) do f
        for l in eachline(f)
            if l in collision_types
                collision_type = l
                collision_target = readline(f)
                if !(collision_type in ["ELASTIC", "EFFECTIVE"])
                    DE = parse(Float64, readline(f))
                else
                    DE = 0.
                end

                readuntil(f, "-----"*nl)
                data_string = readuntil(f, "-----"*nl)
                data_numbers = map(split(data_string, nl)[1:end-1]) do line
                    map(x->parse(Float64, x), split(line, "\t"))
                end
                CS = hcat(data_numbers...)'
                push!(result, Interaction(collision_type*" "*collision_target,
                        CS[:,2], CS[:,1], DE, species1, species2
                    ))
            end
        end
    end
    result
end

mutable struct Interactions
    # structure representing all interactions of Species spec1 with
    # partners from list spec2_list
    spec1::Particles
    list::Vector{Vector{Interaction}} # list[i] contains list of interactions with spec2_list[1]
    spec2_names::Vector{String}
    spec2_list::Vector{Neutrals} # list of interacting species
    svmax_list::Vector{Float64}  # maximum of total sigma*v_rel for each spec2
    prob_list::Vector{Float64}   # interaction rates of spec2 divided by total rate
    rate::Float64 # (s-1) total interaction rate of spec1
end

Interactions(species) = Interactions(species, [], [], [], [], [], 0.)

function add_interaction!(inters, interaction)
    i = findfirst(isequal(interaction.species2.name), inters.spec2_names)
    if isnothing(i)
        push!(inters.list, [interaction])
        push!(inters.spec2_names, interaction.species2.name)
        push!(inters.spec2_list, interaction.species2)
    else
        push!(inters.list[i], interaction)
    end
end

function make_interactions(species, interaction_list, Emax=20.)
    interactions = Interactions(species)
    for interaction in interaction_list
        add_interaction!(interactions, interaction)
    end
    svmax_find!(interactions, Emax)
    init_rates!(interactions)
    interactions
end

function svmax_find!(inters, Emax)
    inters.svmax_list = zeros(length(inters.list))
    inters.prob_list = zeros(length(inters.list))
    for (i, inter) in enumerate(inters.list)
        vmax = sqrt(2*Emax*q_e/inter[1].mu)
        sumv_neg(v) = -sum([sigmav(I, v) for I in inter])
        res = Optim.optimize(sumv_neg, 0, vmax)
        inters.svmax_list[i] = -res.minimum
    end
end

function init_rates!(inters)
    inters.rate = 0
    for (i, svmax) in enumerate(inters.svmax_list)
        inters.prob_list[i] = inters.spec2_list[i].n * svmax
        inters.rate += inters.prob_list[i]
    end
    inters.prob_list ./= inters.rate
end

@inline function scatter(v1, v2, m1, m2, v_rel_norm, DE)
    # transform v1 into center of mass system
    v_cm = (v1*m1 + v2*m2)/(m1 + m2)

    # account for energy loss
    if DE != 0
        mu = m1*m2/(m1+m2)
        E = 0.5*mu*v_rel_norm^2 - DE*q_e
        #if self.collision_type == "ionization": # only the primary electrons are tracked in case of ionization
        #    pE = E*np.random.uniform(0, 1, len(v1))
        #end
        v_rel_norm = sqrt(2*E/mu)
    end

    # random rotation
    v_rel = random_direction(v_rel_norm)
    # reverse transformation
    v1_cm = m2/(m1 + m2) * v_rel
    v1 = v1_cm .+ v_cm
    v1
end

function scatter(v::SVector{3, Float64}, m, inters::Interactions)

    # select colliding species
    i_species = choose(inters.prob_list)
    spec2 = inters.spec2_list[i_species]
    svmax = inters.svmax_list[i_species]
    interlist = inters.list[i_species]

    v_partner = random_sample(spec2)
    v_rel = norm(v-v_partner)

    # select interaction type
    U_interactions = rand()*svmax
    for inter in  interlist
        U_interactions -= sigmav(inter, v_rel)
        if U_interactions < 0
            v = scatter(v, v_partner, m, spec2.m, v_rel, inter.DE)
            break
        end
    end
    v
end

function advance!(particles::Particles, interactions::Interactions, E, tmax)
    q = particles.q
    m = particles.m
    accel = E.*q/m
    tau_mean = 1/interactions.rate

    Threads.@threads for particle in particles.list
        t = particle.t
        v = particle.v
        tau = particle.tau

        while true
            tau = randexp()*tau_mean
            if t+tau > tmax
                tau = max(tmax - t, 0.)
                v, t = advance(v, t, accel, tau)
                break
            end

            v, t = advance(v, t, accel, tau)
            v = scatter(v, m , interactions,)

        end
        particle.v = v
        particle.t = t
    end
end

@inline function advance(v, t, a, tau)
    v+a*tau, t+tau
end

end


text:
type RegisterGradient
    libeval::Function

    dx::Float64
    funcs::Vector{Function}
end

function unarylibrary(model::Model, components::Vector{Symbol}, names::Vector{Symbol})
    function evaluator(xx::Vector, funcs::Vector{Function})
        setparameters(model, components, names, xx)
        run(model)

        map(func -> func(model), funcs)
    end

    evaluator
end

RegisterGradient(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, dx::Float64) =
    RegisterGradient(unarylibrary(model, components, names), dx, Function[])

function addfunction!(rg::RegisterGradient, func::Function)
    rg.funcs = [rg.funcs; func]
end

function getgradients(rg::RegisterGradient, xx::Vector{Float64}; verbose=false, usesparse=true)
    baseline = rg.libeval(xx, rg.funcs)
    if any(isnan(baseline))
        error("Cannot compute the baseline for gradients.")
    end

    if usesparse
        gradients = spzeros(length(rg.funcs), length(xx))
    else
        gradients = Array{Float64}(length(rg.funcs), length(xx))
    end

    for ii in 1:length(xx)
        println(ii / length(xx))
        # Adjust a single value
        xx[ii] = xx[ii] + rg.dx

        newvals = rg.libeval(xx, rg.funcs)
        if any(isnan(newvals))
            error("Cannot calculate the gradient for constraint $ii")
        end

        if usesparse
            for jj in 1:length(rg.funcs)
                if newvals[jj] != baseline[jj]
                    gradients[jj, ii] = (newvals[jj] - baseline[jj]) / rg.dx
                end
            end
        else
            gradients[:, ii] = (newvals - baseline) / rg.dx
        end

        xx[ii] = xx[ii] - rg.dx
    end

    baseline, gradients
end


text:
@testset "Obstacles" begin
    @test typeof(randpoly(rand(2), n=3, scale=rand())) <: ConvexPolygon{2, 3}
    @test typeof(randpoly(rand(2), n=8, scale=π)) <: ConvexPolygon{2, 8}
    @test typeof(randpoly(rand(2), n=15, scale=rand())) <: ConvexPolygon{2, 15}
    @test typeof(@point rand(2)) <: ConvexPolygon{2, 1}
    @test typeof(@point rand(3)) <: ConvexPolygon{3, 1}
    @test typeof(@line rand(2), rand(2)) <: ConvexPolygon{2, 2}
    @test typeof(@line rand(3), rand(3)) <: ConvexPolygon{3, 2}
    @test typeof(@rect rand(2), rand(2)) <: ConvexPolygon{2, 4}
    @test typeof(@rect rand(3), rand(3)) <: ConvexPolygon{3, 8}
    @test typeof(@square rand(2), π) <: ConvexPolygon{2, 4}
    @test typeof(@square rand(3), rand()) <: ConvexPolygon{3, 8}
end


text:
@testset "GriddedTests" begin

include("gridded.jl")
include("mixed.jl")

end


text:
# Module Definition
# #+HTML: Source
# #+HTML:

# [[file:~/.julia/dev/MatrixProductStates/README.org::*Module%20Definition][Module Definition:1]]
module MatrixProductStates

using LinearAlgebra, TensorOperations, TensorCast, LowRankApprox, Arpack, Strided, SparseArrays
#using ProgressMeter

export *, /, ==, ≈, isequal, adjoint, getindex, randn
export MPS, MPO, left, right, compress, imag_time_evolution, real_time_evolution, real_timedep_time_evolution, rightcanonical, leftcanonical
export ground_state, two_point_correlator, realize

include("utils.jl")
include("MPS.jl")
include("MPO.jl")
include("compression.jl")
include("contraction.jl")
include("groundstate.jl")
include("correlation.jl")
include("timeevolution.jl")

end
# Module Definition:1 ends here


text:
# Integration test combining AtomsBase and NBodySimulator

@testset "Integration: AB & NBS" begin
    # TODO
end


text:
Base.show(io::IO, metric::T) where T<:AbstractMetric = print(io, string(T))

# AverageLoss
mutable struct AverageLoss <: AbstractMetric
    loss
    last
    count
    AverageLoss() = new(nothing, nothing, nothing)
end

function on(::EpochBegin, ::AbstractFittingPhase, metric::AverageLoss, learner)
    metric.loss = 0.
    metric.count = 0
end

function on(::BatchEnd, phase::AbstractFittingPhase, metric::AverageLoss, learner)
    metric.loss += learner.batch.loss
    metric.last = learner.batch.loss
    metric.count += 1
end

stateaccess(::AverageLoss) = (batch = (loss = Read(),),)

stepvalue(metric::AverageLoss) = metric.last
epochvalue(metric::AverageLoss) = metric.loss / metric.count

Base.show(io::IO, loss::AverageLoss) = print(io, "Loss()")

# OnlineMetrics metric

mutable struct Metric{T} <: AbstractMetric
    metricfactory
    fn
    name
    metric::OnlineStat{T}
    last::T
    device
end

stepvalue(metric::Metric) = metric.last
epochvalue(metric::Metric) = OnlineStats.value(metric.metric)

"""
    Metric(fn, name = string(fn); device = cpu)

`fn(y_pred, y)`
"""
function Metric(fn, name = string(fn); device = cpu, metric = Mean)
    Metric(() -> metric(), fn, name, metric(), Inf, device)
end

Base.show(io::IO, metric::Metric) = print(io, metric.name)

function on(::EpochBegin, ::AbstractFittingPhase, metric::Metric, learner)
    metric.metric = metric.metricfactory()
end

function on(::BatchEnd, ::AbstractFittingPhase, metric::Metric, learner)
    metric.last = metric.fn(
            metric.device(learner.batch.ŷs),
            metric.device(learner.batch.ys),
        )
    OnlineStats.fit!(
        metric.metric,
        metric.last,
    )
end

stateaccess(::Metric) = (batch = (ŷs = Read(), ys = Read()),)

Accuracy() = Metric(accuracy)


text:
module Samplers

using DocStringExtensions

import ..Bincounts
using ..Configuration
import ..Delayed
import ..GMM
import ..IndexObjects
import ..ForgettingWeights
import ..LogGMM

include("adaptive_parzen.jl")
include("normal.jl")
include("quantnormal.jl")
include("lognormal.jl")
include("logquantnormal.jl")
include("randindex.jl")
include("uniform.jl")
include("quantuniform.jl")
include("categoricalindex.jl")
include("logquantuniform.jl")
include("loguniform.jl")

end # module Samplers