text:
"""
    global TODOs:
        . ensure no function failes

"""

"""
    Computes several groebner bases of the given G
    while substituting coeffs of the underlying ideal with random points
    Allows to grasp the general structure of the groebner basis

    Returns the number of polynomials
    and the number of field generators (i.e coeficients)
    present in the Groebner basis
"""
function discover_groebner_structure(G::GroebnerEvaluator)
    # generate two substitution points for coeffs of G
    p1, p2 = generate_point(G), generate_point(G)

    # compute two groebner bases at that point
    gb1, gb2 = evaluate(G, p1), evaluate(G, p2)
    # obtain general info
    structure1 = groebner_structure(gb1)
    structure2 = groebner_structure(gb2)

    # if groebner structures does not match
    if structure1 != structure2
        # :D
        @debug "discovering failed.. recursive call"
        return discover_groebner_structure(G)
    end

    npolys, ncoeffs = structure1
    return npolys, ncoeffs
end

"""
    Computes the Grobner basis of the given ideal at several random points to
    estimate the variable-wise degrees of the coefficiet polynomials by applying
    univariate interpolation over each variable

    Returns an array of degrees expected to occure in Groebner basis,
    i-th array entry corresponds to coeff degree with respect to i-th variable
"""
function discover_groebner_degrees(G::GroebnerEvaluator)
    # staring from small amount of interpolation points
    initial_n = 8
    n = initial_n

    nxs = nvars(G)
    ground = base_ring(G)
    npolys, ncoeffs = discover_groebner_structure(G)
    ncoeffs = sum(ncoeffs)

    lightring, = Nemo.PolynomialRing(ground, "x")

    # idx of var --> max degree
    answer = zeros(Int, nxs)

    for varidx in 1:nxs
        n = initial_n
        predicted_degrees = [ n for _ in 1:ncoeffs ]
        all_success = false

        # @debug "handling variable " varidx

        while !all_success

            # criterion of interpolation success:
            # interpolated degree is far less
            # than the amount of interpolation points
            check_success = (deg) -> 2*(deg + 2) < n

            p = generate_point(G)
            xs = [ deepcopy(p) for _ in 1:n ]

            # TODO: sep into a function
            # Gleb: how about broadcast: xs[k] = sequence_gen.^(1:n)  ?
            # Alex: this won't work..
            for k in 1:n
                xs[k][varidx] = rand(ground)
            end

            evaluated = map(G, xs)
            generators = map(obtain_generators, evaluated)

            ys = map(c -> ground.(Int.(c)), generators)

            all_success = true
            for (j, deg) in enumerate(predicted_degrees)
                xi = [ x[varidx] for x in xs ]
                yi = [ y[j] for y in ys ]

                # TODO: fix
                interpolated = lightring(0)

                try
                    interpolated = interpolate_rational_function(
                                        lightring, xi, yi
                    )
                catch
                    @warn "rational interpolation failed =(\nincrease sample"
                    all_success = false
                    break
                end

                num, den = unpackfrac(interpolated)
                predicted_degrees[j] = maximum([degree(num), degree(den)])

                all_success = all_success && check_success(predicted_degrees[j])
            end

            # @debug "for $n points and $varidx variable degrees are " predicted_degrees

            n = n * 2
        end

        answer[varidx] = maximum(predicted_degrees)

    end

    answer
end

###############################################################################

"""
    Returns a plain array of coefficients of the given groebner basis
"""
### AA.coefficients would be better
function obtain_generators(groebner)
    generators = []
    for poly in groebner
        push!(generators, coefficients(poly)...)
    end
    generators
end

"""
    Returns an ideal of polynomials in ys of form
        < Fyi * Qx - Fxi * Qy >
    for each Fi in genset
    Here Q stands for the lcm of the denominators occuring in genset.

    So that generators of the resulting ideal are polynomials,
    but not fractions

    Saturates the ideal with Q by adding 1 - Qt to its generators

"""
function generators_to_saturated_ideal(genset)

    basepolyring = parent(numerator( first(genset) ))
    nvariables = length(gens(basepolyring))
    ground = base_ring(basepolyring)

    dens = map(denominator, genset)
    Q = dens[1]
    for d in dens
        Q = lcm(Q, d)
    end

    Fs = map(numerator ∘ (g -> g * Q), genset)

    ystrings = ["y$i" for i in 1:nvariables]
    yoverx, yoverxvars = Nemo.PolynomialRing(
                                                 basepolyring,
                                                 ystrings,
                                                 ordering=:lex
                                                 )

    Fx = Fs
    Fy = map(F -> change_base_ring(basepolyring, F, parent=yoverx), Fs)
    Qx = Q
    Qy = change_base_ring(basepolyring, Q, parent=yoverx)

    # Gleb: to cancel Fyi and Qy by their gcdi
    # oops
    I = [
        Fyi * Qx - Fxi * Qy
        for (Fyi, Fxi) in zip(Fy, Fx)
    ]

    I, t = saturate(I, Q)

    (I=I, yoverx=yoverx, basepolyring=basepolyring,
     nvariables=nvariables, ground=ground, ystrings=ystrings,
     Q=Q, t=t)
end

# TODO: to delete
# Gleb: As far as I see, this function is always followed by saturation (and this makes perfect sense)
# I suggest to do the saturation right here
function generators_to_ideal(genset)
    basepolyring = parent(numerator( first(genset) ))
    nvariables = length(gens(basepolyring))
    ground = base_ring(basepolyring)

    dens = map(denominator, genset)
    Q = dens[1]
    for d in dens
        Q = lcm(Q, d)
    end

    Fs = map(numerator ∘ (g -> g * Q), genset)

    ystrings = ["y$i" for i in 1:nvariables]
    yoverx, yoverxvars = Nemo.PolynomialRing(
                                                 basepolyring,
                                                 ystrings,
                                                 ordering=:lex
                                                 )

    Fx = Fs
    Fy = map(F -> change_base_ring(basepolyring, F, parent=yoverx), Fs)
    Qx = Q
    Qy = change_base_ring(basepolyring, Q, parent=yoverx)

    # Gleb: to cancel Fyi and Qy by their gcd
    I = [
        Fyi * Qx - Fxi * Qy
        for (Fyi, Fxi) in zip(Fy, Fx)
    ]

    # Gleb: how about a named tuple here? Impossible to remember the order
    I, yoverx, basepolyring, nvariables, ground, ystrings, Q
end

"""
    Saturates the given ideal by adding 1 - Q*t to its generators
    where t is the new introduced variable

    Returns the saturated ideal and the t variable
"""
function saturate(I, Q)
    R = parent(I[1])
    base = base_ring(R)
    strings = map(String, symbols(R))
    parentring, vs = Nemo.PolynomialRing(base, [strings..., "t"], ordering=:lex)

    t = last(vs)
    It = [
        change_parent_ring(f, parentring)[1]
        for f in I
    ]
    sat = 1 - Q*t
    push!(It, sat)
    It, t
end

###############################################################################

"""
    Returns the set of new generators for the field genereted by the given genset
    with no interpolation involved in computing
"""
function naive_groebner_ideal(genset)

    I, yoverx, basepolyring, nvariables, ground, ystrings, Q = generators_to_ideal(genset)

    It, t = saturate(I, Q)

    basepolyrings, = Singular.AsEquivalentSingularPolynomialRing(basepolyring)
    yoverxs,  = Singular.PolynomialRing(basepolyrings, [ystrings..., "t"])

    Is = [
          map_coefficients(
              c -> change_base_ring(tosingular(ground), c, parent=basepolyrings),
              f
          )
          for f in It
    ]

    Is = map(F -> change_base_ring(base_ring(yoverxs), F, parent=yoverxs), Is)

    Ideal = Singular.Ideal(yoverxs, Is...)

    gb = collect(gens(Singular.std(Ideal, complete_reduction=true)))

    gb = [
        double_singular2aa(f, base=basepolyring, new_ring=parent(t))
        for f in gb
    ]

    gb_no_t = []
    for f in gb
        if degree(f, t) == 0
            push!(gb_no_t, f)
        end
    end

    gb = [
        change_parent_ring(f, yoverx)
        for f in gb_no_t
    ]

    return gb
end

function idealize_and_eval_element(elem, eval_ring, point)
    A, B = numerator(elem), denominator(elem)

    A, _ = change_parent_ring(A, eval_ring)
    B, _ = change_parent_ring(B, eval_ring)

    Ae = evaluate(A, point)
    Be = evaluate(B, point)

    f = Ae * B - A * Be

    f
end

