text:
const default_accuracy = 1.0e-12
const default_max_iterations = 1000
const near_int_threshold = 1.0e-6
const series_transition_threshold = 0.25

struct Diagnostics # this immutable type has no fields, so constructing it is essentially free
end

"""
    polylog(s, z)

Calculates the Polylogarithm function ``{Li}_s(z)`` defined by

``{Li}_s = \\sum_{n=1}^{\\infty} \\frac{z^n}{n^s}``

Uses double precision complex numbers (not arbitrary precision).
It's goal is an relative error bound 10^{-12}.

## Input Arguments
* ``s`` `::Complex`: the 'fractional' parameter
* ``z`` `::Complex`: the point at which to calculate it

It should also accept input arguments as `Real` or `Rational` or `Integer` but these aren't completely tested.

There are additional keywords, but these are currently intended for testing not use.

## Output Arguments
* ``Li_s(z)``: The result

## Examples
```jldoctest; setup = :(using Polylogarithms)
julia> polylog(0.35, 0.2)
0.23803890574407033
```
"""
function polylog(s::Number, z::Number;
                 level=1, # keep track of recursion
                 accuracy::Float64=default_accuracy,
                 min_iterations::Integer=0,
                 max_iterations::Integer=default_max_iterations)
    polylog(s, z, Diagnostics();
            level=level, accuracy=accuracy, min_iterations=min_iterations, max_iterations=max_iterations)[1] # just output the result
end

# this is the main version, but outputs diagnosstics, which I guess most people won't want
"""
    polylog(s, z, Diagnostics())

Calculates the Polylogarithm function ``{Li}_s(z)`` defined by

``{Li}_s = \\sum_{n=1}^{\\infty} \\frac{z^n}{n^s}``

Uses double precision complex numbers (not arbitrary precision).
It's goal is an relative error bound 10^{-12}.

## Input Arguments
* ``s::`` `Complex`: the 'fractional' parameter
* ``z`` `::Complex`: the point at which to calculate it
* `::Diagnostics`: use this to indicate that the output should include extra information

It should also accept input arguments as `Real` or `Rational` or `Integer` but these aren't completely tested.

There are additional keywords, but these are currently intended for testing not use.

## Output Arguments
* ``Li_s(z)``: The result
* ``n``:       The number of elements used in each series
* `series`:    The series used to compute results (4 = reciprocal)

## Examples
```jldoctest; setup = :(using Polylogarithms)
julia> polylog(0.35, 0.2, Diagnostics() )
(0.23803890574407033, 17, 1)
```
"""
function polylog(s::Number, z::Number, ::Diagnostics;
                 level=1, # keep track of recursion
                 accuracy::Float64=default_accuracy,
                 min_iterations::Integer=0,
                 max_iterations::Integer=default_max_iterations)
    tau_threshold = 1.0e-3
    μ = log(convert(Complex{Float64}, z)) # input z could be an integer or anything
    t = abs(μ/twoπ)
    T = 0.512 # the duplication formula seems to work if T=0.5 in the small associated wedge, but wyh risk it?
    if abs(μ) < 1.0e-14
        # deal with this case separately or it causes pain in Series 2 and 3, which should be resolvable, but its just easier here
        # there are lots of other special cases to add here eventually
        if real(s) > 1
            return SpecialFunctions.zeta(s), 0, 0
        else
            return typeof(z)(Inf), 0 , 0
        end
    elseif abs(z) <= 0.5 && abs(z) < t
        return polylog_series_1(s, z; accuracy=accuracy, min_iterations=min_iterations, max_iterations=max_iterations)
    elseif t <= T && ( abs(round(real(s))-s) > tau_threshold || real(s)<= 0 )
        return polylog_series_2(s, z; accuracy=accuracy, min_iterations=min_iterations, max_iterations=max_iterations)
    elseif t <= T
        return polylog_series_3(s, z; accuracy=accuracy, min_iterations=min_iterations, max_iterations=max_iterations)
#    elseif t <= 2.0
    else
        # println("  main level $level, z=$z,  abs(μ)/2π = ", abs(log(z))/twoπ  )
        return polylog_duplication(s, z; level=level, accuracy=accuracy, min_iterations=min_iterations, max_iterations=max_iterations)
    # elseif abs(z) > 1
    #     return polylog_reciprocal(s, z; accuracy=accuracy, min_iterations=min_iterations, max_iterations=max_iterations)
    # else
    #     throw(DomainError(z, "bad z value $z with |z|=$(abs(z)) and |log(z)/(twoπ)|=$(abs(μ/twoπ))"))
    end
    # we could have a lot more special cases here, particularly for integer input
    # to make the code faster for these cases, but at the moment may use such values for testing
end

# old version without square root
    # if abs(z) <= 0.5 && abs(z) < abs(μ/twoπ)
    #     return polylog_series_1(s, z; accuracy=accuracy, min_iterations=min_iterations, max_iterations=max_iterations)
    # elseif abs(μ/twoπ) < series2_threshold && ( abs(round(real(s))-s) > tau_threshold || real(s)<= 0 )
    #     return polylog_series_2(s, z; accuracy=accuracy, min_iterations=min_iterations, max_iterations=max_iterations)
    # elseif abs(μ/twoπ) < series2_threshold
    #     return polylog_series_3(s, z; accuracy=accuracy, min_iterations=min_iterations, max_iterations=max_iterations)
    # elseif abs(z) > 1
    #     return polylog_reciprocal(s, z; accuracy=accuracy, min_iterations=min_iterations, max_iterations=max_iterations)
    # else
    #     throw(DomainError(z, "bad z value $z with |z|=$(abs(z)) and |log(z)/(twoπ)|=$(abs(μ/twoπ))"))
    # end

# function polylog(s::Number, z::Number, accuracy::Real=default_accuracy)
#     if z ≈ 1.0
#         if real(s) > 1
#             return zeta(s)
#         else
#             return Inf
#         end
#     elseif z ≈ -1.0
#         return -eta(s)
#     elseif s ≈ 0.0
#         return z ./ (1-z)
#     elseif abs(z) <= T
#         ifconfig
#         return polylog_direct(s, z, accuracy)
#     elseif abs(z) >= 1/T && isinteger(s) && real(s) < 0
#         # use reciprocal formula to calculate in terms of Li_n(1/z)
#         # but note for negative integer s, it collapses to something small
#         return -(-1.0)^s .*polylog_direct(s, 1/z, accuracy)
#     elseif  abs(z) >= 1/T
#         # use reciprocal formula to calculate in terms of Li_s(1/z)
#         twopi = 2π
#         z = convert(Complex{Float64}, z)
#         G = (twopi*im)^s * zeta( 1-s, 0.5 + log(-z)/(twopi*im) ) /  gamma(s)
#         F = complex(-1.0)^s * polylog_direct(s, 1/z, accuracy)

#         A = twopi*im*log(z)^(s-1)/(gamma(s))
#         if ( isreal(z) && real(z)>=1 )
#             Θ = 1
#         else
#             Θ = 0
#         end
#         # println("G = $G, F=$F, Θ=$Θ, A=$A")
#         return ( G - F - Θ*A )
#     else
#         # power series around mu=0, for z = e^mu
#         polylog_series_mu(s, z, accuracy)
#     end
# end

####################################################################
#### these are component functions and aren't exported at this point
#### note that for consistency they all have keywords arguments like "accuracy" but
#### these aren't intended for general use, just for testing (at the moment)

# # calculate using the relationship to the Hurwitz zeta function
# function polylog_zeta(s::Number, z::Number, accuracy=default_accuracy)
#     # compute using the Hurwitz-zeta function identity
#     #   N.B. this doesn't seem to work as well as you might think
#     x = im * (log(convert(Complex{Float64}, -z)) / twoπ)
#     ss = 1-s
#     ip = im^ss
#     return ( SpecialFunctions.gamma(ss)/twoπ^(ss) ) * (ip * SpecialFunctions.zeta(ss, 0.5+x) + conj(ip) * SpecialFunctions.zeta(ss, 0.5-x))
# end

# calculate using the duplication formula
function polylog_duplication(s::Number, z::Number;
                             level=0, # keep track of recursion
                             accuracy::Float64=default_accuracy,
                             min_iterations::Integer=0,
                             max_iterations::Integer=default_max_iterations)
    z = convert(Complex{Float64}, z)
    f = min( 0.5, 2.0^(1-real(s)) )
    # println("  dup level $level, z=$z,  abs(μ)/2π = ", abs(log(z))/twoπ  )
    (Li1, k1, series1) = polylog(s,  sqrt(z), Diagnostics();
                                 level=level+1, accuracy=f*accuracy,
                                 min_iterations=min_iterations, max_iterations=max_iterations)
    (Li2, k2, series2) = polylog(s, -sqrt(z), Diagnostics();
                                 level=level+1, accuracy=f*accuracy,
                                 min_iterations=min_iterations, max_iterations=max_iterations)
    if typeof(s) <: Real
        s = convert(Float64, s) # convert s into a double
    elseif typeof(s) <: Complex
        s = convert(Complex{Float64}, s) # convert s into doubles
    end
    return (2^(s-1) * ( Li1 + Li2 ), k1 + k2, 10+series1+series2)
end

# # calculate using the reciprocal formula
# function polylog_reciprocal(s::Number, z::Number;
#                             accuracy::Float64=default_accuracy,
#                             min_iterations::Integer=0,
#                             max_iterations::Integer=default_max_iterations)
#     # z = convert(Complex{Float64}, z)
#     if abs(z) <= 1
#         throw(DomainError(z, "only use this function for |z|>1, and pref |z| > 2"))
#     end
#     # if abs(z) < 2
#     #     warn("Slow convergence for  |z| < 2")
#     # end
#     if abs(s) < 0.1*accuracy
#         return (z/(1-z), 0) # use the identity
#     elseif real(s) < 0 &&  abs(imag(s)) < 0.1*accuracy && abs( round(s) - s ) < 0.1*accuracy
#         G = 0.0 # pole of the Gamma function
#         A = 0.0
#     else
#         # G = (twoπ*im)^s * SpecialFunctions.zeta( 1-s, 0.5 + log(complex(-z))/(twoπ*im) ) /  SpecialFunctions.gamma(s)
#         # A = twoπ*im*log(z)^(s-1) / SpecialFunctions.gamma(s)
#         tmp = exp( s*log(twoπ*im) - SpecialFunctions.loggamma(complex(s)) )  # (twoπ*im)^s /  SpecialFunctions.gamma(s)
#         G = tmp * SpecialFunctions.zeta( 1-s, 0.5 + log(complex(-z))/(twoπ*im) )
#         A = twoπ*im*log(z)^(s-1) / SpecialFunctions.gamma(s)
#     end
#     # accuracy of overall result depends on size of total, which includes these other parts
#     (Li, k, series) = polylog_series_1(s, 1/z; accuracy=0.1*accuracy,
#                                min_iterations=min_iterations, max_iterations=max_iterations, existing_total=G)
#     F = complex(-1.0)^s * Li
#     if ( imag(z) == 0 ) &&  ( real(z) >= 1 )
#         Θ = 1.0
#     else
#         Θ = 0.0
#     end
#     # println("G = $G, F=$F, Θ=$Θ, A=$A")
#     return ( G - F - Θ*A, k, 3+series )
# end

# calculate using direct definition
function polylog_series_1(s::Number, z::Number;
                          accuracy::Float64=default_accuracy,
                          min_iterations::Integer=0,
                          max_iterations::Integer=default_max_iterations,
                          existing_total::Number=0.0)
    # prolly should convert z to a double or complex-double
    if abs(z) > 1 || ( abs(z) ≈ 1  && real(s) <= 2)
        throw(DomainError(z))
    end
    if abs(z) > 1/2
        throw(DomainError(z, "Slow convergence for  |z| > 1/2"))
    end
    total = 0.0
    converged = false
    a = z
    k = 0
    if real(s) < 0
        min_iterations = ceil( real(s) / log(abs(z)) )
    end
    while k<=max_iterations && ~converged
        k = k+1
        total += a
        a *= z * ( k/(k+1.0) )^s
        # println("   total = $total")
        if k > min_iterations && abs(a)/abs(total) < 0.5*accuracy
            converged = true
        end
    end
    return (total, k, 1)
end

