text:
"""
Returns one default FUND model and one model with additional emissions of the specified gas in the specified year.
"""
function getmarginalmodels(; gas = :C, emissionyear = 2010, parameters = nothing, yearstorun = 1050)

    # Get default FUND model
    m1 = MimiFUND.get_model(nsteps = yearstorun, params = parameters)

    # Get model to add marginal emissions to
    m2 = MimiFUND.get_model(nsteps = yearstorun, params = parameters)
    add_comp!(m2, Mimi.adder, :marginalemission, before = :climateco2cycle, first = 1951)
    addem = zeros(yearstorun)
    addem[getindexfromyear(emissionyear)-1:getindexfromyear(emissionyear) + 8] .= 1.0
    set_param!(m2, :marginalemission, :add, addem)

    # Reconnect the appropriate emissions in the marginal model
    if gas == :C
        connect_param!(m2, :marginalemission, :input, :emissions, :mco2)
        connect_param!(m2, :climateco2cycle, :mco2, :marginalemission, :output, repeat([missing], yearstorun + 1))
    elseif gas == :CH4
        connect_param!(m2, :marginalemission, :input, :emissions, :globch4)
        connect_param!(m2, :climatech4cycle, :globch4, :marginalemission, :output, repeat([missing], yearstorun + 1))
    elseif gas == :N2O
        connect_param!(m2, :marginalemission, :input, :emissions, :globn2o)
        connect_param!(m2, :climaten2ocycle, :globn2o, :marginalemission, :output, repeat([missing], yearstorun + 1))
    elseif gas == :SF6
        connect_param!(m2, :marginalemission, :input, :emissions,:globsf6)
        connect_param!(m2, :climatesf6cycle, :globsf6, :marginalemission, :output, repeat([missing], yearstorun + 1))
    else
        error("Unknown gas.")
    end

    # Run each model
    run(m1)
    run(m2)

    return m1, m2
end

"""
Returns the social cost per one ton of additional emissions of the specified gas in the specified year.
Uses the specified eta and prtp for discounting, with the option to use equity weights.
"""
function marginaldamage3(; emissionyear = 2010, parameters = nothing, yearstoaggregate = 1000, gas = :C, useequityweights = false, eta = 1.0, prtp = 0.001)
    yearstorun = min(1050, getindexfromyear(emissionyear) + yearstoaggregate)

    m1, m2 = getmarginalmodels(emissionyear = emissionyear, parameters = parameters, yearstorun = yearstorun, gas = gas)

    damage1 = m1[:impactaggregation, :loss]
    # Take out growth effect effect of run 2 by transforming
    # the damage from run 2 into % of GDP of run 2, and then
    # multiplying that with GDP of run 1
    damage2 = m2[:impactaggregation, :loss] ./ m2[:socioeconomic, :income] .* m1[:socioeconomic, :income]

    # Calculate the marginal damage between run 1 and 2 for each
    # year/region
    marginaldamage = (damage2 .- damage1) / 10000000.0

    ypc = m1[:socioeconomic, :ypc]

    df = zeros(yearstorun + 1, 16)
    if !useequityweights
        for r = 1:16
            x = 1.
            for t = getindexfromyear(emissionyear):yearstorun
                df[t, r] = x
                gr = (ypc[t, r] - ypc[t - 1, r]) / ypc[t - 1,r]
                x = x / (1. + prtp + eta * gr)
            end
        end
    else
        globalypc = m1[:socioeconomic, :globalypc]
        df = Float64[t >= getindexfromyear(emissionyear) ? (globalypc[getindexfromyear(emissionyear)] / ypc[t, r]) ^ eta / (1.0 + prtp) ^ (t - getindexfromyear(emissionyear)) : 0.0 for t = 1:yearstorun + 1, r = 1:16]
    end

    scc = sum(marginaldamage[2:end, :] .* df[2:end, :])

    return scc
end


text:
#negative Bernoulli log likelihood
type NBLL <: Functor{2} end
evaluate(::NBLL, pr, y) = y == one(y)? -log(pr):
                          y == zero(y)? -log(one(y) - pr):
                          -(y * log(pr) + (one(y) - y) * log(one(pr) - pr))

# Gradient helper for quantile regression
type Qgrad <: Functor{2} end
evaluate(::Qgrad, r, tau) = r == zero(r)? zero(r) :
                                 r < zero(r)? oftype(r, -tau):
                                 oftype(r, 1.0-tau)

emptyvector{T}(::Type{T}) = Array(T, 0)

Base.At_mul_B!{T<:Base.LinAlg.BlasFloat}(α::Number,A::StridedMatrix{T},x::StridedVector{T},β::Number,y::StridedVector{T}) =
  BLAS.gemv!('T', oftype(T, α), A, x, oftype(T, β), y)

ext: jl


text:
# Nq is the number of positions used to test convergence.
function fitminimax( wfunc::Function,
                    dfunc::Function,
                    𝓧0::Vector{T},
                    fig_num,
                    config) where T

    ##### parse.
    tol_converged = config.tol_converged
    tol_𝓧_spacing = config.tol_𝓧_spacing
    tol_derivative_zero = config.tol_derivative_zero
    tol_no_update = config.tol_no_update
    N_samples_interval = config.N_samples_interval
    N_test_positions = config.N_test_positions
    max_iters = config.max_iters
    L_interval = config.L_interval

    ##### set up.
    𝓧 = copy(𝓧0)
    L = length(𝓧) - 2

    # diagnostics.
    level_error_history = Vector{T}(undef,max_iters)
    xq = convert(Vector{BigFloat},collect(LinRange(-one(T),one(T),N_test_positions)))

    # pre-allocate.
    𝑤 = get𝑤scaled(𝓧)
    h = geth(𝑤, wfunc, dfunc, 𝓧)
    𝑓 = collect( dfunc(𝓧[i]) - (-1)^(i-1)*h/wfunc(𝓧[i]) for i = 1:length(𝓧))
    int_storage_L∞ = SharedArray{Float64}(length(xq))

    f = xx->wfunc(xx)*dfunc(xx)
    f_itp = xx->getp(xx, 𝑓, 𝑤, h, wfunc, dfunc, 𝓧)
    𝑒 = xx->(f(xx)-f_itp(xx))

    sample_position_template = collect( Chebyshev2ndnode(k,config.N_samples_interval) for k = 0:config.N_samples_interval )

    ##### run algorithm 2.
    iter = 1

    # repeat until convergence.
    while iter <= max_iters

        ##### update interpolant.
        get𝑤scaled!(𝑤, 𝓧)
        h = geth(𝑤, wfunc, dfunc, 𝓧)

        for i = 1:length(𝓧)
            𝑓[i] = dfunc(𝓧[i]) - (-1)^(i-1)*h/wfunc(𝓧[i])
        end

        ##### diagnostics
        level_error_history[iter] = h

        ### brute-force uniform sampling method to find the position in [-1,1]
        #       with the largest absolute error.
        L∞_value, L∞_ind, 𝑒_at_L∞_position = discretizeL∞normparallel!(int_storage_L∞, f,f_itp,xq)

        sign_𝑒_at_L∞_position = sign(𝑒_at_L∞_position) > 0

        L∞_position = xq[L∞_ind]
        score = (L∞_value - abs(h))/L∞_value

        ### outputs.
        if config.verbose_flag
            println("iter: ", iter)
            println("h: ", h)
            println("score: ", score)
        end

        if score <= tol_converged
            resize!(level_error_history, iter)

            status_msg = "Converged"
            return 𝓧, level_error_history, iter, status_msg, fig_num
        end

        ### find next set of reference positions.
        fig_num, proceed_flag, status_msg = innerstep(𝓧, h,
                                                        f, f_itp, 𝑒,
                                                        wfunc, dfunc, L_interval,
                                                        L,
                                                        xq,
                                                        L∞_value,
                                                        L∞_position,
                                                        sign_𝑒_at_L∞_position,
                                                        sample_position_template,
                                                        iter,
                                                        fig_num,
                                                        config)

        if config.verbose_flag
            println()
        end

        if !proceed_flag
            resize!(level_error_history, iter)
            println(status_msg)
            return 𝓧, level_error_history, iter, status_msg, fig_num
        end

        iter += 1
    end

    status_msg = "Max. iterations reached."
    println(status_msg)
    return 𝓧, level_error_history, iter, status_msg, fig_num
end

# interval extrema search.
function innerstep(𝓧::Vector{T},
                    h::T,
                    f::Function,
                    f_itp::Function,
                    𝑒::Function,
                    wfunc::Function,
                    dfunc::Function,
                    L_interval::Int,
                    L::Int,
                    xq,
                    L∞_value::Float64,
                    L∞_position::T,
                    sign_𝑒_at_L∞_position::Bool,
                    sample_position_template,
                    iter,
                    fig_num,
                    config) where T

    tol_𝓧_spacing = config.tol_𝓧_spacing
    tol_derivative_zero = config.tol_derivative_zero
    tol_no_update = config.tol_no_update

    if config.plot_flag
        # PyPlot.figure(fig_num)
        # fig_num += 1
    end

    if config.plot_flag
        min_dist, j = findmin(abs.( 𝓧 .- L∞_position ))
        x_j = 𝓧[j]
    end

    ##### find error extrema positions.
    𝓟, 𝑒_𝓟, d_set, X_set = getcandidatesfromextremas( 𝑒,
                    L_interval, 𝓧,
                    sample_position_template,
                    tol_𝓧_spacing, tol_derivative_zero, h)

    if config.plot_flag
        𝓟_before_processing = copy(𝓟)
    end

    𝓧_new, status_flag = processreferenceset(  𝓧,
                                                h,
                                                𝑒,
                                                𝓟,
                                                𝑒_𝓟,
                                                L∞_value,
                                                L∞_position,
                                                sign_𝑒_at_L∞_position,
                                                tol_𝓧_spacing,
                                                #config.tol_h,
                                                iter,
                                                config.verbose_flag)

    if config.plot_flag
        # xq = convert(Vector{BigFloat},collect(LinRange(-1,1,20*length(𝓧))))
        # boundary_positions = collect( X_set[i][1] for i = 1:length(X_set))
        # push!(boundary_positions, X_set[end][end])
        #
        # d, ν = getChebyshevinterpolator(𝑒, L)
        # 𝑒_Chebyitp_interval = xx->evalChebyshevpolynomialinterval(xx, d_set, X_set)
        #
        # PyPlot.plot(𝓧, 𝑒.(𝓧), "+", label = "e(X)")
        #
        # PyPlot.plot(xq, 𝑒.(xq), label = "e()")
        #
        # PyPlot.plot(𝓧_new, 𝑒.(𝓧_new), "x", label = "X_new")
        #
        # println("x_j is ", x_j)
        # println("L∞_position is ", L∞_position)
        #
        #
        # title_string = Printf.@sprintf("𝑒_Chebyitp_interval: boundary and extrema, iter = %d", iter)
        # PyPlot.title(title_string)
        # PyPlot.legend()
    end

    ### update.

    # check if we're stuck.
    if sum(abs.(𝓧-𝓧_new)) < tol_no_update
        println("Reference set no longer changing.")
        return fig_num, false, "Reference set no longer changing."
    end

    𝓧[:] = 𝓧_new

    return fig_num, true, ""