function idealize_element(elem, basepolyring, basepolyrings, yoverx, yoverxs)
    A, B = numerator(elem), denominator(elem)

    ground = base_ring(basepolyring)
    varstrings = string.(symbols(yoverx))

    tmp, = Nemo.PolynomialRing(basepolyring, varstrings)

    Ay = change_base_ring(basepolyring, A, parent=tmp)
    By = change_base_ring(basepolyring, B, parent=tmp)

    Ay = map_coefficients(c -> c // basepolyring(1), Ay)
    By = map_coefficients(c -> c // basepolyring(1), By)

    f = A*By - Ay*B

    f = map_coefficients(
             c -> change_base_ring(ground, numerator(c), parent=basepolyrings) //
                    change_base_ring(ground, denominator(c), parent=basepolyrings),
             f
    )

    f = change_base_ring(
              base_ring(yoverxs),
              f,
              parent=yoverxs)

    f
end

function aa_ideal_to_singular(I)
    basepolyring = parent(first(I))
    nvs = nvars(basepolyring)

    xstrings = ["x$i" for i in 1:nvs]
    ystrings = ["y$i" for i in 1:nvs]
    # Gleb: do you assume that the input variables are always x1, x2, ... ?

    F = base_ring(basepolyring)

    yoverx, yoverxvars = Nemo.PolynomialRing(basepolyring, ystrings, ordering=:lex)

    basepolyrings, = Singular.AsEquivalentSingularPolynomialRing(basepolyring)
    fracbases = FractionField(basepolyrings)
    yoverxs, = Singular.PolynomialRing(fracbases,
                             ystrings,
                             ordering=ordering_lp(nvs)*ordering_c())

    Fs = base_ring(basepolyrings)

    Is = [
          map_coefficients(c -> change_base_ring(Fs, numerator(c), parent=basepolyrings) // change_base_ring(Fs, denominator(c), parent=basepolyrings), f) for f in I
    ]

    Is = map(c -> change_base_ring(base_ring(yoverxs), c, parent=yoverxs), Is)

    Is, yoverx, basepolyring, yoverxs, basepolyrings, fracbases
end

function groebner_structure(gb)
    l = length(gb)
    a = map(f -> length(coefficients(f)), gb)
    l, a
end

function groebner_ideal_to_singular(polys)

    yoverx = parent(polys[1])
    varstrings = string.(Singular.symbols(yoverx))
    frac_ring = base_ring(polys[1])
    basepolyring = base_ring(frac_ring)
    ground = base_ring(basepolyring)

    s_basepolyring = tosingular(basepolyring)
    s_frac = FractionField(s_basepolyring)
    s_yoverx,  = Singular.PolynomialRing(s_frac, varstrings)

    polys = [
        map_coefficients(c ->
           change_base_ring(ground, numerator(c), parent=s_basepolyring) //
                change_base_ring(ground, denominator(c), parent=s_basepolyring),
           f
        )
        for f in polys
    ]

    polys = [
        change_base_ring(base_ring(s_yoverx), f, parent=s_yoverx)
        for f in polys
    ]

    polys, yoverx, basepolyring, s_yoverx, s_basepolyring
end

"""
    Computes the Groebner basis of the "good" ideal of the given genset
    interpolating coefficients as rational functions at several numerical points
"""
function groebner_ideal_backend(genset)
    # Generating "good" ideal and saturating it

    #=
    ideal_info = generators_to_saturated_ideal(genset)
    It, yoverx, basepolyring, nvariables, ground, ystrings, Q, t = ideal_info.I, ideal_info.yoverx, ideal_info.basepolyring, ideal_info.nvariables, ideal_info.ground, ideal_info.ystrings, ideal_info.Q, ideal_info.t
    =#

    I = idealize(genset, saturated=true)
    ground = I.ground

    # Estimating the largest degree of a coeff in the Groebner basis

    G = GroebnerEvaluator(I)

    true_structure = discover_groebner_structure(G)
    npolys, ncoeffs = true_structure

    @info "Groebner structure: $npolys polys, $(sum(ncoeffs)) coeffs"

    exponents = discover_groebner_degrees(G)
    maxexp = maximum(exponents) + 1

    @info "Maximal exponents per variables: $exponents"
    @info "The largest + 1: $maxexp"

    # Building a Kronecker substitution and generating substitution points

    xs, points = generate_good_kronecker_points(ground, exponents, I.nvariables)
    npoints = length(points)

    @info "The total number of sampling points: $npoints"

    @info "Evaluating Groebner bases..."

    # Evaluating Groebner bases of specializations
    gbs = []
    for (i, point) in enumerate(points)
        push!(gbs, evaluate(G, point))
        if i % 1000 == 0
            @info "Computed $i / $npoints Groebner bases.."
            @info last(gbs)
        end
    end

    # Ensure bases are of the same shape
    for (gb_i, gb) in enumerate(gbs)
        if groebner_structure(gb) != true_structure
            @warn "Beda: GB of $(groebner_structure(gb)) against $true_structure"
            @warn "Aborting.."
            @assert false
        end
    end

    # Performing univariate interpolation
    @info "Started interpolation.."
    uni, = Nemo.PolynomialRing(ground, "x")

    # Extract coefficients
    ys = [ [ zeros(ground, 0) for _ in 1:ncoeffs[k] ] for k in 1:npolys  ]
    gbcoeffs = [
        [
            map(c -> ground(Int(c)), coefficients(gb[i]))
            for i in 1:npolys
        ]
        for gb in gbs
    ]
    for poly_j in 1:npolys
        for coef_i in 1:ncoeffs[poly_j]
            for gb_k in 1:length(gbs)
                push!(
                    ys[poly_j][coef_i],
                    gbcoeffs[gb_k][poly_j][coef_i]
                )
            end
        end
    end

    answer_e = [ Dict() for _ in 1:npolys ]
    for (poly_j, poly) in enumerate(gbs[1])

        @info "Interpolating $poly_j th out of $npolys polynomials.."
        for (c_i, ev) in enumerate(exponent_vectors(poly))

            yscoeff = ys[poly_j][c_i]

            f = interpolate_rational_function(
                    uni,
                    xs,
                    yscoeff
            )

            answer_e[poly_j][ev] = backward_good_kronecker(
                                    f, I.basepolyring, exponents
            )
        end
    end

    # Reconstruct the final Groebner basis and new generators

    @info "Reconstructing polynomials.."

    gb = []
    R, = Nemo.PolynomialRing(
                                FractionField(I.basepolyring),
                                I.ystrings,
                                ordering=:lex
                           )

    for etoc in answer_e
        polybuilder = MPolyBuildCtx(R)
        for e in keys(etoc)
            push_term!(polybuilder, etoc[e], e[1:end-1])
        end
        push!(gb, finish(polybuilder))
    end

    gb
end

"""
    Checks if the good ideal composed of the given generators
    contains in the given groebner_basis ideal
"""
function check_ideal_inclusion(groebner_basis, generators)
    Is, yoverx, basepolyring, yoverxs, basepolyrings = groebner_ideal_to_singular(groebner_basis)
    Is = Ideal(yoverxs, Is)

    for f in generators
        f = idealize_element(f, basepolyring, basepolyrings, yoverx, yoverxs)

        fs = Ideal(yoverxs, f)

        if !Singular.contains(Is, fs)
            return false
        end
    end

    return true
end

"""
    Wrapper for the groebner basis computation function

    Performs Modular reduction and Reconstruction,
    implements CRT to work with big numbers,
    and conducts correctnes checks

    Currently supported input
    .    Fractions of Nemo polynomials with Nemo QQ coefficients
"""
function groebner_ideal(
                   initial_genset::Array{T};
                   modular=true,
                   check=false) where {T<:Frac{Nemo.fmpq_mpoly}}

    #=
        TODO:
            . ensure evaluated generators are of the same structure
            .
    =#

    # reduction moduli
    moduli = [ 2^30+3 ]

    groebner_bases = []

    while true

        modulo = last(moduli)
        FF = Nemo.GF(modulo)

        genset = initial_genset

        if modular
            @info "Modular reduction.."
            genset = modular_reduction(initial_genset, FF)
        end

        gb = groebner_ideal_backend(genset)
        push!(groebner_bases, gb)

        if modular
            @info "CRT reconstruction.."

            # let us just hope no coeff vanishes under a small field
            gb = CRT_reconstruction(groebner_bases, BigInt.(moduli))

            @info "Rational reconstruction.."
            gb = rational_reconstruction(gb, prod(BigInt.(moduli)))

        end

        if check
            @info "Running correctness checks.."
            isideal = check_ideal_inclusion(gb, initial_genset)
            isfield = check_field_inclusion(initial_genset, obtain_generators(gb))
            @info "Done.\n Is correct ideal: $isideal. Is correct field: $isfield"

            if isideal && isfield ### ???
                return gb
            end

            @warn "Computed *incorrect* \n$gb \nmodulo $modulo"

            push!(moduli, Primes.nextprime(BigInt(modulo + 1)))

        else
            return gb
        end

    end

    @error "how did we end up like this?.."
end

function groebner_ideal(
                   initial_genset::Array{T};
                   modular=true,
                   check=false) where {T<:Frac{MPoly{Singular.n_Q}}}

    # ?

    strings = string.(gens(parent(numerator(first(initial_genset)))))

    subpolyring, = Nemo.PolynomialRing(Nemo.QQ, strings)

    initial_genset = map(
            f -> change_base_ring(Nemo.QQ, numerator(f), parent=subpolyring) //
                change_base_ring(Nemo.QQ, denominator(f), parent=subpolyring),
            initial_genset
    )

    groebner_ideal(initial_genset, modular=modular, check=check)
end

## TODO: dup
function contains_randomized(
            generators, elem;
            p=0.99)

    # Theorem 3.3 from
    # Structural identifiability via input-output projections, Pogudin et alia

    # choose appropriate to the given 'p' random upper bound
    n = length(gens(parent(numerator(elem))))
    d = maximum([
            maximum([degrees(numerator(f))..., degrees(denominator(f))...])
            for f in generators
    ])
    M = round(Int, 6d^(n+3) / (1 - p))

    @info "Random values bound M = $M"

    # idealize the generators and the elem in the same rings
    I = idealize(generators)
    i = idealize(elem, I)
    # evaluate at a random point
    G = GroebnerEvaluator(I)
    g = GroebnerEvaluator(i)

    p = generate_point(G, M=M)

    gb = evaluate(G, p)

    elem = evaluate(g, p)

    @debug "Evaluated eLement is $elem"
    @debug "Evaluated groebner basis is $gb"

    Is = Singular.Ideal(I.evalring, gb)

    # check ideal membership
    return iszero(Singular.reduce(elem..., GroebnerBasis.f4(Is)))
end

function check_field_inclusion(generators, elements)
    for elem in elements
        a, b = unpackfrac(elem)
        if Nemo.isconstant(a) && Nemo.isconstant(b)
            continue
        end
        if !contains_randomized(generators, elem, p=0.999)
            return false
        end
    end
    true
end


text:
function encode_string_literal(buf::IOBuffer, data::Array{UInt8}; huffman=true)
    mask::UInt8 = 0x0
    if huffman
        mask = 0x80
        data = huffman_encode_bytes(data)
    end

    h = mask | UInt8(length(data))
    write(buf, h)
    write(buf, data)
end

function encode_literal(buf::IOBuffer, header::HeaderBinary; options...)
    mask::UInt8 = 0x0
    write(buf, mask)

    encode_string_literal(buf, header[1]; options...)
    encode_string_literal(buf, header[2]; options...)
end

function encode_indexed(buf::IOBuffer, index::UInt8, header::HeaderBinary; options...)
    mask::UInt8 = 0x40
    write(buf, mask | index)

    encode_string_literal(buf, header[2]; options...)
end

function encode_totally_indexed(buf::IOBuffer, index::UInt8; options...)
    mask::UInt8 = 0x80
    write(buf, mask | index)
end

find_totally_indexed(name::Vector{UInt8}, value::Vector{UInt8}) = findfirst(x->x==(name,value), STATIC_TABLE)

function encode(table::DynamicTable, headers; options...)
    buf = IOBuffer()
    for header in headers
        index = find_totally_indexed(bytearr(header[1]), bytearr(header[2]))
        if index === nothing
            encode_literal(buf, (bytearr(header[1]), bytearr(header[2])); options...)
        else
            encode_totally_indexed(buf, UInt8(index); options...)
        end
    end

    return take!(buf)
end

ext: jl


text:
#=
A perfect number is a number for which the sum of its proper divisors is exactly equal to the number. For example, the sum of the proper divisors of 28 would be 1 + 2 + 4 + 7 + 14 = 28, which means that 28 is a perfect number.

A number n is called deficient if the sum of its proper divisors is less than n and it is called abundant if this sum exceeds n.

As 12 is the smallest abundant number, 1 + 2 + 3 + 4 + 6 = 16, the smallest number that can be written as the sum of two abundant numbers is 24. By mathematical analysis, it can be shown that all integers greater than 28123 can be written as the sum of two abundant numbers. However, this upper limit cannot be reduced any further by analysis even though it is known that the greatest number that cannot be expressed as the sum of two abundant numbers is less than this limit.

Find the sum of all the positive integers which cannot be written as the sum of two abundant numbers.
=#

function sum_proper_divisors(n)
  prod = 1
  k = 2
  while k*k <= n
    p = 1
    while n%k==0
      p=p*k+1
      n = div(n,k)
    end
    prod*=p
    k += 1
  end
  if n>1
    prod*=1+n
  end
  prod
end

function abundant_numbers_less_than(n)
  arr = Int64[]
  for i in 1:n-1
    if sum_proper_divisors(i) > i+i
      push!(arr, i)
    end
  end
  arr
end

function calc(n)
  abundant_nums = abundant_numbers_less_than(n)
  canbewrittenasabundant = Bool[false for x in 1:n]
  for i in 1:length(abundant_nums)
    for j in 1:length(abundant_nums)
      if abundant_nums[i] + abundant_nums[j] <= n
        canbewrittenasabundant[abundant_nums[i]+abundant_nums[j]] = true
      else
        break
      end
    end
  end

  total = 0
  for i in 1:n
    total += canbewrittenasabundant[i] ? 0 : i
  end
  total
end
@time println(calc(28_124))

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
data = randn(100)
hist(data; horiz=true)


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

JLLWrappers.@generate_main_file_header("ERFA")
JLLWrappers.@generate_main_file("ERFA", UUID("a00367d1-88f1-52c5-9c6b-193aa5f849e1"))
end  # module ERFA_jll


text:
using DrWatson
@quickactivate "INN_Velocity-Migration"

using JLD, Statistics, PyPlot, ImageFiltering, InvertibleNetworks, LinearAlgebra, Distributed, SharedArrays
addprocs(20)
@everywhere include("../src/FWI_64x200.jl")
@everywhere using .FWI_64x200, Statistics, ImageFiltering, Dates

X_orig = load(datadir("vel_4k_samples_64x200_lin_vel.jld"), "m_all")
X_orig = Float32.(X_orig)
n1, n2, nc, nsamples = size(X_orig)
ntrain = Int(nsamples*.9)
X_train_orig=X_orig[:, :, :, 1:ntrain]

AN = ActNorm(ntrain)
X_train_orig = AN.forward(X_train_orig) # zero mean and unit std of the training data X

# define data and figure path
datapath = plotsdir("chint/test3_6000_16_16")
figfolder = "fwi/test3_6000_16_16/posterior_samples2"
mkpath(plotsdir(figfolder))

X_fixed = load(joinpath(datapath, "posterior_samples2.jld"), "X_fixed")
X_post  = load(joinpath(datapath, "posterior_samples2.jld"), "X_post")
X_fixed = AN.inverse(X_fixed) # unnormarlize based on the same parameter of the training data X
X_post = AN.inverse(X_post)   # unnormarlize based on the same parameter of the training data X
X_post_mean = mean(X_post; dims=4)

m = X_fixed[:, :, 1, 1] # true velocity in squared slowness
idx_wb = 9 # index of the water bottom
niterations = 15

####################################################################################################
## Scenario 1: initial velocity m0 = all posterior samples
println(string("Scenario 1: initial velocity m0 = posterior samples"))

test_size = size(X_post)[4] #100
m0_all = SharedArray{Float32,3}(n1, n2, test_size)
fwi_all = SharedArray{Float32,3}(n1, n2, test_size)
obj_all = SharedArray{Float32,2}(niterations, test_size)

t1 = now()
println(string("FWI parallel computation starts at ", t1))

@sync @distributed for i = 1:test_size

    println(string("FWI using the ", i, "th sample of X posterior starts at ", now() ))

    m0 = X_post[:, :, 1, i]
    m0[1:idx_wb, :] .= m[1,1]
    m0[idx_wb+1:end, :] = imfilter(m0[idx_wb+1:end, :], Kernel.gaussian(3)) # smoothed velocity
    fwi_all[:,:,i], obj_all[:,i] = fwi(m, m0, niterations)

    m0_all[:,:,i] = m0

end

t2 = now()
println(string("FWI parallel computation finishes after ", Dates.value.(t2-t1)/3600000, " hours"))

m0_all = Array(m0_all)
fwi_all = Array(fwi_all)
obj_all = Array(obj_all)
save(plotsdir(figfolder, "fwi_post_samples.jld"), "m0_all", m0_all, "fwi_all", fwi_all, "obj_all", obj_all)

i=1
v = sqrt.(1f0 ./ m)
v0 = sqrt.(1f0 ./ m0_all[:,:,i])
fwi_v = sqrt.(1f0 ./ fwi_all[:,:,i])
obj = obj_all[:,i]

figure(figsize=[16,8])
ax1 = subplot(2, 2, 1); imshow(v, cmap="jet", aspect = "auto"); plt.colorbar(); title("True velocity (km/s)")
ax2 = subplot(2, 2, 2); imshow(v0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("Initial model, L2=", round(norm(v-v0,2), sigdigits=4)))
ax3 = subplot(2, 2, 3); imshow(fwi_v, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi_v,2), sigdigits=4)))
ax4 = subplot(2, 2, 4); plot(obj); title("Objective function value")
savefig(plotsdir(figfolder, "fwi_post_samples.png"))