# calculate using power series around μ = log(z) = 0
# this should not be used near positive integer values of s, but we allow it here in order to test
function polylog_series_2(s::Number, z::Number;
                          accuracy::Float64=default_accuracy,
                          min_iterations::Integer=0,
                          max_iterations::Integer=default_max_iterations )
    μ = log(convert(Complex{Float64}, z)) # input z could be an integer or anything
    if typeof(s) <: Real
        s = convert(Float64, s) # convert s into a double
    elseif typeof(s) <: Complex
        s = convert(Complex{Float64}, s) # convert s into doubles
    end
    # println("μ = $μ")
    if abs(μ) > twoπ
        throw(DomainError(z, "we need |log(z)|< 2π for this series"))
    end
    # if real(s) > 0
        # min_iterations = ceil( real(s) ) + 1
    # else
        # min_iterations = ceil( -real(s) ) + 1
        # min_iterations = ceil( real(s) / log(abs( log(z)/twoπ ) ) )
    # end
    oneminuss = 1.0 - s
    total = SpecialFunctions.gamma(oneminuss) * (-μ)^(- oneminuss)
    # total = exp( SpecialFunctions.loggamma(oneminuss) - oneminuss*log(-μ) )
    converged = false
    tmp = 1
    k = 0
    a = Inf
    a_2 = Inf
    # A = abs( 2.0*twoπ^real(s) * exp(abs(imag(π*s)))  )
    # this doesn't work if z=1.0, and hence μ=0, even when that converges, but should already be delt with
    while k<=max_iterations && ~converged
        a_3 = a_2
        a_2 = a
        a = tmp * SpecialFunctions.zeta(s - k)
        total += a
        tmp *= μ/(k+1)
        if k > min_iterations &&
            abs(a)/abs(total) < 0.5*accuracy &&
            abs(a_2)/abs(total) < 0.5*accuracy &&
            abs(a_2) > abs(a)
            # abs( A * (k-real(s))^real(-s) * (μ/twoπ)^k )/abs(total) < accuracy
            # && abs( 2*twoπ^real(s) * (μ/twoπ)^k )/abs(total) < accuracy
            # the stopping rule should be implemented more efficiently as part of the calculation above
            converged = true
        end
        k = k + 1
    end
    # get correct value along the branch
    if isreal(z) && real(z)>=1
        # total -= 2*π*im*μ^(s-1)/SpecialFunctions.gamma(s)
        total -= exp( log(twoπ*im) + (s-1)*log(μ) - SpecialFunctions.loggamma(s) )
    end
    return (total, k, 2)
end

function c_closed(n::Integer, j::Integer,  ℒ::Number)
    d2 = SpecialFunctions.digamma(n+1) - ℒ
    if j==0
        return harmonic(n) - ℒ
    elseif j==1
        # Wood:+stieltjes(1) - d2^2/2 + π^2/6 + SpecialFunctions.polygamma(1,n+1)/2
        return -stieltjes(1) - d2^2/2 - π^2/6 + SpecialFunctions.polygamma(1,n+1)/2
    elseif j==2
        # Wood:stieltjes(2)   + d2^3/6 + d2*( π^2/6 + SpecialFunctions.polygamma(1,n+1)/2 ) + SpecialFunctions.polygamma(2,n+1)/6
        return stieltjes(2)/2 + d2^3/6 + d2*( π^2/6 - SpecialFunctions.polygamma(1,n+1)/2 ) + SpecialFunctions.polygamma(2,n+1)/6
   end
end

function Q_closed(n::Integer, τ::Number, ℒ::Number; n_terms::Integer=3)
    # τ is the distance from the pole s=n>0, ℒ = log(-μ) = log(-log( z ))
    max_n_terms = 3
    if n_terms < 1 || n_terms > max_n_terms
        throw(DomainError(n_terms))
    end
    return sum( c_closed.(n, 0:n_terms-1,  ℒ) .* τ.^(0:n_terms-1) )
end

function Q(n::Integer, τ::Number, ℒ::Number; n_terms::Integer=5) # Crandall,2012, p.35
    # τ is the distance from the pole s=n>0, ℒ = log(-μ) = log(-log( z ))
    if abs(τ) <= 1.0e-14
        # if really close to the integer, then ignore the extra terms
        return c_closed(n, 0,  ℒ)
    else
        max_n_terms = 7
        if n_terms < 1 || n_terms > max_n_terms
            throw(DomainError(n_terms))
        end
        if n_terms <= 3
            # use the direct method in this case
            return Q_closed(n, τ, ℒ; n_terms=n_terms)
        end
        return sum( c_crandall.(n, 0:n_terms-1,  ℒ) .* τ.^(0:n_terms-1) )
    end
end

function c_crandall(k::Integer, j::Integer,  ℒ) # Crandall,2012, p.35
    return (-1)^j*stieltjes(j)/SpecialFunctions.gamma(j+1) - b_crandall(k,j+1,ℒ)
end

function b_crandall(k::Integer, j::Integer,  ℒ) # Crandall,2012, p.36
    total = 0
    for q=0:j
        for t=0:j-q
            p = j-q-t
            a1 = ℒ^p / SpecialFunctions.gamma(p+1)
            a2 = (-1)^t * f_crandall(k,q) # Bailey and Borwein, 2015 correct Crandall (t+q - > t)
            a3 = g_crandall(t)/SpecialFunctions.gamma(t+1)
            total += a1 * a2 * a3
        end
    end
    return total
end

const gamma_t = [1.0,-0.5772156649015315,1.9781119906559432,-5.44487445648531,23.561474084025583,-117.83940826837748,715.0673625273184,-5019.848872629852,40243.62157333573,-362526.2891146549,3.627042412756892e6,-3.990708415143132e7,4.7894329176518273e8,-6.226641351546061e9,8.717563381070836e10]
function g_crandall(t::Integer) # Crandall,2012, p.17
    # t derivate of Gamma function at 1
    # see "gamma_derivatives.jl" for derivations of these numbers
    if t<0
        throw(DomainError(t))
    elseif t>14
        throw(DomainError(t, "only calculate the 1st 14 derivatives"))
    else
        return gamma_t[t+1]
    end
end

function f_crandall(k::Integer, q::Integer) # Crandall,2012, p.36
    # Crandall omits the k=0 case. but this is given in Bailey and Borwein and matches other text
    if k == 0 && q == 0
        return 1
    elseif k == 0
        return 0
    elseif q == 0
        return 1
    # elseif q == 1
    #     return -harmonic(k)
    # elseif q == 2
    #     return (harmonic(k)^2 + harmonic(k,2))/2
    else
        h = 0:q
        return sum( (-1).^h .* f_crandall.(k-1, q .- h)  ./ k.^h )
    end
end

# For the special case that s is near a postive integer n>0
# Calculate in a power series around z=1, and s=n
function polylog_series_3(s::Number, z::Number;
                          accuracy::Float64=default_accuracy,
                          min_iterations::Integer=0,
                          max_iterations::Integer=default_max_iterations,
                          n_terms::Integer=5 )
    μ = log(convert(Complex{Float64}, z))
    if abs(μ) > twoπ
        throw(DomainError(z, "does not converge for abs(ln(z)) > twoπ"))
    end
    if real(s)<=0.5
        throw(DomainError(s, "for this function s should be near a positive integer"))
    end
    # this series assumes s is near a positive integer
    n = Int(round(real(s)))
    τ = s - n
    # if real(s) > 0
    #     min_iterations = ceil( real(s) )
    # end
    ℒ = log(complex(-μ))  # '\u2112'
    # total = μ^(n-1)*Q(n-1, τ, ℒ; n_terms=n_terms)/SpecialFunctions.gamma(n)
    total = exp( (n-1)*log(μ) + log(Q(n-1, τ, ℒ; n_terms=n_terms))  - SpecialFunctions.loggamma(n) )
    converged = false
    a = Inf
    a_2 = Inf
    tmp = 1
    k = 0
    while k<=max_iterations && ~converged
        if n - k != 1
            a_2 = a
            a = tmp * SpecialFunctions.zeta(s - k)
            total += a
        end
        tmp *= μ/(k+1)
        if k > min_iterations &&
            abs(a)/abs(total) < 0.5*accuracy &&
            abs(a_2)/abs(total) < 0.5*accuracy &&
            abs(a_2) > abs(a)
            # abs( (μ/twoπ)^k )/abs(total) < 0.05*accuracy
            converged = true
        end
        k = k + 1
    end
    # get correct value along the branch
    if isreal(z) && real(z)>=1
        total -= 2*π*im*μ^(s-1)/SpecialFunctions.gamma(s)
    end
    return (total, k, 3)
end

# # asymptotic expansion
# function polylog_series_4(s::Number, z::Number;
#                           accuracy::Float64=default_accuracy,
#                           min_iterations::Integer=0,
#                           max_iterations::Integer=default_max_iterations,
#                           existing_total::Number=0.0)
#     x = log(-convert(Complex{Float64}, z))
#     total = 0.0
#     converged = false
#     k = 0
#     a = Inf
#     while k<=max_iterations && ~converged
#         old_a = a
#         a = (-1)^k * (1.0 - 2.0^(1-2*k)) * twoπ^(2*k) * bernoulli(2*k,0.0) * x^(s-2*k) / ( SpecialFunctions.gamma(2*k+1) * SpecialFunctions.gamma(s+1-2*k) )
#         if abs(a) < abs(old_a)
#             total += a
#         else
#             converged = true
#         end
#         k = k+1
#     end
#     return (total, k, 4)
# end

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

# # Special case s=integer for use in testing etc.
#     function polylog(n::Integer, z::Number;
#                       accuracy::Float64=default_accuracy, max_iterations::Integer=default_max_iterations )
#     L = Int(ceil(-log10(accuracy)*log2(10))) # revisit this limit
#     k = 0
#     if n>1
#         μ = log(convert(Complex{Float64}, z))
#         if abs(μ) > twoπ
#             throw(DomainError(z))
#         end
#         total = μ^(n-1)*(harmonic(n-1) - log(-μ))/SpecialFunctions.gamma(n)
#         tmp = 1
#         for m=0:L
#             if n - m != 1
#                 total += tmp * zeta(n - m)
#             end
#             tmp *= μ/(m+1)
#             if abs(tmp)/abs(total) < 1.0e-30
#                 break
#             end
#         end
#         if  isreal(z) && real(z)>=1
#             total -= 2*pi*im*μ^(s-1)/SpecialFunctions.gamma(n)
#         end
#     elseif n==1
#         total = -log(complex(1-z))
#     elseif n==0
#         total = z / (1-z)
#     elseif n==-1
#         total = z / (1-z)^2
#     elseif n<-1
#         # Crandall's 1.5 for s integer
#         total = factorial(-n) * (-μ)^(n-1)
#         tmp = 1
#         for k=0:L
#             # total -= μ^k * bernoulli(k-n+1, 0.0) / ( gamma(k+1)*(k-n+1) )
#             total -= tmp * bernoulli(k-n+1, 0.0) / (k-n+1)
#             tmp *= μ/(k+1)
#             if abs(tmp)/abs(total) < 1.0e-30
#                     break
#             end
#         end
#     else
#         error("Should not get this case")
#     end

#     if isreal(z) && real(z)>=1
#         total -= 2*pi*im*μ^(s-1)/gamma(s)
#     end
#     return (total,k)
# end


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

IndexStyle(A::AbstractQuasiArray) = IndexStyle(typeof(A))
IndexStyle(::Type{<:AbstractQuasiArray}) = IndexCartesian()

IndexStyle(A::AbstractQuasiArray, B::AbstractQuasiArray) = IndexStyle(IndexStyle(A), IndexStyle(B))
IndexStyle(A::AbstractArray, B::AbstractQuasiArray) = IndexStyle(IndexStyle(A), IndexStyle(B))
IndexStyle(A::AbstractQuasiArray, B::AbstractArray) = IndexStyle(IndexStyle(A), IndexStyle(B))
IndexStyle(A::AbstractQuasiArray, B::AbstractQuasiArray...) = IndexStyle(IndexStyle(A), IndexStyle(B...))
IndexStyle(A::AbstractQuasiArray, B::AbstractArray...) = IndexStyle(IndexStyle(A), IndexStyle(B...))

function promote_shape(a::AbstractQuasiArray, b::AbstractQuasiArray)
    promote_shape(axes(a), axes(b))
end

const QuasiIndices{N} = NTuple{N,Union{AbstractQuasiVector{<:Number},AbstractVector{<:Number}}}
function promote_shape(a::QuasiIndices, b::QuasiIndices)
    if length(a) < length(b)
        return promote_shape(b, a)
    end
    for i=1:length(b)
        if a[i] != b[i]
            throw(DimensionMismatch("dimensions must match"))
        end
    end
    for i=length(b)+1:length(a)
        if a[i] != 1:1
            throw(DimensionMismatch("dimensions must match"))
        end
    end
    return a
end