end

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
# convert binary into HDF5 data

using HDF5

datasets = [("train", ["data_batch_$i.bin" for i in 1:5]),
            ("test", ["test_batch.bin"])]

const width      = 32
const height     = 32
const channels   = 3
const batch_size = 10000

mean_model = zeros(Float32, width, height, channels, 1)

for (key, sources) in datasets
  h5open("$key.hdf5", "w") do h5
    dset_data = d_create(h5, "data", datatype(Float32),
        dataspace(width, height, channels, batch_size * length(sources)))
    dset_label = d_create(h5, "label", datatype(Float32),
        dataspace(1, batch_size * length(sources)))

    for n = 1:length(sources)
      open("cifar-10-batches-bin/$(sources[n])") do f
        println("Processing $(sources[n])...")
        mat = readbytes(f, (1 + width*height*channels) * batch_size)
        mat = reshape(mat, 1+width*height*channels, batch_size)

        # random shuffle within batch
        rp = randperm(batch_size)

        label = convert(Array{Float32},mat[1, rp])
        # If I divide by 256 as in the MNIST example, then
        # training on the giving DNN gives me random
        # performance: objective function not changing,
        # and test performance is always 10%...
        # The same results could be observed when
        # running Caffe, as our HDF5 dataset is
        # compatible with Caffe.
        img = convert(Array{Float32},mat[2:end, rp])
        img = reshape(img, width, height, channels, batch_size)

        if key == "train"
          # only accumulate mean from the training data
          global mean_model
          mean_model = (batch_size*mean_model + sum(img, 4)) / (n*batch_size)
        end

        index = (n-1)*batch_size+1:n*batch_size
        dset_data[:,:,:,index] = img
        dset_label[:,index] = label
      end
    end

    # but apply mean subtraction for both training and testing data
    println("Subtracting the mean...")
    for n = 1:length(sources)
      index = (n-1)*batch_size+1:n*batch_size
      dset_data[:,:,:,index] = dset_data[:,:,:,index] .- mean_model
    end
  end
end


text:
# Copyright 2019 Tobias Frilling
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
    ABCcurve(data::Vector{<:Real})

Construct an ABCcurve from the given data.

# Fields
- `effort::AbstractRange`: The effort as equidistant steps between 0 and 1.
- `yield::Vector`: The cumulative relative yield.
- `interpolation::AbstractInterpolation`: A cubic spline interpolation of the curve.
"""
struct ABCcurve
    effort::AbstractRange{Float64}
    yield::Vector{<:Real}
    interpolation::AbstractInterpolation{<:Real,1,BSpline{Cubic{Line{OnGrid}}}}

    function ABCcurve(data::Vector{<:Real})
        @assert all(x->x >= 0, data)
        sorted_data = sort(data, rev = true)
        n = length(sorted_data)

        # (0,0) needs to be first point
        effort = (0:n) / n
        yield = pushfirst!(cumsum(sorted_data) / sum(sorted_data), 0)
        interpolation = CubicSplineInterpolation(effort, yield)

        new(effort, yield, interpolation)
    end
end

# text representation for REPL or Jupyter
function Base.show(io::IO, ::MIME"text/plain", curve::ABCcurve)
    print(io, curve.effort |> length, "-element ABCcurve")
    limit::Bool = get(io, :limit, false)
    if limit
        rows = displaysize(io)[1] - 3
        rows < 2 && (print(io, " …"); return)
        rows -= 2
    else
        rows = typemax(Int)
    end

    show_part(iter) = for v in iter
        print(io, "\n ")
        show(io, v)
    end

    zipped_curve = collect(zip(curve.effort, curve.yield))

    if rows < length(zipped_curve)
        show_part(zipped_curve[1:div(rows-1,2)])
        print(io, "\n ⋮")
        show_part(zipped_curve[end-div(rows-1,2):end])
    else
        show_part(zipped_curve)
    end
end

# plotting
@recipe function f(curve::ABCcurve; comparison=true)
    xlims --> (0,1)
    ylims --> (0,1)
    legend --> :bottomright
    ratio --> 1

    seriestype := :path
    markershape := :none

    @series begin
        label := "Data"
        linecolor := :blue
        linewidth --> 2
        curve.effort, curve.yield
    end
    if comparison
        @series begin
            label := "Uniform"
            linecolor := :green
            linewidth --> 1
            linealpha := 0.3
            curve.effort, map(x -> -x^2 + 2x, curve.effort)
        end
        @series begin
            label := "Identity"
            linecolor := :magenta
            linealpha := 0.3
            x -> x
        end
    end
end


text:
using NonUniformRandomVariateGeneration
using Test

import Random.seed!

include("goodness_of_fit.jl")

@testset "Binomial tests" begin
  include("binomial_test.jl")
end

@testset "Multinomial tests" begin
  include("multinomial_test.jl")
end

@testset "Poisson tests" begin
  include("poisson_test.jl")
end

@testset "Gamma tests" begin
  include("gamma_test.jl")
end

@testset "Beta tests" begin
  include("beta_test.jl")
end

@testset "Uniform tests" begin
  include("uniform_test.jl")
end

@testset "Categorical tests" begin
  include("categorical_test.jl")
end


text:
# Model
include_model("rocket")
model_nom = free_time_model(model_nominal)

# Horizon
T = 41

# Bounds
ul, uu = control_bounds(model_nom, T, [-5.0; 0.0; 0.01], [5.0; 100.0; 1.0])

# Initial and final states
x1 = [15.0;
	  model_nom.l1 + 10.0;
	  -45.0 * pi / 180.0;
	  -10.0;
	  -10.0;
	  -1.0 * pi / 180.0]

xT = [0.0;
	  model_nom.l1;
	  0.0;
	  0.0;
	  0.0;
	  0.0]

# Radius of landing pad
r_pad = 0.25

# xl <= x <= xl
_xl = -Inf * ones(model_nom.n)
_xl[2] = model_nom.l1
_xu = Inf * ones(model_nom.n)

xl, xu = state_bounds(model_nom, T, _xl, _xu, x1 = x1, xT = xT)

xl[T][1] = -1.0 * r_pad
xu[T][1] = r_pad
xl[T][2] = xT[2] - 0.001
xu[T][2] = xT[2] + 0.001
xl[T][3] = -1.0 * pi / 180.0
xu[T][3] = 1.0 * pi / 180.0
xl[T][4] = -0.001
xu[T][4] = 0.001
xl[T][5] = -0.001
xu[T][5] = 0.001
xl[T][6] = -0.01 * pi / 180.0
xu[T][6] = 0.01 * pi / 180.0

# Objective
obj = quadratic_time_tracking_objective(
	[(t != T ? Diagonal([1.0; 10.0; 1.0; 1.0; 10.0; 1.0])
		: Diagonal([10.0; 100.0; 10.0; 10.0; 100.0; 10.0])) for t = 1:T],
	[Diagonal([1.0e-1*ones(2); 0.0]) for t = 1:T-1],
    [xT for t = 1:T],
	[zeros(model_nom.m) for t = 1:T-1],
	1.0)

# Time step constraints
include_constraints("free_time")
con_free_time = free_time_constraints(T)

# Problem
prob_nominal = trajectory_optimization_problem(model_nom,
					obj,
					T,
                    xl = xl,
                    xu = xu,
                    ul = ul,
                    uu = uu,
                    con = con_free_time)

# Trajectory initialization
x0 = linear_interpolation(x1, xT, T) # linear interpolation on state
u0 = [ones(model_nom.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob_nominal)

# Solve
if false # set to true to reoptimize
	include_snopt()
    @time z̄_nom , info = solve(prob_nominal, copy(z0),
		nlp = :SNOPT7,
		time_limit = 60 * 10)
    @save joinpath(@__DIR__, "sol_nom.jld2") z̄_nom
else
    println("Loading solution...")
    @load joinpath(@__DIR__, "sol_nom.jld2") z̄_nom
end

x̄_nom, ū_nom = unpack(z̄_nom, prob_nominal.prob)
@show sum([ū_nom[t][end] for t = 1:T-1]) # ~2.72


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct NLSolversFG!{F<:GradientFunction}
    gradfunc::F
end

NLSolversBase.only_fg!(f::GradientFunction) = NLSolversBase.only_fg!(NLSolversFG!(f))

function (gf::NLSolversFG!)(val_f::Real, grad_f::Any, v::Any)
    return gf.gradfunc(!, grad_f, v)
end

function (gf::NLSolversFG!)(val_f::Nothing, grad_f::Any, v::Any)
    gf.gradfunc(!, grad_f, v)
    return Nothing
end

function _bat_findmode_impl_optim(rng::AbstractRNG, target::AnySampleable, algorithm::AbstractModeEstimator)
    transformed_density, trafo = transform_and_unshape(algorithm.trafo, target)

    initalg = apply_trafo_to_init(trafo, algorithm.init)
    x_init = collect(bat_initval(rng, transformed_density, initalg).result)

    f = negative(logdensityof(transformed_density))
    optim_result = _run_optim(f, x_init, algorithm)
    r_optim = Optim.MaximizationWrapper(optim_result)
    transformed_mode = Optim.minimizer(r_optim.res)
    result_mode = inverse(trafo)(transformed_mode)

    dummy_f_x = f(x_init) # ToDo: Avoid recomputation
    trace_trafo = StructArray(;_neg_opt_trace(optim_result, x_init, dummy_f_x) ...)

    (result = result_mode, result_trafo = transformed_mode, trafo = trafo, trace_trafo = trace_trafo, info = r_optim)
end

# Wrapper for type stability of optimize result (why does this work?):
function _optim_optimize(f, x0::AbstractArray, method::Optim.AbstractOptimizer, options = Optim.Options(); kwargs...)
    Optim.optimize(f, x0, method, options; kwargs...)
end

function _neg_opt_trace(
    @nospecialize(optim_result::Optim.MultivariateOptimizationResults),
    dummy_x::AbstractVector{<:Real}, dummy_f_x::Real
)
    trc = Optim.trace(optim_result)
    tr_len = length(eachindex(trc))
    nd = length(eachindex(dummy_x))

    v = nestedview(similar(dummy_x, nd, tr_len))
    foreach((a,b) -> a[:] = b, v, Optim.x_trace(optim_result))

    logd = similar(dummy_x, typeof(dummy_f_x), tr_len)
    logd[:] = - Optim.f_trace(optim_result)

    if optim_result isa Optim.MultivariateOptimizationResults{<:Optim.ZerothOrderOptimizer}
        (v = v, logd = logd)
    else
        grad_logd = nestedview(similar(dummy_x, nd, tr_len))
        foreach((a,b) -> a[:] = -b.metadata["g(x)"], grad_logd, trc)
        (v = v, logd = logd, grad_logd = grad_logd)
    end
end

function _neg_opt_trace(
    @nospecialize(optim_result::Optim.MultivariateOptimizationResults{<:Optim.NelderMead}),
    dummy_x::AbstractVector{<:Real}, dummy_f_x::Real
)
    trc = Optim.trace(optim_result)
    tr_len = length(eachindex(trc))
    nd = length(eachindex(dummy_x))

    v = nestedview(similar(dummy_x, nd, tr_len))
    foreach((a,b) -> a[:] = b, v, Optim.centroid_trace(optim_result))

    (;v = v)
end

"""
    struct MaxDensityNelderMead <: AbstractModeEstimator