####################################################################################################
## Scenario 2: initial velocity m0 = Mean of the FWI results in Scenario 1
println(string("Scenario 2: initial velocity m0 = fwi_mean"))

fwi_mean =  dropdims( mean(fwi_all, dims=3), dims=3 )
fwi_std =  dropdims( std(fwi_all, dims=3), dims=3 )
figure(figsize=[16,8])
ax1 = subplot(1, 2, 1); imshow(fwi_mean, cmap="jet", aspect = "auto"); plt.colorbar(); title("Mean of FWI")
ax2 = subplot(1, 2, 2); imshow(fwi_std, cmap="jet", aspect = "auto"); plt.colorbar(); title("Standard deviation of FWI")
savefig(plotsdir(figfolder, "fwi_mean&std.png"))

fwi_mean[1:idx_wb, :] .= m[1,1]
fwi1, obj1 = fwi(m, fwi_mean, niterations)
save(plotsdir(figfolder, "fwi_fwi_mean.jld"), "m", m, "m0", fwi_mean, "fwi", fwi1, "obj", obj1)

v = sqrt.(1f0 ./ m)
v0 = sqrt.(1f0 ./ fwi_mean)
fwi_v = sqrt.(1f0 ./ fwi3)
obj = obj3

figure(figsize=[16,8])
ax1 = subplot(2, 2, 1); imshow(v, cmap="jet", aspect = "auto"); plt.colorbar(); title("True velocity (km/s)")
ax2 = subplot(2, 2, 2); imshow(v0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("Initial model, L2=", round(norm(v-v0,2), sigdigits=4)))
ax3 = subplot(2, 2, 3); imshow(fwi_v, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi_v,2), sigdigits=4)))
ax4 = subplot(2, 2, 4); plot(obj); title("Objective function value")
savefig(plotsdir(figfolder, "fwi_fwi_mean.png"))

####################################################################################################
## Scenario 3: initial velocity m0 = posterior mean
println(string("Scenario 3: initial velocity m0 = posterior mean"))

m0 = X_post_mean[:, :, 1, 1]
m0[1:idx_wb, :] .= m[1,1]
m0[idx_wb+1:end, :] = imfilter(m0[idx_wb+1:end, :], Kernel.gaussian(3)) # smoothed velocity

fwi2, obj2 = fwi(m, m0, niterations)
save(plotsdir(figfolder, "fwi_post_mean.jld"), "m", m, "m0", m0, "fwi", fwi2, "obj", obj2)

v = sqrt.(1f0 ./ m)
v0 = sqrt.(1f0 ./ m0)
fwi_v = sqrt.(1f0 ./ fwi2)
obj = obj2

figure(figsize=[16,8])
ax1 = subplot(2, 2, 1); imshow(v, cmap="jet", aspect = "auto"); plt.colorbar(); title("True velocity (km/s)")
ax2 = subplot(2, 2, 2); imshow(v0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("Initial model, L2=", round(norm(v-v0,2), sigdigits=4)))
ax3 = subplot(2, 2, 3); imshow(fwi_v, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi_v,2), sigdigits=4)))
ax4 = subplot(2, 2, 4); plot(obj); title("Objective function value")
savefig(plotsdir(figfolder, "fwi_post_mean.png"))

####################################################################################################
## Scenario 4: initial velocity m0 = migration velocity used in rtm
println(string("Scenario 4: initial velocity m0 = migration velocity used in rtm"))

# build a migration velocity m0 whose value increases linearly with depth
v0 = ones(n1, 1)
v0[idx_wb+1:end] = collect(range(3.5, 5.5, length = n1-idx_wb))
v0 = (1f0 ./ v0).^2 # convert to slowness squared
v0[1:idx_wb] .= m[1,1] # water layer velocity
m_lin = repeat(v0, 1, n2)

fwi3, obj3 = fwi(m, m_lin, niterations)
save(plotsdir(figfolder, "fwi_lin_vel.jld"), "m", m, "m0", m_lin, "fwi", fwi3, "obj", obj3)

v = sqrt.(1f0 ./ m)
v0 = sqrt.(1f0 ./ m_lin)
fwi_v = sqrt.(1f0 ./ fwi3)
obj = obj3

figure(figsize=[16,8])
ax1 = subplot(2, 2, 1); imshow(v, cmap="jet", aspect = "auto"); plt.colorbar(); title("True velocity (km/s)")
ax2 = subplot(2, 2, 2); imshow(v0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("Initial model, L2=", round(norm(v-v0,2), sigdigits=4)))
ax3 = subplot(2, 2, 3); imshow(fwi_v, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi_v,2), sigdigits=4)))
ax4 = subplot(2, 2, 4); plot(obj); title("Objective function value")
savefig(plotsdir(figfolder, "fwi_lin_vel.png"))

####################################################################################################
## Scenario 5: initial velocity m0 = smoothed true velocity
println(string("Scenario 5: initial velocity m0 = smoothed true velocity"))

m_smooth = m[1,1] .* ones(Float32, size(m)) # smoothed velocity
m_smooth[idx_wb+1:end, :] = imfilter(m[idx_wb+1:end, :], Kernel.gaussian(10)) # smoothed velocity

fwi4, obj4 = fwi(m, m_smooth, niterations)
save(plotsdir(figfolder, "fwi_smooth_vel.jld"), "m", m, "m0", m_smooth, "fwi", fwi4, "obj", obj4)

v = sqrt.(1f0 ./ m)
v0 = sqrt.(1f0 ./ m_smooth)
fwi_v = sqrt.(1f0 ./ fwi4)
obj = obj4

figure(figsize=[16,8])
ax1 = subplot(2, 2, 1); imshow(v, cmap="jet", aspect = "auto"); plt.colorbar(); title("True velocity (km/s)")
ax2 = subplot(2, 2, 2); imshow(v0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("Initial model, L2=", round(norm(v-v0,2), sigdigits=4)))
ax3 = subplot(2, 2, 3); imshow(fwi_v, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi_v,2), sigdigits=4)))
ax4 = subplot(2, 2, 4); plot(obj); title("Objective function value")
savefig(plotsdir(figfolder, "fwi_smooth_vel.png"))


text:
#==========================================================================================#

# ONE SAMPLE: INTERFACE

function hausman_test(obj₁::ParModel, obj₂::ParModel)
    hausman_test(obj₁, obj₂, intersect(coefnames(obj₁), coefnames(obj₂)))
end

function hausman_test(obj₁::ParModel, obj₂::ParModel, name::String)
    hausman_test(obj₁, obj₂, [name])
end

function hausman_test(obj₁::ParModel, obj₂::ParModel, names::Vector{String})

    i₁    = findall((in)(names), coefnames(obj₁))
    i₂    = findall((in)(names), coefnames(obj₂))
    w₁    = getweights(obj₁)
    w₂    = getweights(obj₂)
    corr₁ = getcorr(obj₁)
    corr₂ = getcorr(obj₂)

    (iszero(i₁) | iszero(i₂)) && throw("missing coefficients in at least one model")
    isequal(corr₁, corr₂) || throw("different correlation structures")
    isequal(w₁, w₂) || throw("different weighting schemes")

    output       = ParEstimate()
    output.names = copy(names)

    _hausman_test!(output, obj₁, i₁, obj₂, i₂, w₁, corr₁)

    return output
end

# ONE SAMPLE: ESTIMATION

function _hausman_test!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        obj₂::ParModel,
        i₂::Vector{Int},
        w::UnitWeights,
        corr::Heteroscedastic
    )

    Ψ₁  = influence(obj₁, w)
    ψ₁  = view(Ψ₁, :, i₁)
    Ψ₂  = influence(obj₂, w)
    ψ₂  = view(Ψ₂, :, i₂)
    V₁₂ = ψ₁' * ψ₂

    corr.corrected && lmul!(varcorrection(corr, w), V₁₂)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_test!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        obj₂::ParModel,
        i₂::Vector{Int},
        w::FrequencyWeights,
        corr::Heteroscedastic
    )

    Ψ₁  = influence(obj₁, w)
    ψ₁  = view(Ψ₁, :, i₁)
    Ψ₂  = influence(obj₂, w)
    ψ₂  = view(Ψ₂, :, i₂) ; lmul!(Diagonal(w), ψ₂)
    V₁₂ = ψ₁' * ψ₂

    corr.corrected && lmul!(varcorrection(corr, w), V₁₂)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_test!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        obj₂::ParModel,
        i₂::Vector{Int},
        w::AbstractWeights,
        corr::Heteroscedastic
    )

    Ψ₁  = influence(obj₁, w)
    ψ₁  = view(Ψ₁, :, i₁) ; lmul!(Diagonal(w), ψ₁)
    Ψ₂  = influence(obj₂, w)
    ψ₂  = view(Ψ₂, :, i₂) ; lmul!(Diagonal(w), ψ₂)
    V₁₂ = ψ₁' * ψ₂

    corr.corrected && lmul!(varcorrection(corr, w), V₁₂)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_test!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        obj₂::ParModel,
        i₂::Vector{Int},
        w::UnitWeights,
        corr::Clustered
    )

    Ψ₁  = influence(obj₁, w)
    ψ₁  = corr.mat * view(Ψ₁, :, i₁)
    Ψ₂  = influence(obj₂, w)
    ψ₂  = corr.mat * view(Ψ₂, :, i₂)
    V₁₂ = ψ₁' * ψ₂

    corr.corrected && lmul!(varcorrection(corr, w), V₁₂)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_test!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        obj₂::ParModel,
        i₂::Vector{Int},
        w::AbstractWeights,
        corr::Clustered
    )

    Ψ₁  = influence(obj₁, w)
    ψ₁  = corr.mat * view(Ψ₁, :, i₁) ; lmul!(Diagonal(w), ψ₁)
    Ψ₂  = influence(obj₂, w)
    ψ₂  = corr.mat * view(Ψ₂, :, i₂) ; lmul!(Diagonal(w), ψ₂)
    V₁₂ = ψ₁' * ψ₂

    corr.corrected && lmul!(varcorrection(corr, w), V₁₂)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_test!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        obj₂::ParModel,
        i₂::Vector{Int},
        w::UnitWeights,
        corr::CrossCorrelated
    )

    Ψ₁  = influence(obj₁, w)
    ψ₁  = view(Ψ₁, :, i₁)
    Ψ₂  = influence(obj₂, w)
    ψ₂  = view(Ψ₂, :, i₂)
    V₁₂ = ψ₁' * corr.mat * ψ₂

    corr.corrected && lmul!(varcorrection(corr, w), V₁₂)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_test!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        obj₂::ParModel,
        i₂::Vector{Int},
        w::AbstractWeights,
        corr::CrossCorrelated
    )

    Ψ₁  = influence(obj₁, w)
    ψ₁  = view(Ψ₁, :, i₁) ; lmul!(Diagonal(w), ψ₁)
    Ψ₂  = influence(obj₂, w)
    ψ₂  = view(Ψ₂, :, i₂) ; lmul!(Diagonal(w), ψ₂)
    V₁₂ = ψ₁' * corr.mat * ψ₂

    corr.corrected && lmul!(varcorrection(corr, w), V₁₂)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