# check for valid sizes in A[I...] = X where X <: AbstractQuasiArray
# we want to allow dimensions that are equal up to permutation, but only
# for permutations that leave array elements in the same linear order.
# those are the permutations that preserve the order of the non-singleton
# dimensions.
function setindex_shape_check(X::AbstractQuasiArray, I::Integer...)
    li = ndims(X)
    lj = length(I)
    i = j = 1
    while true
        ii = length(axes(X,i))
        jj = I[j]
        if i == li || j == lj
            while i < li
                i += 1
                ii *= length(axes(X,i))
            end
            while j < lj
                j += 1
                jj *= I[j]
            end
            if ii != jj
                throw_setindex_mismatch(X, I)
            end
            return
        end
        if ii == jj
            i += 1
            j += 1
        elseif ii == 1
            i += 1
        elseif jj == 1
            j += 1
        else
            throw_setindex_mismatch(X, I)
        end
    end
end

setindex_shape_check(X::AbstractQuasiArray) =
    (length(X)==1 || throw_setindex_mismatch(X,()))

setindex_shape_check(X::AbstractQuasiArray, i::Integer) =
    (length(X)==i || throw_setindex_mismatch(X, (i,)))

setindex_shape_check(X::AbstractQuasiArray{<:Any,1}, i::Integer) =
    (length(X)==i || throw_setindex_mismatch(X, (i,)))

setindex_shape_check(X::AbstractQuasiArray{<:Any,1}, i::Integer, j::Integer) =
    (length(X)==i*j || throw_setindex_mismatch(X, (i,j)))

function setindex_shape_check(X::AbstractQuasiArray{<:Any,2}, i::Integer, j::Integer)
    if length(X) != i*j
        throw_setindex_mismatch(X, (i,j))
    end
    sx1 = length(axes(X,1))
    if !(i == 1 || i == sx1 || sx1 == 1)
        throw_setindex_mismatch(X, (i,j))
    end
end

to_index(I::AbstractQuasiArray{Bool}) = LogicalIndex(I)
to_index(I::AbstractQuasiArray) = I
to_index(I::AbstractQuasiArray{<:Union{AbstractArray, Colon}}) =
    throw(ArgumentError("invalid index: $I of type $(typeof(I))"))

to_quasi_index(i::Number) = i
to_quasi_index(i) = Base.to_index(i)
to_index(A::AbstractQuasiArray, i) = to_quasi_index(i)

LinearIndices(A::AbstractQuasiArray) = LinearIndices(axes(A))

"""
   Inclusion(domain)

Represents the inclusion operator of a domain (that is, a type that overrides in)
as an AbstractQuasiVector. That is, if `v = Inclusion(domain)`, then
`v[x] == x` if `x in domain`, otherwise it throws a `DomainError`.

Inclusions are useful for turning domains into axes. They also serve the same
role as `Slice` does for offset arrays.
"""
struct Inclusion{T,AX} <: AbstractQuasiVector{T}
    domain::AX
end
Inclusion{T}(domain) where T = Inclusion{T,typeof(domain)}(domain)
Inclusion{T}(S::Inclusion) where T = Inclusion{T}(S.domain)
Inclusion{T}(S::Slice) where T = Inclusion{T}(S.indices)
Inclusion(domain) = Inclusion{eltype(domain)}(domain)
Inclusion(S::Inclusion) = S
Inclusion(S::Slice) = Inclusion(S.indices)

convert(::Type{Inclusion}, d::Inclusion) = d
convert(::Type{Inclusion{T}}, d::Inclusion) where T = Inclusion{T}(d)
convert(::Type{AbstractVector}, d::Inclusion{<:Any,<:AbstractVector}) =
    convert(AbstractVector, d.domain)
convert(::Type{AbstractArray}, d::Inclusion{<:Any,<:AbstractVector}) =
    convert(AbstractArray, d.domain)
Vector(d::Inclusion{<:Any,<:AbstractVector}) = Vector(d.domain)
Array(d::Inclusion{<:Any,<:AbstractVector}) = Array(d.domain)

copy(d::Inclusion) = d

==(A::Inclusion, B::Inclusion) = A.domain == B.domain
domain(A::Inclusion) = A.domain
domain(A::AbstractUnitRange) = A
axes(S::Inclusion) = (S,)
unsafe_indices(S::Inclusion) = (S,)
axes1(S::Inclusion) = S
axes(S::Inclusion{<:Any,<:OneTo}) = (S.domain,)
unsafe_indices(S::Inclusion{<:Any,<:OneTo}) = (S.domain,)
axes1(S::Inclusion{<:Any,<:OneTo}) = S.domain

first(S::Inclusion) = first(S.domain)
last(S::Inclusion) = last(S.domain)
size(S::Inclusion) = (cardinality(S.domain),)
length(S::Inclusion) = cardinality(S.domain)
unsafe_length(S::Inclusion) = cardinality(S.domain)
cardinality(S::Inclusion) = cardinality(S.domain)
getindex(S::Inclusion{T}, i::Number) where T =
    (@_inline_meta; @boundscheck checkbounds(S, i); convert(T,i))
getindex(S::Inclusion{T}, i::AbstractVector{<:Number}) where T =
    (@_inline_meta; @boundscheck checkbounds(S, i); convert(AbstractVector{T},i))
getindex(S::Inclusion, i::Inclusion) =
    (@_inline_meta; @boundscheck checkbounds(S, i); copy(S))
getindex(S::Inclusion, ::Colon) = copy(S)
show(io::IO, r::Inclusion) = print(io, "Inclusion(", r.domain, ")")
iterate(S::Inclusion, s...) = iterate(S.domain, s...)

in(x, S::Inclusion) = x in S.domain

checkindex(::Type{Bool}, inds::Inclusion, i::Number) = i ∈ inds.domain
checkindex(::Type{Bool}, inds::Inclusion, ::Colon) = true
checkindex(::Type{Bool}, inds::Inclusion, ::Inclusion) = true
function checkindex(::Type{Bool}, inds::Inclusion, I::AbstractArray)
    @_inline_meta
    b = true
    for i in I
        b &= checkindex(Bool, inds, i)
    end
    b
end

function checkindex(::Type{Bool}, inds::Inclusion, r::AbstractRange)
    @_propagate_inbounds_meta
    isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
end
checkindex(::Type{Bool}, indx::Inclusion, I::AbstractVector{Bool}) = indx == axes1(I)
checkindex(::Type{Bool}, indx::Inclusion, I::AbstractArray{Bool}) = false

ext: jl


text:
export PDataMA57
export preprocessMA57, AInv, reconstructH
export preprocessMA57_2

mutable struct PDataMA57{T} <: PDataFact{T}
    L :: SparseMatrixCSC{T, Int64} # sparse L
    D :: SparseMatrixCSC{T, Int64} # block diagonal 1X1 and 2X2
    #P::Array{T,2}  # permutation matrix: should be vector
                         # LDL' = PHP'
    pp :: Array{Int64, 1}   # permutation vector LDL' = H[pp,pp]
                         # P*v = v[pp]
                         # P'*v = v[invperm(pp)]
    s :: Array{T, 1}  # Scaling matrix S=diagm(s)
    Δ :: Array{T, 1}  # diagonal, eigenvalues of D
    Q :: SparseMatrixCSC{T, Int64}  # orthogonal matrix, eigenvectors of D:  should be sparse
                         # QΔQ'  =  D
    g̃ :: Array{T, 1}  # transformed gradient
    λ :: T
    success :: Bool        # previous iteration was successfull
    OK :: Bool             # preprocess success

    PDataMA57() = new{Nothing}()
    PDataMA57(L, D, pp, s, Δ, Q, g̃, l, success, OK) = new{eltype(g̃)}(L, D, pp, s, Δ, Q, g̃, l, success, OK)
end

function preprocessMA57(H, g, params::Tparams, n1, n2)
    M     = Ma57
    L     = SparseMatrixCSC{Float64, Int64}
    D57   = SparseMatrixCSC{Float64, Int64}
    pp    = Array{Int64, 1}
    s     = Array{Float64}
    ρ     = Float64
    ncomp = Int64

    H57 = convert(SparseMatrixCSC{Float64, Int64}, H)  #  Hard coded Int
    try
        M = Ma57(H57)#, print_level = -1)
        ma57_factorize(M)
    catch
 	println("*******   Problem in MA57_0")
        M = Ma57(H57,print_level = -1)
        ma57_factorize(M)
        res = PDataMA57()
        res.OK = false
        return res
    end
    try
        (L, D57, s, pp) = ma57_get_factors(M)
    catch
        println("*******   Problem after MA57_1")
        # println(" Cond(H) = $(cond(full(H)))")
        res = PDataMA57_0()
        res.OK = false
        return res
    end

    #################  Future object BlockDiag operator?
    vD1 = diag(D57)       # create internal representation for block diagonal D
    vD2 = diag(D57, 1)     #
    vQ1 = ones(length(vD1))       # vector representation of orthogonal matrix Q
    vQ2 = zeros(length(vD2))      #
    vQ2m = zeros(length(vD2))     #
    veig = copy(vD1)      # vector of eigenvalues of D, initialized to diagonal of D
                          # if D diagonal, nothing more will be computed

    i = 1;
    while i < length(vD1)
        if vD2[i] == 0.0
            i += 1
        else
            mA = [vD1[i] vD2[i]; vD2[i] vD1[i + 1]] #  2X2 submatrix
            # DiagmA, Qma = eig(mA)                   #  spectral decomposition of mA
            X = eigen(mA)
            DiagmA = X.values
            Qma = X.vectors
            veig[i] = DiagmA[1]
            vQ1[i] = Qma[1, 1]
            vQ2[i] = Qma[1, 2]
            vQ2m[i] = Qma[2, 1]
            vQ1[i + 1] = Qma[2, 2]
            veig[i + 1] = DiagmA[2]
            i += 2
        end
    end

    Q = sparse(SparseArrays.spdiagm(0 => vQ1, 1 => vQ2m, -1 => vQ2))
    Δ = veig

    l_m, = findmin(Δ)
    sg = s .* g
    L = SparseMatrixCSC{Float64,Int64}(L)  #### very important, the \ command doesn't work with Int32
    ĝ = L \ (sg[pp])
    g̃ = Q' * ĝ

    n_g = norm(g)
    λ =  max(-l_m, 0.0)
    return  PDataMA57(L, D57, pp, s, Δ, Q, g̃, λ, true, true)
end

function preprocessMA57_2(H, g, params::Tparams, n1, n2)
    M     = Ma57
    L     = SparseMatrixCSC{Float64, Int64}
    D57   = SparseMatrixCSC{Float64, Int64}
    pp    = Array{Int64, 1}
    s     = Array{Float64}
    ρ     = Float64
    ncomp = Int64

    H57 = convert(SparseMatrixCSC{Float64, Int64}, H)  #  Hard coded Int
    try
        M = Ma57(H57)#, print_level = -1)
        ma57_factorize(M)
    catch
 	println("*******   Problem in MA57_0")
        M = Ma57(H57,print_level = -1)
        ma57_factorize(M)
        res = PDataMA57()
        res.OK = false
        return res
    end
    M.control.cntl[1] = 1.0
    M.control.icntl[7] = 0
    # @show M.control.cntl[1]
    try
        (L, D57, s, pp) = ma57_get_factors(M)
    catch
        println("*******   Problem after MA57_1")
        # println(" Cond(H) = $(cond(full(H)))")
        res = PDataMA57_0()
        res.OK = false
        return res
    end

    #################  Future object BlockDiag operator?
    vD1 = diag(D57)       # create internal representation for block diagonal D
    vD2 = diag(D57, 1)     #
    vQ1 = ones(length(vD1))       # vector representation of orthogonal matrix Q
    vQ2 = zeros(length(vD2))      #
    vQ2m = zeros(length(vD2))     #
    veig = copy(vD1)      # vector of eigenvalues of D, initialized to diagonal of D
                          # if D diagonal, nothing more will be computed

    i = 1;
    while i < length(vD1)
        if vD2[i] == 0.0
            i += 1
        else
            mA = [vD1[i] vD2[i]; vD2[i] vD1[i + 1]] #  2X2 submatrix
            # DiagmA, Qma = eig(mA)                   #  spectral decomposition of mA
            X = eigen(mA)
            DiagmA = X.values
            Qma = X.vectors
            veig[i] = DiagmA[1]
            vQ1[i] = Qma[1, 1]
            vQ2[i] = Qma[1, 2]
            vQ2m[i] = Qma[2, 1]
            vQ1[i + 1] = Qma[2, 2]
            veig[i + 1] = DiagmA[2]
            i += 2
        end
    end

    Q = sparse(SparseArrays.spdiagm(0 => vQ1, 1 => vQ2m, -1 => vQ2))
    Δ = veig

    l_m, = findmin(Δ)
    sg = s .* g
    L = SparseMatrixCSC{Float64,Int64}(L)  #### very important, the \ command doesn't work with Int32
    ĝ = L \ (sg[pp])
    g̃ = Q' * ĝ

    n_g = norm(g)
    λ =  max(-l_m, 0.0)
    return  PDataMA57(L, D57, pp, s, Δ, Q, g̃, λ, true, true)