Estimates the mode of a probability density using Nelder-Mead optimization
(currently via [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl),
subject to change).

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct MaxDensityNelderMead{
    TR<:AbstractDensityTransformTarget,
    IA<:InitvalAlgorithm
} <: AbstractModeEstimator
    trafo::TR = NoDensityTransform()
    init::IA = InitFromTarget()
end
export MaxDensityNelderMead

function _run_optim(f::Function, x_init::AbstractArray{<:Real}, algorithm::MaxDensityNelderMead)
    opts = Optim.Options(store_trace = true, extended_trace=true)
    _optim_optimize(f, x_init, Optim.NelderMead(), opts)
end

bat_findmode_impl(rng::AbstractRNG, target::AnySampleable, algorithm::MaxDensityNelderMead) = _bat_findmode_impl_optim(rng, target, algorithm)

"""
    struct MaxDensityLBFGS <: AbstractModeEstimator

Estimates the mode of a probability density using LBFGS optimization
(currently via [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl),
subject to change).

The gradient of the target density is computed via auto-differentiation.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct MaxDensityLBFGS{
    TR<:AbstractDensityTransformTarget,
    IA<:InitvalAlgorithm
} <: AbstractModeEstimator
    trafo::TR = PriorToGaussian()
    init::IA = InitFromTarget()
end
export MaxDensityLBFGS

bat_findmode_impl(rng::AbstractRNG, target::AnySampleable, algorithm::MaxDensityLBFGS) = _bat_findmode_impl_optim(rng, target, algorithm)

function _run_optim(f::Function, x_init::AbstractArray{<:Real}, algorithm::MaxDensityLBFGS)
    fg = valgradof(f)
    opts = Optim.Options(store_trace = true, extended_trace=true)
    _optim_optimize(Optim.only_fg!(fg), x_init, Optim.LBFGS(), opts)
end


text:
module Agent

using JLD, App
include("flex.jl")
include("util.jl")

global vocab = load(joinpath(App.config.server_document_root, "models/vocab_sail.jld"), "vocab")

function loadmodel(fname; flex=true)
    w = Dict()
    d = load(fname, "weights")
    for k in keys(d)
        if flex && startswith(k, "filter")
            w[k] = d[k]
        else
            w[k] = d[k]
        end
    end
    return w
end

global models = [loadmodel(joinpath(App.config.server_document_root, "models/sail_vtest_jelly2_"*"$i"*"_vtest_grid_l.jld")) for i=1:10]
#global models = [loadmodel(joinpath(App.config.server_document_root, "models/sailx_base_"*"$i"*".jld")) for i=1:10]

function predict(instruction, startpos, navimap)
    global models
    global vocab

    ins_text = split(strip(instruction), "_")
    words = ins_arr(vocab, ins_text)
    args = Dict()
    args["hidden"] = 120
    args["embed"] = 120
    args["limactions"] = 35
    args["bs"] = 1
    args["wvecs"] = false
    args["greedy"] = true
    args["percp"] = true
    args["preva"] = true
    args["worldatt"] = 100
    args["attinwatt"] = 0
    args["att"] = false
    args["inpout"] = true
    args["prevaout"] = false
    args["beam"] = true
    args["attout"] = false
    args["beamsize"] = 10
    args["encoding"] = "grid"
    args["atype"] = Array{Float32}

    predict_beam(models, words, navimap, startpos; args=args)
end

end

max_stars_repo_name: ozanarkancan/navi-app

max_issues_repo_name: ozanarkancan/navi-app

max_forks_repo_name: ozanarkancan/navi-app


text:
using CSV
using DataFrames
using Test

using Stonks.Stores: WriteOperation, WriteTransaction, list_partition_nesting
using Stonks.Stores: backup, cleanup, commit, execute, rollback, write, write_tmp

include("test_utils.jl")

@testset "Store writes" begin
  symbols = ["AAPL", "IBM"]
  data = fake_stock_data(7, last_workday(), symbols)
  dest = joinpath(@__DIR__, "data/test_stonks")
  tr = @chain symbols begin
    map(s -> WriteOperation(filter(:symbol => ==(s), data), "csv", "symbol=$s"), _)
    WriteTransaction(_, "csv", dest)
  end
  writer_csv(df::AbstractDataFrame, path::String) = CSV.write(path, df)
  get_type_param(::Vector{S}) where {S} = S

  @testset "Write using a WriteOperation" begin
    path = joinpath(@__DIR__, "data/test_stonks")
    @chain begin
      WriteOperation(data, "csv")
      write(_, path, writer_csv)
    end
    @test isfile("$path/data.csv")
    @test filesize("$path/data.csv") > 0
    rm(path; recursive=true, force=true)
  end

  @testset "Write transaction individual functions" begin
    write_tmp(tr, writer_csv)
    @test sort(readdir("$dest/.tmp")) == ["symbol=AAPL", "symbol=IBM"]
    @test filesize("$dest/.tmp/symbol=IBM") > 0
    @test filesize("$dest/.tmp/symbol=AAPL") > 0

    commit(tr)
    dirs = @chain readdir(dest) filter(x -> !startswith(x, "."), _) sort
    @test list_partition_nesting(dest) == ["symbol"]
    @test dirs == ["symbol=AAPL", "symbol=IBM"]
    @test filesize("$dest/symbol=IBM") > 0
    @test filesize("$dest/symbol=AAPL") > 0

    has_backup, bkps = backup(tr)
    @test sort(readdir("$dest/.bkp")) == ["symbol=AAPL", "symbol=IBM"]
    @test sort(readdir("$dest/.bkp")) == ["symbol=AAPL", "symbol=IBM"]

    foreach(x -> rm("$dest/symbol=$x"; recursive=true, force=true), ["AAPL", "IBM"])
    rollback(bkps)
    dirs = @chain readdir(dest) filter(x -> !startswith(x, "."), _) sort
    @test dirs == ["symbol=AAPL", "symbol=IBM"]

    cleanup(tr)
    @test isdir("$dest/.tmp") == false
    @test isdir("$dest/.bkp") == false

    rm(dest; recursive=true, force=true)
  end

  @testset "Write transaction execute method" begin
    status, err = execute(tr, writer_csv)
    @test isnothing(err)
    @test isdir(dest)
    @test isdir("$dest/.bkp") == false
    @test isdir("$dest/.tmp") == false
    rm(dest; recursive=true, force=true)
  end
end


text:
function create_trigonometric_blocks(cmp_bus, ref_bus, t_buses, va_fr, va_to; extra_block_needed::Bool=false)

    vat = []
    vaf_ref = [0.0, -2.0943951023931953, 2.0943951023931953] # angle of reference buses
    vaf_block = cmp_bus == ref_bus ? (x, fc, fd) -> vaf_ref[fc] - vaf_ref[fd] : (x, fc, fd) -> x[va_fr["$fc"]] - x[va_fr["$fd"]]

    if cmp_bus == ref_bus
        vat_block = (x, fc, br, td, d) -> vaf_ref[fc] - x[va_to[br]["$td"]]
        if extra_block_needed  extra_block = (x, fc) -> vaf_ref[fc] end
    elseif ref_bus ∈ t_buses
        for bus in t_buses
            bus != ref_bus ? push!(vat, (x, br, td) -> [x[va_to[br]["$td"]], x[va_to[br]["$td"]], x[va_to[br]["$td"]] ]) : push!(vat, (x, br, td) -> [0.0, -2.9044, 2.0944])
        end
        vat_block = (x, fc, br, td, d) -> x[va_fr["$fc"]] - vat[br](x, br, td)[d]
        if extra_block_needed  extra_block = (x, fc) -> x[va_fr["$fc"]] end
    else
        vat_block = (x, fc, br, td, d) -> x[va_fr["$fc"]] - x[va_to[br]["$td"]]
        if extra_block_needed  extra_block = (x, fc) -> x[va_fr["$fc"]] end
    end
    extra_block_needed ? (return vaf_block, vat_block, extra_block) : (return vaf_block, vat_block)
end

function build_measurement_function_array(data::Dict, variable_dict::Dict)

    functions = []
    ref_bus = [bus["index"] for (_, bus) in data["bus"] if bus["bus_type"]==3][1]

    for (m, meas) in data["meas"]
        add_h_function!(meas["var"], m, data, ref_bus, variable_dict, functions)
    end
    return functions
end

function add_h_function!(qmeas::Symbol, m::String, data::Dict, ref_bus::Int64, variable_dict::Dict, functions::Array)

    meas = data["meas"][m]

    if qmeas ∈ [:pd, :qd, :pg, :qg, :p, :q]
        input = qmeas ∈ [:pd, :qd, :pg, :qg] ? build_multibranch_input(qmeas, meas, data, variable_dict) : build_singlebranch_input(meas, data, variable_dict)
        eval(Meta.parse(string(qmeas)[1]*"inj_function"))(functions, ref_bus, input[1], input[2], input[3], input[4], input[5], input[6], input[7], input[8], input[9], input[10], input[11], input[12], input[13])
    elseif qmeas ∈ [:cmd, :cmg, :cm]
        input = qmeas ∈ [:cmd, :cmg] ? build_multibranch_input(qmeas, meas, data, variable_dict) : build_singlebranch_input(meas, data, variable_dict)
        cm_function(functions, ref_bus, input[1], input[2], input[3], input[4], input[5], input[6], input[7], input[8], input[9], input[10], input[11], input[12], input[13])
    elseif qmeas ∈ [:crd, :cid, :crg, :cig, :cr, :ci]
        input = qmeas ∈ [:crd, :cid, :crg, :cig] ? build_multibranch_input(qmeas, meas, data, variable_dict) : build_singlebranch_input(meas, data, variable_dict)
        eval(Meta.parse(string(qmeas)[1:2]*"_function"))(functions, ref_bus, input[1], input[2], input[3], input[4], input[5], input[6], input[7], input[8], input[9], input[10], input[11], input[12], input[13])
    elseif qmeas ∈ [:vm, :va]
        v_bus = variable_dict[string(qmeas)]["$(meas["cmp_id"])"]
        for i in collect(keys(v_bus))
            push!(functions, x->x[v_bus["$i"]])
        end
    else
        error("Measured quantity $qmeas of measurement $m not recognized.")
    end
end
"""
Measurement function for both (active) power flow measurement and (reactive) power injections measurement
"""
function pinj_function(functions, ref_bus, cmp_bus, vm_fr, va_fr, conn_branches, f_conn, t_conn, g, b, G_fr, B_fr, t_buses, vm_to, va_to)

    vaf_block, vat_block = create_trigonometric_blocks(cmp_bus, ref_bus, t_buses, va_fr, va_to)

    br_func = Matrix{Any}(undef, (length(conn_branches), 3))
    for br in 1:length(conn_branches)
        for (c, (fc, _)) in enumerate(zip(f_conn[br], t_conn[br]))
            br_func[br, fc] = x-> (g[br][c,c]+G_fr[br][c,c])*x[vm_fr["$fc"]]^2 +
                                   sum( (g[br][c,d]+G_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*cos(vaf_block(x, fc, fd)) +
                                        (b[br][c,d]+B_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*sin(vaf_block(x, fc, td))
                                        for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br])) if d!=c; init=0 ) +
                                  sum( -g[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*cos(vat_block(x, fc, br, td, d))
                                       -b[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*sin(vat_block(x, fc, br, td, d))
                                       for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br])))
        end
    end
    for i in 1:size(br_func)[1]
        for j in 1:size(br_func)[2]
            if !isassigned(br_func, i, j) br_func[i,j] = x->0.0 end
        end
    end
    for c in 1:maximum([length(t_conn[br]) for br in 1:length(conn_branches)])
        if maximum([length(t_conn[br]) for br in 1:length(conn_branches)]) < 3
            fc = f_conn[1][c]
            push!(functions, x->sum([br_func[i,fc](x) for i in 1:size(br_func)[1]]))
        else
            push!(functions, x->sum([br_func[i,c](x) for i in 1:size(br_func)[1]]))
        end
    end
end
"""
Measurement function for both (reactive) power flow measurement and (reactive) power injections measurement
"""
function qinj_function(functions, ref_bus, cmp_bus, vm_fr, va_fr, conn_branches, f_conn, t_conn, g, b, G_fr, B_fr, t_buses, vm_to, va_to)

    vaf_block, vat_block = create_trigonometric_blocks(cmp_bus, ref_bus, t_buses, va_fr, va_to)

    br_func = Matrix{Any}(undef, (length(conn_branches), 3))
    for br in 1:length(conn_branches)
        for (c, (fc, _)) in enumerate(zip(f_conn[br], t_conn[br]))
            br_func[br, fc] = x-> ( -(b[br][c,c]+B_fr[br][c,c])*x[vm_fr["$fc"]]^2
                                -sum( (b[br][c,d]+B_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*cos(vaf_block(x, fc, fd))
                                        -(g[br][c,d]+G_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*sin(vaf_block(x, fc, fd))
                                        for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br])) if d!=c; init=0 )
                                -sum( -b[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*cos(vat_block(x, fc, br, td, d))
                                        +g[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*sin(vat_block(x, fc, br, td, d))
                                        for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br])))
            )# closepush
        end
    end
    for i in 1:size(br_func)[1]
        for j in 1:size(br_func)[2]
            if !isassigned(br_func, i, j) br_func[i,j] = x->0.0 end
        end
    end
    for c in 1:maximum([length(t_conn[br]) for br in 1:length(conn_branches)])
        if maximum([length(t_conn[br]) for br in 1:length(conn_branches)]) < 3
            fc = f_conn[1][c]
            push!(functions, x->sum([br_func[i,fc](x) for i in 1:size(br_func)[1]]))
        else
            push!(functions, x->sum([br_func[i,c](x) for i in 1:size(br_func)[1]]))
        end
    end
end
"""
Current magnitude function
"""
function cm_function(functions, ref_bus, cmp_bus, vm_fr, va_fr, conn_branches, f_conn, t_conn, g, b, G_fr, B_fr, t_buses, vm_to, va_to)

    vaf_block, vat_block = create_trigonometric_blocks(cmp_bus, ref_bus, t_buses, va_fr, va_to)

    br_func = Matrix{Any}(undef, (length(conn_branches), 3))
    for br in 1:length(conn_branches)
        for (c, (fc, _)) in enumerate(zip(f_conn[1], t_conn[1]))
            br_func[br, fc] = x-> (sqrt(         # P^2
                                        ( -(b[br][c,c]+B_fr[br][c,c])*x[vm_fr["$fc"]]^2
                                            -sum( (b[br][c,d]+B_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*cos(vaf_block(x, fc, fd))
                                            -(g[br][c,d]+G_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*sin(vaf_block(x, fc, fd))
                                                for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br])) if d!=c; init=0 )
                                            -sum( -b[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*cos(vat_block(x, fc, br, td, d))
                                            +g[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*sin(vat_block(x, fc, br, td, d))
                                                for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br]))) )^2 +
                                                # Q^2
                                        ( (g[br][c,c]+G_fr[br][c,c])*x[vm_fr["$fc"]]^2 +
                                            sum( (g[br][c,d]+G_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*cos(vaf_block(x, fc, fd)) +
                                                (b[br][c,d]+B_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*sin(vaf_block(x, fc, fd))
                                                for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br])) if d!=c; init=0 ) +
                                            sum( -g[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*cos(vat_block(x, fc, br, td, d))
                                                -b[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*sin(vat_block(x, fc, br, td, d))
                                                for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br]))) )^2 )
                                                # / Vi
                                                / x[vm_fr["$fc"]] )
        end
    end
    for i in 1:size(br_func)[1]
        for j in 1:size(br_func)[2]
            if !isassigned(br_func, i, j) br_func[i,j] = x->0.0 end
        end
    end
    for c in 1:maximum([length(t_conn[br]) for br in 1:length(conn_branches)])
        if maximum([length(t_conn[br]) for br in 1:length(conn_branches)]) < 3
            fc = f_conn[1][c]
            push!(functions, x->sum([br_func[i,fc](x) for i in 1:size(br_func)[1]]))
        else
            push!(functions, x->sum([br_func[i,c](x) for i in 1:size(br_func)[1]]))
        end
    end
end
"""
Reactive current function
"""
function cr_function(functions, ref_bus, cmp_bus, vm_fr, va_fr, conn_branches, f_conn, t_conn, g, b, G_fr, B_fr, t_buses, vm_to, va_to)

    vaf_block, vat_block, extra_block = create_trigonometric_blocks(cmp_bus, ref_bus, t_buses, va_fr, va_to, extra_block_needed = true)

    br_func = Matrix{Any}(undef, (length(conn_branches), 3))
    for br in 1:length(conn_branches)
        for (c, (fc, tc)) in enumerate(zip(f_conn[br], t_conn[br]))
            br_func[br,fc] = x-> (          #  Q * sin(Θ)
                                        ( -(b[br][c,c]+B_fr[br][c,c])*x[vm_fr["$fc"]]^2
                                        -sum( (b[br][c,d]+B_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*cos(vaf_block(x, fc, fd))
                                        -(g[br][c,d]+G_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*sin(vaf_block(x, fc, fd))
                                            for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br])) if d!=c; init=0 )
                                        -sum( -b[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*cos(vat_block(x, fc, br, td, d))
                                        +g[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*sin(vat_block(x, fc, br, td, d))
                                            for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br]))))*sin(extra_block(x, fc))+
                                            # P * cos(Θ)
                                        ( (g[br][c,c]+G_fr[br][c,c])*x[vm_fr["$fc"]]^2 +
                                        sum( (g[br][c,d]+G_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*cos(vaf_block(x, fc, fd)) +
                                            (b[br][c,d]+B_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*sin(vaf_block(x, fc, fd))
                                            for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br])) if d!=c; init=0 ) +
                                        sum( -g[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*cos(vat_block(x, fc, br, td, d))
                                            -b[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*sin(vat_block(x, fc, br, td, d))
                                            for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br]))))*cos(extra_block(x, fc))
                                            # / |U|
                                            / x[vm_fr["$fc"]]
                )
        end
    end
    for i in 1:size(br_func)[1]
        for j in 1:size(br_func)[2]
            if !isassigned(br_func, i, j) br_func[i,j] = x->0.0 end
        end
    end
    for c in 1:maximum([length(t_conn[br]) for br in 1:length(conn_branches)])
        if maximum([length(t_conn[br]) for br in 1:length(conn_branches)]) < 3
            fc = f_conn[1][c]
            push!(functions, x->sum([br_func[i,fc](x) for i in 1:size(br_func)[1]]))
        else
            push!(functions, x->sum([br_func[i,c](x) for i in 1:size(br_func)[1]]))
        end
    end
end

function ci_function(functions, ref_bus, cmp_bus, vm_fr, va_fr, conn_branches, f_conn, t_conn, g, b, G_fr, B_fr, t_buses, vm_to, va_to)

    vaf_block, vat_block, extra_block = create_trigonometric_blocks(cmp_bus, ref_bus, t_buses, va_fr, va_to, extra_block_needed = true)

    br_func = Matrix{Any}(undef, (length(conn_branches), 3))
    for br in 1:length(conn_branches)
        for (c, (fc, tc)) in enumerate(zip(f_conn[br], t_conn[br]))
            br_func[br,fc] = x-> (          #  - Q * cos(Θ)
                                - ( -(b[br][c,c]+B_fr[br][c,c])*x[vm_fr["$fc"]]^2
                                -sum( (b[br][c,d]+B_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*cos(vaf_block(x, fc, fd))
                                -(g[br][c,d]+G_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*sin(vaf_block(x, fc, fd))
                                    for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br])) if d!=c; init=0 )
                                -sum( -b[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*cos(vat_block(x, fc, br, td, d))
                                +g[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*sin(vat_block(x, fc, br, td, d))
                                    for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br]))))*cos(extra_block(x, fc))+
                                    # + P * sin(Θ)
                                ( (g[br][c,c]+G_fr[br][c,c])*x[vm_fr["$fc"]]^2 +
                                sum( (g[br][c,d]+G_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*cos(vaf_block(x, fc, fd)) +
                                    (b[br][c,d]+B_fr[br][c,d])*x[vm_fr["$fc"]]*x[vm_fr["$fd"]]*sin(vaf_block(x, fc, fd))
                                    for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br])) if d!=c; init=0 ) +
                                sum( -g[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*cos(vat_block(x, fc, br, td, d))
                                    -b[br][c,d]*x[vm_fr["$fc"]]*x[vm_to[br]["$td"]]*sin(vat_block(x, fc, br, td, d))
                                    for (d, (fd, td)) in enumerate(zip(f_conn[br], t_conn[br]))))*sin(extra_block(x, fc))
                                    # / |U|
                                    / x[vm_fr["$fc"]]
                )
        end
    end
    for i in 1:size(br_func)[1]
        for j in 1:size(br_func)[2]
            if !isassigned(br_func, i, j) br_func[i,j] = x->0.0 end
        end
    end
    for c in 1:maximum([length(t_conn[br]) for br in 1:length(conn_branches)])
        if maximum([length(t_conn[br]) for br in 1:length(conn_branches)]) < 3
            fc = f_conn[1][c]
            push!(functions, x->sum([br_func[i,fc](x) for i in 1:size(br_func)[1]]))
        else
            push!(functions, x->sum([br_func[i,c](x) for i in 1:size(br_func)[1]]))
        end
    end
end


text:
# HACK: work around Pkg.jl#2500
test_project = Base.active_project()
preferences_file = joinpath(dirname(@__DIR__), "LocalPreferences.toml")
test_preferences_file = joinpath(dirname(test_project), "LocalPreferences.toml")
if isfile(preferences_file) && !isfile(test_preferences_file)
    cp(preferences_file, test_preferences_file)
end

using Enzyme
using Test
using FiniteDifferences
using ForwardDiff
using Statistics

using Enzyme_jll
@info "Testing against" Enzyme_jll.libEnzyme

# Test against FiniteDifferences
function test_scalar(f, x; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)
    ∂x, = autodiff(f, Active, Active(x))
    if typeof(x) <: Complex
    else
      @test isapprox(∂x, fdm(f, x); rtol=rtol, atol=atol, kwargs...)
    end

    rm = ∂x
    if typeof(x) <: Integer
        x = Float64(x)
    end
    ∂x, = fwddiff(f, Duplicated(x, one(typeof(x))))
    if typeof(x) <: Complex
      @test ∂x ≈ rm
    else
      @test isapprox(∂x, fdm(f, x); rtol=rtol, atol=atol, kwargs...)
    end
end

include("abi.jl")
include("typetree.jl")

@testset "Internal tests" begin
    f(x) = 1.0 + x
    thunk_a = Enzyme.Compiler.thunk(f, nothing, Active, Tuple{Active{Float64}})
    thunk_b = Enzyme.Compiler.thunk(f, nothing, Const, Tuple{Const{Float64}})
    thunk_c = Enzyme.Compiler.thunk(f, nothing, Active{Float64}, Tuple{Active{Float64}})
    @test thunk_a.adjoint !== thunk_b.adjoint
    @test thunk_c.adjoint === thunk_a.adjoint

    @test thunk_a(Active(2.0), 1.0) == (1.0,)
    @test thunk_a(Active(2.0), 2.0) == (2.0,)
    @test thunk_b(Const(2.0)) === ()

    forward, pullback = Enzyme.Compiler.thunk(f, nothing, Active, Tuple{Active{Float64}}, Val(Enzyme.API.DEM_ReverseModeGradient))
    # @test thunk_split.primal !== C_NULL
    # @test thunk_split.primal !== thunk_split.adjoint
    # @test thunk_a.adjoint !== thunk_split.adjoint
end

@testset "Reflection" begin
    Enzyme.Compiler.enzyme_code_typed(Active, Tuple{Active{Float64}}) do x
        x ^ 2
    end
    f(x) = 1.0 + x
    sprint() do io
        Enzyme.Compiler.enzyme_code_native(io, f, Active, Tuple{Active{Float64}})
    end

    sprint() do io
        Enzyme.Compiler.enzyme_code_llvm(io, f, Active, Tuple{Active{Float64}})
    end
end

# @testset "Split Tape" begin
#     f(x) = x[1] * x[1]

#     thunk_split = Enzyme.Compiler.thunk(f, Tuple{Duplicated{Array{Float64,1}}}, Val(Enzyme.API.DEM_ReverseModeGradient))
#     @test thunk_split.primal !== C_NULL
#     @test thunk_split.primal !== thunk_split.adjoint
# end

@testset "Simple tests" begin
    f1(x) = 1.0 + x
    f2(x) = x*x
    @test autodiff(f1, Active, Active(1.0))[1] ≈ 1.0
    @test fwddiff(f1, DuplicatedNoNeed, Duplicated(1.0, 1.0))[1] ≈ 1.0
    @test fwddiff(f1, Duplicated, Duplicated(1.0, 1.0))[2] ≈ 1.0
    @test autodiff(f2, Active, Active(1.0))[1] ≈ 2.0
    @test fwddiff(f2, Duplicated(1.0, 1.0))[1] ≈ 2.0
    @test autodiff(tanh, Active, Active(1.0))[1] ≈ 0.41997434161402606939
    @test fwddiff(tanh, Duplicated(1.0, 1.0))[1] ≈ 0.41997434161402606939
    @test autodiff(tanh, Active, Active(1.0f0))[1] ≈ Float32(0.41997434161402606939)
    @test fwddiff(tanh, Duplicated(1.0f0, 1.0f0))[1] ≈ Float32(0.41997434161402606939)
    test_scalar(f1, 1.0)
    test_scalar(f2, 1.0)

    @test autodiff((x)->log(x), Active(2.0)) == (0.5,)
end

@testset "Duplicated" begin
    x = Ref(1.0)
    y = Ref(2.0)

    ∇x = Ref(0.0)
    ∇y = Ref(0.0)

    autodiff((a,b)->a[]*b[], Active, Duplicated(x, ∇x), Duplicated(y, ∇y))

    @test ∇y[] == 1.0
    @test ∇x[] == 2.0
end

@testset "Simple tests" begin
    g(x) = real((x + im)*(1 - im*x))
    @test first(autodiff(g, Active, Active(2.0))) ≈ 2.0
    @test first(fwddiff(g, Duplicated(2.0, 1.0))) ≈ 2.0
    @test first(autodiff(g, Active, Active(3.0))) ≈ 2.0
    @test first(fwddiff(g, Duplicated(3.0, 1.0))) ≈ 2.0
    test_scalar(g, 2.0)
    test_scalar(g, 3.0)
end

@testset "Taylor series tests" begin

# Taylor series for `-log(1-x)`
# eval at -log(1-1/2) = -log(1/2)
function euroad(f::T) where T
    g = zero(T)
    for i in 1:10^7
        g += f^i / i
    end
    return g
end

euroad′(x) = first(autodiff(euroad, Active, Active(x)))

@test euroad(0.5) ≈ -log(0.5) # -log(1-x)
@test euroad′(0.5) ≈ 2.0 # d/dx -log(1-x) = 1/(1-x)
test_scalar(euroad, 0.5)
end

@testset "Array tests" begin

    function arsum(f::Array{T}) where T
        g = zero(T)
        for elem in f
            g += elem
        end
        return g
    end

    inp = Float64[1.0, 2.0]
    dinp = Float64[0.0, 0.0]
    autodiff(arsum, Active, Duplicated(inp, dinp))
    @test inp ≈ Float64[1.0, 2.0]
    @test dinp ≈ Float64[1.0, 1.0]

    @test fwddiff(arsum, Duplicated(inp, dinp))[1] ≈ 2.0
end

@testset "Advanced array tests" begin
    function arsum2(f::Array{T}) where T
        return sum(f)
    end
    inp = Float64[1.0, 2.0]
    dinp = Float64[0.0, 0.0]
    autodiff(arsum2, Active, Duplicated(inp, dinp))
    @test inp ≈ Float64[1.0, 2.0]
    @test dinp ≈ Float64[1.0, 1.0]

    @test fwddiff(arsum2, Duplicated(inp, dinp))[1] ≈ 2.0
end

function grad_closure(f, x)
    function noretval(x,res)
        y = f(x)
        copyto!(res,y)
        return nothing
    end
    n = length(x)
    dx = zeros(n)
    y  = zeros(n)
    dy = zeros(n)
    dy[1] = 1.0

    autodiff(noretval, Duplicated(x,dx), Duplicated(y, dy))
    return dx
end

@testset "Closure" begin
    x = [2.0,6.0]
    dx = grad_closure(x->[x[1], x[2]], x)
    @test dx == [1.0, 0.0]
end

@testset "Bithacks" begin
    function fneg(x::Float64)
        xptr = reinterpret(Int64, x)
        y = Int64(-9223372036854775808)
        out = y ⊻ xptr;
        return reinterpret(Float64, out)
    end
    @test autodiff(fneg, Active, Active(2.0))[1] ≈ -1.0
    @test fwddiff(fneg, Duplicated(2.0, 1.0))[1] ≈ -1.0
    function expor(x::Float64)
        xptr = reinterpret(Int64, x)
        y = UInt64(4607182418800017408)
        out = y | xptr;
        return reinterpret(Float64, out)
    end
    @test autodiff(expor, Active, Active(0.42))[1] ≈ 4.0
    @test fwddiff(expor, Duplicated(0.42, 1.0))[1] ≈ 4.0
end

@testset "GC" begin
    function gc_alloc(x)  # Basically g(x) = x^2
        a = Array{Float64, 1}(undef, 10)
        for n in 1:length(a)
            a[n] = x^2
        end
        return mean(a)
    end
    @test autodiff(gc_alloc, Active, Active(5.0))[1] ≈ 10
    @test fwddiff(gc_alloc, Duplicated(5.0, 1.0))[1] ≈ 10

    # TODO (after BLAS)
    # A = Vector[2.0, 3.0]
    # B = Vector[4.0, 5.0]
    # dB = Vector[0.0, 0.0]
    # f = (X, Y) -> sum(X .* Y)
    # Enzyme.autodiff(f, Active, A, Duplicated(B, dB))

    function gc_copy(x)  # Basically g(x) = x^2
        a = x * ones(10)
        for n in 1:length(a)
            a[n] = x^2
        end
        return mean(a)
    end
    # TODO(wsmoses): Assertion failed: (pp->getNumUses() == 0), function eraseFictiousPHIs
    # @test Enzyme.autodiff(gc_copy, Active, Active(5.0))[1] ≈ 10
    # TODO: https://github.com/wsmoses/Enzyme/issues/393
    # @test Enzyme.fwddiff(gc_copy, Duplicated(5.0, 1.0))[1] ≈ 10
end

@testset "Compare against" begin
    x = 3.0
    fd = central_fdm(5, 1)(sin, x)

    @test fd ≈ ForwardDiff.derivative(sin, x)
    @test fd ≈ first(autodiff(sin, Active, Active(x)))
    @test fd ≈ first(fwddiff(sin, Duplicated(x, 1.0)))

    x = 0.2 + sin(3.0)
    fd = central_fdm(5, 1)(asin, x)

    @test fd ≈ ForwardDiff.derivative(asin, x)
    @test fd ≈ first(autodiff(asin, Active, Active(x)))
    @test fd ≈ first(fwddiff(asin, Duplicated(x, 1.0)))
    test_scalar(asin, x)

    function foo(x)
        a = sin(x)
        b = 0.2 + a
        c = asin(b)
        return c
    end

    x = 3.0
    fd = central_fdm(5, 1)(foo, x)

    @test fd ≈ ForwardDiff.derivative(foo, x)
    @test fd ≈ first(autodiff(foo, Active, Active(x)))
    @test fd ≈ first(fwddiff(foo, Duplicated(x, 1.0)))
    test_scalar(foo, x)

    # Input type shouldn't matter
    x = 3
    @test fd ≈ ForwardDiff.derivative(foo, x)
    @test fd ≈ first(autodiff(foo, Active, Active(x)))
    # They do matter for duplicated, which can't be auto promoted
    # @test fd ≈ first(fwddiff(foo, Duplicated(x, 1)))

    f74(a, c) = a * √c
    @test √3 ≈ first(autodiff(f74, Active, Active(2), 3))
    @test √3 ≈ first(fwddiff(f74, Duplicated(2.0, 1.0), 3))
end

"""
    J(ν, z) := ∑ (−1)^k / Γ(k+1) / Γ(k+ν+1) * (z/2)^(ν+2k)