#==========================================================================================#

# TWO INDEPENDENT SAMPLES

function chow_test(obj₁::ParModel, obj₂::ParModel)
    chow_test(obj₁, obj₂, intersect(coefnames(obj₁), coefnames(obj₂)))
end

function chow_test(obj₁::ParModel, obj₂::ParModel, name::String)
    chow_test(obj₁, obj₂, [name])
end

function chow_test(obj₁::ParModel, obj₂::ParModel, names::Vector{String})

    i₁ = findall((in)(names), coefnames(obj₁))
    i₂ = findall((in)(names), coefnames(obj₂))

    (iszero(i₁) | iszero(i₂)) && throw("missing coefficients in at least one model")

    output       = ParEstimate()
    output.names = copy(names)
    output.β     = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V     = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂)

    return output
end

#==========================================================================================#

# TWO DEPENDENT SAMPLES: INTERFACE

function chow_test(obj₁::ParModel, obj₂::ParModel, corr::CorrStructure)
    chow_test(obj₁, obj₂, corr, intersect(coefnames(obj₁), coefnames(obj₂)))
end

function chow_test(obj₁::ParModel, obj₂::ParModel, corr::CorrStructure, name::String)
    chow_test(obj₁, obj₂, corr, [name])
end

function chow_test(
        obj₁::ParModel,
        obj₂::ParModel,
        corr::CorrStructure,
        names::Vector{String}
    )

    i₁ = findall((in)(names), coefnames(obj₁))
    i₂ = findall((in)(names), coefnames(obj₂))
    w₁ = getweights(obj₁)
    w₂ = getweights(obj₂)

    (iszero(i₁) | iszero(i₂)) && throw("missing coefficients in at least one model")

    if !(typeof(corr) == typeof(getcorr(obj₁)) == typeof(getcorr(obj₂)))
        throw("incompatible correlation structures")
    end

    output       = ParEstimate()
    output.names = copy(names)

    _chow_test!(output, obj₁, i₁, w₁, obj₂, i₂, w₂, corr)

    return output
end

# TWO DEPENDENT SAMPLES: ESTIMATION

function _chow_test!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        w₁::UnitWeights,
        obj₂::ParModel,
        i₂::Vector{Int},
        w₂::UnitWeights,
        corr::Heteroscedastic
    )

    nonmissing₁ = getnonmissing(obj₁)
    nonmissing₂ = getnonmissing(obj₂)
    touse       = nonmissing₁ .* nonmissing₂
    touse₁      = touse[nonmissing₁]
    touse₂      = touse[nonmissing₂]

    Ψ₁  = influence(obj₁, w₁)
    ψ₁  = view(Ψ₁, touse₁, i₁)
    Ψ₂  = influence(obj₂, w₂)
    ψ₂  = view(Ψ₂, touse₂, i₂)
    V₁₂ = ψ₁' * ψ₂

    corr.corrected && lmul!(varcorrection(corr, UnitWeights(sum(touse))), V₁₂)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _chow_test!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        w₁::W,
        obj₂::ParModel,
        i₂::Vector{Int},
        w₂::W,
        corr::Heteroscedastic
    ) where W <: AbstractWeights

    nonmissing₁ = getnonmissing(obj₁)
    nonmissing₂ = getnonmissing(obj₂)
    touse       = nonmissing₁ .* nonmissing₂
    touse₁      = touse[nonmissing₁]
    touse₂      = touse[nonmissing₂]

    Ψ₁  = influence(obj₁, w₁)
    ψ₁  = view(Ψ₁, touse₁, i₁) ; lmul!(Diagonal(w₁), ψ₁)
    Ψ₂  = influence(obj₂, w₂)
    ψ₂  = view(Ψ₂, touse₂, i₂) ; lmul!(Diagonal(w₂), ψ₂)
    V₁₂ = ψ₁' * ψ₂

    corr.corrected && lmul!(varcorrection(corr, UnitWeights(sum(touse))), V₁₂)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _chow_test!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        w₁::UnitWeights,
        obj₂::ParModel,
        i₂::Vector{Int},
        w₂::UnitWeights,
        corr::Clustered
    )

    nonmissing₁ = getnonmissing(obj₁)
    nonmissing₂ = getnonmissing(obj₂)

    Ψ₁  = influence(obj₁, w₁)
    ψ₁  = view(corr.mat, :, nonmissing₁) * view(Ψ₁, :, i₁)
    Ψ₂  = influence(obj₂, w₂)
    ψ₂  = view(corr.mat, :, nonmissing₂) * view(Ψ₂, :, i₂)
    V₁₂ = ψ₁' * ψ₂

    corr.corrected && lmul!(varcorrection(corr, UnitWeights(sum(touse))), V₁₂)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _chow_test!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        w₁::W,
        obj₂::ParModel,
        i₂::Vector{Int},
        w₂::W,
        corr::Clustered
    ) where W <: AbstractWeights

    nonmissing₁ = getnonmissing(obj₁)
    nonmissing₂ = getnonmissing(obj₂)

    Ψ₁  = influence(obj₁, w₁)
    ψ₁  = view(corr.mat, :, nonmissing₁) * view(Ψ₁, :, i₁) ; lmul!(Diagonal(w₁), ψ₁)
    Ψ₂  = influence(obj₂, w₂)
    ψ₂  = view(corr.mat, :, nonmissing₂) * view(Ψ₂, :, i₂) ; lmul!(Diagonal(w₂), ψ₂)
    V₁₂ = ψ₁' * ψ₂

    corr.corrected && lmul!(varcorrection(corr, UnitWeights(sum(touse))), V₁₂)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _chow_test!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        w₁::UnitWeights,
        obj₂::ParModel,
        i₂::Vector{Int},
        w₂::UnitWeights,
        corr::CorrStructure
    )

    nonmissing₁ = getnonmissing(obj₁)
    nonmissing₂ = getnonmissing(obj₂)

    Ψ₁  = influence(obj₁, w₁)
    ψ₁  = view(Ψ₁, :, i₁)
    Ψ₂  = influence(obj₂, w₂)
    ψ₂  = view(Ψ₂, :, i₂)
    V₁₂ = ψ₁' * view(corr.mat, nonmissing₁, nonmissing₂) * ψ₂

    corr.corrected && lmul!(varcorrection(corr, UnitWeights(sum(touse))), V₁₂)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _chow_test!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        w₁::W,
        obj₂::ParModel,
        i₂::Vector{Int},
        w₂::W,
        corr::CorrStructure
    ) where W <: AbstractWeights

    nonmissing₁ = getnonmissing(obj₁)
    nonmissing₂ = getnonmissing(obj₂)

    Ψ₁  = influence(obj₁, w₁)
    ψ₁  = view(Ψ₁, :, i₁) ; lmul!(Diagonal(w₁), ψ₁)
    Ψ₂  = influence(obj₂, w₂)
    ψ₂  = view(Ψ₂, :, i₂) ; lmul!(Diagonal(w₂), ψ₂)
    V₁₂ = ψ₁' * view(corr.mat, nonmissing₁, nonmissing₂) * ψ₂

    corr.corrected && lmul!(varcorrection(corr, UnitWeights(sum(touse))), V₁₂)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

max_stars_repo_name: lbittarello/Microeconometrics

max_issues_repo_name: lbittarello/Microeconometrics

max_forks_repo_name: lbittarello/Microeconometrics


text:
@testset "4.6.4.2 (a+b csc)^m (d csc)^n (A+B csc+C csc^2)" begin
    (A, B, C, a, b, x, ) = @variables A B C a b x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(a+b*csc(e+f*x))^m*(d*csc(e+f*x))^n*(A+B*csc(e+f*x)+C*csc(e+f*x)^2)=#

    #= ::Section:: =#
    #=Integrands*of*the*form*(d*csc(e+f*x))^n*(a+b*csc(e+f*x))^m*(A+C*csc(e+f*x)^2)=#

    #= ::Section:: =#
    #=Integrands*of*the*form*(d*csc(e+f*x))^n*(a+b*csc(e+f*x))^m*(B*csc(e+f*x)+C*csc(e+f*x)^2)=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(d*csc(e+f*x))^n*(a+b*csc(e+f*x))^m*(A+B*csc(e+f*x)+C*csc(e+f*x)^2)=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(d*csc(e+f*x))^(n/2)*(a+b*csc(e+f*x))^m*(A+B*csc(e+f*x)+C*csc(e+f*x)^2)=#

    @test_int [((a + b*csc(x))*(A + B*csc(x) + C*csc(x)^2))/csc(x)^(1/2), x, 7, -2*(b*B + a*C)*cos(x)*sqrt(csc(x)) - (2/3)*b*C*cos(x)*csc(x)^(3/2) + 2*(b*B - a*(A - C))*sqrt(csc(x))*Elliptic.E(pi/4 - x/2, 2)*sqrt(sin(x)) - (2/3)*(3*A*b + 3*a*B + b*C)*sqrt(csc(x))*Elliptic.F(pi/4 - x/2, 2)*sqrt(sin(x))]
    end

max_stars_repo_path: test/4 Trig functions/4.6 Cosecant/4.6.4.2 (a+b csc)^m (d csc)^n (A+B csc+C csc^2).jl


text:
struct TrackedArray{T, N, S <: AbstractArray{T, N}} <: AbstractArray{T, N}
    value::S
    record::Record{S}
end

struct TrackedRealArray{T <: Real, N, S <: AbstractArray{T, N}} <: AbstractArray{TrackedReal{T}, N}
    value::S
    record::Record{S}
end

struct TrackedComplexArray{T <: Complex, N, S <: AbstractArray{T, N}} <: AbstractArray{TrackedComplex{T}, N}
    value::S
    record::Record{S}
end

const TrackedArrayType = Union{TrackedArray, TrackedRealArray, TrackedComplexArray}

function Base.show(io::IO, mime::MIME"text/plain", x::TrackedArrayType)
    print(io, "tracked ")
    show(io, mime, x.value)
end

function track(A::AbstractArray, record::Record=leaf(A))
    return if eltype(A) <: Real
        TrackedRealArray(A, record)
    elseif eltype(A) <: Complex
        TrackedComplexArray(A, record)
    else
        TrackedArray(A, record)
    end
end