end

function AInv(X :: PDataMA57, d̃ ::  Array{T, 1}) where T
    d̂ = X.Q * d̃
    u = X.L' \ d̂
    return u[invperm(X.pp)] .* X.s
end

function reconstructH(X :: PDataMA57)
    A = X.L * X.D * X.L'
    return A[X.pp.X.pp]
end


text:
module Interior
"""
I wrote this when, for some reason, I thought I needed to consider thermal
gradients within a chondrule. In fact, the thermal diffusivity of most minerals
is about 1 mm^2/s (e.g. Drury,Geothermics Volume 16, Issue 2, 1987, Pages 105-115),
which means that any mm-sized chondrule will thermally equilibrate in about a
second (!!!) so, everything in here is complete overkill and totally unecessary.

I've kept this here in case I find a use for it at some later time. It's really
just 1D Crank-Nicolson in spherical geometry based on FV formalism.
"""

export T, mote, heat!, cool!, cp, Tμ, speck, coolrad!, evaporate!#, hμ
# models the interior of a chondrule grain (thermal transport, melting, evaporation)
include("Consts.jl")
using .Consts: σSB, ρₛ, cₚ, κ

include("DragAndHeat.jl")
using .DragAndHeat: evap

using LinearAlgebra
const f64 = Float64

J = 1.0e7
kJ = 1.0e3 * J
mol = 24.305 + 28.0855 + 3 * 15.999 # obvs, mol wt may differ a bit bc isotope ratio not same as Earth
ε = 0.9 # emissivity of chondrule (total wild-ass guess)
cp0 = 1.2e7 # will override something?
Δh = 5.0e9 # heat of fusion
Tℓ = 273.15 + 0.5*(1475+1600) # melting pt; just averaging range in Lofgren's old report
h1 = Tℓ * cp0
h2 = Tℓ * cp0 + Δh
ΔH0 = 2812 * kJ / mol # heat of vaporization (estimate)

function cp(T)
    #= specific heat capacity (erg/g/K)

     This is the heat capacity for enstatite (MgSiO₃), from 398K-1000K, from:
        Kenneth M. Krupka
        Bruce S. Hemingway
        Richard A. Robie
        Derrill M. Kerrick
        American Mineralogist (1985) 70 (3-4): 261–271.

    If you look at the curve, it starts to look like nonsense below 150K but it
    at least "looks" reasonable up to 2000K. Not that you can trust it, but you
    can use it for development purposes anyway. In the future, I'll need a better fit.
    =#
    # in J/mol/K
    T = max(min(1000,T),398) # just artificially flattening the curve, for now
    cpJmol = 350.7 - 0.1473*T + 1.679e6*T^(-2) - 4296 * T^(-0.5) + 5.826e-5*T^2
    # cpJmol is about 127.5 J/mol/K at 1000K, which is equivalent to 1.27e7 erg/g/K
    # (note that the molecular weight is 100.3875, so don't be startled by the
    # coincidence in the first 3 significant figures )
    return (cpJmol * J / mol, cpJmol)
end

"""
    Tμ(h): (T,μ)

For working purposes, I need a very crude measure of the temperature T and the
melt fraction μ as a function of specific enthalpy h (or specific internal energy -
same thing, since pressure is basically zero).

I've done a bit of looking into this and the literature is on the one hand dense and
on the other hand hard to generalize from for the specific conditions of interest.

Melt temperature ("liquidus") is 1965C for pure enstatite, but a bit lower for
chondrules. Lofgren, LSPC XXI (old reference), gives 1475-1600C for melt.
Heat of fusion hard to find. Basalt heat of fusion is about 500J/g, equivalent to
about 50kJ/mol (mol of enstatite =100g about). Heating from 0K to 2000K takes,
within a factor of 2x, about 250kJ/mol. So, the heat of fusion isn't huge but it's
not negligible either.

Heat of vaporization IS huge, being about 2812 kJ/mol (minus heat of fusion and heat to get
to whatever temp the thing is). That means that evaporative
cooling will be substantial.

Here, assume a constant cp for solid and liquid of about 1.2e7 erg/g, and assume
heat of fusion of 500J/g = 5.0e9 erg/g.
"""
function Tμ(h)
    if h ≤ h1 # solid
        return (h/cp0,0.0)
    elseif h ≤ h2 # partially melted
        return (Tℓ,(h-h1)/Δh)
    elseif h > h2 # fully melted
        return (Tℓ + (h-h2)/cp0,1.0)
    else # should not get here
        return (NaN,NaN)
    end
end

function hμ(T) # inverse of the above - only needed for initialization
    if T ≤ Tℓ
        return (T*cp0,0.0)
    elseif T > Tℓ
        return (h2 + (T-Tℓ)*cp0,1.0)
    else
        return (NaN,NaN)
    end
end

function T(sie)
    #=
    Not using this yet. Plan is to eventually use a caloric equation of state
    that's a bit more sophisticated than constant cp and E = cp*T.
    =#
    cps = 1.0
    cpℓ = 0.5
    hℓ = 10.0
    he = 30.0
    Tℓ = 1000.0
    Te = 2000.0
    #
    sie0 = cps * Tℓ
    sie1 = sie0 + hℓ
    sie2 = sie1 + cpℓ*(Te-Tℓ)
    sie3 = sie2 + he
    if (sie < sie0)
        return sie / cps
    elseif (sie < sie1)
        return Tℓ
    elseif (sie < sie2)
        return Tℓ + (sie-sie1)/cpℓ
    elseif (sie < sie3)
        return Te
    else
        return NaN # needs some thought...
    end
end

mutable struct speck # physical state of a bulked (unresolved) grain
    D::f64 # Diameter
    m::f64 # mass
    T::f64 # Temperature
    Tmax::f64 # max historical temperature
    E::f64 # (internal) Energy
    μ::f64 # melt fraction
    μmax::f64 # max historical melt fraction
    dead::Bool
end

function speck(D,T)
    m = ρₛ * (π/6) * D^3
    (e,μ) = hμ(T)
    E = e * m
    speck(D,m,T,T,E,μ,μ,false)
end

function heat!(s::speck,ΔE) # Q is a heating *rate*; total energy delivered is Q * Δt
    s.E += ΔE
    (T,μ) = Tμ(s.E/s.m)
    s.T = T
    s.Tmax = max(T,s.Tmax)
    s.μ = μ
    s.μmax = max(μ,s.μmax)
    return nothing
end

heat!(s::speck,Q,Δt) = heat!(s,Q*Δt) # Q is a heating *rate*; total energy delivered is Q * Δt

function coolrad!(s::speck,Δt;Tx=0.0)
    T0 = s.T
    Q = - ε * σSB * π * s.D^2 * (T0^4 - Tx^4)
    heat!(s,Q,Δt)
    T1 = s.T
    # @assert abs((T1-T0)/T0) < 0.02 # if temperature changes by more than 2%, time step too big.
    return nothing
end

#function evaporate!(s::speck,Δm) #Note: a *positive* Δm is a mass *loss*.
#    @assert (Δm < s.m) | (s.m < 1.70e-9) # This is the mass of a 10μm enstatite grain
function evaporate!(s::speck,Δt)
    Δm = evap(s.T) * Δt
    T0 = s.T
    if Δm > s.m
        s.dead = true
        s.m = s.E = s.T = s.D = 0.0
        return nothing
    else
        e1 = s.E / s.m # mean specific energy
        s.m -= Δm
        s.D = (s.m * (6/π) / ρₛ)^(1/3)
        if s.μ ≈ 1.0
            e = e1
        elseif (s.μ > 0.0) & (Δm ≤ s.m * s.μ)
            e = h2
        else
#            @warn "Evaporating solid grain."
            e = e1
        end
        ΔE = Δm * (ΔH0-e) # Heat of vaporiation is relative to STP, so cooling is a bit less effective if you're already hot, obvs
    end
    heat!(s,-ΔE)
    T1 = s.T
#    @assert abs((T1-T0)/T0) < 0.10 # if temperature changes by more than 10%, time step (or mass loss) too big.
end

#function heatadv(s::speck,)
#    Qdot = π * D * κ * ΔT * Nu
#end

mutable struct mote
    D::f64 # Diameter (plan is to update this due to evaporation...)
    m::f64 # mass (plan is to update this due to evaporation...)
    N::Integer # number of cells
    Δt::f64 # suggested time step
    m_::Array{f64,1} # mass interior
    E_::Array{f64,1} # internal energy
    T_::Array{f64,1} # Temperature
    Tmax_::Array{f64,1} # historical max temperature
    Cᵐ::Any #Diagonal{f64,Array{f64,1}} # heat capacity matrix
    Cⁱᵐ::Any# Diagonal{f64,Array{f64,1}} # inverse heat capacity matrix
    Fᵐ::Any#Tridiagonal{f64,Array{f64,1}} # LHS tridiagonal matrix
    #Bᵐ::Any#Tridiagonal{f64,Array{f64,1}} # RHS tridiagonal matrix
    Dᵐ::Any# differencing matrix (is this really needed? probably not...)
    #phase_::Array{f64,1} # phase: 0 = solid, 1 = liquid, 2 = gas
end
function mote(D,N=5,T0=0.0)
    # We could imagine a non-uniform grid, but even though using Crank-Nicholson,
    # time-step is set by von Neumann stability criterion (to ensure accuracy),
    # so there's no real advantage to the added complication.
    # All of the matrices for time-stepping are set here, and "carried" with the
    # struct, to make time-stepping itself easier.
    # NB: R_, M_ = vertex-centered; r_, m_ = cell=centered. (T_, E_ are cell-centered)
    m = ρₛ * (π/6)D^3
    R_ = collect(range(0.0,D/2,length=N+1))[2:end]
    ΔR = (D/2.0)/N
    ΔtVN = 0.5 * ΔR^2 / (κ/ρₛ/cₚ) # von Neumann stability max Δt
    @assert ΔR ≈ R_[2]-R_[1]
    R0_ = [0.0, R_...] # used to construct r_
    r_ = sqrt.((R0_[2:end].^2 + R0_[1:end-1].^2)/2) # too clever by half?
    M_ = m * R_.^3 / (D/2)^3
    m_ = [M_[1],diff(M_)...]
    Cᵐ  = Diagonal(cₚ * m_) # heat capacity matrix
    Cⁱᵐ = Diagonal(1.0./(cₚ * m_)) # inverse heat capacity matrix
    T_ = T0 * ones(N)
    E_ = Cᵐ * T_
    # E_ = e_ .* m_
    # T_ = T.(e_)
    Δt = ΔtVN #* 2
    # constructing k-matrix
    kdiag_ = [R_[1:end-1].^2 ./ diff(r_) ..., 0.0]
    ksuper_ = -kdiag_[1:end-1]
    ksub_ = zeros(N-1) # 0.0 * ksuper_
    kᵐ = 4π * κ * Tridiagonal(ksub_,kdiag_,ksuper_)
    # constructing D-matrix
    Dᵐ = Tridiagonal(ones(N-1),-ones(N),zeros(N-1))
    #
    #α = 0.5
    #Aᵐ = I(N) -    α  * Δt * Dᵐ * kᵐ * Cⁱᵐ
    #Bᵐ = I(N) + (1-α) * Δt * Dᵐ * kᵐ * Cⁱᵐ
    Fᵐ = Dᵐ * kᵐ * Cⁱᵐ
    return mote(D,m,N,Δt,m_,E_,T_,0.0*T_,Cᵐ,Cⁱᵐ,Fᵐ,Dᵐ)
end