"""
function mybesselj(ν, z, atol=1e-8)
    k = 0
    s = (z/2)^ν / factorial(ν)
    out = s
    while abs(s) > atol
        k += 1
        s *= (-1) / k / (k+ν) * (z/2)^2
        out += s
    end
    out
end
mybesselj0(z) = mybesselj(0, z)
mybesselj1(z) = mybesselj(1, z)

@testset "Bessel" begin
    autodiff(mybesselj, Active, Const(0), Active(1.0))
    autodiff(mybesselj, Active, 0, Active(1.0))
    fwddiff(mybesselj, Const(0), Duplicated(1.0, 1.0))
    fwddiff(mybesselj, 0, Duplicated(1.0, 1.0))
    @testset "besselj0/besselj1" for x in (1.0, -1.0, 0.0, 0.5, 10, -17.1,) # 1.5 + 0.7im)
        test_scalar(mybesselj0, x, rtol=1e-5, atol=1e-5)
        test_scalar(mybesselj1, x, rtol=1e-5, atol=1e-5)
    end
end

## https://github.com/JuliaDiff/ChainRules.jl/tree/master/test/rulesets
if !Sys.iswindows()
    include("packages/specialfunctions.jl")
end

@testset "Threads" begin
    cmd = `$(Base.julia_cmd()) --threads=1 --startup-file=no threads.jl`
   	@test success(pipeline(cmd, stderr=stderr, stdout=stdout))
    cmd = `$(Base.julia_cmd()) --threads=2 --startup-file=no threads.jl`
   	@test success(pipeline(cmd, stderr=stderr, stdout=stdout))
end

@testset "DiffTest" begin
    include("DiffTests.jl")

    n = rand()
    x, y = rand(5, 5), rand(26)
    A, B = rand(5, 5), rand(5, 5)

    # f returns Number
    @testset "Number to Number" for f in DiffTests.NUMBER_TO_NUMBER_FUNCS
        test_scalar(f, n)
    end

    # TODO(vchuravy/wsmoses): Enable these tests
    # for f in DiffTests.VECTOR_TO_NUMBER_FUNCS
    #     @test isa(f(y), Number)
    # end

    # for f in DiffTests.MATRIX_TO_NUMBER_FUNCS
    #     @test isa(f(x), Number)
    # end

    # for f in DiffTests.TERNARY_MATRIX_TO_NUMBER_FUNCS
    #     @test isa(f(A, B, x), Number)
    # end

    # # f returns Array

    # for f in DiffTests.NUMBER_TO_ARRAY_FUNCS
    #     @test isa(f(n), Array)
    # end

    # for f in DiffTests.ARRAY_TO_ARRAY_FUNCS
    #     @test isa(f(A), Array)
    #     @test isa(f(y), Array)
    # end

    # for f in DiffTests.MATRIX_TO_MATRIX_FUNCS
    #     @test isa(f(A), Array)
    # end

    # for f in DiffTests.BINARY_MATRIX_TO_MATRIX_FUNCS
    #     @test isa(f(A, B), Array)
    # end

    # # f! returns Nothing

    # for f! in DiffTests.INPLACE_ARRAY_TO_ARRAY_FUNCS
    #     @test isa(f!(y, x), Nothing)
    # end

    # for f! in DiffTests.INPLACE_NUMBER_TO_ARRAY_FUNCS
    #     @test isa(f!(y, n), Nothing)
    # end

end

@testset "IO" begin

    function printsq(x)
        println(x)
        x*x
    end

    @test 4.6 ≈ first(autodiff(printsq, Active, Active(2.3)))
    @test 4.6 ≈ first(fwddiff(printsq, Duplicated(2.3, 1.0)))

    function tostring(x)
        string(x)
        x*x
    end

    @test 4.6 ≈ first(autodiff(tostring, Active, Active(2.3)))
    @test 4.6 ≈ first(fwddiff(tostring, Duplicated(2.3, 1.0)))
end

@testset "hmlstm" begin
    sigm(x)  = @fastmath 1 / (1 + exp(-x))
    @fastmath function hmlstm_update_c_scalar(z, zb, c, f, i, g)
        if z == 1.0f0 # FLUSH
            return sigm(i) * tanh(g)
        elseif zb == 0.0f0 # COPY
            return c
        else # UPDATE
            return sigm(f) * c + sigm(i) * tanh(g)
        end
    end

    N = 64
    Z = round.(rand(Float32, N))
    Zb = round.(rand(Float32, N))
    C = rand(Float32, N, N)
    F = rand(Float32, N, N)
    I = rand(Float32, N, N)
    G = rand(Float32, N, N)

    function broadcast_hmlstm(out, Z, Zb, C, F, I, G)
        out .= hmlstm_update_c_scalar.(Z, Zb, C, F, I, G)
        return nothing
    end

    ∇C = zeros(Float32, N, N)
    ∇F = zeros(Float32, N, N)
    ∇I = zeros(Float32, N, N)
    ∇G = zeros(Float32, N, N)

    # TODO(wsmoses): Check after updating Enzyme_jll
    # autodiff(broadcast_hmlstm, Const,
    #          Const(zeros(Float32, N, N)), Const(Z), Const(Zb),
    #          Duplicated(C, ∇C), Duplicated(F, ∇F), Duplicated(I, ∇I), Duplicated(G, ∇G))
    # fwddiff(broadcast_hmlstm, Const,
    #          Const(zeros(Float32, N, N)), Const(Z), Const(Zb),
    #          Duplicated(C, ∇C), Duplicated(F, ∇F), Duplicated(I, ∇I), Duplicated(G, ∇G))
end

genlatestsin(x)::Float64 = Base.invokelatest(sin, x)
function genlatestsinx(xp)
    x = @inbounds xp[1]
    @inbounds xp[1] = 0.0
    Base.invokelatest(sin, x)::Float64 + 1
end

function loadsin(xp)
    x = @inbounds xp[1]
    @inbounds xp[1] = 0.0
    sin(x)
end
function invsin(xp)
    xp = Base.invokelatest(convert, Vector{Float64}, xp)
    loadsin(xp)
end

@testset "generic" begin
    @test -0.4161468365471424 ≈ Enzyme.autodiff(genlatestsin, Active, Active(2.0))[1]
    @test -0.4161468365471424 ≈ Enzyme.fwddiff(genlatestsin, Duplicated(2.0, 1.0))[1]

    x = [2.0]
    dx = [0.0]
    Enzyme.autodiff(genlatestsinx, Active, Duplicated(x, dx))
    @test 0 ≈ x[1]
    @test -0.4161468365471424 ≈ dx[1]

    x = [2.0]
    dx = [0.0]
    Enzyme.autodiff(invsin, Active, Duplicated(x, dx))
    @test 0 ≈ x[1]
    @test -0.4161468365471424 ≈ dx[1]
end

@testset "invoke" begin
    @noinline apply(@nospecialize(func)) = func()

    function invtest(arr)
        function f()
           arr[1] *= 5.0
           nothing
        end
        apply(f)
    end

    x  = [2.0]
    dx = [1.0]

    Enzyme.autodiff(invtest, Duplicated(x, dx))

    @test 10.0 ≈ x[1]
    @test 5.0 ≈ dx[1]
end

@testset "broadcast" begin
    A = rand(10); B = rand(10); R = similar(A)
    dA = zero(A); dB = zero(B); dR = fill!(similar(R), 1)

    function foo_bc!(R, A, B)
        R .= A .+ B
        return nothing
    end

    autodiff(foo_bc!, Const, Duplicated(R, dR), Duplicated(A, dA), Duplicated(B, dB))

    # works since aliasing is "simple"
    autodiff(foo_bc!, Const, Duplicated(R, dR), Duplicated(R, dR), Duplicated(B, dB))

    A = rand(10,10); B = rand(10, 10)
    dA = zero(A); dB = zero(B); dR = fill!(similar(A), 1)

    autodiff(foo_bc!, Const, Duplicated(A, dR), Duplicated(transpose(A), transpose(dA)), Duplicated(B, dB))
end

@testset "DuplicatedReturn" begin
    moo(x) = fill(x, 10)

    @test_throws ErrorException autodiff(moo, Active(2.1))
    fo, = fwddiff(moo, Duplicated(2.1, 1.0))
    for i in 1:10
        @test 1.0 ≈ fo[i]
    end
end

@testset "GCPreserve" begin
    function f(x, y)
        GC.@preserve x y begin
            ccall(:memcpy, Cvoid,
                (Ptr{Float64},Ptr{Float64},Csize_t), x, y, 8)
        end
        nothing
    end
    autodiff(f, Duplicated([1.0], [0.0]), Duplicated([1.0], [0.0]))
    fwddiff(f, Duplicated([1.0], [0.0]), Duplicated([1.0], [0.0]))
end

@testset "GCPreserve2" begin
    function f!(a_out, a_in)
           a_out[1:end-1] .= a_in[2:end]
           return nothing
    end
    a_in = rand(4)
    a_out = a_in

    shadow_a_out = ones(4)
    shadow_a_in = shadow_a_out

    autodiff(f!, Const, Duplicated(a_out, shadow_a_out), Duplicated(a_in, shadow_a_in))

    @test shadow_a_in ≈ Float64[0.0, 1.0, 1.0, 2.0]
    @test shadow_a_out ≈ Float64[0.0, 1.0, 1.0, 2.0]

    fwddiff(f!, Const, Duplicated(a_out, shadow_a_out), Duplicated(a_in, shadow_a_in))

    @test shadow_a_in ≈ Float64[1.0, 1.0, 2.0, 2.0]
    @test shadow_a_out ≈ Float64[1.0, 1.0, 2.0, 2.0]
end

@testset "UndefVar" begin
    function f(x, y)
        if x
            undefinedfnthowmagic()
        end
        y
    end
    @test 1.0 ≈ autodiff(f, false, Active(2.14))[1]
    @test_throws Base.UndefVarError autodiff(f, true, Active(2.14))

    @test 1.0 ≈ fwddiff(f, false, Duplicated(2.14, 1.0))[1]
    @test_throws Base.UndefVarError fwddiff(f, true, Duplicated(2.14, 1.0))

    function foo(x, y)
        if x
            Threads.@threads for N in 1:5:20
                println("The number of this iteration is $N")
            end
        end
        y
    end
    @test 1.0 ≈ autodiff(foo, false, Active(2.14))[1]
    @test 1.0 ≈ fwddiff(foo, false, Duplicated(2.14, 1.0))[1]
end

@testset "Return GC error" begin
	t = 0.0

	function tobedifferentiated(cond, a)::Float64
		if cond
			t + t
		else
			0.0
		end
	end

	@test 0.0 ≈ autodiff(tobedifferentiated, true, Active(2.1))[1]
	@test 0.0 ≈ fwddiff(tobedifferentiated, true, Duplicated(2.1, 1.0))[1]

	function tobedifferentiated2(cond, a)::Float64
		if cond
			a + t
		else
			0.0
		end
	end

	@test 1.0 ≈ autodiff(tobedifferentiated2, true, Active(2.1))[1]
	@test 1.0 ≈ fwddiff(tobedifferentiated2, true, Duplicated(2.1, 1.0))[1]

    @noinline function copy(dest, p1, cond)
        bc = convert(Broadcast.Broadcasted{Nothing}, Broadcast.instantiate(p1))

        if cond
            return nothing
        end

        bc2 = Broadcast.preprocess(dest, bc)
        @inbounds    dest[1] = bc2[1]

        nothing
    end

    function mer(F, F_H, cond)
        p1 = Base.broadcasted(Base.identity, F_H)
        copy(F, p1, cond)

        # Force an apply generic
        flush(stdout)
        nothing
    end

    L_H = Array{Float64, 1}(undef, 2)
    L = Array{Float64, 1}(undef, 2)

    F_H = [1.0, 0.0]
    F = [1.0, 0.0]

    autodiff(mer, Duplicated(F, L), Duplicated(F_H, L_H), true)
    fwddiff(mer, Duplicated(F, L), Duplicated(F_H, L_H), true)
end

@testset "Split GC" begin
    @noinline function bmat(x)
        data = [x]
        return data
    end

    function f(x::Float64)
        @inbounds return bmat(x)[1]
    end
    @test 1.0 ≈ autodiff(f, Active(0.1))[1]
    @test 1.0 ≈ fwddiff(f, Duplicated(0.1, 1.0))[1]
end

@testset "Array Copy" begin
	F = [2.0, 3.0]

	dF = [0.0, 0.0]

	function copytest(F)
		F2 = copy(F)
		@inbounds F[1] = 1.234
		@inbounds F[2] = 5.678
		@inbounds F2[1] * F2[2]
	end
	autodiff(copytest, Duplicated(F, dF))
	@test F ≈ [1.234, 5.678]
	@test dF ≈ [3.0, 2.0]

    @test 31.0 ≈ fwddiff(copytest, Duplicated([2.0, 3.0], [7.0, 5.0]))[1]
end

@testset "No inference" begin
    c = 5.0
    @test 5.0 ≈ autodiff((A,)->c * A, Active, Active(2.0))[1]
    @test 5.0 ≈ fwddiff((A,)->c * A, Duplicated(2.0, 1.0))[1]
end

@testset "Type-instable capture" begin
    L = Array{Float64, 1}(undef, 2)

    F = [1.0, 0.0]

    function main()
        t = 0.0

        function cap(m)
            t = m
        end

        @noinline function inner(F, cond)
            if cond
                genericcall(F)
            end
        end

        function tobedifferentiated(F, cond)
            inner(F, cond)
            # Force an apply generic
            -t
            nothing
        end
        autodiff(tobedifferentiated, Duplicated(F, L), false)
        fwddiff(tobedifferentiated, Duplicated(F, L), false)
    end

    main()
end

@testset "Arrays are double pointers" begin
    @noinline function func_scalar(X)
        return X
    end

    function timsteploop_scalar(FH1)
        G = Float64[FH1]
        k1 = @inbounds func_scalar(G[1])
        return k1
    end
    @test Enzyme.autodiff(timsteploop_scalar, Active(2.0))[1] ≈ 1.0
    @test Enzyme.fwddiff(timsteploop_scalar, Duplicated(2.0, 1.0))[1] ≈ 1.0

    @noinline function func(X)
        return @inbounds X[1]
    end
    function timsteploop(FH1)
        G = Float64[FH1]
        k1 = func(G)
        return k1
    end
    @test Enzyme.autodiff(timsteploop, Active(2.0))[1] ≈ 1.0
    @test Enzyme.fwddiff(timsteploop, Duplicated(2.0, 1.0))[1] ≈ 1.0
end

@testset "Type" begin
    function foo(in::Ptr{Cvoid}, out::Ptr{Cvoid})
        markType(Float64, in)
        ccall(:memcpy,Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), out, in, 8)
    end

    x = [2.0]
    y = [3.0]
    dx = [5.0]
    dy = [7.0]

    @test markType(x) === nothing
    @test markType(zeros(Float32, 64)) === nothing
    @test markType(view(zeros(64), 16:32)) === nothing

    GC.@preserve x y begin
        foo(Base.unsafe_convert(Ptr{Cvoid}, x), Base.unsafe_convert(Ptr{Cvoid}, y))
    end

    GC.@preserve x y dx dy begin
      autodiff(foo,
                Duplicated(Base.unsafe_convert(Ptr{Cvoid}, x), Base.unsafe_convert(Ptr{Cvoid}, dx)),
                Duplicated(Base.unsafe_convert(Ptr{Cvoid}, y), Base.unsafe_convert(Ptr{Cvoid}, dy)))
    end
end

@testset "Exception" begin
    f_exc(x) = x'*x
    y = [1.0, 2.0]
    f_x = zero.(y)
    @test_throws Enzyme.CompilationException autodiff(f_exc, Duplicated(y, f_x))
end


text:
#!/usr/bin/env julia
using ArgParse
using CSV
using JSON
using Markdown

# define the command line interface
function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--root", "-r"
            help = "root of the v3 repository"
            arg_type = String
            default = joinpath("..", "..")
        "--track", "-t"
            help = "track slug"
            arg_type = String
            default = "julia"
        "--no-config-check"
            help = "disable config.json checks"
            action = :store_true
        "--no-directory-check"
            help = "disable exercise directory checks"
            action = :store_true
        "--no-extraction-check"
            help = "disable concept extraction checks"
            action = :store_true
    end

    return parse_args(ARGS, s)
end

"""
Check if all concepts used by exercises in config.json are defined in concepts.csv