is_tracked(::TrackedArrayType) = true
Base.IndexStyle(X::TrackedArrayType) = IndexStyle(untrack(X))
Base.size(X::TrackedArrayType, idx::Int...) = size(untrack(X), idx...)
Base.length(X::TrackedArrayType) = length(untrack(X))


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   J4 orbit propagator algorithm.
#
#   This algorithm propagates the orbit considering the secular perturbations of
#   central body zonal harmonics as presented in [1, p. 647-654] and in [2].
#   Notice that only the terms J2, J2², and J4 are considered, i.e. the J6 is
#   assumed to be 0. The effect of the drag is also taken into account. This
#   can be used as a propagator of mean elements for mission analysis in which
#   the satellite orbit is maintained.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [2] Merson, R. H (1961). The motion of a satellite in an axi-symmetric
#       gravitational field. Geophysical journal international, Vol. 4(1),
#       p. 17-52.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
export j4_gc_egm08, j4_gc_egm96, j4_gc_jgm02, j4_gc_jgm03
export j4_gc_egm08_f32, j4_gc_egm96_f32, j4_gc_jgm02_f32, j4_gc_jgm03_f32
export j4_init, j4!

################################################################################
#                                  Constants
################################################################################

# These constants were obtained from the GFC files. Remember that:
#
#   J_n = -C_n,0 * sqrt(2n+1)
#

# EGM-08 Gravitational constants.
const j4_gc_egm08 = J4_GravCte(
    6378137.0,
    sqrt(3.986004415e14 / 6378137.0^3),
    0.0010826261738522227,
    -1.6198975999169731e-6
)

const j4_gc_egm08_f32 = J4_GravCte{Float32}(
    6378137.0,
    sqrt(3.986004415e14 / 6378137.0^3),
    0.0010826261738522227,
    -1.6198975999169731e-6
)

# EGM-96 Gravitational constants.
const j4_gc_egm96 = J4_GravCte(
    6378136.3,
    sqrt(3.986004415e14/6378136.3^3),
    0.0010826266835531513,
    -1.619621591367e-6
)

const j4_gc_egm96_f32 = J4_GravCte{Float32}(
    6378136.3,
    sqrt(3.986004415e14/6378136.3^3),
    0.0010826266835531513,
    -1.619621591367e-6
)

# JGM-02 Gravitational constants.
const j4_gc_jgm02 = J4_GravCte(
    6378136.3,
    sqrt(3.986004415e14/6378136.3^3),
    0.0010826269256388149,
    -1.62042999e-6
)

const j4_gc_jgm02_f32 = J4_GravCte{Float32}(
    6378136.3,
    sqrt(3.986004415e14/6378136.3^3),
    0.0010826269256388149,
    -1.62042999e-6
)

# JGM-03 Gravitational constants.
const j4_gc_jgm03 = J4_GravCte(
    6378136.3,
    sqrt(3.986004415e14/6378136.3^3),
    0.0010826360229829945,
    -1.619331205071e-6
)

const j4_gc_jgm03_f32 = J4_GravCte{Float32}(
    6378136.3,
    sqrt(3.986004415e14/6378136.3^3),
    0.0010826360229829945,
    -1.619331205071e-6
)

################################################################################
#                                  Functions
################################################################################

"""
    j4_init(epoch::Number, a_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, f_0::Number, dn_o2::Number = 0, ddn_o6::Number = 0; j4_gc::J4_GravCte{T} = j4_gc_egm08) where T

Initialize the data structure of J4 orbit propagator algorithm.

!!! note
    The type used in the propagation will be the same as used to define the
    gravitational constants in the structure `j4_gc`.

# Args

- `epoch::Number`: Epoch of the initial mean orbital elements [Julian Day].
- `a_0::Number`: Initial mean semi-major axis [m].
- `e_0::Number`: Initial mean eccentricity.
- `i_0::Number`: Initial mean inclination [rad].
- `Ω_0::Number`: Initial mean right ascension of the ascending node [rad].
- `ω_0::Number`: Initial mean argument of perigee [rad].
- `f_0::Number`: Initial mean true anomaly [rad].
- `dn_o2::Number`: First time derivative of the mean motion divided by two
    [rad/s^2].
- `ddn_o6::Number`: Second time derivative of the mean motion divided by six
    [rad/s^3].

# Keywords

* `j4_gc::J4_GravCtr`: J4 orbit propagator gravitational constants (see
    [`J4_GravCte`](@ref)). (**Default** = `j4_gc_egm08`)

# Returns

The structure [`J4_Structure`](@ref) with the initialized parameters.
"""
function j4_init(
    epoch::Tepoch,
    a_0::Number,
    e_0::Number,
    i_0::Number,
    Ω_0::Number,
    ω_0::Number,
    f_0::Number,
    dn_o2::Number = 0,
    ddn_o6::Number = 0;
    j4_gc::J4_GravCte{T} = j4_gc_egm08
) where {Tepoch, T}
    # Unpack the gravitational constants to improve code readability.
    @unpack_J4_GravCte j4_gc

    # Convert all inputs to the correct type.
    a_0    = T(a_0)
    e_0    = T(e_0)
    i_0    = T(i_0)
    Ω_0    = T(Ω_0)
    ω_0    = T(ω_0)
    f_0    = T(f_0)
    dn_o2  = T(dn_o2)
    ddn_o6 = T(ddn_o6)

    # Initial values.
    al_0 = a_0 / R0               # ............ Normalized semi-major axis [er]
    n_0  = μm / al_0^(T(3 / 2))   # ............ Unperturbed mean motion [rad/s]
    p_0  = al_0 * (1 - T(e_0)^2)  # ..................... Semi-latus rectum [er]
    M_0  = f_to_M(e_0, f_0)       # ................. Initial mean anomaly [rad]

    # Auxiliary variables.
    dn               = 2 * dn_o2  # .... Time-deriv. of the mean motion [rad/s²]
    e_0²             = e_0^2
    e_0⁴             = e_0^4
    sin_i_0, cos_i_0 = sincos(i_0)
    sin_i_0          = sin(i_0)
    sin_i_0²         = sin_i_0^2
    sin_i_0⁴         = sin_i_0^4
    aux              = (1 - e_0²)
    saux             = sqrt(aux)
    p_0²             = p_0^2
    p_0⁴             = p_0^4
    J2²              = J2^2

    # First-order time-derivative of the orbital elements.
    #
    # See [1, p 692].

    δa   = -T( 2 /  3) * al_0 * dn / n_0
    δe   = -T( 2 /  3) * (1 - e_0) * dn / n_0
    δΩ   = -T( 3 /  2) * n_0 * J2 /p_0² * cos_i_0 +
            T( 3 / 32) * n_0 * J2²/ p_0⁴ * cos_i_0 * (12 -  4e_0² - (80 +  5e_0²) * sin_i_0²) +
            T(15 / 32) * n_0 * J4 / p_0⁴ * cos_i_0 * ( 8 + 12e_0² - (14 + 21e_0²) * sin_i_0²)

    δω   = T( 3 /   4) * n_0 * J2  / p_0² * (4 - 5sin_i_0²) +
           T( 9 / 384) * n_0 * J2² / p_0⁴ * (     56e_0² + (760 -  36e_0²) * sin_i_0² - (890 +  45e_0²) * sin_i_0⁴) -
           T(15 / 128) * n_0 * J4  / p_0⁴ * (64 + 72e_0² - (248 + 252e_0²) * sin_i_0² + (196 + 189e_0²) * sin_i_0⁴)

    δM_0 = T( 3 /   4) * n_0 * J2 / p_0² * saux * (2 - 3sin_i_0²) +
           T( 3 / 512) * n_0 * J2²/ p_0⁴ / saux * (
               (         320e_0² - 280e_0⁴) +
               ( 1600 - 1568e_0² + 328e_0⁴) * sin_i_0² +
               (-2096 + 1072e_0² +  79e_0⁴) * sin_i_0⁴
           ) - T(45 / 128) * n_0 * J4 / p_0⁴ * saux * e_0² * (-8 + 40sin_i_0 - 35sin_i_0²)

    # Create the output structure with the data.
    J4_Structure{Tepoch, T}(
        epoch  = epoch,
        al_0   = al_0,
        n_0    = n_0,
        e_0    = e_0,
        i_0    = i_0,
        Ω_0    = Ω_0,
        ω_0    = ω_0,
        f_0    = f_0,
        M_0    = M_0,
        dn_o2  = dn_o2,
        ddn_o6 = ddn_o6,
        j4_gc  = j4_gc,
        Δt     = 0,
        al_k   = al_0,
        e_k    = e_0,
        i_k    = i_0,
        Ω_k    = Ω_0,
        ω_k    = ω_0,
        f_k    = f_0,
        M_k    = M_0,
        δa     = δa,
        δe     = δe,
        δΩ     = δΩ,
        δω     = δω,
        δM_0   = δM_0,
    )
end

"""
    j4!(j4d::J4_Structure{T}, t::Number) where T

Propagate the orbit defined in `j4d` (see [`J4_Structure`](@ref)) until the time
`t` [s].

!!! note
    The internal values in `j4d` will be modified.

# Returns

- The position vector represented in the inertial frame at time `t` [m].
- The velocity vector represented in the inertial frame at time `t` [m/s]

# Remarks

The inertial frame in which the output is represented depends on which frame it
was used to generate the orbit parameters. Notice that the perturbation theory
requires an inertial frame with true equator.
"""
function j4!(j4d::J4_Structure{Tepoch, T}, t::Number) where {Tepoch, T}
    # Unpack the variables.
    @unpack epoch, al_0, n_0, e_0, i_0, Ω_0, ω_0, f_0, M_0 = j4d
    @unpack dn_o2, ddn_o6, j4_gc, δa, δe, δΩ, δω, δM_0 = j4d
    @unpack R0, μm, J2, J4 = j4_gc

    # Time elapsed since epoch.
    Δt = T(t)

    # Propagate the orbital elements.
    al_k = al_0 + δa * Δt
    e_k  = e_0  + δe * Δt
    i_k  = i_0
    Ω_k  = mod(Ω_0 + δΩ * Δt, T(2π))
    ω_k  = mod(ω_0 + δω * Δt, T(2π))
    M_k  = mod(M_0 + (δM_0 + n_0) * Δt, T(2π))
    f_k  = M_to_f(e_k, M_k)

    # Make sure that eccentricity is not lower than 0.
    e_k < 0 && (e_k = T(0))

    # Compute the position and velocity vectors given the orbital elements.
    r_i_k, v_i_k = kepler_to_rv(al_k * R0, e_k, i_k, Ω_k, ω_k, f_k)

    # Update the J4 orbit propagator structure.
    @pack! j4d = Δt, al_k, e_k, i_k, Ω_k, ω_k, M_k, f_k

    # Return the position and velocity vector represented in the inertial
    # reference frame.
    return r_i_k, v_i_k
end


text:
"""
MOI constraints
"""

"""
Single variable bound constraints
"""
MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{MOI.EqualTo{Float64}}) = true

"""
Binary/Integer variable support
"""
MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{<:VAR_TYPES}) = true

"""
Linear constraints
"""
MOI.supports_constraint(::Optimizer, ::Type{SAF}, ::Type{<:BOUNDS}) = true

"""
Quadratic constraints (scalar i.e., vectorized constraints are not supported)
"""
MOI.supports_constraint(::Optimizer, ::Type{SQF}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{SQF}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{SQF}, ::Type{MOI.EqualTo{Float64}}) = true

"""
``define_add_constraint`` macro - from Ipopt.jl
"""
macro define_add_constraint(function_type, set_type, array_name)
    quote
        function MOI.add_constraint(model::Optimizer, func::$function_type, set::$set_type)
            check_inbounds(model, func)
            push!(model.$(array_name), (func, set))
            return MOI.ConstraintIndex{$function_type, $set_type}(length(model.$(array_name)))
        end
    end
end