"""
    heat!(grain,Q;Δt,α): heats (or cools) grain
    -  Q is the heat power into the grain, integrated over the surface (erg/s)
    -  Δt is the time step (default: 1x von Neumann for grain); can be set arbitrarily
            large - if so, heat!() just iterates until specified Δt is reached. The
            actual Δt in the loop is either 1x the grain's von Neumann Δt, or slightly less.
    - α is weighting for implicit/explicit. α=0 is explicit; α=1 is fully implicit.
"""
function heat!(mot::mote,Q=0.0;Δt=mot.Δt,α=0.5)
    # Δt may be set larger than mote's internal Δt, but not less (i.e., that is overridden).
    # If more, then heat equation is iterated at mote's Δt until specified Δt is reached.
    Ni = ceil(Δt/mot.Δt) # number of iterations
    Δtℓ = Δt / Ni # "loop" Δt
    #E_ = mote.E_;
    Fᵐ = mot.Fᵐ #;  Bᵐ = mot.Bᵐ
    N = mot.N;    Cⁱᵐ = mot.Cⁱᵐ;  Dᵐ = mot.Dᵐ
    #
    Qdot_ = [zeros(N-1)..., -Q]
    P_ = Δtℓ * Dᵐ * Qdot_  # inhomogeneous source vector on RHS
    #= could also use "diff()", in which case:
    Qdot_ = [zeros(N)..., -Q]
    P_ = -Δt * diff(Qdot_) =#
    Aᵐ = I(N) -   α   * Δtℓ * Fᵐ
    Bᵐ = I(N) + (1-α) * Δtℓ * Fᵐ
    E0 = sum(mot.E_)
    for i in 1:Ni
        E_ = mot.E_
        mot.E_ = Aᵐ \ (Bᵐ * E_ + P_)
        mot.T_ = Cⁱᵐ * mot.E_
        mot.Tmax_ = max.(mot.Tmax_,mot.T_)
    end
    E1 = sum(mot.E_)
    ΔE = E1-E0;  ΔQ = Q*Δt
    @assert ΔE ≈ ΔQ
    return nothing #ΔE ≈ ΔQ #nothing
end

function cool!(mot::mote;Δt=mot.Δt,α=0.5)
    Q = -σSB * π * mot.D^2 * mot.T_[end]^4
    heat!(mot,Q;Δt,α)
end

end#module

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
module TestShow
using Test
using Dates
using CategoricalArrays

@testset "show() for CategoricalPool{String} and its values" begin
    pool = CategoricalPool(["c", "b", "a"])

    opool = CategoricalPool(["c", "b", "a"], ["a", "b", "c"], true)

    nv1 = CategoricalValue(1, pool)
    nv2 = CategoricalValue(2, pool)
    nv3 = CategoricalValue(3, pool)

    ov1 = CategoricalValue(1, opool)
    ov2 = CategoricalValue(2, opool)
    ov3 = CategoricalValue(3, opool)

    @test sprint(show, pool) == "$CategoricalPool{String,UInt32}([\"c\",\"b\",\"a\"])"
    @test sprint(show, opool) == "$CategoricalPool{String,UInt32}([\"a\",\"b\",\"c\"]) with ordered levels"

    @test sprint(show, nv1) == "$CategoricalValue{String,UInt32} \"c\""
    @test sprint(show, nv2) == "$CategoricalValue{String,UInt32} \"b\""
    @test sprint(show, nv3) == "$CategoricalValue{String,UInt32} \"a\""

    @test sprint(show, ov1) == "$CategoricalValue{String,UInt32} \"c\" (3/3)"
    @test sprint(show, ov2) == "$CategoricalValue{String,UInt32} \"b\" (2/3)"
    @test sprint(show, ov3) == "$CategoricalValue{String,UInt32} \"a\" (1/3)"

    @test sprint(show, nv1, context=:typeinfo=>typeof(nv1)) == "\"c\""
    @test sprint(show, nv2, context=:typeinfo=>typeof(nv2)) == "\"b\""
    @test sprint(show, nv3, context=:typeinfo=>typeof(nv3)) == "\"a\""

    @test sprint(show, ov1, context=:typeinfo=>typeof(ov1)) == "\"c\""
    @test sprint(show, ov2, context=:typeinfo=>typeof(ov2)) == "\"b\""
    @test sprint(show, ov3, context=:typeinfo=>typeof(ov3)) == "\"a\""

    @test sprint(print, nv1) == sprint(print, ov1) == "c"
    @test sprint(print, nv2) == sprint(print, ov2) == "b"
    @test sprint(print, nv3) == sprint(print, ov3) == "a"

    @test string(nv1) == string(ov1) == "c"
    @test string(nv2) == string(ov2) == "b"
    @test string(nv3) == string(ov3) == "a"

    @test String(nv1) == String(ov1) == "c"
    @test String(nv2) == String(ov2) == "b"
    @test String(nv3) == String(ov3) == "a"

    @test repr(nv1) == repr(ov1) == "\"c\""
    @test repr(nv2) == repr(ov2) == "\"b\""
    @test repr(nv3) == repr(ov3) == "\"a\""

    b = IOBuffer()
    @test write(b, nv1) == 1
    @test String(take!(b)) == "c"
    @test write(b, nv2) == 1
    @test String(take!(b)) == "b"
    @test write(b, nv3) == 1
    @test String(take!(b)) == "a"
    @test write(b, ov1) == 1
    @test String(take!(b)) == "c"
    @test write(b, ov2) == 1
    @test String(take!(b)) == "b"
    @test write(b, ov3) == 1
    @test String(take!(b)) == "a"
end

@testset "show() for CategoricalPool{Date} and its values" begin
    pool = CategoricalPool([Date(1999, 12), Date(1991, 8), Date(1993, 10)])

    opool = CategoricalPool([Date(1999, 12), Date(1991, 8), Date(1993, 10)],
                            [Date(1991, 8), Date(1993, 10), Date(1999, 12)], true)

    nv1 = CategoricalValue(1, pool)
    nv2 = CategoricalValue(2, pool)
    nv3 = CategoricalValue(3, pool)

    ov1 = CategoricalValue(1, opool)
    ov2 = CategoricalValue(2, opool)
    ov3 = CategoricalValue(3, opool)

    @test sprint(show, pool) == "$CategoricalPool{Dates.Date,UInt32}([1999-12-01,1991-08-01,1993-10-01])"
    @test sprint(show, opool) == "$CategoricalPool{Dates.Date,UInt32}([1991-08-01,1993-10-01,1999-12-01]) with ordered levels"

    @test sprint(show, nv1) == "$CategoricalValue{Dates.Date,UInt32} 1999-12-01"
    @test sprint(show, nv2) == "$CategoricalValue{Dates.Date,UInt32} 1991-08-01"
    @test sprint(show, nv3) == "$CategoricalValue{Dates.Date,UInt32} 1993-10-01"

    @test sprint(show, ov1) == "$CategoricalValue{Dates.Date,UInt32} 1999-12-01 (3/3)"
    @test sprint(show, ov2) == "$CategoricalValue{Dates.Date,UInt32} 1991-08-01 (1/3)"
    @test sprint(show, ov3) == "$CategoricalValue{Dates.Date,UInt32} 1993-10-01 (2/3)"

    @test sprint(show, nv1, context=:typeinfo=>typeof(nv1)) == "1999-12-01"
    @test sprint(show, nv2, context=:typeinfo=>typeof(nv2)) == "1991-08-01"
    @test sprint(show, nv3, context=:typeinfo=>typeof(nv3)) == "1993-10-01"

    @test sprint(show, ov1, context=:typeinfo=>typeof(ov1)) == "1999-12-01"
    @test sprint(show, ov2, context=:typeinfo=>typeof(ov2)) == "1991-08-01"
    @test sprint(show, ov3, context=:typeinfo=>typeof(ov3)) == "1993-10-01"

    @test sprint(print, nv1) == sprint(print, ov1) == "1999-12-01"
    @test sprint(print, nv2) == sprint(print, ov2) == "1991-08-01"
    @test sprint(print, nv3) == sprint(print, ov3) == "1993-10-01"

    @test string(nv1) == string(ov1) == "1999-12-01"
    @test string(nv2) == string(ov2) == "1991-08-01"
    @test string(nv3) == string(ov3) == "1993-10-01"

    @test repr(nv1) == repr(ov1) == "1999-12-01"
    @test repr(nv2) == repr(ov2) == "1991-08-01"
    @test repr(nv3) == repr(ov3) == "1993-10-01"
end

using JSON
@testset "JSON.lower" for pool in (CategoricalPool(["a"]),
                                   CategoricalPool([1]),
                                   CategoricalPool([1.0]))
    v = CategoricalValue(1, pool)
    @test JSON.lower(v) == JSON.lower(get(v))
    @test typeof(JSON.lower(v)) == typeof(JSON.lower(get(v)))
end

end


text:
module IFGF

const PROJECT_ROOT =  pkgdir(IFGF)

using LinearAlgebra
using StaticArrays
using TimerOutputs
using Printf
using FFTW

import AbstractTrees
import WavePropBase
import LinearAlgebra: mul!

import WavePropBase:
    HyperRectangle,
    ClusterTree,
    UniformCartesianMesh,
    center,
    CardinalitySplitter,
    return_type,
    assert_concrete_type,
    partition_by_depth,
    ambient_dimension,
    container,
    radius,
    distance,
    isleaf,
    isroot,
    svector,
    children,
    parent,
    elements,
    low_corner,
    high_corner,
    ElementIterator,
    root_elements,
    index_range,
    coords,
    loc2glob

include("utils.jl")
include("targettree.jl")
include("sourcetree.jl")
include("chebinterp.jl")
include("ifgfoperator.jl")

export
    # re-export WavePropBase
    WavePropBase,
    # types
    IFGFOp,
    ClusterTree,
    # methods
    assemble_ifgf,
    # macros
    @hprofile
end

max_stars_repo_name: WaveProp/IFGF

max_issues_repo_name: WaveProp/IFGF

max_forks_repo_name: WaveProp/IFGF


text:
@enum(MatchingStrategy::Int,COMPLETE=1,COMPLETE2x2=2,CONSTRAINTS=3)

@kwdef mutable struct Options <: AbstractOptions
    pardiso_matching_strategy::MatchingStrategy = COMPLETE2x2
    pardiso_max_inner_refinement_steps::Int = 1
    pardiso_msglvl::Int = 0
    pardiso_order::Int = 2
    pardiso_log_level::String = ""
end

mutable struct Solver <: AbstractLinearSolver
    pt::Vector{Int}
    iparm::Vector{Int32}
    dparm::Vector{Float64}
    perm::Vector{Int32}
    msglvl::Ref{Int32}
    err::Ref{Int32}
    csc::SparseMatrixCSC{Float64,Int32}
    w::Vector{Float64}
    opt::Options
    logger::Logger
end

_pardisoinit(
    pt::Vector{Int},mtype::Ref{Cint},solver::Ref{Cint},iparm::Vector{Cint},
    dparm::Vector{Cdouble},err::Ref{Cint}) =
        ccall(
            (:pardisoinit,libpardiso),
            Cvoid,
            (Ptr{Int},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},Ptr{Cint}),
            pt,mtype,solver,iparm,dparm,err)
_pardiso(
    pt::Vector{Int},maxfct::Ref{Cint},mnum::Ref{Cint},mtype::Ref{Cint},
    phase::Ref{Cint},n::Ref{Cint},a::Vector{Cdouble},ia::Vector{Cint},ja::Vector{Cint},
    perm::Vector{Cint},nrhs::Ref{Cint},iparm::Vector{Cint},msglvl::Ref{Cint},
    b::StrideOneVector{Cdouble},x::StrideOneVector{Cdouble},err::Ref{Cint},dparm::Vector{Cdouble}) =
        ccall(
            (:pardiso,libpardiso),
            Cvoid,
            (Ptr{Int},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},Ptr{Cint},Ptr{Cint},Ptr{Cint},
             Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cint},Ptr{Cdouble}),
            pt,maxfct,mnum,mtype,phase,n,a,ia,ja,perm,nrhs,iparm,msglvl,b,x,err,dparm)

function Solver(csc::SparseMatrixCSC{Float64,Int32};
                opt=Options(),logger=Logger(),
                option_dict::Dict{Symbol,Any}=Dict{Symbol,Any}(),
                kwargs...)
    !isempty(kwargs) && (for (key,val) in kwargs; option_dict[key]=val; end)
    set_options!(opt,option_dict)

    w   = Vector{Float64}(undef,csc.n)

    pt = Vector{Int}(undef,64)
    iparm = Vector{Int32}(undef,64)
    dparm = Vector{Float64}(undef,64)
    perm = Vector{Int32}(undef,csc.n)
    msglvl = Ref{Int32}(0)
    err = Ref{Int32}(0)

    pt.=0
    iparm[1]=0

    _pardisoinit(pt,Ref{Int32}(-2),Ref{Int32}(0),iparm,dparm,err)
    err.x < 0  && throw(SymbolicException())

    iparm[1]=1
    iparm[2]=opt.pardiso_order # METIS
    iparm[3]=haskey(ENV,"OMP_NUM_THREADS") ? parse(Int32,ENV["OMP_NUM_THREADS"]) : 1
    iparm[6]=1
    iparm[8]=opt.pardiso_max_inner_refinement_steps
    iparm[10]=12 # pivot perturbation
    iparm[11]=0 # disable scaling
    iparm[13]=1 # matching strategy
    iparm[21]=3 # bunch-kaufman pivotin
    iparm[24]=1 # parallel factorization
    iparm[25]=1 # parallel solv

    _pardiso(pt,Ref{Int32}(1),Ref{Int32}(1),Ref{Int32}(-2),Ref{Int32}(11),
             Ref{Int32}(csc.n),csc.nzval,csc.colptr,csc.rowval,perm,
             Ref{Int32}(1),iparm,msglvl,Float64[],Float64[],err,dparm)
    err.x < 0  && throw(SymbolicException())

    M = Solver(pt,iparm,dparm,perm,msglvl,err,csc,w,opt,logger)

    finalizer(finalize,M)

    return M