Returns a dictionary of exercises containing undefined concepts: exercise slug ↦ undefined concepts
"""
function check_config_json(concepts, config)
    undefined_concepts = Dict{String,Vector{String}}() # exercise slug ↦ undefined concepts

    # union of concept and practice exercises
    exercises = config["exercises"]["concept"] ∪ config["exercises"]["practice"]

    for exercise in exercises
        # union of concepts and prerequisites
        ex_concepts = exercise["concepts"] ∪ exercise["prerequisites"]

        # vector to store all undefined concepts
        undefined = String[]
        for c in ex_concepts
            # if the concept is not in concepts.csv, push it to undefined
            if c ∉ concepts[!, :concept]
                push!(undefined, c)
            end
        end

        # if undefined is empty, all concepts are defined
        # if not, add it to the returned dictionary
        if !isempty(undefined)
            undefined_concepts[exercise["slug"]] = undefined
        end
    end

    return undefined_concepts
end

"""
Check if a concept extraction doc only contains valid concepts.

**Assumptions**:
- The first list is ignored, as it is assumed to contain the example implementations
- In the rest of the doc, all lists are assumed to have the following format: `- <concept>: <explanation>`. Entries that aren't lists are ignored

Returns a dictionary of concepts that don't exist in concepts mapping to the reason for why it's used.
"""
function check_concept_extraction_doc(file, concepts)
    # ignore meta information and the first two blocks
    md_content = Markdown.parse_file(file, flavor = :github).content

    undefined = Dict{String,String}()

    skippedfirst = false

    for md_block in md_content
        # if it's not a list, ignore it
        typeof(md_block) == Markdown.List || continue

        # if it's the very first list, it's the list of example implementations
        skippedfirst || (skippedfirst = true; continue)

        for entry in md_block.items
            entry_str = strip(Markdown.plain(entry))
            concept, explanation = split(entry_str, ": ")

            # if the concept is not in concepts.csv, push it to undefined
            if concept ∉ concepts[!, :concept]
                undefined[concept] = explanation
            end
        end
    end

    undefined
end

function check_concept_extraction_docs(concepts, rootpath, track)
    base_path = joinpath(rootpath, "languages", track, "reference", "exercise-concepts")

    undefined = Dict{String,Dict{String,String}}()

    # If the directory doesn't exist, we can't check anything, so everything is correct
    isdir(base_path) || return undefined

    for ex in readdir(base_path)
        # ignore README.md and _sidebar.md
        ex ∉ ("README.md", "_sidebar.md") || continue

        d = check_concept_extraction_doc(joinpath(base_path, ex), concepts)
        if !isempty(d)
            undefined[ex] = d
        end
    end

    undefined
end

"""
Check if all exercise directory names are valid concepts.