"""
``MOI.add_constraint()`` overloads for all the supported constraint types
"""
@define_add_constraint(SAF, MOI.LessThan{Float64}, linear_le_constraints)
@define_add_constraint(SAF, MOI.GreaterThan{Float64}, linear_ge_constraints)
@define_add_constraint(SAF, MOI.EqualTo{Float64}, linear_eq_constraints)
@define_add_constraint(SQF, MOI.LessThan{Float64}, quadratic_le_constraints)
@define_add_constraint(SQF, MOI.GreaterThan{Float64}, quadratic_ge_constraints)
@define_add_constraint(SQF, MOI.EqualTo{Float64}, quadratic_eq_constraints)

"""
Binary variable support
"""
function MOI.add_constraint(model::Optimizer, v::SVF, ::MOI.ZeroOne)
    println("==============flag! MOI_wrapper>constraints>add_constraint===========================")
    vi = v.variable
    println(model)
    check_inbounds(model, vi)
    # the bounds are set using info_array_of_variables
    # according to mlubin the bounds should not be set here as the bounds should stay when
    # the binary constraint is deleted
    model.variable_info[vi.value].is_binary = true

    return MOI.ConstraintIndex{SVF, MOI.ZeroOne}(vi.value)
end

"""
Integer variable support
"""
function MOI.add_constraint(model::Optimizer, v::SVF, ::MOI.Integer)
    vi = v.variable
	check_inbounds(model, vi)
    model.variable_info[vi.value].is_integer = true

    return MOI.ConstraintIndex{SVF, MOI.Integer}(vi.value)
end

"""
ConstraintIndex support
"""
MOI.supports(::Optimizer, ::MOI.ConstraintName, ::Type{CI}) = true


text:
function tmp(posterior, t, numsamples)

    μ, Σ = mean(posterior), cov(posterior)

    a_idx, b_idx = [1;2], [3;4]

    μ_a = μ[a_idx]
    μ_b = μ[b_idx]

    Σ_aa = Σ[a_idx, a_idx]

    Σ_ab = Σ[a_idx, b_idx]

    Σ_ba = Σ[b_idx, a_idx]

    Σ_bb = Σ[b_idx, b_idx]

    samples = [rand(posterior) for i in 1:numsamples]

    aux(x) = x[a_idx]*t + x[b_idx]

    # mean(map(aux, samples)), μ_a*t + μ_b

    Σ_a_given_b = Σ_aa - Σ_ab*(Σ_bb\Σ_ba)

    K = t*Σ_ab/Σ_bb + I

    cov(map(aux, samples)), t*Σ_a_given_b*t + K*Σ_bb*K'

end


text:
using ConstrainedDynamics

# Parameters
ex = [1.;0.;0.]
ey = [0.;1.;0.]
ez = [0.;0.;1.]

l1 = 1.0
x, y = .1, .1
b1 = Box(x, l1, y, l1, color = RGBA(1., 1., 0.))

vert11 = [0.;l1 / 2;0.]
vert12 = -vert11

# Initial orientation
phi1, phi2, phi3 = 0., 2pi / 3, 4pi / 3
q1, q2, q3 = UnitQuaternion(RotZ(phi1)), UnitQuaternion(RotZ(phi2)), UnitQuaternion(RotZ(phi3))

# Links
origin = Origin{Float64}()
link1 = Body(b1)
link2 = deepcopy(link1, b1)

# Constraints
joint0to123 = EqualityConstraint(PlanarFree(origin, link1, ez; p2=vert12), PlanarFree(origin, link2, ez; p2=vert12))
joint1to23 = EqualityConstraint(Spherical(link1, link2; p1=vert11, p2=vert11))

links = [link1; link2]
constraints = [joint0to123;joint1to23]
shapes = [b1]

mech = Mechanism(origin, links, constraints, shapes = shapes)
setPosition!(origin,link1,p2 = vert11,Δq = q1)
setPosition!(origin,link2,p2 = vert11,Δq = q2)


text:
# -------------------- RANDOM EFFECTS ------------------------------

function rkey(key::AbstractString, ilevel::AbstractString="") ilevel != "" ? key*" ("*ilevel*")" : key  end
function rkey(row::DataFrameRow) return rkey(row[:class], row[:level]) end

abstract RanEffect <: ModelEffect

function getRandomFormula(rlist::Vector{Symbol}, excludeCol::Symbol=:NA, df_name::AbstractString="df_data")
    fmula=""
    for s in setdiff(rlist, [excludeCol])
        fmula=fmula*"+"*df_name*"[:"*string(s)*"]"
    end
    return fmula[2:end]
end

function genRandomTotals(reff::RanEffect, df_data::DataFrame)
    model=reff.v_model
    idf=reff.src
    rcols=Symbol[]
    for ranfx in unique(idf[:class])    # -- Loop through each random effect :: create random df_data columns - and add list to reff
        df_data[Symbol(ranfx)] = map(x->string(x), df_data[Symbol(ranfx)])    # Convert all ranfx to string ... since they are categorical anyway
        dest_colname=Symbol("reff_"*lowercase(ranfx)*"_"*model)
        push!(rcols,dest_colname)
        df_data[dest_colname]=0.0
        for row in eachrow(idf[idf[:class].==ranfx,:])   # Loop through each level for outter random effect
            level=string(row[:level])
            println("populating : ",dest_colname," : for :: (",level,")")
            df_data[ (df_data[Symbol(ranfx)].==level) & (df_data[:group].==0),dest_colname]=row[:B0]
            df_data[ (df_data[Symbol(ranfx)].==level) & (df_data[:group].==1),dest_colname]=row[:B1]
        end
    end
    reff.rcols=rcols
    reff.fmula=getRandomFormula(rcols)

    function calc_reff(row::DataFrameRow)
        tot=0
        for x in 1:length(row)
            tot=tot+row[x]
        end
        return tot
    end
    df_data[Symbol("reff_"*model)] = map(x->calc_reff(x),   eachrow(df_data[rcols]) )
    df_data[Symbol("reff_"*model)] = map(Float64,df_data[Symbol("reff_"*model)])      #convert to float

    # --- Random Effect Mean Scor cols ---
    for row = eachrow(reff.src)
        k = row[:key]
        #needs review
        dest_col0=Symbol(k*"_"*model*"_0")
        dest_col1=Symbol(k*"_"*model*"_1")
        rcols = setdiff(reff.rcols, [Symbol("reff_"*lowercase(row[:class])*"_"*model)])
        println("generating : ",k*"_"*model*"_X := sum(",rcols,", :pre_"*model*"_scoreX,+ BX) ")
        function tot_reff(row::DataFrameRow, B::Float64)
            tot=0
            for x in 1:length(row)
                tot=tot+row[x]
            end
            tot = exp(tot+B)
            if model== "pen"
                return tot = tot / (1+tot)  #exp(eval(parse(fmula0))) ./ (1+exp(eval(parse(fmula0))))
            end
            return tot
        end
        #model
        df_data[dest_col0] = map(x->tot_reff(x, Float64(row[:B0]) ),   eachrow(df_data[vcat(rcols,Symbol("pre_"*model*"_score0"))]) )
        df_data[dest_col1] = map(x->tot_reff(x, row[:B1] ),   eachrow(df_data[vcat(rcols,Symbol("pre_"*model*"_score1"))]) )
        df_data[dest_col0] = map(Float64,df_data[dest_col0])      #convert to float
        df_data[dest_col1] = map(Float64,df_data[dest_col1])      #convert to float
    end

end

type ROcc <: RanEffect
    sdf::DataFrame
    src::DataFrame
    rcols::Vector{Symbol}
    fmula::AbstractString
    v_model::AbstractString
    function ROcc(df_data::DataFrame, idf::DataFrame)  this=new(); this.v_model="occ"; this.sdf=deepcopy(SDF);
         this.src = idf
         for row in eachrow(idf)
             pushSDFrow!(this.sdf,OrderedDict(:key => row[:key], :model => this.v_model))
         end
        delete!(this.sdf, :B)
        delete!(this.sdf, :SE)
        delete!(this.sdf, :P)
        if length(this.sdf[1]) > 0 genRandomTotals(this,df_data)  end
        return this
    end
end

type RDolOcc <: RanEffect
    sdf::DataFrame
    src::DataFrame
    rcols::Vector{Symbol}
    fmula::AbstractString
    v_model::AbstractString
    function RDolOcc(df_data::DataFrame, idf::DataFrame)  this=new(); this.v_model="dolocc"; this.sdf=deepcopy(SDF);
         this.src = idf
         for row in eachrow(idf)
             pushSDFrow!(this.sdf,OrderedDict(:key => row[:key], :model => this.v_model))
         end
        delete!(this.sdf, :B)
        delete!(this.sdf, :SE)
        delete!(this.sdf, :P)
        if length(this.sdf[1]) > 0  genRandomTotals(this,df_data) end
        return this
    end
end

type RPen <: RanEffect
    sdf::DataFrame
    src::DataFrame
    rcols::Vector{Symbol}
    fmula::AbstractString
    v_model::AbstractString
    function RPen(df_data::DataFrame, idf::DataFrame)  this=new(); this.v_model="pen"; this.sdf=deepcopy(SDF);
         this.src = idf
         for row in eachrow(idf)
             pushSDFrow!(this.sdf,OrderedDict(:key => row[:key], :model => this.v_model))
         end
        delete!(this.sdf, :B)
        delete!(this.sdf, :SE)
        delete!(this.sdf, :P)
        if length(this.sdf[1]) > 0 genRandomTotals(this,df_data)  end
        return this
    end
end

function Base.show(io::IO, raneffect::RanEffect)
    show(io,raneffect.sdf)
end

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
using SparseArrays
using ArgCheck

"""
    Permute{A} <: Bijector{1}

A bijector implementation of a permutation. The permutation is performed
using a matrix of type `A`. There are a couple of different ways to construct `Permute`:

    Permute([0 1; 1 0])          # will map [1, 2] => [2, 1]
    Permute([2, 1])              # will map [1, 2] => [2, 1]
    Permute(2, 2 => 1, 1 => 2)   # will map [1, 2] => [2, 1]
    Permute(2, [1, 2] => [2, 1]) # will map [1, 2] => [2, 1]

If this is not clear, the examples might be of help.

# Examples
A simple example is permuting a vector of size 3.

```julia-repl
julia> b1 = Permute([
           0 1 0;
           1 0 0;
           0 0 1
       ])
Permute{Array{Int64,2}}([0 1 0; 1 0 0; 0 0 1])

julia> b2 = Permute([2, 1, 3])           # specify all elements at once
Permute{SparseArrays.SparseMatrixCSC{Float64,Int64}}(

  [2, 1]  =  1.0
  [1, 2]  =  1.0
  [3, 3]  =  1.0)

julia> b3 = Permute(3, 2 => 1, 1 => 2)    # element-wise
Permute{SparseArrays.SparseMatrixCSC{Float64,Int64}}(
  [2, 1]  =  1.0
  [1, 2]  =  1.0
  [3, 3]  =  1.0)

julia> b4 = Permute(3, [1, 2] => [2, 1])  # block-wise
Permute{SparseArrays.SparseMatrixCSC{Float64,Int64}}(
  [2, 1]  =  1.0
  [1, 2]  =  1.0
  [3, 3]  =  1.0)

julia> b1.A == b2.A == b3.A == b4.A
true

julia> b1([1., 2., 3.])
3-element Array{Float64,1}:
 2.0
 1.0
 3.0

julia> b2([1., 2., 3.])
3-element Array{Float64,1}:
 2.0
 1.0
 3.0

julia> b3([1., 2., 3.])
3-element Array{Float64,1}:
 2.0
 1.0
 3.0

julia> b4([1., 2., 3.])
3-element Array{Float64,1}:
 2.0
 1.0
 3.0

julia> inv(b1)
Permute{LinearAlgebra.Transpose{Int64,Array{Int64,2}}}([0 1 0; 1 0 0; 0 0 1])

julia> inv(b1)(b1([1., 2., 3.]))
3-element Array{Float64,1}:
 1.0
 2.0
 3.0
```
"""
struct Permute{A} <: Bijector{1}
    A::A
end

Base.:(==)(b1::Permute, b2::Permute) = b1.A == b2.A