end
function factorize!(M::Solver)
    _pardiso(M.pt,Ref{Int32}(1),Ref{Int32}(1),Ref{Int32}(-2),Ref{Int32}(22),
             Ref{Int32}(M.csc.n),M.csc.nzval,M.csc.colptr,M.csc.rowval,M.perm,
             Ref{Int32}(1),M.iparm,M.msglvl,Float64[],Float64[],M.err,M.dparm)
    M.err.x < 0  && throw(FactorizationException())
    M.iparm[14] == 0 && return M
    _pardiso(M.pt,Ref{Int32}(1),Ref{Int32}(1),Ref{Int32}(-2),Ref{Int32}(12),
             Ref{Int32}(M.csc.n),M.csc.nzval,M.csc.colptr,M.csc.rowval,M.perm,
             Ref{Int32}(1),M.iparm,M.msglvl,Float64[],Float64[],M.err,M.dparm)
    M.err.x < 0  && throw(FactorizationException())
    return M
end
function solve!(M::Solver,rhs::StrideOneVector{Float64})
    _pardiso(M.pt,Ref{Int32}(1),Ref{Int32}(1),Ref{Int32}(-2),Ref{Int32}(33),
             Ref{Int32}(M.csc.n),M.csc.nzval,M.csc.colptr,M.csc.rowval,M.perm,
             Ref{Int32}(1),M.iparm,M.msglvl,rhs,M.w,M.err,M.dparm)
    M.err.x < 0  && throw(SolveException())
    return rhs
end
function finalize(M::Solver)
    _pardiso(M.pt,Ref{Int32}(1),Ref{Int32}(1),Ref{Int32}(-2),Ref{Int32}(-1),
             Ref{Int32}(M.csc.n),M.csc.nzval,M.csc.colptr,M.csc.rowval,M.perm,
             Ref{Int32}(1),M.iparm,M.msglvl,Float64[],M.w,M.err,M.dparm)
end
is_inertia(::Solver)=true
function inertia(M::Solver)
    pos = M.iparm[22]
    neg = M.iparm[23]
    return (pos,M.csc.n-pos-neg,neg)
end
function improve!(M::Solver)
    @debug(M.logger,"improve quality failed.")
    return false
end

introduce(::Solver)="pardiso"


text:
### A Pluto.jl notebook ###
# v0.17.4

using Markdown
using InteractiveUtils

# ╔═╡ 303488a5-8102-4077-8cad-12655d8c5777
begin
	using ScikitLearn
	using Statistics
	using StatsBase: sample
	using DataFrames
	using CSV
	using Bootstrap
end

# ╔═╡ 6790c593-f5d3-4e45-be83-c974e581bbf2
md"""
# Bootstrap Notes pt. 2

This document will work through Jason Brownlee's article [A Gentle Introduction to the Bootstrap Method](https://machinelearningmastery.com/a-gentle-introduction-to-the-bootstrap-method/). This is my second document on this method, so I'll try to be more concise. I'll be translating the python in the article to Julia.

Further down in this document, I'll work on another of Brownlee's articles, [How to Calculate Bootstrap Confidence Intervals For Machine Learning Results in Python
](https://machinelearningmastery.com/calculate-bootstrap-confidence-intervals-machine-learning-results-python/).

Date: 2022 April 23
"""

# ╔═╡ 8b12df14-0a1f-4fc4-8677-02f26f16e1d6
md"""
---
The bootstrap method is a **resampling technique** used to estimate statistics on a population. It can be used to estimate **summary statistics** on a population such as the mean or standard deviation. It is used to estimate the skill of a model when making predictions on data not in the training set.

This method provides us confidence intervals which is something cross-validation does not.

This article will teach us:

 - The bootstrap method involves iteratively resampling a dataset with replacement.
 - The bootstrap method has two parameters, the size of the sample and the number of repeats.
 - Scikit-learn provides a resampling funciton. (Note, Julia provides a wrapper for ScikitLearn)
"""

# ╔═╡ 1023b0ba-fc8e-4d49-b795-ed3d75004b69
md"""
## Bootstrap Method

First, get one sample from the data set. This is done by choosing a random set of observations to add to the sample. I.e. this is done *without* replacement.

Then, we draw bootstrap samples from this original sample, but this time *with* replacement. We calculate the statistic on each sample, then calculate the mean of those statistics.

How would we use this method in the case of statistical inference?

For each bootstrap sample:
1. Fit a model on this sample
2. Estimate the skill on the model on the out-of-bag sample

Then calculate the mean of the sample of model skill estimates.

*Note*: The model is evaluated against the samples that weren't selected, *not* against the test dataset, which we can't use until we want to fully evluate the model. The data in the original sample but not the bootstrap sample are called out-of-bag (oob).

*Note*: A useful feature of the bootstrap method is that the resulting sample often forms a normal distribution.

"""

# ╔═╡ 66d8365c-cb3a-4a85-8da3-00eab0e1d1eb
md"""
## Example

In this example, we have a dataset with 6 elements. We will be using the entire dataset to create our EDF, meaning we can sample from the whole thing to create our bootstrap samples.

"""

# ╔═╡ 1a0f56fa-94a9-420f-b318-b6a0f197b140
dataset = 0.1:0.1:0.6

# ╔═╡ 8b712ac8-37fe-4c0e-98c1-495f632ed305
begin
	sampleset = [0.2, 0.1, 0.2, 0.6] # we can have repeats
	oob    = [0.3, 0.4, 0.5] # everything not in sample
end

# ╔═╡ ef6a30bf-2f90-486d-879d-a1a089a89c84

# ╔═╡ 68ad9943-a5a8-492f-8c5b-93bc7b49d36f
md"""
Note that above, we have a single sample (that we made by hand instead of actually using random sampling). This would be sufficient for a single iteration of the bootstrap procedure.

The rest of this iteration would look as follows:
1. Fit the model on `sample`
2. Evaluate the model on `oob`

## ScikitLearn

We can use ScikitLearn's `resample` method to generate the sample for us, but I'm not quite sure how to use the interface so I'll leave this one here.
"""

# ╔═╡ 4006c742-dc11-403d-b921-459b64ddd3c2
md"""
# Bootstrap Confidence Intervals

Now we'll get to the second article, as promised.

Step 1 is to do our bootstrap sampling:

```julia
B = 1000
stats = [calculate_stat(sample(data,N,replace=True)) for _ in 1:B]
```

We can get an **Empirical Confidence Interval** by simply sorting this list and taking the values at the percentile we are interested in. For intance if we want a 95% confidence interval, we take our interval to be the range starting at the 2.5 percentil and ending at the 97.5 percentile.

```julia
alpha=0.95
ordered = sort(stats)
lower = quantile(ordered, (1-alpha)/2)
upper = quantile9ordered, (alpha+((1-alpha)/2))
```
"""

# ╔═╡ e74d5d8f-7d69-4e32-a036-c711194d216e
md"""
## Calculate Classification Accuracy Confidence Interval

We'll read our dataset into a dataframe first.
"""

# ╔═╡ 54890bac-8050-4141-8e7e-2fb603bfe3b8
data = CSV.read("./datasets/pima-indians-diabetes.data.csv", DataFrame;
				header=false)

# ╔═╡ b8ccebb7-b303-4315-94fb-8ff253a90271
begin
	n_iterations = 1
	n_size = trunc(Int64, nrow(data)*0.5)
	stats = Array{Int64, n_iterations}
end

# ╔═╡ bde70fbd-d415-4d6b-b1d0-6eaba93012b9
#for i in 1:1
#	train = bootstrap(data, )
#	train = sample(std, data, BasicSampling(n_size))
#	#test = setdiff(data, train)
#end

# ╔═╡ 0148021c-b8a8-4c47-bcbe-3eae9ef6e75d
bootstrap(std, data[!, "Column6"], BasicSampling(n_iterations))

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Bootstrap = "e28b5b4c-05e8-5b66-bc03-6f0c0a0a06e0"
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
ScikitLearn = "3646fa90-6ef7-5e7e-9f22-8aca16db6324"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"

[compat]
Bootstrap = "~2.3.3"
CSV = "~0.10.4"
DataFrames = "~1.3.3"
ScikitLearn = "~0.6.4"
StatsBase = "~0.33.16"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Bootstrap]]
deps = ["DataFrames", "Distributions", "Random", "Statistics", "StatsBase", "StatsModels"]
git-tree-sha1 = "fb4bdaa6c24fea362cc54965c3f645f7cfac659f"
uuid = "e28b5b4c-05e8-5b66-bc03-6f0c0a0a06e0"
version = "2.3.3"

[[deps.CSV]]
deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
git-tree-sha1 = "873fb188a4b9d76549b81465b1f75c82aaf59238"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.10.4"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.14.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "b153278a25dd42c65abbf4e62344f9d22e59191b"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.43.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Conda]]
deps = ["Downloads", "JSON", "VersionParsing"]
git-tree-sha1 = "6e47d11ea2776bc5627421d59cdcc1296c058071"
uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
version = "1.7.0"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "6c19003824cbebd804a51211fd3bbd81bf1ecad5"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.3.3"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "5a4168170ede913a2cd679e53c2123cb4b889795"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.53"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FilePathsBase]]
deps = ["Compat", "Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "129b104185df66e408edd6625d480b7f9e9823a0"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.18"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "246621d23d1f43e3b9c368bf3b72b2331a27c286"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.13.2"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "61feba885fac3a407465726d0c330b3055df897f"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.1.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.3"

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.IterTools]]
git-tree-sha1 = "fa6287a4469f5e048d763df38279ee729fbd44e5"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.4.0"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "a970d55c2ad8084ca317a4658ba6ce99b7523571"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.12"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "e8185b83b9fc56eb6456200e873ce598ebc7f262"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.7"

[[deps.Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.3"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "3b429f37de37f1fc603cc1de4a799dc7fbe4c0b6"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.3.0"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "28ef6c7ce353f0b35d0df0d5930e0d072c1f5b9b"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.1"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.3.0"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "dfb54c4e414caa595a1f2ed759b160f5a3ddcba5"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.3.1"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.PyCall]]
deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
git-tree-sha1 = "1fc929f47d7c151c839c5fc1375929766fb8edcc"
uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
version = "1.93.1"

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.ScikitLearn]]
deps = ["Compat", "Conda", "DataFrames", "Distributed", "IterTools", "LinearAlgebra", "MacroTools", "Parameters", "Printf", "PyCall", "Random", "ScikitLearnBase", "SparseArrays", "StatsBase", "VersionParsing"]
git-tree-sha1 = "ccb822ff4222fcf6ff43bbdbd7b80332690f168e"
uuid = "3646fa90-6ef7-5e7e-9f22-8aca16db6324"
version = "0.6.4"