Returns an array of exercises named after undefined concepts.
"""
function check_exercise_directories(concepts, rootpath, track)
    # vector to store all undefined exercise directory names
    undefined = String[]

    base_path = joinpath(rootpath, "languages", track, "exercises", "concept")

    for ex in readdir(base_path)
        # ignore README.md, _sidebar.md and other files
        isdir(joinpath(base_path, ex)) || continue

        # if the concept is not in concepts.csv, push it to undefined
        if ex ∉ concepts[!, :concept]
            push!(undefined, ex)
        end
    end

    undefined
end

function main()
    anyerror = false

    @info "Parsing commandline arguments..."
    args = parse_commandline()
    conceptspath = joinpath(args["root"], "languages", args["track"], "reference", "concepts.csv")
    configpath = joinpath(args["root"], "languages", args["track"], "config.json")

    @info "Reading concepts from $conceptspath..."
    concepts = CSV.read(conceptspath)

    @info "Reading config from $configpath..."
    config = JSON.parsefile(configpath)

    if !args["no-config-check"]
        @info "Checking $configpath..."
        undefined_concepts = check_config_json(concepts, config)
        if !isempty(undefined_concepts)
            for (exercise, concepts) in undefined_concepts
                @error "Found undefined concepts in $exercise: $(join(concepts, ", "))"
            end
            anyerror = true
        end
    else
        @info "Skipping $configpath check..."
    end

    if !args["no-directory-check"]
        @info "Checking exercise directory names..."
        wrong_dirs = check_exercise_directories(concepts, args["root"], args["track"])
        if !isempty(wrong_dirs)
            @error "Found exercises named after undefined concepts: $(join(wrong_dirs, ", "))"
            anyerror = true
        end
    else
        @info "Skipping directory names check..."
    end

    if !args["no-extraction-check"]
        @info "Checking extracted concepts from v2..."
        undefined_concepts = check_concept_extraction_docs(concepts, args["root"], args["track"])
        if !isempty(undefined_concepts)
            for (ex, errors) in undefined_concepts
                pretty_errors = join(collect("- $(rpad("$k:", 20)) $v" for (k, v) in errors), '\n')
                @error "Found undefined concepts in $ex:\n$pretty_errors"
            end
            anyerror = true
        end
    else
        @info "Skipping extracted concepts check..."
    end
    anyerror
end

if abspath(PROGRAM_FILE) == @__FILE__
    main() ? exit(1) : exit(0)
end


text:
push!(LOAD_PATH,"./src")
push!(LOAD_PATH,"./app/resources/insurancecontracts")
using Documenter
using Pkg
Pkg.add(url="https://github.com/michaelfliegner/BitemporalPostgres.jl")
import BitempWebApp
using InsuranceContracts, InsurancePartners, InsuranceTariffs

makedocs(
    sitename = "BitempWebApp",
    format = Documenter.HTML(),
    modules = [BitempWebApp,InsuranceContracts, InsurancePartners,InsuranceTariffs],
    pages = [
        "Home" => "index.md",
        "BitempWebApp API" => [
            "BitempWebApp" => "api/BitempWebApp.md",
            "InsuranceContractsController" => "api/InsuranceContractsController.md"
        ]
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/michaelfliegner/BitempWebApp.jl"
)

max_stars_repo_name: michaelfliegner/BitempWebApp

max_issues_repo_name: michaelfliegner/BitempWebApp


text:
using DelimitedFiles

puzzleinput = readdlm("inputs/day10.txt", ' ', Int, '\n')
puzzleinput = sort(vec(puzzleinput))
# Sort joltages

# First Half

diffs = vcat(puzzleinput, last(puzzleinput) + 3) - vcat(0, puzzleinput)
# The charging outlet has an effective rating of 0 jolts and
# your device's built-in adapter is always 3 higher than the highest adapter

diff1jolt = count(d -> d == 1, diffs)
# How many differences of 1 jolt

diff3jolt = count(d -> d == 3, diffs)
# How many differences of 3 jolts

diff1jolt * diff3jolt

# Second Half

puzzleinput = append!([0], puzzleinput)
puzzleinput = append!(puzzleinput, [last(puzzleinput) + 3])
# Add lowest (0) and highest (max jolt in puzzle input + 3) jolts

len = length(puzzleinput)
paths = BitArray(
    abs(puzzleinput[i] - puzzleinput[j]) < 4 &&
    abs(puzzleinput[i] - puzzleinput[j]) != 0 &&
    j > i for i = 1:len, j = 1:len
)
# Create a diagonal matrix for the cartesian product of joltages
# but where (i, j) value is:
#   1 if joltage(i) - joltage(j) < 3, there is a possible electrical connection
#   0 otherwise, no possible connection
# Lower part of the diagonal matrix (j > i) is set to 0

sums = sum(paths, dims=2)[:, 1]
# For every joltage calculates how many possible connections,
# of course the min is 1 and the max is 3 (condition of max 3 jolt o diff)

pop!(sums)
# Removes last value, always 0

sums = prod(string.(sums))
# Convert the vector of numbers in a string for pattern search

# Because the nodes was ordered the string repesents a sequence
# of the number of paths starting from every node.

sums = replace(sums, "332" => "7")
# The sequence of nodes with a pattern "332" is a graph with 7 possible paths

sums = replace(sums, "32" => "4")
# The pattern "32" is a graph with 4 possible paths

# Other patterns (eg, 121 or 111) are not a problem because the preceding
# and the following has only 1 possible path.
# Is not possible to have the "31" pattern, why? think it about it :)

sums = collect.(sums)
sums = map(x -> Int(x) - 48, sums)
# Convert strings in array of integers

prod(sums)
# From the combinatorial calculus, the number of possible paths is
# the product of the possible paths between the subseguent joltages