function Permute(indices::AbstractVector{Int})
    # construct a sparse-matrix for use in the multiplication
    n = length(indices)
    A = spzeros(n, n)

    for (i, idx) in enumerate(indices)
        A[idx, i] = 1.0
    end

    return Permute(A)
end

function Permute(n::Int, indices::Pair{Int, Int}...)
    A = spdiagm(0 => ones(n))

    dests = Set{Int}()
    sources = Set{Int}()

    for (src, dst) in indices
        @argcheck dst ∉ dests
        @argcheck src ∉ sources

        push!(dests, dst)
        push!(sources, src)

        A[dst, src] = 1.0
        A[src, src] = 0.0  # <= remove `src => src`
    end

    @argcheck (sources ∩ dests) == (sources ∪ dests) "$sources ∩ $dests ≠ $sources ∪ $dests"

    dropzeros!(A)
    return Permute(A)
end

function Permute(n::Int, indices::Pair{Vector{Int}, Vector{Int}}...)
    A = spdiagm(0 => ones(n))

    dests = Set{Int}()
    sources = Set{Int}()

    for (srcs, dsts) in indices
        @argcheck length(srcs) == length(dsts) "$srcs => $dsts is not bijective"

        for (src, dst) in zip(srcs, dsts)
            @argcheck dst ∉ dests "$dst used more than once"
            @argcheck src ∉ sources "$src used more than once"

            push!(dests, dst)
            push!(sources, src)

            A[dst, src] = 1.0
            A[src, src] = 0.0  # <= remove `src => src`
        end
    end

    @argcheck (sources ∩ dests) == (sources ∪ dests) "$sources ∩ $dests ≠ $sources ∪ $dests"

    dropzeros!(A)
    return Permute(A)
end

@inline (b::Permute)(x::AbstractVecOrMat) = b.A * x
@inline inv(b::Permute) = Permute(transpose(b.A))

logabsdetjac(b::Permute, x::AbstractVector) = zero(eltype(x))
logabsdetjac(b::Permute, x::AbstractMatrix) = zero(eltype(x), size(x, 2))


text:
# // process synthes1
# // mixed-integer nonlinear program
#
# // Source: Test problem 1 (Synthesis of processing system) in
# // M. Duran & I.E. Grossmann,
# // "An outer approximation algorithm for a class of mixed integer nonlinear
# //  programs", Mathematical Programming 36, pp. 307-339, 1986.
#
# // Number of variables:   6 (3 binary variables)
# // Number of constraints: 6
# // Objective nonlinear
# // Nonlinear constraints

function minlp_demo(gm::Bool = false)
    m = JuMP.Model()
    @variable(m, x[1:6] >= 0)
    JuMP.set_binary.(x[4:6])
    JuMP.set_upper_bound(x[1], 2)
    JuMP.set_upper_bound(x[2], 2)
    JuMP.set_upper_bound(x[3], 1)
    JuMP.set_upper_bound(x[4], 1)
    JuMP.set_upper_bound(x[5], 1)
    JuMP.set_upper_bound(x[6], 1)

    @objective(m, Min, 10x[1] - 17x[3] -5x[4] + 6x[5] + 8x[6])

    @constraint(m, c3, x[2] - x[1] <= 0)
    @constraint(m, c4, x[2] - 2*x[4] <= 0)
    @constraint(m, c5, x[1] - x[2] - 2*x[5] <= 0)
    @constraint(m, c6, x[4] + x[5] <= 1)

    if !gm
        @NLconstraint(m, c1, 0.8*log(x[2] + 1) + 0.96*log(x[1] - x[2] + 1) - 0.8*x[3] >= 0)
        @NLconstraint(m, c2, log(x[2] + 1) + 1.2*log(x[1] - x[2] + 1) - x[3] - 2*x[6] >= -2)
        set_optimizer(m, BARON_SILENT)
        return m
    else
        gm = GlobalModel(model = m, name = "minlp")
        set_optimizer(gm, CPLEX_SILENT)
        add_nonlinear_constraint(gm, :(x -> 0.8*log(x[2] + 1) + 0.96*log(x[1] - x[2] + 1) - 0.8*x[3]),
                                vars =[x[1], x[2], x[3]], name = "c1")
        add_nonlinear_constraint(gm, :(x -> log(x[2] + 1) + 1.2*log(x[1] - x[2] + 1) - x[3] - 2*x[6] + 2),
                                vars =[x[1], x[2], x[3], x[6]], name = "c2")
        return gm
    end
end


text:
# Implementation of 'isprime' from the Julia standard library: (uses Miller Rabin)
witnesses(n::Union(Uint8,Int8,Uint16,Int16)) = (2, 3)
witnesses(n::Union(Uint32,Int32)) = n < 1373653 ? (2, 3) : (2, 7, 61)
witnesses(n::Union(Uint64,Int64)) =
        n < 1373653         ? (2, 3) :
        n < 4759123141      ? (2, 7, 61) :
        n < 2152302898747   ? (2, 3, 5, 7, 11) :
        n < 3474749660383   ? (2, 3, 5, 7, 11, 13) :
                              (2, 325, 9375, 28178, 450775, 9780504, 1795265022)

function isPrime(n::Integer)
    n == 2 && return true
    (n < 2) | iseven(n) && return false
    s = trailing_zeros(n - 1)
    d = (n - 1) >>> s
    for a in witnesses(n)
        a < n || break
        x = powermod(a,d,n)
        x == 1 && continue
        t = s
        while x != n-1
            (t -= 1) <= 0 && return false
            x = oftype(n, Base.widemul(x, x) % n)
            x == 1 && return false
        end
    end
    return true
end


text:
# You come across a field of hydrothermal vents on the ocean floor! These vents
# constantly produce large, opaque clouds, so it would be best to avoid them if
# possible.
#
# They tend to form in lines; the submarine helpfully produces a list of nearby
# lines of vents (your puzzle input) for you to review. For example:
#
# 0,9 -> 5,9
# 8,0 -> 0,8
# 9,4 -> 3,4
# 2,2 -> 2,1
# 7,0 -> 7,4
# 6,4 -> 2,0
# 0,9 -> 2,9
# 3,4 -> 1,4
# 0,0 -> 8,8
# 5,5 -> 8,2
#
# Each line of vents is given as a line segment in the format x1,y1 -> x2,y2
# where x1,y1 are the coordinates of one end the line segment and x2,y2 are the
# coordinates of the other end. These line segments include the points at both
# ends. In other words:
#
# An entry like 1,1 -> 1,3 covers points 1,1, 1,2, and 1,3.
# An entry like 9,7 -> 7,7 covers points 9,7, 8,7, and 7,7.
#
# For now, only consider horizontal and vertical lines: lines where either x1 =
# x2 or y1 = y2.
#
# So, the horizontal and vertical lines from the above list would produce the
# following diagram:
#
# .......1..
# ..1....1..
# ..1....1..
# .......1..
# .112111211
# ..........
# ..........
# ..........
# ..........
# 222111....
#
# In this diagram, the top left corner is 0,0 and the bottom right corner is
# 9,9. Each position is shown as the number of lines which cover that point or
# . if no line covers that point. The top-left pair of 1s, for example, comes
# from 2,2 -> 2,1; the very bottom row is formed by the overlapping lines 0,9
# -> 5,9 and 0,9 -> 2,9.
#
# To avoid the most dangerous areas, you need to determine the number of points
# where at least two lines overlap. In the above example, this is anywhere in
# the diagram with a 2 or larger - a total of 5 points.
#
# Consider only horizontal and vertical lines. At how many points do at least
# two lines overlap?

map = zeros(1000, 1000)
for l in eachline("input.txt")
  # parse the line and skip the input if it's not horizontal or vertical
  (x1, y1, x2, y2) = parse.(Int, split(l, r",| -> "))
  x1 != x2 && y1 != y2 && continue

  # make sure x1<x2, and y1<y2 - it simplifies the loop below
  x1 > x2 && ((x1,x2) = (x2,x1))
  y1 > y2 && ((y1,y2) = (y2,y1))

  # add one to each point
  for x in x1:x2, y in y1:y2
    map[x, y] += 1
  end
end

# count points where at least two lines overlap
println("Result: ", count(>(1), map))


text:
type ExampleIndividual <: Individual
  genome::Array{Int64, 1}
  fitness::Float64
  code

  function ExampleIndividual(size::Int64, max_value::Int64)
    genome = rand(1:max_value, size)
    return new(genome, -1.0, nothing)
  end

  ExampleIndividual(genome::Array{Int64, 1}) = new(genome, -1.0, nothing)
end

type ExamplePopulation <: Population
  individuals::Array{ExampleIndividual, 1}

  function ExamplePopulation(population_size::Int64, genome_size::Int64)
    individuals = Array(ExampleIndividual, 0)
    for i=1:population_size
      push!(individuals, ExampleIndividual(genome_size, 1000))
    end

    return new(individuals)
  end

  function ExamplePopulation(individuals::Array{ExampleIndividual, 1})
    return new(copy(individuals))
  end
end

max_stars_repo_name: abeschneider/GrammaticalEvolution

max_issues_repo_name: abeschneider/GrammaticalEvolution

max_forks_repo_name: abeschneider/GrammaticalEvolution


text:
using Plots

function plot_results(sp, pv, w, d; hd = 0, s=1, stage_1=[:gci, :gco], stage_2=[:gci2, :gco2], inv_dict = 0)
    if inv_dict != 0
        u_pv = inv_dict[:u_pv]
        u_wind = inv_dict[:u_wind]
        u_storage = inv_dict[:u_storage]
        u_heatpump = inv_dict[:u_heatpump]
        u_heat_storage = inv_dict[:u_heat_storage]
    else
        u_pv = value.(sp[1, :u_pv])
        u_wind = value.(sp[1, :u_wind])
        u_storage = value.(sp[1, :u_storage])
        u_heatpump = value.(sp[1, :u_heatpump])
        u_heat_storage = value.(sp[1, :u_heat_storage])
    end
    plt_sto = plot()
    plt_invest = plot()
    plt = plot()
    plot!(plt_invest, pv .* u_pv, label="pv")
    plot!(plt_invest, w .* u_wind, label="wind")
    plot!(plt_invest, d, label="demand")
    t_xi = scenarios(sp)[s].data.t_xi
    recovery_time = sp.stages[2].parameters[:recovery_time]

    stor_charge = value.(sp[1, :sto_soc])
    plot!(plt_sto, 1:length(pv), stor_charge, label="global storage charge")
    plot!(plt_sto, (t_xi):(t_xi+recovery_time-1), value.(sp[2, :sto_soc2],s), label=string("stochastic storage charge")*string(s), linestyle=:dash, linewidth=2)

    for var in stage_1
        plot!(plt, 1:length(pv), value.(sp[1, var]), label=string(var))
    end

    for var in stage_2
        plot!(plt, (t_xi+1):(t_xi+recovery_time), value.(sp[2, var], s), label=string(var)*string(s), linestyle=:dash, linewidth=2)
    end
    if hd != 0.
        plt_heat = plot()
        COP = sp.stages[2].parameters[:COP]

        plot!(plt_heat, hd, label = "heat demand")
        plot!(plt_heat, COP*value.(sp[1, :heatpumpflow]), label = "heatpump")
        plot!(plt_heat, value.(sp[1, :heat_sto_in])-value.(sp[1, :heat_sto_out]), label = "heat storage use")
        plot!(plt_heat, (t_xi+1):(t_xi+recovery_time), COP*value.(sp[2, :heatpumpflow2], s), label = "heatpump$s", linestyle=:dash, linewidth=2)
        plot!(plt_heat, (t_xi+1):(t_xi+recovery_time), value.(sp[2, :heat_sto_in2], s)-value.(sp[2, :heat_sto_out2], s), label = "heat storage use $s", linestyle=:dash, linewidth=2)

        display(plot(plt_invest, plt, plt_sto, plt_heat, layout = (4,1)))
    else
        display(plot(plt_invest, plt, plt_sto, layout = (3,1)))
    end

end