[[deps.ScikitLearnBase]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "7877e55c1523a4b336b433da39c8e8c08d2f221f"
uuid = "6e75b9c4-186b-50bd-896f-2d2496a4843e"
version = "0.5.0"

[[deps.SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "6a2f7d70512d205ca8c7ee31bfa9f142fe74310c"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.12"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.ShiftedArrays]]
git-tree-sha1 = "22395afdcf37d6709a5a0766cc4a5ca52cb85ea0"
uuid = "1277b4bf-5013-50f5-be3d-901d8477a67a"
version = "1.0.0"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.4"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "c82aaa13b44ea00134f8c9c89819477bd3986ecd"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.3.0"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8977b17906b0a1cc74ab2e3a05faa16cf08a8291"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.16"

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "5950925ff997ed6fb3e985dcce8eb1ba42a0bbe7"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.18"

[[deps.StatsModels]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Printf", "REPL", "ShiftedArrays", "SparseArrays", "StatsBase", "StatsFuns", "Tables"]
git-tree-sha1 = "03c99c7ef267c8526953cafe3c4239656693b8ab"
uuid = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
version = "0.6.29"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
git-tree-sha1 = "5ce79ce186cc678bbb5c5681ca3379d1ddae11a1"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.7.0"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.VersionParsing]]
git-tree-sha1 = "58d6e80b4ee071f5efd07fda82cb9fbe17200868"
uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
version = "1.3.0"

[[deps.WeakRefStrings]]
deps = ["DataAPI", "InlineStrings", "Parsers"]
git-tree-sha1 = "b1be2855ed9ed8eac54e5caff2afcdb442d52c23"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.4.2"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═303488a5-8102-4077-8cad-12655d8c5777
# ╟─6790c593-f5d3-4e45-be83-c974e581bbf2
# ╟─8b12df14-0a1f-4fc4-8677-02f26f16e1d6
# ╟─1023b0ba-fc8e-4d49-b795-ed3d75004b69
# ╟─66d8365c-cb3a-4a85-8da3-00eab0e1d1eb
# ╠═1a0f56fa-94a9-420f-b318-b6a0f197b140
# ╠═8b712ac8-37fe-4c0e-98c1-495f632ed305
# ╠═ef6a30bf-2f90-486d-879d-a1a089a89c84
# ╟─68ad9943-a5a8-492f-8c5b-93bc7b49d36f
# ╟─4006c742-dc11-403d-b921-459b64ddd3c2
# ╠═e74d5d8f-7d69-4e32-a036-c711194d216e
# ╠═54890bac-8050-4141-8e7e-2fb603bfe3b8
# ╠═b8ccebb7-b303-4315-94fb-8ff253a90271
# ╠═bde70fbd-d415-4d6b-b1d0-6eaba93012b9
# ╠═0148021c-b8a8-4c47-bcbe-3eae9ef6e75d
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# Mixed precision monetary arithmetic

# Promote to larger type if precision same
function Base.promote_rule(
        ::Type{Monetary{A,B,D}},
        ::Type{Monetary{A,C,D}}) where {A,B,C,D}
    Monetary{A, promote_type(B, C), D}
end

# Promote to BigInt if precision different
function Base.promote_rule(
        ::Type{Monetary{A,B,D}},
        ::Type{Monetary{A,C,E}}) where {A,B,C,D,E}
    Monetary{A, BigInt, max(D, E)}
end

# Convert with same kind of currency
function Base.convert(::Type{Monetary{A,B,C}}, m::Monetary{A,D,E}) where {A,B,C,D,E}
    Monetary{A,B,C}(m.val)
end

Base.isless(m::Monetary{T}, n::Monetary{T}) where {T} = isless(promote(m, n)...)
+(m::Monetary{T}, n::Monetary{T}) where {T} = +(promote(m, n)...)
/(m::Monetary{T}, n::Monetary{T}) where {T} = /(promote(m, n)...)

for fns in DIVS
    for fn in fns
        @eval function Base.$(fn)(m::Monetary{T}, n::Monetary{T}) where T
            $(fn)(promote(m, n)...)
        end
    end
end


text:
module StackTraces

import Base: ==

VERSION >= v"0.5.0-dev+2420" && warn("StackTraces.jl is deprecated. Please use Base.stacktrace() and Base.catch_stacktrace() instead.")

export StackFrame, StackTrace
export stacktrace, catch_stacktrace, format_stacktrace, format_stackframe, show_stacktrace, ==

immutable StackFrame
    func::Symbol
    file::Symbol
    line::Int
    inlined_file::Symbol
    inlined_line::Int
    from_c::Bool
    pointer::Int64
end

typealias StackTrace Vector{StackFrame}

const UNKNOWN = StackFrame(Symbol("???"), Symbol("???"), -1, Symbol(""), -1, true, 0)

function ==(a::StackFrame, b::StackFrame)
    a.line == b.line && a.from_c == b.from_c && a.func == b.func && a.file == b.file
end

function hash(frame::StackFrame, h::UInt)
    h += 0xf4fbda67fe20ce88 % UInt
    h = hash(frame.line, h)
    h = hash(frame.file, h)
    h = hash(frame.func, h)
end

"""
Given a pointer, looks up stack frame context information.
"""
function lookup(pointer::Ptr{Void})
    frame_info = ccall(:jl_lookup_code_address, Any, (Ptr{Void}, Cint), pointer - 1, 0)
    return (length(frame_info) == 7) ? StackFrame(frame_info...) : UNKNOWN
end

lookup(pointer::UInt) = lookup(convert(Ptr{Void}, pointer))

"""
Returns a stack trace in the form of a vector of StackFrames. Each StackFrame contains a
function name, a file name, a line number, and a flag indicating whether it's a C function.
(By default `stacktrace` doesn't return C functions, but this can be enabled.)
"""
function stacktrace(trace::Vector{Ptr{Void}}, c_funcs::Bool=false)
    stack = map(lookup, trace)

    # Remove frames that come from C calls.
    if !c_funcs
        filter!(frame -> !frame.from_c, stack)
    end

    # Remove frame for this function (and any functions called by this function).
    remove_frames!(stack, :stacktrace)
end

stacktrace(c_funcs::Bool=false) = stacktrace(backtrace(), c_funcs)

"""
Returns the stack trace for the most recent error thrown, rather than the current context.
"""
catch_stacktrace(c_funcs::Bool=false) = stacktrace(catch_backtrace(), c_funcs)

"""
Takes a StackTrace (a vector of StackFrames) and a function name (a Symbol) and removes the
StackFrame specified by the function name from the StackTrace (also removing all functions
above the specified function). Primarily used to remove StackTraces functions from the Stack
prior to returning it.
"""
function remove_frames!(stack::StackTrace, name::Symbol)
    splice!(stack, 1:findlast(frame -> frame.func == name, stack))
    return stack
end

function remove_frames!(stack::StackTrace, names::Vector{Symbol})
    splice!(stack, 1:findlast(frame -> in(frame.func, names), stack))
    return stack
end

function format_stackframe(frame::StackFrame; full_path::Bool=false)
    file_info = "$(full_path ? frame.file : basename(string(frame.file))):$(frame.line)"

    if frame.inlined_file != Symbol("")
        inline_info = "[inlined code from $file_info] "
        file_info = string(
            full_path ? frame.inlined_file : basename(string(frame.inlined_file)),
            ":", frame.inlined_line
        )
    else
        inline_info = ""
    end

    "$inline_info$(frame.func != "" ? frame.func : "?") at $file_info"
end

function format_stacktrace(
    stack::StackTrace, separator::AbstractString, start::AbstractString="",
    finish::AbstractString=""; full_path::Bool=false
)
    if isempty(stack)
        return ""
    end

    string(
        start,
        join(map(f -> format_stackframe(f, full_path=full_path), stack), separator),
        finish
    )
end

function show_stacktrace(io::IO, stack::StackTrace; full_path::Bool=false)
    println(
        io, "StackTrace with $(length(stack)) StackFrames$(isempty(stack) ? "" : ":")",
        format_stacktrace(stack, "\n  ", "\n  "; full_path=full_path)
    )
end

show_stacktrace(; full_path::Bool=false) = show_stacktrace(STDOUT; full_path=full_path)

function show_stacktrace(io::IO; full_path::Bool=false)
    show_stacktrace(io, remove_frames!(stacktrace(), :show_stacktrace); full_path=full_path)
end

function show_stacktrace(stack::StackTrace; full_path::Bool=false)
    show_stacktrace(STDOUT, stack; full_path=full_path)
end

end


text:
import Laplacian
import ConjugateGradient
import p1

# Run question 1 with grid size N
function q1(N::Integer)

    # Boundary conditions as matrices
    Vb, B = p1.boundaryconditions1(N)
    r = p1.chargedensity1(N)

    # Reshape conditions into vectors
    ϕb = reshape(Vb, (N+1)^2)
    β  = reshape(B,  (N+1)^2)
    ρ  = reshape(r,  (N+1)^2)

    # Laplacian matrix
    Δ = Laplacian.lap(N)
    Δf, Δb = Laplacian.opsplit(Δ, β)
    Δe = Laplacian.deltaeff(Δf, β)

    # Effective charge density
    ρe = Laplacian.rhoeff(ρ, Δb, ϕb, β)

    # Grip it and rip it; calculate ϕfree, use it to find total ϕ
    ϕf = ConjugateGradient.conjgrad(Δe, -ρe)
    ϕ = β.*ϕb + (1.-β).*ϕf

    # Reshape ϕ to give the physical grid and return that value
    return reshape(ϕ, N+1, N+1), r, reshape(ρe, N+1, N+1)

end

# Return X, Y, and Z coordinates of ϕ gridi matrix, assuming 1×1 grid size
function xyzcoords(ϕ)

    N = size(ϕ)[1] - 1
    a = 1/N

    ii = [0:N]  .* ones(N+1)'
    jj = [0:N]' .* ones(N+1)

    xs = a .* ii
    ys = a .* jj
    zs = ϕ

    X = reshape(xs, (N+1)^2)
    Y = reshape(ys, (N+1)^2)
    Z = reshape(zs, (N+1)^2)

    return X, Y, Z

end

function writexyxcoords(prefix,ϕ)

    X, Y, Z = xyzcoords(ϕ)

    writedlm("$prefix-phi-X.dat", X)
    writedlm("$prefix-phi-Y.dat", Y)
    writedlm("$prefix-phi-Z.dat", Z)

end


text:
using Parameters, ForwardDiff

struct DSGE_Model{data <: AbstractArray, alfa <: Float64}
    "Data"
    data::data
    "Alpha"
    alfa::alfa
end

function (problem::DSGE_Model)(pars)
    @unpack alfa, data = problem
    @unpack bet, epsilon, theta, sig, s2, phi, phi_pi, phi_y, rho_v = pars
    log_p_bet(bet) = logpdf(prior_bet,bet)
    log_p_epsilon(epsilon) = logpdf(prior_epsilon,epsilon)
    log_p_theta(theta) = logpdf(prior_theta,theta)
    log_p_sig(sig) = logpdf(prior_sig,sig)
    log_p_s2(s2) = logpdf(prior_s2,s2)
    log_p_phi(phi) = logpdf(prior_phi,phi)
    log_p_phi_pi(phi_pi) = logpdf(prior_phi_pi,phi_pi)
    log_p_phi_y(phi_y) = logpdf(prior_phi_y,phi_y)
    log_p_rho_v(rho_v) = logpdf(prior_rho_v,rho_v)
    llh,dll = log_like_dsge([alfa,bet,epsilon,theta,sig,s2,phi,phi_pi,phi_y,rho_v],data)
    llh = llh + log_p_bet(bet) + log_p_epsilon(epsilon) + log_p_theta(theta) + log_p_sig(sig) + log_p_s2(s2) + log_p_phi(phi) + log_p_phi_pi(phi_pi) + log_p_phi_y(phi_y) + log_p_rho_v(rho_v)
    dll = dll + [0 ForwardDiff.derivative(log_p_bet,bet) ForwardDiff.derivative(log_p_epsilon,epsilon) ForwardDiff.derivative(log_p_theta,theta) ForwardDiff.derivative(log_p_sig,sig) ForwardDiff.derivative(log_p_s2,s2) ForwardDiff.derivative(log_p_phi,phi) ForwardDiff.derivative(log_p_phi_pi,phi_pi) ForwardDiff.derivative(log_p_phi_y,phi_y) ForwardDiff.derivative(log_p_rho_v,rho_v)]

    return llh,dll
end

function problem_transform(p::DSGE_Model)
    as((bet = as_unit_interval, epsilon=as_positive_real,theta = as_unit_interval, sig = as_positive_real, s2 = as_positive_real, phi = as_positive_real, phi_pi = as_positive_real, phi_y = as_positive_real, rho_v = as_unit_interval))
end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
function included_file_path(path::String)
    normpath(joinpath(@__DIR__, "..", path))
end

function included_file(path::String)
  read(included_file_path(path), String)
end


text:
module test_h3_api_unidirectional_edges

using Test
using H3.API # h3IndexesAreNeighbors getH3UnidirectionalEdge h3UnidirectionalEdgeIsValid getOriginH3IndexFromUnidirectionalEdge getDestinationH3IndexFromUnidirectionalEdge getH3IndexesFromUnidirectionalEdge getH3UnidirectionalEdgesFromHexagon getH3UnidirectionalEdgeBoundary ≈

sfGeo = GeoCoord(0.659966917655, -2.1364398519396)
sf = geoToH3(sfGeo, 9)
ring = hexRing(sf, 1)
@test ring == [0x089283080ddbffff, 0x089283080c37ffff, 0x089283080c27ffff, 0x089283080d53ffff, 0x089283080dcfffff, 0x089283080dc3ffff]
@test !h3IndexesAreNeighbors(sf, sf)

sf2 = ring[1]
edge = getH3UnidirectionalEdge(sf, sf2)
@test edge == 0x149283080dcbffff

@test h3UnidirectionalEdgeIsValid(edge)

@test sf == getOriginH3IndexFromUnidirectionalEdge(edge)

@test sf2 == getDestinationH3IndexFromUnidirectionalEdge(edge)

originDestination = getH3IndexesFromUnidirectionalEdge(edge)
@test originDestination == (0x089283080dcbffff, 0x089283080ddbffff)

edges = getH3UnidirectionalEdgesFromHexagon(sf)
@test edges == [0x119283080dcbffff, 0x129283080dcbffff, 0x139283080dcbffff, 0x149283080dcbffff, 0x159283080dcbffff, 0x169283080dcbffff]

edgeBoundary = getH3UnidirectionalEdgeBoundary(edges[1])
@test edgeBoundary ≈ [GeoCoord(0.6600077078134922, -2.1364818365456384), GeoCoord(0.6599783466210994, -2.136500544207272)]

end # module test_h3_api_unidirectional_edges


text:
"""

A structure containing the
"""
mutable struct SubSolvers{Q<:MOI.AbstractOptimizer, S<:MOI.AbstractOptimizer, T<:ExtensionType}
    relaxed_optimizer::Q
    upper_optimizer::S
    ext_typ::T
end
SubSolvers{Q,S,T}(; r::Q = Cbc.Optimizer(), u::S = Ipopt.Optimizer(), t::T = DefaultExt()) where {Q,S,T} = SubSolvers{Q,S,T}(r, u, t)
SubSolvers(; r::Q = Cbc.Optimizer(), u::S = Ipopt.Optimizer(), t::T = DefaultExt()) where {Q,S,T} = SubSolvers{Q,S,T}(r,u,t)

function _relaxed_optimizer(d::SubSolvers{Q,S,T}) where {Q <: MOI.AbstractOptimizer,
                                                         S <: MOI.AbstractOptimizer,
                                                         T <: ExtensionType}
    return d.relaxed_optimizer
end

function _upper_optimizer(d::SubSolvers{Q,S,T}) where {Q <: MOI.AbstractOptimizer,
                                                       S <: MOI.AbstractOptimizer,
                                                       T <: ExtensionType}
    return d.upper_optimizer
end

function _ext_type(d::SubSolvers{Q,S,T}) where {Q <: MOI.AbstractOptimizer,
                                                S <: MOI.AbstractOptimizer,
                                                T <: ExtensionType}
    return d.ext_typ
end

function isempty(d::SubSolvers{Q,S,T}) where {Q,S,T}
    MOI.is_empty(d.relaxed_optimizer) && MOI.is_empty(d.upper_optimizer)
end
function MOI.empty!(d::SubSolvers{Q,S,T}) where {Q,S,T}
    MOI.empty!(d.relaxed_optimizer)
    MOI.empty!(d.upper_optimizer)
    return
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@testset "EnvelopeMatrices" begin
    trials = 100

    for m = 1:100, h = 1:4
        n = m
        A = dropzeros!(sprand(UInt, m, n, 0.5))
        env = rowenvelope(A)
        @testset "h = $h, m = $m" begin
            for _ = 1:trials
                j = rand(1:n + 1)
                j′ = rand(j:n + 1)
                if length(A.rowval[A.colptr[j]:A.colptr[j′] - 1]) > 0
                    @test env[j, j′] == extrema(A.rowval[A.colptr[j]:A.colptr[j′] - 1])
                else
                    @test env[j, j′] == (m + 1, 0)
                end
            end
        end
    end
end


text:
# # Immersed layer caches

#md # ```@meta
#md # CurrentModule = ImmersedLayers
#md # ```

#=
This package uses precomputed caches to efficiently implement the immersed layer
operators. The starting point for a cache is the specification of the
(discretized) body shape and the grid. Let's use an example.
=#

# ## Setting up a cache

using ImmersedLayers
using Plots
using LinearAlgebra

#=
### Set up a grid
First, we will set up a grid for performing the operations. We use the `PhysicalGrid`
constructor of the `CartesianGrids.jl` package to create this.
=#
Δx = 0.01
Lx = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx)

#=
### Set the shape
Now let's set a shape to immerse into the grid. We will use a circle, but
there are a variety of other shapes available. Many of these are in the
`RigidBodyTools.jl` package. Note that we set the spacing between the points
on this shape equal to 1.4 times the grid spacing. This is not critical, but it
is generally best to set it to a value between 1 and 2.
=#
RadC = 1.0
Δs = 1.4*cellsize(g)
body = Circle(RadC,Δs)

#=
### Create the cache
After setting the grid and the surface to immerse, the next step for using the
immersed layer tools is to set up a *surface cache*. This allocates a set of
data structures, as well as the critical *regularization* and *interpolation*
operators that will get used.

There are a few choices to make when setting this up
* **What kind of data (scalar or vector) are we dealing with?**

We will demonstrate with scalar data, which means we use the [`SurfaceScalarCache`](@ref)
function. For vector data, use [`SurfaceVectorCache`](@ref).

* **What type of scaling (grid or index) do we wish to apply to the operators?**

Grid scaling, set with `scaling = GridScaling`, means that the various operators are scaled with the physical grid spacing
and/or surface point spacing so that they approximate the continuous operators. This
means that regularization and interpolation are transposes with respect to
inner products that incorporate these physical spacings, rather than the usual
linear algebra inner products. Also, differential operations on the grid are true
approximations of their continuous counterparts. This choice of scaling is usually the best, and
the [`dot`](@ref) operator is extended in this package to implement the physically-
scaled inner products.

On the other hand, `scaling = IndexScaling` does not scale these, but rather, uses
pure differencing for the grid differential operators, and regularization is
the simple matrix transpose of interpolation.

* **What discrete Dirac delta function (DDF) do we wish to use?**

This is specified with the `ddftype` keyword argument. The default is
`ddftype=CartesianGrids.Yang3` [^1]. However, there are
other choices, such as `CartesianGrids.Roma`, `CartesianGrids.Goza`, `CartesianGrids.Witchhat`,
`CartesianGrids.M3`, `CartesianGrids.M4prime`.

* **Are there other data or functions you need to pass along?**

In more complicated cases, you might want to pass along some physical
parameters, in which case you can supply these with the `phys_params`
keyword. Boundary conditions can be provided with the `bc` keyword,
and forcing functions with the `f_funcs` keyword. These are all optional,
of course. We will see examples later in [Problems and the system](@ref)
=#

cache = SurfaceScalarCache(body,g,scaling=GridScaling)

#=
## Plotting the immersed points
We can plot the immersed points with the `plot` function of the `Plots.jl`
package, using a special plot recipe. This accepts all of the
usual attribute keywords of the basic plot function.
=#
plot(cache,xlims=(-2,2),ylims=(-2,2))

#=
In this plotting demonstration, we have used the `points` function
to obtain the coordinates of the immersed points. Other useful
utilities are `normals` and `areas`, e.g.
=#
normals(cache)

#=
## Some basic utilities
We will see deeper uses of this cache in [Surface-grid operations](@ref).
However, for now we can learn how to get basic copies of the grid
and surface data. For example, a copy of the grid data, all initialized to zero:
=#
zeros_grid(cache)

#=
or similarly on the surface
=#
zeros_surface(cache)

#=
We also might need to initialize data on the grid to accept the curl of
a vector field. This is used in conjunction with [`surface_curl!`](@ref),
for example.
=#
zeros_gridcurl(cache)

#=
If we want them initialized to unity, then use
=#
ones_grid(cache)

# and
ones_surface(cache)

#=
To evaluate functions on the grid, it is useful to be able to
fill grid data with the x and y coordinates. For this, we use
=#
x_grid(cache)
#-
y_grid(cache)

#=
## Using norms and inner products
It is useful to compute norms and inner products on grid and surface data.
These tools are easily accessible, e.g., `dot(u,v,cache)` and `norm(u,cache)`,
and they respect the scaling associated with the cache. For example,
the following gives an approximation of the circumference of the circle:
=#
os = ones_surface(cache)
dot(os,os,cache)

#md # ## Cache types and constructors

#md # ```@docs
#md # BasicILMCache
#md # SurfaceScalarCache
#md # SurfaceVectorCache
#md # AbstractExtraILMCache
#md # ```

#md # ## Utilities for creating instances of data

#md # ```@docs
#md # similar_grid
#md # similar_gridgrad
#md # similar_gridcurl
#md # similar_surface
#md # zeros_grid
#md # zeros_gridgrad
#md # zeros_gridcurl
#md # zeros_surface
#md # ones_grid
#md # ones_gridgrad
#md # ones_gridcurl
#md # ones_surface
#md # x_grid
#md # y_grid
#md # x_gridcurl
#md # y_gridcurl
#md # ```

#md # ## Utilities for accessing surface information

#md # ```@docs
#md # areas(::BasicILMCache)
#md # normals(::BasicILMCache)
#md # points(::BasicILMCache)
#md # ```

#md # ## Inner products, norms, and integrals

#md # ```@docs
#md # dot(::GridData,::GridData,::BasicILMCache)
#md # dot(::PointData,::PointData,::BasicILMCache)
#md # norm(::GridData,::BasicILMCache)
#md # norm(::PointData,::BasicILMCache)
#md # integrate(::ScalarData{N},::BasicILMCache{N}) where {N}
#md # dot(::PointData,::PointData,::BasicILMCache,::Int)
#md # norm(::PointData,::BasicILMCache,::Int)
#md # integrate(::ScalarData{N},::BasicILMCache{N},::Int) where {N}
#md # ```

#md # ## Other cache utilities
#md #
#md # ```@docs
#md # view(::PointData,::BasicILMCache,::Int)
#md # RegularizationMatrix(::BasicILMCache,::PointData,::GridData)
#md # InterpolationMatrix(::BasicILMCache,::GridData,::PointData)
#md # ```

#md # [^1]: Yang, X., et al., (2009) "A smoothing technique for discrete delta functions with application to immersed boundary method in moving boundary simulations," J. Comput. Phys., 228, 7821--7836.


text:
# This file roughly corresponds to functions documented in the

"""
    Canvas.groups(co::Course; kwargs...) -> ????

"""
function groups(co::Course; kwargs...)
    json, page_data = paged_request("GET", "/api/v1$(Internals.pid(co))/groups"; kwargs...)
    return Group.(json), page_data
end


text: Dict{String, Any}("arrays-and-ints" => Any[Dict{String, Any}("value" => "1", "type" => "integer"), Any[Dict{String, Any}("value" => "Arrays are not integers.", "type" => "string")]])


text:
is_gurobi(m) = solver_name(m) == "Gurobi"

"""set the MIPGap parameter of the passed CPLEX model to the passed tolerance"""
function set_tolerance!(model, tolerance)
	if !is_gurobi(model)
		set_optimizer_attribute(model, "CPXPARAM_MIP_Tolerances_MIPGap", tolerance)
	else
		set_optimizer_attribute(model, "MIPGap", tolerance)
	end
end

function set_time!(model, time)
	if !is_gurobi(model)
		set_optimizer_attribute(model, "CPXPARAM_TimeLimit", time)
	else
		set_optimizer_attribute(model, "TimeLimit", time)
	end
end

function set_threads!(model, num_threads)
	if !is_gurobi(model)
		set_optimizer_attribute(model, "CPXPARAM_Threads", num_threads)
	else
		set_optimizer_attribute(model, "Threads", num_threads)
	end
end

"""accept a struct with a bitlist attribute and set the model to have the same
bitlist as its start value"""
function set_bitlist!(model, sol)
	set_start_value.(model[:x], convert.(Float64, sol.bitlist))
end

"""Run CPLEX optimization without printing to the console"""
function silent_optimize_slow!(m)
	tempout = stdout # save stream
	start_time, end_time = 0, 0
	try
		redirect_stdout(nothing) # redirect to null
		start_time = time()
		optimize!(m)
		end_time = time()
		GC.gc()
		redirect_stdout(tempout)
	catch e
		#restore stream if user interrupts process, or other error
		redirect_stdout(tempout)
		return e
	end
	end_time - start_time
end

function silent_optimize!(m)
	start_time = time()
	optimize!(m)
	end_time = time()
	GC.gc()
	end_time - start_time
end

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
# Program to find GCD using Euclidean Algorithm

function gcd(x,y)
	if x==0
		return y
	else
		return gcd(y%x,x)
	end
end

println("x = ")
x = parse(Int64,readline())
println("\ny = ")
y = parse(Int64,readline())

println("gcd of ($x , $y) = ",gcd(x,y))

max_stars_repo_name: cashewcodes/Basic-Julia-Projects

max_issues_repo_name: cashewcodes/Basic-Julia-Projects

max_forks_repo_name: cashewcodes/Basic-Julia-Projects