function plot_difference(sp; s = 1)
    plt_gc = plot()
    plt_sto = plot()
    t_xi = scenarios(sp)[s].data.t_xi
    recovery_time = sp.stages[2].parameters[:recovery_time]
    plot!(plt_gc, value.(sp[1, :heatpumpflow])[t_xi+1:t_xi+recovery_time], label = "global heatpumpflow")
    plot!(plt_gc, value.(sp[2, :heatpumpflow2], s), label = "stochastic heatpumpflow")
    plot!(plt_gc, -value.(sp[1, :gci])[t_xi+1:t_xi+recovery_time]+value.(sp[2, :gci2], s), label = "global grid connection use")
    plot!(plt_gc, value.(sp[2, :gco2], s) - value.(sp[1, :gco])[t_xi+1:t_xi+recovery_time], xlabel = "time after the event, h", label = "stochastic grid connection use")
    #plot!(plt_gc, -value.(sp[1, :gco])[t_xi+1:t_xi+recovery_time]+value.(sp[2, :gco2], s), xlabel = "time after the event, h", label = "gco2-gco")

    stor_charge = value.(sp[1, :sto_soc])[t_xi+1:t_xi+recovery_time]
    stor_charge2 = value.(sp[2, :sto_soc2], s)

    plot!(plt_sto, stor_charge2-stor_charge, label = "soc2-soc")

    display(plot(plt_gc, plt_sto, layout = (2, 1)))
end


text:
using BeerTest
using Test, SQLite

@testset "BeerTest.jl" begin
    @testset "database.jl" begin
        @test isfile(BeerTest.path_beers)
        @test isfile(BeerTest.path_breweries)
        @test isfile(BeerTest.path_geocodes)

        @test isfile(BeerTest.path_database)
        @test issetequal(first(SQLite.tables(SQLite.DB(BeerTest.path_database))), ["beers", "locations"])
    end

    @testset "dbloader.jl" begin
        coord, locinfo = BeerTest.get_mapinfo()
        @test size(coord.longitudes) == size(coord.latitudes)
        @test size(coord.longitudes, 1) == coord.dim
        @test size(coord.pairwise, 1) == coord.dim
        @test size(coord.pairwise, 2) == coord.dim
        @test length(locinfo.beers) == coord.dim
        @test size(locinfo.locations, 1) == coord.dim
        @test size(locinfo.locations[!, :name], 1) == coord.dim
    end

    @testset "features.jl" begin
        coord, locinfo = BeerTest.get_mapinfo()
        feat = BeerTest.Features(coord, locinfo)
        @test size(coord.pairwise) == size(feat.pairwise_distance)
        @test coord.dim == size(feat.home_distance, 1)
        @test coord.dim == size(feat.beer_count, 1)
        @test coord.dim == size(feat.scores, 1)

        @test sum(feat.pairwise_distance[:, 1].^2) ≈ 1.
        @test sum(feat.beer_count.^2) ≈ 1.
    end

    @testset "main.jl" begin
        state = initialize()
        output = false
        beer_counts = collect(1:3) .|> _ -> main(;state, output)
        @test all(first(beer_counts) .== beer_counts)
    end
end


text:
using NestedEnvironments
using DifferentialEquations
using Transducers
using Test

# `BaseEnv` is a kind of syntax sugar for definition of environment; provided by `src/zoo.jl`.
# Note: `NestedEnvironments.initial_condition(env::BaseEnv) = env.initial_state`
struct Env2 <: AbstractEnv
    env21::BaseEnv
    env22::BaseEnv
end
struct Env <: AbstractEnv
    env1::BaseEnv
    env2::Env2
    gain::Float64
end

# differential equations are regarded as nested envs (NamedTuple)
function dynamics(env::Env)
    return function (x, p, t)
        x1 = x.env1
        x21 = x.env2.env21
        x22 = x.env2.env22
        ẋ1 = -x1 - env.gain*sum(x21 + x22)
        ẋ21 = -x21
        ẋ22 = -x22
        (; env1 = ẋ1, env2 = (; env21 = ẋ21, env22 = ẋ22))
    end
end

# for convenience
function make_env()
    env21, env22 = BaseEnv(reshape(collect(1:8), 2, 4)), BaseEnv(reshape(collect(9:16), 2, 4))
    env1 = BaseEnv(-1)
    env2 = Env2(env21, env22)
    gain = 2.0
    env = Env(env1, env2, gain)
    env
end
# register env; do it in global scope
__env = make_env()
__x0 = NestedEnvironments.initial_condition(__env)
@reg_env __env __x0

# test
function test()
    env = make_env()
    # initial condition
    # if you extend `NestedEnvironments.initial_condition` for all sub environments, then `NestedEnvironments.initial_condition(env::Env)` will automatically complete a readable initial condition as NamedTuple.
    x0 = NestedEnvironments.initial_condition(env)  # auto-completion of initial condition
    @show x0  # x0 = (env1 = -1, env2 = (env21 = [1 3 5 7; 2 4 6 8], env22 = [9 11 13 15; 10 12 14 16]))
    t0 = 0.0
    tf = 10.0
    tspan = (t0, tf)
    Δt = 0.01  # saveat; not numerical integration
    ts = t0:Δt:tf
    prob = ODEProblem(env, dynamics(env), x0, tspan)
    @time sol = solve(prob, Tsit5(), saveat=ts)
    # readable
    xs = sol.u |> Map(_x -> @readable _x) |> collect  # nested states
    @test xs[1].env1 == x0.env1
    @test xs[1].env2.env21 == x0.env2.env21
    @test xs[1].env2.env22 == x0.env2.env22
    # raw
    _x0 = @raw x0
    @show _x0  # _x0 = [-1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
    @test _x0 == sol.u[1]
end


text:
struct CachedFunction{F,I,O,S}
    f!::F
    input::I #prototype info for input
    output::O #prototype info for input
    x_cache::IdDict{DataType, S}
    closures::IdDict{DataType, Any}
    f_calls::Vector{Int64}
    lock::ReentrantLock
end

function CachedFunction(f!,_in::SIZE{N1},_out::SIZE{N2}) where {N1,N2}
    input = Prototype(Array{Float64,N1},(Float64,_in))
    output = Prototype(Array{Float64,N2},(Float64,_out))
    x_cache =  IdDict{DataType,Array{Any,N1}}()
    closures = IdDict{DataType, Any}()
    f_calls = [0]
    lock = ReentrantLock()
    return CachedFunction(f!,input,output,x_cache,closures,f_calls,lock)
end

function CachedFunction(f!,_in::T1,out::T2) where {T1,T2}
    input = prototype(_in)
    output = prototype(out)
    x_cache =  IdDict{DataType,T1.name.wrapper}()
    closures = IdDict{DataType, Any}()
    f_calls = [0]
    lock = ReentrantLock()
    return CachedFunction(f!,input,output,x_cache,closures,f_calls,lock)
end

function CachedFunction(f!,_in::Int,_out::Int)
    return CachedFunction(f!,(_in,),(_out,))
end

#creates an instance of the cache, using the storage_type and the size. uses type constructor

#generates an array of _eltype in

function make_closure!(f!,prototype,elem_type)
    out_cache = make_cache!(f!,prototype,elem_type)
    return x -> f!(out_cache, x)
end

function make_closure!(f!,prototype,elem_type,t)
    out_cache = make_cache!(f!,prototype,elem_type)
    return (x,t) -> f!(out_cache, x,t)
end

function make_closure!(f!,prototype,elem_type,t1,t2)
    out_cache = make_cache!(f!,prototype,elem_type)
    return (x,t1,t2) -> f!(out_cache, x,t1,t2)
end

eval_f(f::F, x) where {F} = f(x)
eval_f(f::F, x,t) where {F} = f(x,t)
eval_f(f::F, x,t1,t2) where {F} = f(x,t1,t2)

function (cfn::CachedFunction{F})(x) where {F}
    X = valtype(x)
    f = get!(() -> make_closure!(cfn.f!, cfn.output, X),cfn.closures,X)
    cfn.f_calls[] +=1
    return eval_f(f, x)
end

function (cfn::CachedFunction{F})(x,t) where {F}
    X = promote_type(valtype(x),eltype(t))
    f = get!(() -> make_closure!(cfn.f!, cfn.output, X,t),cfn.closures,X)
    cfn.f_calls[] +=1
    return eval_f(f, x, t)
end

function (cfn::CachedFunction{F})(x,t1,t2) where {F}
    X = promote_type(valtype(x),eltype(t1),eltype(t2))
    f = get!(() -> make_closure!(cfn.f!, cfn.output, X,t1,t2),cfn.closures,X)
    cfn.f_calls[] +=1
    return eval_f(f, x, t)
end

function allocate!(cfn::CachedFunction{F},x)  where {F}
    X = valtype(x)
    get!(() -> x,cfn.x_cache,X)
    get!(() -> make_closure!(cfn.f!, cfn.output, X),cfn.closures,X)
    return nothing
end

function allocate!(cfn::CachedFunction{F},x::Type{X})  where {F,X}
    get!(() -> make_cache!(cfn.f!, cfn.input, X),cfn.x_cache,X)
    get!(() -> make_closure!(cfn.f!, cfn.output, X),cfn.closures,X)
    return nothing
end
#@btime c($[1.0,2.0])

###calls###
"""
returns how many times has ´f´ been called.
"""
calls(f::CachedFunction) = f.f_calls[]

###output####

"""
returns the output cache stored in ´f´
Returns `f(x)`, it doesn' t store the input.
if your function has only one type cached, you can get the output cache without adding the type.
"""
function output(cfn::CachedFunction)
    if length(cfn.closures) == 1
        return last(first(cfn.closures)).out_cache
    elseif length(cfn.closures) == 0
        error("there aren't any output caches allocated in cached function")
    else
        error("there are more than one type of output cache stored in cached function.")
    end
end

function output(cfn::CachedFunction,x::Type{T}) where T
    if length(cfn.closures) == 0
        error("there aren't any output caches allocated in cached function.")
    else
        res = get(cfn.closures,T,nothing)
        if isnothing(res)
            TT = T.name.wrapper
        error("there aren't any output caches of type $TT in cached function.")
        else
            return res.out_cache
        end
    end
end

output(f::CachedFunction,x) = output(f,valtype(x))

###input####

function input(cfn::CachedFunction,x::Type{T}) where T
    if length(cfn.x_cache) == 0
        error("there aren't any input caches allocated in cached function")
    else
        res = get(cfn.x_cache,T,nothing)
        if isnothing(res)
        TT = T.name.wrapper
        error("there aren't any input caches of type $TT in cached function.")
        else
            return res
        end
    end
end

input(f::CachedFunction,x) = input(f,valtype(x))

####evaluate#####
"""
Evaluates the cached function at `x`.
Returns `f(x)`, it doesn' t store the input
"""
evaluate(f::CachedFunction,x) = f(x)

####evaluate!#####
##evaluates cached_f and stores x

"""
    evaluate!(f,x)

Force (re-)evaluation of the cached function at `x`.
Returns `f(x)` and stores the input value.

"""
function evaluate!(cfn::CachedFunction,x)
    X = valtype(x)
    out = cfn.out_size
    if haskey(cfn.x_cache,X)
        xx = get(cfn.x_cache,X,nothing)
        copyto!(xx,x)
    else
        xx = get!(cfn.x_cache,X,x)
    end
    f = get!(() -> make_closure!(cfn.f!, cfn.output, X),cfn.closures,X)
    cfn.f_calls[1] +=1
    return eval_f(f, x)
end

function Base.show(io::IO, fn::CachedFunction)
    n_methods = length(fn.closures)
    _methods = n_methods == 1 ? "method" : "methods"
    name = string(fn.f!)
    print(io, "cached version of $name (function with $n_methods cached $_methods)")
end

cached_methods(f::CachedFunction) = f.closures

#lock methods
Base.islocked(cfn::CachedFunction) = islocked(cfn.lock)
Base.lock(f, cfn::CachedFunction) = lock(f, cfn.lock)
Base.trylock(f, cfn::CachedFunction) = trylock(f, cfn.lock)