text:
import Codex
import Codex.TransformExport
import Codex.Questionnaires

using DataFrames
using Test

@testset "personality" begin
    path = joinpath(export_dir, "responses", "personality.csv")
    simple = Codex.TransformExport.read_csv(path, delim=';')
    scores = Codex.Questionnaires.personality2scores(simple)
    @test typeof(first(scores)[:A1]) == Int
end


text:
# FIXME: Remove once fma is implemented in julia
using ErrorFreeTransforms.fma

facts("Fused multiply-add") do
  context("Float32") do
    a = float32(0x1000002)
    b = 3.0f0
    c = -float32(0x300000a)
    @fact fma(a, b, c) => big(a)*big(b)+big(c)
  end

  context("Float64") do
    a = float64(0x10000000000001)
    b = 3.0e0
    c = -float64(0x30000000000004)
    @fact fma(a, b, c) => big(a)*big(b)+big(c)
  end
end

facts("Error-free fma") do
  a = ldexp(float64(0x10000000000001), -52)
  b = ldexp(float64(0x10000000000001), -51)
  c = ldexp(float64(0x10000000000001), -48)

  context("Standard") do
    x, y, z = err_fma(a, b, c)
    @fact x => fma(a, b, c)
    @fact big(x)+big(y)+big(z) => big(a)*big(b)+big(c)
  end

  context("Compressed error") do
    x, y, z = err_fma_compr(a, b, c)
    @fact x => fma(a, b, c)
    @fact big(x)+big(y)+big(z) => big(a)*big(b)+big(c)
    # assuming non-zero y, z
    @fact exponent(y)-exponent(z) => greater_than_or_equal(53)
  end

  context("Nearest error") do
    x, y, z = err_fma(a, b, c)
    @fact err_fma_nearest(a, b, c) => x, y+z
  end

  context("Approximated error") do
    x, y = err_fma_approx(a, b, c)
    @fact big(a)*big(b)+big(c) => roughly(big(x)+big(y), ldexp(14*abs(x),-53))
  end
end

ext: jl


text:
function iszero{T}(x::T)
    return x == zero(T)
end

function Maybe(T::Type)
    return Union(T, Nothing)
end

# Cycle-zip. Zip two or more arrays, cycling the short ones.
function cyclezip(xs::AbstractArray...)
    if any(map(isempty, xs))
        return {}
    end
    n = maximum([length(x) for x in xs])
    return takestrict(zip([cycle(x) for x in xs]...), n)
end


text:
#=

    plotting.jl

    Contains code to produce inference related plots.

    Author:     Alexander P. Browning
                ======================
                School of Mathematical Sciences
                Queensland University of Technology
                ======================
                ap.browning@icloud.com
                alexbrowning.me

=#

"""
    plot_posteriors(Θ)

Produce density plot of posteriors. Optionally provide names of variables (a vector of `String`s or `Symbol`s) and
`Matrix{Symbol}` to specify the plot porder and porder.

"""
function plot_posteriors(Θ::Matrix,W=ones(size(Θ,2));
        idx = trues(size(Θ,2)),
        param_names = ["p$i" for i = 1:size(Θ,1)],
        porder = nothing,
        lims = nothing,
        nbins = nothing,
        layout=nothing,
        kwargs...
    )

    n = length(param_names)
    param_names = Symbol.(param_names)
    porderv = porder === nothing ? param_names : filter(x -> x != :_,permutedims(porder)[:])

    # Convert to named tuple for plotting
    Θ2 = NamedTuple{Tuple(param_names)}([Θ[i,:] for i = 1:n])
    lims2 = lims === nothing ? nothing : NamedTuple{Tuple(param_names)}(lims)

    # Loop through porder and create plots
    plts = [stephist(Θ2[porderv[i]][idx],weights=W[idx],xlabel=porderv[i],
        xlim = lims === nothing ? :auto : lims2[porderv[i]],
        frange = 0.0,
        c = :grey,
        bins = (lims === nothing) || (nbins === nothing) ? :auto : range(lims2[porderv[i]]...,length=nbins),
        kwargs...
    ) for i = 1:length(porderv)]

    # Put together
    Plots.plot(plts...,
        layout=layout === nothing ? length(porderv) : layout,
        xtickfontsize=6,
        yticks=[],
        size=(700,400),
        legend=:none,
        widen=false
    )

end

plot_posteriors(P::Vector{Particle};kwargs...) = plot_posteriors(locations(P),weights(P);kwargs...)
plot_posteriors(C::Chains;burnin=0,skip=1,kwargs...) = plot_posteriors(Matrix(C[burnin+1:skip:end]);param_names=C.name_map.parameters,kwargs...)

function plot_chains(C::Chains;
        porder = nothing,
        lims = nothing
    )

    n = size(C,2) - 1
    param_names = C.name_map.parameters
    porderv = porder === nothing ? param_names : permutedims(porder)[:]

    # Convert to named tuple for plotting
    lims2 = lims === nothing ? nothing : NamedTuple{Tuple(param_names)}(lims)

    # Loop through porder and create plots
    plts = [plot(C[:,porderv[i],1],xlabel=porderv[i],
        ylim = lims === nothing ? :auto : lims2[porderv[i]]) for i = 1:n]

    # Put together
    Plots.plot(plts...,
        porder=porder===nothing ? length(param_names) : grid(size(porder)...),
        xtickfontsize=6,
        ytickfontsize=6,
        size=(700,400),
        legend=:none
    )

end

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
"""
    Hermite( nodes::NodesType, epsilon:Real, gamma:Real )

Hermite interpolation

"""
mutable struct Hermite <: InterpolationType

    nnodes     :: Int64
    epsilon    :: Float64
    gamma      :: Float64
    nodes      :: Array{Float64,1}
    colloc_mat :: Array{Float64,2}

    function Hermite( nodes_t::NodesType, epsilon::Float64, gamma::Float64 )

        nodes  = nodes_t.xk
        nnodes = nodes_t.nx
        colloc_mat = evaluate_hermite(nodes, nnodes, epsilon, gamma)
        new( nnodes, epsilon, gamma, nodes, colloc_mat )

    end

end

"""
    Hermite( xe )

This function returns the matrix He of the values of the HermiteGF basis
functions. The computation is done via three term recurrence for Hermite
functions with an argument gamma*x and then appropriate exponential scaling.

More details can be found in Section 5.1 of the paper

*STABLE EVALUATION OF GAUSSIAN RADIAL BASIS
 FUNCTIONS USING HERMITE POLYNOMIALS*
 by Anna Yurova and Katharina Kormann.

"""
function (interp::Hermite)( x::Array{Float64,1} )

   nx      = interp.nnodes
   epsilon = interp.epsilon
   gamma   = interp.gamma

   evaluate_hermite(x, nx, epsilon, gamma) / interp.colloc_mat

end

"""

    evaluate_hermite(xk, n, epsilon, gamma)

This function returns the matrix He of the values of the HermiteGF basis
functions. The computation is done via three term recurrence for Hermite
functions with an argument gamma*x and then appropriate exponential scaling.

More details can be found in Section 5.1 of the paper

*STABLE EVALUATION OF GAUSSIAN RADIAL BASIS
 FUNCTIONS USING HERMITE POLYNOMIALS*
 by Anna Yurova and Katharina Kormann.

This function should be used both for computation of the
collocation and evaluation matrices.

"""
function evaluate_hermite(xk, n, epsilon, gamma)

  # Initialize the result matrix
  result = zeros(Int(length(xk)), n)

  # Write the values of the first two Hermite functions to initialize the
  # three-term recurrence
  result[:,1] = exp.(-0.5.*((gamma*xk).^2))
  result[:,2] = gamma*xk.*sqrt(2).*exp.(-0.5.*((gamma*xk).^2))

  # Three term recurrence for the Hermite functions with argument gamma*x
  for i = 3:n
    result[:,i] = sqrt(2.0/(i-1)) .* (gamma*xk) .* result[:,i-1] - sqrt((i-2)/(i-1)) .* result[:,i-2]
  end

  # Scaling the Hermite functions with the exponential factor
  # exp(-epsilon^2 x^2 + (gamma x)^2/2)

  for i=1:n
      result[:,i] = (pi^(1/4)) * result[:,i] .* exp.((xk.^2) .* (gamma*gamma*0.5 - epsilon^2))
  end

  result

end


text:
"""
    GeoMeanBridge{T, F, G, H}

The `GeometricMeanCone` is `SecondOrderCone` representable; see [1, p. 105].

The reformulation is best described in an example.

Consider the cone of dimension 4:
```math
t \\le \\sqrt[3]{x_1 x_2 x_3}
```
This can be rewritten as ``\\exists x_{21} \\ge 0`` such that:
```math
\\begin{align*}
  t & \\le x_{21},\\\\
  x_{21}^4 & \\le x_1 x_2 x_3 x_{21}.
\\end{align*}
```

Note that we need to create ``x_{21}`` and not use ``t^4`` directly as ``t`` is
allowed to be negative. Now, this is equivalent to:
```math
\\begin{align*}
  t & \\le x_{21}/\\sqrt{4},\\\\
  x_{21}^2 & \\le 2x_{11} x_{12},\\\\
  x_{11}^2 & \\le 2x_1 x_2, & x_{12}^2 & \\le 2x_3(x_{21}/\\sqrt{4}).
\\end{align*}
```

[1] Ben-Tal, Aharon, and Arkadi Nemirovski. *Lectures on modern convex
    optimization: analysis, algorithms, and engineering applications*. Society
    for Industrial and Applied Mathematics, 2001.
"""
struct GeoMeanBridge{T,F,G,H} <: AbstractBridge
    d::Int
    xij::Vector{MOI.VariableIndex}
    t_upper_bound_constraint::MOI.ConstraintIndex{F,MOI.LessThan{T}}
    rsoc_constraints::Vector{MOI.ConstraintIndex{G,MOI.RotatedSecondOrderCone}}
    x_nonnegative_constraint::Union{
        Nothing,
        MOI.ConstraintIndex{H,MOI.Nonnegatives},
    }
    # We use an internal constructor here so that users cannot pass
    # GeoMeanBridge(args..., nothing) can cause a `H not defined` error.
    function GeoMeanBridge{T,F,G,H}(
        d,
        xij,
        t_upper_bound_constraint,
        rsoc,
        x_nonnegative_constraint,
    ) where {T,F,G,H}
        return new{T,F,G,H}(
            d,
            xij,
            t_upper_bound_constraint,
            rsoc,
            x_nonnegative_constraint,
        )
    end
end

_ilog2(n, i = 0) = n <= (1 << i) ? i : _ilog2(n, i + 1)

function bridge_constraint(
    ::Type{GeoMeanBridge{T,F,G,H}},
    model,
    f::MOI.AbstractVectorFunction,
    s::MOI.GeometricMeanCone,
) where {T,F,G,H}
    f_scalars = MOI.Utilities.eachscalar(f)
    t = f_scalars[1]
    if MOI.dimension(s) == 2
        # The case {(t, x) ∈ R²: x₁ >= 0, t <= x₁}
        # The constraint t <= x₁
        t_upper_bound_constraint = MOI.Utilities.normalize_and_add_constraint(
            model,
            MOI.Utilities.operate!(-, T, t, f_scalars[2]),
            MOI.LessThan(zero(T)),
            allow_modify_function = true,
        )
        # The constraint x₁ >= 0
        x_nonnegative_constraint =
            MOI.add_constraint(model, f_scalars[2:2], MOI.Nonnegatives(1))
        return GeoMeanBridge{T,F,G,H}(
            MOI.dimension(s),
            MOI.VariableIndex[],
            t_upper_bound_constraint,
            MOI.ConstraintIndex{G,MOI.RotatedSecondOrderCone}[],
            x_nonnegative_constraint,
        )
    end
    SG = MOIU.scalar_type(G)
    d = MOI.dimension(s)
    n = d - 1
    l = _ilog2(n)
    N = 1 << l
    rsoc_constraints =
        Vector{MOI.ConstraintIndex{G,MOI.RotatedSecondOrderCone}}(undef, N - 1)
    xij = MOI.add_variables(model, N - 1)
    xl1 = xij[1]
    sN = one(T) / sqrt(N)
    _getx(i)::SG = i > n ? sN * xl1 : f_scalars[1+i]

    # With sqrt(2)^l*t - xl1, we should scale both the ConstraintPrimal and
    # ConstraintDual
    t_upper_bound_constraint = MOIU.normalize_and_add_constraint(
        model,
        MOIU.operate!(+, T, t, -sN * xl1),
        MOI.LessThan(zero(T)),
        allow_modify_function = true,
    )
    offset = offset_next = 0
    for i in 1:l
        num_lvars = 1 << (i - 1)
        offset_next = offset + num_lvars
        for j in 1:num_lvars
            if i == l
                a = _getx(2j - 1)
                b = _getx(2j)
            else
                a = convert(SG, xij[offset_next+2j-1])
                b = convert(SG, xij[offset_next+2j])
            end
            c = xij[offset+j]
            rsoc_constraints[offset+j] = MOI.add_constraint(
                model,
                MOIU.operate(vcat, T, a, b, c),
                MOI.RotatedSecondOrderCone(3),
            )
        end
        offset = offset_next
    end
    @assert offset == length(rsoc_constraints)
    return GeoMeanBridge{T,F,G,H}(
        d,
        xij,
        t_upper_bound_constraint,
        rsoc_constraints,
        nothing,
    )
end

function MOI.supports_constraint(
    ::Type{GeoMeanBridge{T}},
    ::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.GeometricMeanCone},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{<:GeoMeanBridge})
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(
    ::Type{<:GeoMeanBridge{T,F,G}},
) where {T,F,G}
    return Tuple{Type,Type}[
        (F, MOI.LessThan{T}),
        (G, MOI.RotatedSecondOrderCone),
        (G, MOI.Nonnegatives),
    ]
end

function concrete_bridge_type(
    ::Type{<:GeoMeanBridge{T}},
    H::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.GeometricMeanCone},
) where {T}
    S = MOIU.scalar_type(H)
    A = MOIU.promote_operation(*, T, T, MOI.VariableIndex)
    F = MOIU.promote_operation(+, T, S, A)
    G = MOIU.promote_operation(vcat, T, A, S, MOI.VariableIndex)
    return GeoMeanBridge{T,F,G,H}
end

# Attributes, Bridge acting as a model
MOI.get(b::GeoMeanBridge, ::MOI.NumberOfVariables)::Int64 = length(b.xij)

MOI.get(b::GeoMeanBridge, ::MOI.ListOfVariableIndices) = copy(b.xij)

function MOI.get(
    ::GeoMeanBridge{T,F},
    ::MOI.NumberOfConstraints{F,MOI.LessThan{T}},
)::Int64 where {T,F}
    return 1
end

function MOI.get(
    b::GeoMeanBridge{T,F,G},
    ::MOI.NumberOfConstraints{G,MOI.RotatedSecondOrderCone},
)::Int64 where {T,F,G}
    return length(b.rsoc_constraints)
end

function MOI.get(
    b::GeoMeanBridge{T,F,G},
    ::MOI.NumberOfConstraints{G,MOI.Nonnegatives},
)::Int64 where {T,F,G}
    return b.d > 2 ? 0 : 1
end

function MOI.get(
    b::GeoMeanBridge{T,F},
    ::MOI.ListOfConstraintIndices{F,MOI.LessThan{T}},
) where {T,F}
    return [b.t_upper_bound_constraint]
end

function MOI.get(
    b::GeoMeanBridge{T,F,G},
    ::MOI.ListOfConstraintIndices{G,MOI.RotatedSecondOrderCone},
) where {T,F,G}
    return copy(b.rsoc_constraints)
end

function MOI.get(
    b::GeoMeanBridge{T,F,G,H},
    ::MOI.ListOfConstraintIndices{H,MOI.Nonnegatives},
) where {T,F,G,H}
    if b.d > 2
        return MOI.ConstraintIndex{H,MOI.Nonnegatives}[]
    end
    return [b.x_nonnegative_constraint]
end

function MOI.delete(model::MOI.ModelLike, bridge::GeoMeanBridge)
    MOI.delete(model, bridge.xij)
    MOI.delete(model, bridge.t_upper_bound_constraint)
    MOI.delete(model, bridge.rsoc_constraints)
    if bridge.d == 2
        MOI.delete(model, bridge.x_nonnegative_constraint)
    end
    return
end

# Attributes, Bridge acting as a constraint
function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, bridge::GeoMeanBridge)
    return MOI.GeometricMeanCone(bridge.d)
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    bridge::GeoMeanBridge{T,F,G,H},
) where {T,F,G,H}
    d = bridge.d
    f_scalars = Vector{MOIU.scalar_type(H)}(undef, bridge.d)
    tub = MOI.get(model, attr, bridge.t_upper_bound_constraint)
    rhs = MOI.constant(
        MOI.get(model, MOI.ConstraintSet(), bridge.t_upper_bound_constraint),
    )
    tub = MOIU.operate(-, T, tub, rhs)
    if d == 2
        x = MOI.get(model, attr, bridge.x_nonnegative_constraint)
        f_scalars[2] = MOIU.eachscalar(x)[1]
        f_scalars[1] = MOIU.operate(+, T, tub, f_scalars[2])
    else
        t = MOIU.remove_variable(tub, bridge.xij[1])
        f_scalars[1] = t
        n = d - 1
        l = _ilog2(n)
        num_lvars = 1 << (l - 1)
        offset = num_lvars - 1
        for j in 1:num_lvars
            if 2j <= bridge.d
                func = MOI.get(model, attr, bridge.rsoc_constraints[offset+j])
                func_scalars = MOIU.eachscalar(func)
                # Numerical issues can arise when a VectorOfVariables function
                # is turned into a VectorAffineFunction, because the RSOC
                # constraints might return something like `0.999999999x` instead
                # of exactly `1.0x`. To counteract this, use `convert_approx`.
                f_scalars[2j] = MOI.Utilities.convert_approx(
                    MOI.Utilities.scalar_type(H),
                    func_scalars[1],
                )
                if 2j + 1 <= bridge.d
                    f_scalars[2j+1] = MOI.Utilities.convert_approx(
                        MOI.Utilities.scalar_type(H),
                        func_scalars[2],
                    )
                end
            end
        end
    end
    return MOIU.vectorize(f_scalars)
end

function _getconstrattr(model, attr, bridge::GeoMeanBridge{T}) where {T}
    output = Vector{T}(undef, bridge.d)
    output[1] = MOI.get(model, attr, bridge.t_upper_bound_constraint)
    if bridge.d == 2
        output[2] = MOI.get(model, attr, bridge.x_nonnegative_constraint)[1]
    else
        N = length(bridge.xij) + 1
        # div(N, 2) gets layer before original problem variables are involved
        offset = div(N, 2) - 1 # 1 + 2 + ... + n/4
        for i in 1:(bridge.d-1)
            j = ((i - 1) >> 1) + 1
            k = i - 2(j - 1)
            output[1+i] =
                MOI.get(model, attr, bridge.rsoc_constraints[offset+j])[k]
        end
    end
    return output
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimal,
    bridge::GeoMeanBridge,
)
    output = _getconstrattr(model, attr, bridge)
    N = length(bridge.xij) + 1
    # the constraint is t - x_l1/sqrt(2^l) ≤ 0, we need to add the value of x_l1
    if bridge.d == 2
        output[1] += MOI.get(model, attr, bridge.x_nonnegative_constraint)[1]
    else
        output[1] +=
            MOI.get(model, MOI.VariablePrimal(), bridge.xij[1]) / sqrt(N)
    end
    return output
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintDual,
    bridge::GeoMeanBridge,
)
    output = _getconstrattr(model, attr, bridge)
    if bridge.d == 2
        # The transformation is:
        # [t]                [1 -1] [t]
        # [x] in GeoMean <=> [0  1] [x] in R- x R+
        # Hence the dual transformation is
        # [ 1 0] [u]                 [u]
        # [-1 1] [v] in GeoMean* <=> [v] in R- x R+
        output[2] -= MOI.get(model, attr, bridge.t_upper_bound_constraint)
    end
    # Otherwise, the transformation is:
    # [t]                                  [t]
    # [x] in GeoMean <=> exists xij s.t. A [x] + B xij in R- x RSOC^(N-1)
    # Hence the dual transformation is
    # A' y in GeoMean* <=>
    # B' y = 0         <=> y in R- x RSOC^(N-1)
    # Here A has at most one nonzero entry by line and column and this entry is 1.
    # Therefore we can A' = inv(A) hence we can use the same `getconstrattr` function
    # than for getting the `ConstraintPrimal`.
    # We don't need to use B' y = 0.
    return output
end


text:
"""
Created on Fri May 04 2020

@author: Yoann Pradat

    CentraleSupelec
    MICS laboratory
    9 rue Juliot Curie, Gif-Sur-Yvette, 91190 France

Algorithm for evaluating factorisation fidelity  in rank selection method in NMF based on average stability of latent
factor across perturbed versions of the input matrix V.
"""

abstract type AbstractRSFidNMFResults <: AbstractNMFResults end

"""
    RSFidNMFResults(list_err, list_fid)

NMFResults returned by `evaluate_fidelity`.
"""
mutable struct RSFidNMFResults{T} <: AbstractRSFidNMFResults
    list_V_err::Vector{Matrix{T}}
    list_fid::Vector{Float64}

    function RSFidNMFResults{T}(
        list_V_err::Vector{Matrix{T}},
        list_fid::Vector{Float64}) where T <: Real

        new{T}(list_V_err, list_fid)
    end
end

"""
    evaluate_fidelity(list_nmf_results, list_V_pert)

Evaluates the fidelity of multiple runs of a NMF factorisation by compute the reconstruction error across
runs. Note that the reconstruction error is different from the cost minimized in the NMF algorithm in case
the latter includes regularization terms.

#Arguments
- list_nmf_results::Vector{NMFResults{T}} : vector containing results of each NMF run
- list_V_pert::Vector{Matrix{T}}       : vector containing the perturbed V matrix of each NMF run
"""
function evaluate_fidelity(list_nmf_results::Vector{NMFResults{T}}, list_V_pert::Vector{Matrix{T}}) where T <: Real
    list_V_err = Vector{Matrix{T}}()
    list_fid   = Vector{Float64}()

    for (nmf_results, V_pert) in zip(list_nmf_results, list_V_pert)
        V_appr = nmf_results.W * nmf_results.H
        V_err  = V_pert - V_appr

        if typeof(nmf_results.local_params) == FIParams
            fid = β_divergence(V_pert, V_appr, nmf_results.local_params.β)
        else
            if nmf_results.local_params.div == :α
                fid = α_divergence(V_pert, V_appr, nmf_results.local_params.α)
            elseif nmf_results.local_params.div == :β
                fid = β_divergence(V_pert, V_appr, nmf_results.local_params.β)
            end
        end

        push!(list_V_err, V_err)
        push!(list_fid, fid)
    end

    RSFidNMFResults{T}(
        list_V_err,
        list_fid
    )
end


text:
#######################
## Unit Conversions ###
#######################

#' Conversion between Latent Heat Flux and Evapotranspiration
#'
#' converts evaporative water flux from mass (ET=evapotranspiration)
#'              to energy (LE=latent heat flux) units, or vice versa.
#'
#' @aliases LE_to_ET ET_to_LE
#'
#' - LE   Latent heat flux (W m-2)
#' - ET   Evapotranspiration (kg m-2 s-1)
#' - Tair Air temperature (deg C)
#'
#' # Details

#' The conversions are given by:
#'
#' ``ET = LE/\\lambda``
#'
#' ``LE = \\lambda ET``
#'
#' where ``\\lambda`` is the latent heat of vaporization (J kg-1) as calculated by
#' [`latent_heat_vaporization`](@ref).
#'
#' ```@example; output = false
#' ```
#' # LE of 200 Wm-2 and air temperature of 25degC
#' LE_to_ET(200,25)
#'
"""
"""
function LE_to_ET(LE,Tair)

  lambda = latent_heat_vaporization(Tair)
  ET     = LE/lambda

  return(ET)
end

#' @rdname LE_to_ET
"""
"""
function ET_to_LE(ET,Tair)

  lambda = latent_heat_vaporization(Tair)
  LE     = ET*lambda

  return(LE)
end

#' Conversion between Conductance Units
#'
#' Converts conductances from mass (m s-1)
#'              to molar units (mol m-2 s-1), or vice versa.
#'
#' @aliases ms_to_mol mol_to_ms
#'
#' - G_ms       Conductance (m s-1)
#' - G_mol      Conductance (mol m-2 s-1)
#' - Tair       Air temperature (deg C)
#' - pressure   Atmospheric pressure (kPa)
#' - constants  Kelvin - conversion degree Celsius to Kelvin
#'                   Rgas - universal gas constant (J mol-1 K-1)
#'                   kPa2Pa - conversion kilopascal (kPa) to pascal (Pa)
#'
#' # Details

#' The conversions are given by:
#'
#' ``G_mol = G_ms * pressure / (Rgas * Tair)``
#'
#' ``G_ms = G_mol * (Rgas * Tair) / pressure``
#'
#' where Tair is in Kelvin and pressure in Pa (converted from kPa internally)
#'
#' #References
#' Jones, H_G. 1992. Plants and microclimate: a quantitative approach to environmental plant physiology.
#'             2nd Edition., Cambridge University Press, Cambridge. 428 p
#'
#' ```@example; output = false
#' ```
#' ms_to_mol(0.005,25,100)
#'
"""
"""
function ms_to_mol(G_ms,Tair,pressure,constants=BigleafConstants())

  Tair     = Tair + constants.Kelvin
  pressure = pressure * constants.kPa2Pa

  G_mol  = G_ms * pressure / (constants.Rgas * Tair)

  return(G_mol)
end

#' @rdname ms_to_mol
"""
"""
function mol_to_ms(G_mol,Tair,pressure,constants=BigleafConstants())

  Tair     = Tair + constants.Kelvin
  pressure = pressure * constants.kPa2Pa

  G_ms  = G_mol * (constants.Rgas * Tair) / (pressure)

  return(G_ms)
end

#' Conversions between Humidity Measures
#'
#' Conversion between vapor pressure (e), vapor pressure deficit (VPD),
#'              specific humidity (q), and relative humidity (rH).
#'
#' - Tair      Air temperature (deg C)
#' - pressure  Atmospheric pressure (kPa)
#' - e         Vapor pressure (kPa)
#' - q         Specific humidity (kg kg-1)
#' - VPD       Vapor pressure deficit (kPa)
#' - rH        Relative humidity (-)
#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
#'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
#'                      See [`Esat_slope`](@ref).
#' - constants eps - ratio of the molecular weight of water vapor to dry air (-)
#'                  Pa2kPa - conversion pascal (Pa) to kilopascal (kPa)
#'
#' @family humidity conversion
#'
#' #References
#' Foken, T, 2008: Micrometeorology. Springer, Berlin, Germany.
#'
"""
"""
function VPD_to_rH(VPD,Tair,Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                      constants=BigleafConstants())

  esat = Esat_slope(Tair,Esat_formula,constants)[,"Esat"]
  rH   = 1 - VPD/esat
  return(rH)
end

#' @rdname VPD_to_rH
#' @family humidity conversion
"""
"""
function rH_to_VPD(rH,Tair,Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                      constants=BigleafConstants())

  if(any(rH > 1 & !ismissing(rH)))
    @warn"relative humidity (rH) has to be between 0 and 1.")
end

  esat = Esat_slope(Tair,Esat_formula,constants)[,"Esat"]
  VPD  = esat - rH*esat
  return(VPD)
end

#' @rdname VPD_to_rH
#' @family humidity conversion
"""
"""
function e_to_rH(e,Tair,Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                    constants=BigleafConstants())

  esat = Esat_slope(Tair,Esat_formula,constants)[,"Esat"]

  if (any(e > esat + .Machine$double_eps^0.5 & !ismissing(e)))
    @warn"Provided vapour pressure that was higher than saturation.
             Returning rH=1 for those cases.")
end

  rH  = pmin(1, e/esat)
  return(rH)
end

#' @rdname VPD_to_rH
#' @family humidity conversion
"""
"""
function VPD_to_e(VPD,Tair,Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                     constants=BigleafConstants())

  esat = Esat_slope(Tair,Esat_formula,constants)[,"Esat"]
  e    = esat - VPD
  return(e)
end

#' @rdname VPD_to_rH
#' @family humidity conversion
"""
"""
function e_to_VPD(e,Tair,Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                     constants=BigleafConstants())

  esat = Esat_slope(Tair,Esat_formula,constants)[,"Esat"]
  VPD  = esat - e
  return(VPD)
end

#' @rdname VPD_to_rH
#' @family humidity conversion
"""
"""
function e_to_q(e,pressure,constants=BigleafConstants())
  q = constants.eps] * e / (pressure - (1-constants[:eps) * e)
  return(q)
end

#' @rdname VPD_to_rH
#' @family humidity conversion
"""
"""
function q_to_e(q,pressure,constants=BigleafConstants())
  e = q * pressure / ((1-constants.eps]) * q + constants[:eps)
  return(e)
end

#' @rdname VPD_to_rH
#' @family humidity conversion
"""
"""
function q_to_VPD(q,Tair,pressure,Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                     constants=BigleafConstants())

  esat = Esat_slope(Tair,Esat_formula,constants)[,"Esat"]
  e    = q_to_e(q,pressure,constants)
  VPD  = esat - e
  return(VPD)
end

#' @rdname VPD_to_rH
#' @family humidity conversion
"""
"""
function VPD_to_q(VPD,Tair,pressure,Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                     constants=BigleafConstants())

  esat = Esat_slope(Tair,Esat_formula,constants)[,"Esat"]
  e    = esat - VPD
  q    = e_to_q(e,pressure,constants)
  return(q)
end

#' Conversions between Global Radiation and Photosynthetic Photon Flux Density
#'
#' Converts radiation from W m-2 to umol m-2 s-1 and vice versa.
#'
#' - Rg       Global radiation = incoming short-wave radiation at the surface (W m-2)
#' - PPFD     Photosynthetic photon flux density (umol m-2 s-1)
#' - J_to_mol Conversion factor from J m-2 s-1 (= W m-2) to umol (quanta) m-2 s-1
#' - frac_PAR Fraction of incoming solar irradiance that is photosynthetically
#'                 active radiation (PAR); defaults to 0.5
#'
#' # Details

#' The conversion is given by:
#'
#'  ``PPFD = Rg * frac_PAR * J_to_mol``
#'
#' by default, the combined conversion factor (`frac_PAR * J_to_mol`) is 2.3
#'
#' ```@example; output = false
#' ```
#' # convert a measured incoming short-wave radiation of 500 Wm-2 to
#' # PPFD in umol m-2 s-1 and backwards
#' Rg_to_PPFD(500)
#' PPFD_to_Rg(1150)
#'
"""
"""
function Rg_to_PPFD(Rg,J_to_mol=4.6,frac_PAR=0.5)
  PPFD = Rg * frac_PAR * J_to_mol
  return(PPFD)
end

#' @rdname Rg_to_PPFD
"""
"""
function PPFD_to_Rg(PPFD,J_to_mol=4.6,frac_PAR=0.5)
  Rg = PPFD / frac_PAR / J_to_mol
  return(Rg)
end

#' Conversion between Mass and Molar Units
#'
#' Converts mass units of a substance to the corresponding molar units
#'              and vice versa.
#'
#' - mass      Numeric vector of mass in kg
#' - molarMass Numeric vector of molar mass of the substance (kg mol-1)
#'                  e.g. as provided by [`bigleaf_constants`](@ref)()$H2Omol
#'                  Default is molar mass of Water.
#'
#' # Value
 Numeric vector of amount of substance in mol.
"""
"""
function kg_to_mol(mass, molarMass=BigleafConstants()$H2Omol)

  moles = mass / molarMass

  return(moles)

end

#' Conversion between Mass and Molar Units of Carbon and CO2
#'
#' Converts CO2 quantities from umol CO2 m-2 s-1 to g C m-2 d-1 and vice versa.
#'
#' - CO2_flux  CO2 flux (umol CO2 m-2 s-1)
#' - C_flux    Carbon (C) flux (gC m-2 d-1)
#' - constants Cmol - molar mass of carbon (kg mol-1)
#'                  umol2mol - conversion micromole (umol) to mol (mol)
#'                  mol2umol - conversion mole (mol) to micromole (umol)
#'                  kg2g - conversion kilogram (kg) to gram (g)
#'                  g2kg - conversion gram (g) to kilogram (kg)
#'                  days2seconds - seconds per day
#'
#' ```@example; output = false
#' ```
#' umolCO2_to_gC(20)  # gC m-2 d-1
#'
"""
"""
function umolCO2_to_gC(CO2_flux,constants=BigleafConstants())

  C_flux = CO2_flux * constants.umol2mol] * constants[:Cmol] * constants[:kg2g] * constants[:days2seconds

  return(C_flux)
end

#' @rdname umolCO2_to_gC
"""
"""
function gC_to_umolCO2(C_flux,constants=BigleafConstants())

  CO2_flux = (C_flux * constants.g2kg] / constants[:days2seconds]) / constants[:Cmol] * constants[:mol2umol

  return(CO2_flux)
end


text:
# With this file you may reproduce the results for cases 1, 2, 3, 4 in Table 2,
# by changing the number of neurons (variable neurons below), and ϵ (tol1).
# The activation function and it's inverse may be changed in the main file ANNBN.jl
# The Random Forests and Gradient Boosting calls utilized in the manuscript are aplso appended bellow

# 1st Step: Read the MNIST Database
using Plots, Statistics, MLDatasets, LinearAlgebra, Clustering, Printf, Dates, Random, MLDatasets, Plots
using Printf, Debugger, IterativeSolvers
train_x, yy_train_all = MLDatasets.MNIST.traindata()
test_x,  yy_test_all  = MLDatasets.MNIST.testdata()
xx_train=MLDatasets.MNIST.convert2features(MLDatasets.MNIST.traintensor())'
xx_test=MLDatasets.MNIST.convert2features(test_x)'
xx_train=convert(Array{Float64,2},xx_train)
xx_test=convert(Array{Float64,2},xx_test)
yy_train_all=convert(Array{Float64,1},yy_train_all)
yy_test_all=convert(Array{Float64,1},yy_test_all)
i_train=size(xx_train,1)
i_test=size(xx_test,1)
vars=size(xx_train,2)

# Import ANNBN module
path1=realpath(dirname(@__FILE__)*"/../..")
include(string(path1,"/src/ANNBN.jl"))

# 2nd Step: Train for a single digit num_to_test. tol1=ϵ normalization in manuscript
num_to_test=0;tol1=0.02;
yy_train=copy(yy_train_all)
i1=abs.(yy_train.-num_to_test).<0.01;i2=abs.(yy_train.-num_to_test).>0.01;yy_train[i1].=1;;yy_train[i2].=0;
mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=1.0-2tol1;yy_train.+=tol1
yy_test=copy(yy_test_all)
i1=abs.(yy_test.-num_to_test).<0.01;i2=abs.(yy_test.-num_to_test).>0.01;yy_test[i1].=1;yy_test[i2].=0;
yy_test.-=mi;yy_test./=ma;yy_test.*=1.0-2tol1;yy_test.+=tol1

neurons=1000; # For higher accuracy, increase the number of neurons. RAM demands will increase as well.
inds_all=1:i_train;items_per_neuron=(Int64(floor(i_train/(neurons))))*ones(Int64,neurons);n_per_part=[0;cumsum(items_per_neuron)];n_per_part[end]=i_train
@time a_all,a_layer1,layer1,mat1=ANNBN.train_layer_1_sigmoid_fast(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train)
predl1=[layer1 ones(i_train)]*a_layer1
maetr=mean(abs.(yy_train-predl1))
predl1_test,layer1_test=ANNBN.predict_new(a_all,a_layer1,xx_test,i_test,neurons);
maete=mean(abs.(yy_test-predl1_test))
cctr=100mean(abs.(round.(predl1).-round.(yy_train)).<0.1)
ccte=100mean(abs.(round.(predl1_test).-round.(yy_test)).<0.1)

# 3rd Step: Train all digits
cctrs=Array{Float64}(undef,0)
cctes=Array{Float64}(undef,0)
dts=Array{Float64}(undef,0)
predl1_train_all=Array{Float64}(undef,i_train,0)
predl1_test_all=Array{Float64}(undef,i_test,0)
neurons=1000;tol1=0.02;
@time for num_to_test=0:9
    global yy_train_all,yy_test_all,predl1_train_all,predl1_test_all,neurons,tol1
    yy_train=copy(yy_train_all)
    i1=abs.(yy_train.-num_to_test).<0.01;i2=abs.(yy_train.-num_to_test).>0.01;yy_train[i1].=1;;yy_train[i2].=0;
    mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=1.0-2tol1;yy_train.+=tol1
    yy_test=copy(yy_test_all)
    i1=abs.(yy_test.-num_to_test).<0.01;i2=abs.(yy_test.-num_to_test).>0.01;yy_test[i1].=1;yy_test[i2].=0;
    yy_test.-=mi;yy_test./=ma;yy_test.*=1.0-2tol1;yy_test.+=tol1
    inds_all=1:i_train;items_per_neuron=(Int64(floor(i_train/(neurons))))*ones(Int64,neurons);n_per_part=[0;cumsum(items_per_neuron)];n_per_part[end]=i_train
    t1=now();a_all,a_layer1,layer1=ANNBN.train_layer_1_sigmoid_fast(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train)
    t2=now();dt=convert(Float64,Dates.value(t2-t1))/1000;push!(dts,dt)
    predl1=[layer1 ones(i_train)]*a_layer1; maetr=mean(abs.(yy_train-predl1))
    predl1_test,layer1_test=ANNBN.predict_new(a_all,a_layer1,xx_test,i_test,neurons); maete=mean(abs.(yy_test-predl1_test))

    predl1_train_all=[predl1_train_all predl1]; predl1_test_all=[predl1_test_all predl1_test]
    cctr=100mean(abs.(round.(predl1).-round.(yy_train)).<0.1)
    ccte=100mean(abs.(round.(predl1_test).-round.(yy_test)).<0.1)
    push!(cctrs,cctr);push!(cctes,ccte)
    println(num_to_test," ",ccte)
end
println(dts)
println(cctes)

yy_train_pred=zeros(i_train)
for i=1:i_train
    yy_train_pred[i]=sortperm(predl1_train_all[i,:],rev=true)[1]-1
end
CC=100-100mean(yy_train_all-yy_train_pred.>0.01)

yy_test_pred=zeros(i_test)
for i=1:i_test
    yy_test_pred[i]=sortperm(predl1_test_all[i,:],rev=true)[1]-1
end
CC=100-100mean(yy_test_all-yy_test_pred.>0.01)

# # # # # # # # # # # Random Forests
using DecisionTree
num_trees=261
regr_3 = RandomForestClassifier(n_trees=num_trees);
t1=now()
DecisionTree.fit!(regr_3, xx_train, yy_train);
t2=now()
convert(Float64,Dates.value(t2-t1))/1000
predl1 = DecisionTree.predict(regr_3, xx_train)
predl1_test = DecisionTree.predict(regr_3, xx_test)
i1=predl1.<=0.5;predl1[i1].=tol1;i1=predl1.>0.5;predl1[i1].=1.0-tol1;i2=abs.(predl1.-yy_train).<0.01;maetr=100sum(i2)/length(i2)
i1=predl1_test.<=0.5;predl1_test[i1].=tol1;i1=predl1_test.>0.5;predl1_test[i1].=1.0-tol1;i2=abs.(predl1_test.-yy_test).<0.01;
maete=100sum(i2)/length(i2)

maetrs=Array{Float64}(undef,0)
maetes=Array{Float64}(undef,0)
dts=Array{Float64}(undef,0)
@time for num_to_test=0:9
    global yy_train_all,yy_test_all

    yy_train=copy(yy_train_all)
    yy_test=copy(yy_test_all)
    tol1=0.0
    yy_train=copy(yy_train_all)
    i1=abs.(yy_train.-num_to_test).<0.01;i2=abs.(yy_train.-num_to_test).>0.01;yy_train[i1].=1;;yy_train[i2].=0;
    mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=1.0-2tol1;yy_train.+=tol1
    yy_test=copy(yy_test_all)
    i1=abs.(yy_test.-num_to_test).<0.01;i2=abs.(yy_test.-num_to_test).>0.01;yy_test[i1].=1;yy_test[i2].=0;
    yy_test.-=mi;yy_test./=ma;yy_test.*=1.0-2tol1;yy_test.+=tol1

    num_trees=261
    regr_3 = RandomForestClassifier(n_trees=num_trees);
    t1=now()
    DecisionTree.fit!(regr_3, xx_train, yy_train);
    t2=now();dt=convert(Float64,Dates.value(t2-t1))/1000;push!(dts,dt)
    predl1 = DecisionTree.predict(regr_3, xx_train)
    predl1_test = DecisionTree.predict(regr_3, xx_test)

    i1=predl1.<=0.5;predl1[i1].=tol1;i1=predl1.>0.5;predl1[i1].=1.0-tol1;i2=abs.(predl1.-yy_train).<0.01;maetr=100sum(i2)/length(i2)
    i1=predl1_test.<=0.5;predl1_test[i1].=tol1;i1=predl1_test.>0.5;predl1_test[i1].=1.0-tol1;i2=abs.(predl1_test.-yy_test).<0.01;maete=100sum(i2)/length(i2)
    push!(maetrs,maetr)
    push!(maetes,maete)
    println(" > > > > > > > > > > > > >",num_to_test," ",maete," ",dt)
end
println(dts)
println(maetes)

# # # # # # # # # # # # # # #  XGBoost
using XGBoost
num_round=200
@time bst = xgboost(xx_train, num_round, label = yy_train, eta = 1, max_depth = 7, verbose=false)
predl1 = XGBoost.predict(bst, copy(xx_train))
predl1_test = XGBoost.predict(bst, copy(xx_test))
i1=predl1.<=0.5;predl1[i1].=tol1;i1=predl1.>0.5;predl1[i1].=1.0-tol1;i2=abs.(predl1.-yy_train).<0.01;maetr=100sum(i2)/length(i2)
i1=predl1_test.<=0.5;predl1_test[i1].=tol1;i1=predl1_test.>0.5;predl1_test[i1].=1.0-tol1;i2=abs.(predl1_test.-yy_test).<0.01;
maete=100sum(i2)/length(i2)

num_round = 200
maetrs=Array{Float64}(undef,0)
maetes=Array{Float64}(undef,0)
dts=Array{Float64}(undef,0)
@time for num_to_test=0:9
    global yy_train_all,yy_test_all,num_round

    yy_train=copy(yy_train_all)
    yy_test=copy(yy_test_all)
    tol1=0.0
    yy_train=copy(yy_train_all)
    i1=abs.(yy_train.-num_to_test).<0.01;i2=abs.(yy_train.-num_to_test).>0.01;yy_train[i1].=1;;yy_train[i2].=0;
    mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=1.0-2tol1;yy_train.+=tol1
    yy_test=copy(yy_test_all)
    i1=abs.(yy_test.-num_to_test).<0.01;i2=abs.(yy_test.-num_to_test).>0.01;yy_test[i1].=1;yy_test[i2].=0;
    yy_test.-=mi;yy_test./=ma;yy_test.*=1.0-2tol1;yy_test.+=tol1
    t1=now()
    bst = xgboost(xx_train, num_round, label = yy_train, eta = 1, max_depth = 7)
    t2=now();dt=convert(Float64,Dates.value(t2-t1))/1000;push!(dts,dt)
    predl1 = XGBoost.predict(bst, xx_train)
    predl1_test = XGBoost.predict(bst, xx_test)

    i1=predl1.<=0.5;predl1[i1].=tol1;i1=predl1.>0.5;predl1[i1].=1.0-tol1;i2=abs.(predl1.-yy_train).<0.01;maetr=100sum(i2)/length(i2)
    i1=predl1_test.<=0.5;predl1_test[i1].=tol1;i1=predl1_test.>0.5;predl1_test[i1].=1.0-tol1;i2=abs.(predl1_test.-yy_test).<0.01;maete=100sum(i2)/length(i2)
    push!(maetrs,maetr)
    push!(maetes,maete)
    println(" > > > > > > > > > > > > >",num_to_test," ",maete," ",dt)
end
println(dts)
println(maetes)


text:
using DrWatson, Test

ps = DrWatson.PATH_SEPARATOR

# Tests for parse_savename

function test_convert(prefix::AbstractString,c;kwargs...)
    name = savename(prefix,c;kwargs...)
    _prefix, _b, _suffix = DrWatson.parse_savename(name)
    dicts_equal(_b,c) && _prefix == prefix && _suffix == ""
end

function test_convert(c,suffix::AbstractString;kwargs...)
    name = savename(c,suffix;kwargs...)
    _prefix, _b, _suffix = DrWatson.parse_savename(name)
    dicts_equal(_b,c) && _suffix == suffix && _prefix == ""
end

function test_convert(prefix::AbstractString,c,suffix::AbstractString;kwargs...)
    name = savename(prefix,c,suffix;kwargs...)
    _prefix, _b, _suffix = DrWatson.parse_savename(name)
    dicts_equal(_b,c) && _suffix == suffix && _prefix == prefix
end

function test_convert(c;kwargs...)
    name = savename(c;kwargs...)
    _prefix, _b, _suffix = DrWatson.parse_savename(name)
    dicts_equal(_b,c) && _prefix == "" && _suffix == ""
end

function dicts_equal(a,b)
    kA,kB = keys(a), keys(b)
    kA != kB && return false
    for k ∈ kA
        a[k] != b[k] && return false
    end
    return true
end

@test test_convert(
    Dict("c" => 0.1534, "u" => 5.1, "r"=>101, "mode" => "double"),
    digits=4)
@test test_convert("prefix",
                   Dict("c" => 0.1534, "u" => 5.1, "r"=>101, "mode" => "double"),
                   digits=4)
@test test_convert(
    Dict("c" => 0.1534, "u" => 5.1, "r"=>101, "mode" => "double"),
    "suffix",
    digits=4)
@test test_convert("prefix",
                   Dict("c" => 0.1534, "u" => 5.1, "r"=>101, "mode" => "double"),
                   "suffix",
                   digits=4)
@test test_convert("prefix",
                   Dict("c" => 0.1534, "u" => 5.1, "r"=>101, "variable_with_underscore" => "double"),
                   "suffix",
                   digits=4)
@test test_convert(joinpath("a=10_mode=double","prefix"),
                   Dict("c" => 0.1534, "u" => 5.1, "r"=>101, "variable_with_underscore" => "double"),
                   "suffix",
                   digits=4)
@test test_convert(Dict("c" => 0.1534, "u" => 5.1),
                   digits=4)
@test test_convert(Dict("never" => "gonna", "give" => "you", "up" => "!"))
@test test_convert(Dict("c" => 0.1534),
                   digits=4)

b = Dict("c" => 0.1534, "u" => 5.1, "r"=>101, "mo_de" => "double")
name = savename("prefix",b,connector="_",digits=4)
_prefix, _b, _suffix = DrWatson.parse_savename(name,connector="_")
@test dicts_equal(_b,b) && _prefix == "prefix" && _suffix == ""

_prefix, _b, _suffix = DrWatson.parse_savename(joinpath("some_random_path_a=10.0","prefix")*"_a=10_just_a_string=I'm not allowed to use underscores here_my_value=10.1.suffix_with_underscore-but-don't-use-dots")
@test _prefix == joinpath("some_random_path_a=10.0","prefix")
@test _suffix == "suffix_with_underscore-but-don't-use-dots"
@test _b["a"] == 10.0
@test _b["just_a_string"] == "I'm not allowed to use underscores here"
@test _b["my_value"] == 10.1

@test_throws ErrorException DrWatson.parse_savename("a=10",connector="__")
@test_throws ErrorException("Savename cannot be parsed. There is a '_' after the last '='. "*
        "Values containing '_' are not allowed when parsing.") DrWatson.parse_savename("a=10_1")

unicode_chars_name="3pb_Gc₁=3_Gc₂=0.1_or=T_α₁=4_α₂=1_β=3.vtu"

@test parse_savename(unicode_chars_name) == ("3pb",Dict(
    "Gc₁"=>3,
    "Gc₂"=>0.1,
    "or"=>"T",
    "α₁"=>4,
    "α₂"=>1,
    "β"=>3),"vtu")


text:
using Test
import Base.//

isapproxstr(s1::AbstractString, s2::AbstractString) =
    isequal(map(s->replace(s, r"\s|\n"=>""), String.((s1, s2)))...)

# stricter than isapproxstr, just strips the outside.
(//)(s1::AbstractString, s2::AbstractString) = strip(s1) == strip(s2)
(//)(o::FP.AbstractSpan, s::AbstractString)  = o.ss // s

# --------------------------------------------------------------------------------------

md_blockifier = s -> FP.default_md_tokenizer(s) |> FP.default_md_blockifier

function check_tokens(tokens, idx, name)
    @test all([t.name for t in tokens[idx]] .== name)
    compl = [i for i in eachindex(tokens) if i ∉ idx]
    @test all([t.name for t in tokens[compl]] .!= name)
end

toks    = FP.default_md_tokenizer
slice   = FP.md_partition
text(b) = FP.prepare_text(b)
ct(b)   = FP.content(b)
ctf(b::FP.Group) = FP.content(first(b.blocks))
grouper = FP.md_grouper ∘ slice
isp(g)  = g.role == :PARAGRAPH
printel(bv) = foreach(e -> println(strip(e.ss)), bv)


text:
using DiffEqBase, DiffEqJump, StaticArrays
using Test

reactstoch = [
    [1 => 1]
]

netstoch = [
    [1 => -1]
]

rates = [1.]
spec_to_dep_jumps = [[1]]
jump_to_dep_specs = [[1]]
dg = [[1]]
majump = MassActionJump(rates, reactstoch, netstoch)
u0 = [10]
dprob = DiscreteProblem(u0,(0.,100.),rates)
algs = DiffEqJump.JUMP_AGGREGATORS

for ssa in algs
    jprob = JumpProblem(dprob, ssa, majump; vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs, save_positions=(false,false))
    sol = solve(jprob, SSAStepper(), saveat=100.)
    @test sol[1,end] == 0
    @test sol.t[end] < Inf
end

u0 = SA[10]
dprob = DiscreteProblem(u0,(0.,100.),rates)

for ssa in algs
    jprob = JumpProblem(dprob, ssa, majump; vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs, save_positions=(false,false))
    sol = solve(jprob, SSAStepper(), saveat=100.)
    @test sol[1,end] == 0
    @test sol.t[end] < Inf
end


text:
include("../spear.jl")
activities = int(readcsv("activities.csv"))

C(x) = x.^.5      # credit scoring function
A = generateAdjacencyMatrix(activities, C)
(expertise, quality) = runSpear(A)

# Display the results
users = readcsv("users.csv")
resources = readcsv("resources.csv")

iSortUsers = reverse(sortperm(vec(expertise)))
sortedPrintU = reduce(*,map((i -> "$(users[i]) =>  $(round(expertise[i],3))\n"),iSortUsers))
println("Users sorted by expertise \n$sortedPrintU")

iSortResources = reverse(sortperm(vec(quality)))
sortedPrintR = reduce(*,map((i -> "$(resources[i]) =>  $(round(quality[i],3))\n"),iSortResources))
println("Resources sorted by quality \n$sortedPrintR")


text:
# A framework for implementing CoR
#
# A chain of functions are executed but we allow early
# exit when a function decides that the request is
# completely handled e.g. A CoR implementation of
# HTTP handlers may skip the rest of the handlers when
# an error is detected early in the chain.

module ChainOfResponsibilityExample

mutable struct DepositRequest
    id::Int
    amount::Float64
end

@enum Status CONTINUE HANDLED

"""
Apply handlers.  Stop the chain when a handler
return `HANDLED`.
"""
function apply(req::DepositRequest, handlers::AbstractVector{Function})
    for f in handlers
        status = f(req)
        status == HANDLED && return nothing
    end
end

# Example usage

function update_account_handler(req::DepositRequest)
    println("Deposited $(req.amount) to account $(req.id)")
    return CONTINUE
end

function send_gift_handler(req::DepositRequest)
    req.amount > 100_000 &&
        println("=> Thank you for your business")
    return CONTINUE
end

function notify_customer(req::DepositRequest)
    println("deposit is finished")
    return HANDLED
end

handlers = [
    update_account_handler,
    send_gift_handler,
    notify_customer
]

function test()
    println("Test: customer depositing a lot of money")
    amount = 300_000
    apply(DepositRequest(1, amount), handlers)

    println("\nTest: regular customer")
    amount = 1000
    apply(DepositRequest(2, amount), handlers)
end

end #module

using .ChainOfResponsibilityExample
ChainOfResponsibilityExample.test()

#=
julia> amount = 100_000
100000

julia> apply(DepositRequest(amount), handlers)
Your new balance is 105000.0
Thank you for your business
finalized request, ok to exit

julia> amount = 1000
1000

julia> apply(DepositRequest(amount), handlers)
Your new balance is 1050.0
finalized request, ok to exit
=#


text:
include("../imports.jl")

@testset ExtendedTestSet "`Recorder`" begin
    learner = testlearner(coeff = 3, Recorder())
    FluxTraining.handle(FluxTraining.Events.Init(), learner, TrainingPhase())
    h = learner.cbstate.history
    @test h.epochs == 0
    @test h.steps == 0
    fit!(learner, 10)
    h = learner.cbstate.history
    @test h.epochs == 10
    @test h.stepsepoch == length(getdataiter(TrainingPhase(), learner))
    @test h.steps == 10 * length(getdataiter(TrainingPhase(), learner))
end


text:
module CMBLensing

using Adapt
using Base.Broadcast: AbstractArrayStyle, ArrayStyle, Broadcasted,
    DefaultArrayStyle, preprocess_args, Style, result_style, Unknown
using Base.Iterators: flatten, product, repeated, cycle, countfrom, peel, partition
using Base.Threads
using Base: @kwdef, @propagate_inbounds, Bottom, OneTo, showarg, show_datatype,
    show_default, show_vector, typed_vcat, typename
using Combinatorics
using CompositeStructs
using DataStructures
using DelimitedFiles
using Distributed: pmap, nworkers, myid, workers, addprocs, @everywhere, remotecall_wait,
    @spawnat, pgenerate, procs, @fetchfrom, default_worker_pool, RemoteChannel
using EllipsisNotation
using FileIO
using FFTW
using InteractiveUtils
using IterTools: flagfirst
using JLD2
using JLD2: jldopen, JLDWriteSession
using KahanSummation
using Loess
using LinearAlgebra
using LinearAlgebra: diagzero, matprod, promote_op
using MacroTools: @capture, combinedef, isdef, isexpr, postwalk, prewalk, rmlines, splitdef
using Match
using Markdown
using Measurements
using Memoization
using NamedTupleTools: select, delete
using Pkg
using Printf
using ProgressMeter
using QuadGK
using Random
using Random: seed!, AbstractRNG
using Roots
using Requires
using Serialization
using Setfield
using SparseArrays
using StaticArrays: @SMatrix, @SVector, SMatrix, StaticArray, StaticArrayStyle,
    StaticMatrix, StaticVector, SVector, SArray, SizedArray
using Statistics
using StatsBase
using TimerOutputs: @timeit, get_defaulttimer, reset_timer!
using Tullio
using UnPack
using Zygote
using Zygote: unbroadcast, Numeric, @adjoint, @nograd
using Zygote.ChainRules: @thunk, NoTangent

import Adapt: adapt_structure
import Base: +, -, *, \, /, ^, ~, ≈, <, <=, |, &, ==, !,
    abs, adjoint, all, any, axes, broadcast, broadcastable, BroadcastStyle, conj, copy, convert,
    copy, copyto!, eltype, eps, fill!, getindex, getproperty, hash, hcat, hvcat, inv, isfinite,
    iterate, keys, lastindex, length, literal_pow, mapreduce, materialize!,
    materialize, one, permutedims, print_array, promote, promote_rule,
    promote_rule, promote_type, propertynames, real, setindex!, setproperty!, show,
    show_datatype, show_vector, similar, size, sqrt, string, sum, summary,
    transpose, zero
import Base.Broadcast: materialize, preprocess, broadcasted
import Zygote.ChainRules: rrule
import LinearAlgebra: checksquare, diag, dot, isnan, ldiv!, logdet, mul!, norm,
    pinv, StructuredMatrixStyle, structured_broadcast_alloc, tr
import Measurements: ±
import Statistics: std

export
    @⌛, @show⌛, @ismain, @namedtuple, @repeated, @unpack, @cpu!, @gpu!, @cu!, animate,
    argmaxf_lnP, AzFourier, BandPassOp, BaseDataSet, batch, batch_index, batch_length,
    BlockDiagEquiRect, beamCℓs, cache, CachedLenseFlow, camb, cov_to_Cℓ, cpu, Cℓ_2D,
    Cℓ_to_Cov, DataSet, DerivBasis, diag, Diagonal, DiagOp, dot, EBFourier, EBMap, expnorm,
    Field, FieldArray, fieldinfo, FieldMatrix, FieldOrOpArray, FieldOrOpMatrix, FieldOrOpRowVector,
    FieldOrOpVector, FieldRowVector, FieldTuple, FieldVector, FieldVector,
    firsthalf, fixed_white_noise, BlockDiagIEB, Fourier, FuncOp, get_max_lensing_step,
    get_Cℓ, get_Cℓ, get_Dℓ, get_αℓⁿCℓ, get_ρℓ, get_ℓ⁴Cℓ, gradhess, gradient, HighPass,
    IEBFourier, IEBMap, InterpolatedCℓs, IQUAzFourier, IQUFourier, IQUMap, kde,
    lasthalf, LazyBinaryOp, LenseBasis, LenseFlow, FieldOp, lnP, load_camb_Cℓs,
    load_chains, load_nolensing_sim, load_sim, LowPass, make_mask, Map, MAP_joint, MAP_marg,
    mean_std_and_errors, MidPass, mix, nan2zero, noiseCℓs,
    ParamDependentOp, pixwin, PowerLens, ProjLambert, ProjEquiRect, ProjHealpix, project,
    QUAzFourier, QUFourier, QUMap, resimulate!, resimulate, RK4Solver, sample_joint, shiftℓ,
    simulate, SymmetricFuncOp, symplectic_integrate, Taylens, toCℓ, toDℓ,
    ud_grade, unbatch, unmix, white_noise, Ð, Ł,
    ℓ², ℓ⁴, ∇, ∇², ∇ᵢ, ∇ⁱ

# bunch of sampling-related exports
export gibbs_initialize_f!, gibbs_initialize_ϕ!, gibbs_initialize_θ!,
    gibbs_sample_f!, gibbs_sample_ϕ!, gibbs_sample_slice_θ!,
    gibbs_mix!, gibbs_unmix!, gibbs_postprocess!,
    once_every, start_after_burnin, mass_matrix_ϕ, hmc_step

# util
include("util.jl")
include("util_fft.jl")
include("util_parallel.jl")
include("numerical_algorithms.jl")

# generic field / operator structure
include("generic.jl")
include("cls.jl")
include("field_tuples.jl")
include("field_vectors.jl")
include("base_fields.jl")
include("specialops.jl")
include("flowops.jl")
include("batched_reals.jl")

# lensing operators
include("lenseflow.jl")
include("powerlens.jl")

# field types
include("proj_cartesian.jl")
include("proj_lambert.jl")
include("proj_equirect.jl")
include("proj_healpix.jl")

# other field-specific stuff
include("masking.jl")
include("taylens.jl")
include("bilinearlens.jl")

# plotting
function animate end
@init @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" include("plotting.jl")

# sampling and maximizing the posteriors
include("dataset.jl")
include("posterior.jl")
include("maximization.jl")
include("sampling.jl")
include("chains.jl")

# other estimates
include("quadratic_estimate.jl")

# AD
include("autodiff.jl")

# gpu
is_gpu_backed(x) = false
@init @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" include("gpu.jl")

# misc init
# see https://github.com/timholy/ProgressMeter.jl/issues/71 and links therein
@init if ProgressMeter.@isdefined ijulia_behavior
    ProgressMeter.ijulia_behavior(:clear)
end

end


text:
# call{P <: Parameter}(param::P, t::Float64) = value(param, t)

mutable struct ConstantParameter{C <: Constraint} <: Parameter
  data::Vector{Float64}
  constraint::C
end

get_data(c::ConstantParameter) = c.data
set_params!(c::ConstantParameter, i::Int, val::Float64) = c.data[i] = val
# call(c::ConstantParameter, t::Float64) = value(c, t)

mutable struct G2FittingParameter{T <: TermStructure} <: Parameter
  a::Float64
  sigma::Float64
  b::Float64
  eta::Float64
  rho::Float64
  ts::T
end

(g::G2FittingParameter)(t::Float64) = value(g, t)

function value(param::G2FittingParameter, t::Float64)
  forward = forward_rate(param.ts, t, t, ContinuousCompounding(), NoFrequency()).rate
  temp1 = param.sigma * (-expm1(-param.a * t)) / param.a
  temp2 = param.eta * (-expm1(-param.b * t)) / param.b

  val = 0.5 * temp1 * temp1 + 0.5 * temp2 * temp2 + param.rho * temp1 * temp2 + forward

  return val
end

mutable struct HullWhiteFittingParameter{T <: TermStructure} <: Parameter
  a::Float64
  sigma::Float64
  ts::T
end

(h::HullWhiteFittingParameter)(t::Float64) = value(h, t)

function value(param::HullWhiteFittingParameter, t::Float64)
  forward = forward_rate(param.ts, t, t, ContinuousCompounding(), NoFrequency()).rate
  temp = param.a < sqrt(eps()) ? param.sigma * t : param.sigma * (-expm1(-param.a * t)) / param.a

  return forward + 0.5 * temp * temp
end

mutable struct TermStructureFittingParameter{T <: TermStructure} <: Parameter
  times::Vector{Float64}
  values::Vector{Float64}
  ts::T
end

TermStructureFittingParameter(ts::T) where {T <: TermStructure} = TermStructureFittingParameter{T}(zeros(0), zeros(0), ts)

(tsp::TermStructureFittingParameter)(t::Float64) = value(tsp, t)

function reset_param_impl!(param::TermStructureFittingParameter)
  param.times = zeros(length(param.times))
  param.values = zeros(length(param.values))

  return param
end

function set_params!(param::TermStructureFittingParameter, tm::Float64, val::Float64)
  push!(param.times, tm)
  push!(param.values, val)

  return param
end

function value(param::TermStructureFittingParameter, t::Float64)
  idx = findfirst(isequal(t), param.times)
  return param.values[idx]
end

mutable struct PiecewiseConstantParameter{C <: Constraint} <: Parameter
  times::Vector{Float64}
  constraint::C

  function PiecewiseConstantParameter{C}(times::Vector{Float64}, constraint::C) where C
    retTimes = push!(times, 0.0)
    return new{C}(retTimes, constraint)
  end
end

PiecewiseConstantParameter(times::Vector{Float64}, constraint::C) where {C <: Constraint} = PiecewiseConstantParameter{C}(times, constraint)

# call(p::PiecewiseConstantParameter, t::Float64) = value(p, t)

set_params!(param::PiecewiseConstantParameter, i::Int, val::Float64) = param.times[i] = val
get_data(param::PiecewiseConstantParameter) = param.times

NullParameter(_type::DataType) = _type([0.0], NoConstraint())

test_params(c::ConstantParameter, params::Vector{Float64}) = QuantLib.Math.test(c.constraint, params)


text:
# Autogenerated wrapper script for nghttp2_jll for x86_64-w64-mingw32
export libnghttp2

JLLWrappers.@generate_wrapper_header("nghttp2")
JLLWrappers.@declare_library_product(libnghttp2, "libnghttp2-14.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libnghttp2,
        "bin\\libnghttp2-14.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text: abstract type Layer{T<:Measurement} end


text:
# Calculation of T-matrix using only real numbers.
# Unlike the complex version, this is compatible with automatic differentiation

# TODO:
# 1- Find a way to enforce concrrete types at the T_matrix function. So far, I have defined abstract types as inputs for all functions. I am not sure how to convert to a concrete type mathcing my computation requirement (e.g., BigFloats)
# 2- Check why trapz_ELZOUKA is not returning exactly the same result as Trapz.trapz
# 3- The functions I have wrote separating real and imaginary parts are slower, figure out why
# 4- I need to get rid of Abstract Types, according to https://docs.julialang.org/en/v1/manual/performance-tips/ . For now, I am not sure how to make my code work with arbitrary float size
# -- I have tried to use r_array::Array{Float64}, but got an error when r_array is ::Vector{Float64}
# 5- Use StaticArrays when possible.

# NOTE: Lessons learned:
# 1- In order to be compatible with auto differentiation, all functions should return only one Array (no Tuples), which contains no nested arrays. This can be done by vcat(array...)

include("utils.jl")
include("geometry.jl")
include("electromagnetics.jl")

using Tmatrix
using ComplexOperations
using LinearAlgebra
using VectorSphericalWaves
using Trapz
using StaticArrays

export calculate_Tmatrix_for_spheroid_SeparateRealImag

"""
    Separate real and imaginary parts. I assume that A is a hcat of real and imag parts
"""
function separate_real_imag(A)
    return A[:,1:Int(end / 2)], A[:,Int(end / 2) + 1:end]
end

"""
    return the complex matrix, given matrix A which is hcat of real and imag parts
"""
function get_complex_matrix_from_concatenated_real_imag(A)
    A_r, A_i = separate_real_imag(A)
    return A_r + im * A_i
end

function J_mn_m_n__integrand_SeparateRealImag(
    m::Int, n::Int, m_::Int, n_::Int,
    k1_r::R, k1_i::R, k2_r::R, k2_i::R,
    r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}, ϕ_array::AbstractVecOrMat{R}, n̂_array::Any,
    kind::String, J_superscript::Int
) where {R <: Real}

    # TODO: do you think I should have done these multiplication outside, rather than repeating them whenever the function is called? or Julia should eliminate these duplicate calculation?
    kr1_r = k1_r .* r_array
    kr1_i = k1_i .* r_array
    kr2_r = k2_r .* r_array
    kr2_i = k2_i .* r_array

    # calculate the integrand
    if J_superscript == 11 # TODO: this if-statement can be done more nicely. We separate J_superscript into two pieces, the number 1 represents M_mn_wave_SeparateRealImag, while number 2 represents N_mn_wave_SeparateRealImag
        first_function = M_mn_wave_SeparateRealImag
        second_function = M_mn_wave_SeparateRealImag
    elseif J_superscript == 12
        first_function = M_mn_wave_SeparateRealImag
        second_function = N_mn_wave_SeparateRealImag
    elseif J_superscript == 21
        first_function = N_mn_wave_SeparateRealImag
        second_function = M_mn_wave_SeparateRealImag
    elseif J_superscript == 22
        first_function = N_mn_wave_SeparateRealImag
        second_function = N_mn_wave_SeparateRealImag
    else
        throw(DomainError("J_superscript must be any of [11,12,21,22]"))
    end

    kind_first_function = "regular"
    if kind == "irregular"
        kind_second_function = "irregular"
    elseif kind == "regular"
        kind_second_function = "regular"
    else
        throw(DomainError("""kind must be any of ["regular", "irregular"]"""))
    end

    # the cross product
    cross_product_MN = complex_vector_cross_product.(
        first_function.(m_, n_, kr2_r, kr2_i, θ_array, ϕ_array, kind_first_function),
        second_function.(-m, n, kr1_r, kr1_i, θ_array, ϕ_array, kind_second_function),
    )
    cross_product_MN_dot_n̂ = (-1).^m .* complex_vector_dot_product.(cross_product_MN, n̂_array)

    J_integrand = Tmatrix.surface_integrand(cross_product_MN_dot_n̂, r_array, θ_array)

    # This line used to have trouble
    # return vcat(J_integrand...) # I had to flatten all nested arrays.
    return hcat([i[1] for i in J_integrand], [i[2] for i in J_integrand])
end

function J_mn_m_n__integrand_SeparateRealImag_SMatrix(
        m::Int, n::Int, m_::Int, n_::Int,
        k1_r::R, k1_i::R, k2_r::R, k2_i::R,
        r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}, ϕ_array::AbstractVecOrMat{R}, n̂_array::Any,
        kind::String, J_superscript::Int
    ) where {R <: Real}

    # TODO: do you think I should have done these multiplication outside, rather than repeating them whenever the function is called? or Julia should eliminate these duplicate calculation?
    kr1_r = k1_r .* r_array
    kr1_i = k1_i .* r_array
    kr2_r = k2_r .* r_array
    kr2_i = k2_i .* r_array

    # calculate the integrand
    if J_superscript == 11 # TODO: this if-statement can be done more nicely. We separate J_superscript into two pieces, the number 1 represents M_mn_wave_SeparateRealImag, while number 2 represents N_mn_wave_SeparateRealImag
        first_function = M_mn_wave_SeparateRealImag_SMatrix
        second_function = M_mn_wave_SeparateRealImag_SMatrix
    elseif J_superscript == 12
        first_function = M_mn_wave_SeparateRealImag_SMatrix
        second_function = N_mn_wave_SeparateRealImag_SMatrix
    elseif J_superscript == 21
        first_function = N_mn_wave_SeparateRealImag_SMatrix
        second_function = M_mn_wave_SeparateRealImag_SMatrix
    elseif J_superscript == 22
        first_function = N_mn_wave_SeparateRealImag_SMatrix
        second_function = N_mn_wave_SeparateRealImag_SMatrix
    else
        throw(DomainError("J_superscript must be any of [11,12,21,22]"))
    end

    kind_first_function = "regular"
    if kind == "irregular"
        kind_second_function = "irregular"
    elseif kind == "regular"
        kind_second_function = "regular"
    else
        throw(DomainError("""kind must be any of ["regular", "irregular"]"""))
    end

    # the cross product
    cross_product_MN = complex_vector_cross_product.(
        first_function.(m_, n_, kr2_r, kr2_i, θ_array, ϕ_array, kind_first_function),
        second_function.(-m, n, kr1_r, kr1_i, θ_array, ϕ_array, kind_second_function),
    )
    cross_product_MN_dot_n̂ = (-1).^m .* complex_vector_dot_product.(cross_product_MN, n̂_array)

    J_integrand = Tmatrix.surface_integrand(cross_product_MN_dot_n̂, r_array, θ_array)

    # This line used to have trouble
    # return vcat(J_integrand...) # I had to flatten all nested arrays.
    return hcat([i[1] for i in J_integrand], [i[2] for i in J_integrand])
end

"""Computer surface integrals in equations 5.184 and 5.190
"""
function J_mn_m_n__SeparateRealImag(
        m::Int, n::Int, m_::Int, n_::Int,
        k1_r::R, k1_i::R, k2_r::R, k2_i::R,
        r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}, ϕ_array::AbstractVecOrMat{R}, n̂_array::Any, # TODO: I don't know why I get an error when I use n̂_array::AbstractVecOrMat{Vector{Float64}}
        kind::String, J_superscript::Int, rotationally_symmetric::Bool,
    ) where {R <: Real}

    if rotationally_symmetric
        # make sure that θ_array is 1D
        if length(size(θ_array)) != 1
            throw(DomainError("Since you have indicated << rotationally_symmetric = true >>, θ_array has to be 1D. Now it is $(length(size(θ_array)))D"))
        end
        ϕ_array = zeros(typeof(θ_array[1]), size(θ_array))
    end

    if rotationally_symmetric && (m != m_)
        # the integral over ϕ is 2π * δ_m_m_, so it is zero if m != m_
        J_integrand_dS_r = zeros(size(θ_array))
        J_integrand_dS_i = zeros(size(θ_array))
        return hcat(typeof(θ_array[1]), typeof(θ_array[1]))

    else
        J_integrand_dS = J_mn_m_n__integrand_SeparateRealImag(
            m, n, m_, n_,
            k1_r, k1_i, k2_r, k2_i,
            r_array, θ_array, ϕ_array, n̂_array,
            kind, J_superscript,
        )
        ## next three lines don't work
        # because J_mn_m_n__integrand_SeparateRealImag returns a flattened J, I have to reshape it to make it work for trapz
        J_integrand_dS_r = reshape(J_integrand_dS[:,1:Int(end/2)], size(θ_array))
        J_integrand_dS_i = reshape(J_integrand_dS[:,Int(end/2+1):end], size(θ_array))
    end

    # surface integral
    if rotationally_symmetric
        # integrate over θ only (i.e., 1D integration or line integral)
        J_r = 2π * trapz_ELZOUKA(θ_array, J_integrand_dS_r)
        J_i = 2π * trapz_ELZOUKA(θ_array, J_integrand_dS_i)

    else
        # integrate over θ and ϕ (i.e., 2D integration or surface integral)
        # TODO: replace this integral with surface mesh quadrature, like this one: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5095443/ , https://github.com/jareeger/Smooth_Closed_Surface_Quadrature_RBF-julia
        # assuming that θ_array, ϕ_array were created with meshgrid function

        # TODO: I need to define trapz_ELZOUKA for 2D integrals. Otherwise, autodiff will not work.
        #J_r = trapz((θ_array[:,1], ϕ_array[1,:]), J_integrand_dS_r)
        #J_i = trapz((θ_array[:,1], ϕ_array[1,:]), J_integrand_dS_i)

        # using the custom function trapz_ELZOUKA
        J_r = trapz_ELZOUKA(θ_array[:,1], ϕ_array[1,:], J_integrand_dS_r)
        J_i = trapz_ELZOUKA(θ_array[:,1], ϕ_array[1,:], J_integrand_dS_i)
    end

    return hcat(J_r, J_i)
end

function J_mn_m_n__SeparateRealImag_SMatrix(
        m::Int, n::Int, m_::Int, n_::Int,
        k1_r::R, k1_i::R, k2_r::R, k2_i::R,
        r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}, ϕ_array::AbstractVecOrMat{R}, n̂_array::Any, # TODO: I don't know why I get an error when I use n̂_array::AbstractVecOrMat{Vector{Float64}}
        kind::String, J_superscript::Int, rotationally_symmetric::Bool,
    ) where {R <: Real}

    if rotationally_symmetric
        # make sure that θ_array is 1D
        if length(size(θ_array)) != 1
            throw(DomainError("Since you have indicated << rotationally_symmetric = true >>, θ_array has to be 1D. Now it is $(length(size(θ_array)))D"))
        end
        ϕ_array = zeros(typeof(θ_array[1]), size(θ_array))
    end

    if rotationally_symmetric && (m != m_)
        # the integral over ϕ is 2π * δ_m_m_, so it is zero if m != m_
        J_integrand_dS_r = zeros(size(θ_array))
        J_integrand_dS_i = zeros(size(θ_array))
        return hcat(typeof(θ_array[1]), typeof(θ_array[1]))

    else
        J_integrand_dS = J_mn_m_n__integrand_SeparateRealImag_SMatrix(
            m, n, m_, n_,
            k1_r, k1_i, k2_r, k2_i,
            r_array, θ_array, ϕ_array, n̂_array,
            kind, J_superscript,
        )
        ## next three lines don't work
        # because J_mn_m_n__integrand_SeparateRealImag_SMatrix returns a flattened J, I have to reshape it to make it work for trapz
        J_integrand_dS_r = reshape(J_integrand_dS[:,1:Int(end / 2)], size(θ_array))
        J_integrand_dS_i = reshape(J_integrand_dS[:,Int(end / 2 + 1):end], size(θ_array))
    end

    # surface integral
    if rotationally_symmetric
        # integrate over θ only
        J_r = 2π * trapz_ELZOUKA(θ_array, J_integrand_dS_r)
        J_i = 2π * trapz_ELZOUKA(θ_array, J_integrand_dS_i)

    else
        # integrate over θ and ϕ
        # TODO: replace this integral with surface mesh quadrature, like this one: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5095443/ , https://github.com/jareeger/Smooth_Closed_Surface_Quadrature_RBF-julia
        # assuming that θ_array, ϕ_array were created with meshgrid function

        # TODO: I need to define trapz_ELZOUKA for 2D integrals. Otherwise, autodiff will not work.
        # J_r = trapz((θ_array[:,1], ϕ_array[1,:]), J_integrand_dS_r)
        # J_i = trapz((θ_array[:,1], ϕ_array[1,:]), J_integrand_dS_i)

        # using the custom function trapz_ELZOUKA
        J_r = trapz_ELZOUKA(θ_array[:,1], ϕ_array[1,:], J_integrand_dS_r)
        J_i = trapz_ELZOUKA(θ_array[:,1], ϕ_array[1,:], J_integrand_dS_i)
    end

    return hcat(J_r, J_i)
end

function Q_mn_m_n_SeparateRealImag(
        m::Int, n::Int, m_::Int, n_::Int,
        k1_r::R, k1_i::R, k2_r::R, k2_i::R,
        r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}, ϕ_array::AbstractVecOrMat{R}, n̂_array::Any, # TODO: how to define a multidimensional array type with Float64 elements?
        kind::String, Q_superscript::Int, rotationally_symmetric::Bool, symmetric_about_plane_perpendicular_z::Bool,
    )  where {R <: Real}

    if Q_superscript == 11; J_superscript_1 = 21 ; J_superscript_2 = 12
    elseif Q_superscript == 12; J_superscript_1 = 11 ; J_superscript_2 = 22
    elseif Q_superscript == 21; J_superscript_1 = 22 ; J_superscript_2 = 11
    elseif Q_superscript == 22; J_superscript_1 = 12 ; J_superscript_2 = 21
    end

    if rotationally_symmetric && (m != m_)
        return hcat(zero(θ_array[1]), zero(θ_array[1])) # TODO: find a better way to get the type of array elements

    # TODO: adding these lines causes Zygote to fail. Why? The code can work without them.
    # elseif rotationally_symmetric && symmetric_about_plane_perpendicular_z && (m == m_) && (Q_superscript == 11 || Q_superscript == 22) && ! iseven(n + n_)
    #    return hcat(0, 0)

    # elseif rotationally_symmetric && symmetric_about_plane_perpendicular_z && (m == m_) && (Q_superscript == 12 || Q_superscript == 21) && ! isodd(n + n_)
    #    return hcat(0, 0)

    else
        Q_r, Q_i = (
            - one(R) .* complex_multiply(
                [zero(R), one(R)]', complex_multiply(
                    [k1_r, k1_i]', complex_multiply(
                        [k2_r, k2_i]', J_mn_m_n__SeparateRealImag(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, kind, J_superscript_1, rotationally_symmetric)
                    )
                )
            )

            - one(R) .* complex_multiply(
                [zero(R), one(R)]', complex_multiply(
                    [k1_r, k1_i]', complex_multiply(
                        [k1_r, k1_i]', J_mn_m_n__SeparateRealImag(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, kind, J_superscript_2, rotationally_symmetric)
                    )
                )
            )
        )
        return hcat(Q_r, Q_i)
    end
end

function Q_mn_m_n_SeparateRealImag_SMatrix(
    m::Int, n::Int, m_::Int, n_::Int,
    k1_r::R, k1_i::R, k2_r::R, k2_i::R,
    r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}, ϕ_array::AbstractVecOrMat{R}, n̂_array::Any, # TODO: how to define a multidimensional array type with Float64 elements?
    kind::String, Q_superscript::Int, rotationally_symmetric::Bool, symmetric_about_plane_perpendicular_z::Bool,
)  where {R <: Real}

    if Q_superscript == 11; J_superscript_1 = 21 ; J_superscript_2 = 12
    elseif Q_superscript == 12; J_superscript_1 = 11 ; J_superscript_2 = 22
elseif Q_superscript == 21; J_superscript_1 = 22 ; J_superscript_2 = 11
elseif Q_superscript == 22; J_superscript_1 = 12 ; J_superscript_2 = 21
    end

    if rotationally_symmetric && (m != m_)
        return hcat(zero(θ_array[1]), zero(θ_array[1])) # TODO: find a better way to get the type of array elements

# TODO: adding these lines causes Zygote to fail. Why? The code can work without them.
# elseif rotationally_symmetric && symmetric_about_plane_perpendicular_z && (m == m_) && (Q_superscript == 11 || Q_superscript == 22) && ! iseven(n + n_)
#    return hcat(0, 0)

# elseif rotationally_symmetric && symmetric_about_plane_perpendicular_z && (m == m_) && (Q_superscript == 12 || Q_superscript == 21) && ! isodd(n + n_)
#    return hcat(0, 0)

    else
        Q_r, Q_i = (
        - one(R) .* complex_multiply_SMatrix(
            SMatrix{1,2}(zero(R), one(R)), complex_multiply_SMatrix(
                SMatrix{1,2}(k1_r, k1_i), complex_multiply_SMatrix(
                    SMatrix{1,2}(k2_r, k2_i), J_mn_m_n__SeparateRealImag(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, kind, J_superscript_1, rotationally_symmetric)
                )
            )
        )

        - one(R) .* complex_multiply_SMatrix(
            SMatrix{1,2}(zero(R), one(R)), complex_multiply_SMatrix(
                SMatrix{1,2}(k1_r, k1_i), complex_multiply_SMatrix(
                    SMatrix{1,2}(k1_r, k1_i), J_mn_m_n__SeparateRealImag(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, kind, J_superscript_2, rotationally_symmetric)
                )
            )
        )
    )
        return hcat(Q_r, Q_i)
    end
end

"""Calculate the entire regular and irregular Q matrix.

Parameters
----------
kind : can be either "regular" or "irregular"

returns
-------
Q_matrix : hcat of real and imaginary parts of Q matrix
"""
function Q_matrix_SeparateRealImag(
        n_max::Int, k1_r::R, k1_i::R, k2_r::R, k2_i::R,
        r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}, ϕ_array::AbstractVecOrMat{R}, n̂_array::Any,
        kind::String, rotationally_symmetric::Bool, symmetric_about_plane_perpendicular_z::Bool,
    ) where {R <: Real}

    m, n, m_, n_ = get_m_n_m__n__matrices_for_T_matrix(n_max)
    Q_mn_m_n_11 = Q_mn_m_n_SeparateRealImag.(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, [r_array], [θ_array], [ϕ_array], [n̂_array], kind, 11, rotationally_symmetric, symmetric_about_plane_perpendicular_z)
    Q_mn_m_n_12 = Q_mn_m_n_SeparateRealImag.(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, [r_array], [θ_array], [ϕ_array], [n̂_array], kind, 12, rotationally_symmetric, symmetric_about_plane_perpendicular_z)
    Q_mn_m_n_21 = Q_mn_m_n_SeparateRealImag.(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, [r_array], [θ_array], [ϕ_array], [n̂_array], kind, 21, rotationally_symmetric, symmetric_about_plane_perpendicular_z)
    Q_mn_m_n_22 = Q_mn_m_n_SeparateRealImag.(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, [r_array], [θ_array], [ϕ_array], [n̂_array], kind, 22, rotationally_symmetric, symmetric_about_plane_perpendicular_z)

    Q = vcat(
        (hcat(Q_mn_m_n_11, Q_mn_m_n_12)),
        (hcat(Q_mn_m_n_21, Q_mn_m_n_22))
    )
    get_real_of_nested(x) = x[1] # TODO: is there a better way to separate the first element of a nested array?
    get_imag_of_nested(x) = x[2] # TODO: is there a better way to separate the first element of a nested array?
    Q_r = get_real_of_nested.(Q) # TODO: is there a better way to separate the first element of a nested array?
    Q_i = get_imag_of_nested.(Q) # TODO: is there a better way to separate the first element of a nested array?

    # TODO: is it better to return a 2D square matrix, each element is a 2-vector, representing real and imag parts? This will not work with Zygote
    # TODO: or concatenate 2 2D square arrays, the first and second represent the real and imaginary parts?
    return hcat(Q_r, Q_i)
end
function Q_matrix_SeparateRealImag_SMatrix(
        n_max::Int, k1_r::R, k1_i::R, k2_r::R, k2_i::R,
        r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}, ϕ_array::AbstractVecOrMat{R}, n̂_array::Any,
        kind::String, rotationally_symmetric::Bool, symmetric_about_plane_perpendicular_z::Bool,
    ) where {R <: Real}

    m, n, m_, n_ = get_m_n_m__n__matrices_for_T_matrix(n_max)
    Q_mn_m_n_11 = Q_mn_m_n_SeparateRealImag_SMatrix.(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, [r_array], [θ_array], [ϕ_array], [n̂_array], kind, 11, rotationally_symmetric, symmetric_about_plane_perpendicular_z)
    Q_mn_m_n_12 = Q_mn_m_n_SeparateRealImag_SMatrix.(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, [r_array], [θ_array], [ϕ_array], [n̂_array], kind, 12, rotationally_symmetric, symmetric_about_plane_perpendicular_z)
    Q_mn_m_n_21 = Q_mn_m_n_SeparateRealImag_SMatrix.(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, [r_array], [θ_array], [ϕ_array], [n̂_array], kind, 21, rotationally_symmetric, symmetric_about_plane_perpendicular_z)
    Q_mn_m_n_22 = Q_mn_m_n_SeparateRealImag_SMatrix.(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, [r_array], [θ_array], [ϕ_array], [n̂_array], kind, 22, rotationally_symmetric, symmetric_about_plane_perpendicular_z)

    Q = vcat(
        (hcat(Q_mn_m_n_11, Q_mn_m_n_12)),
        (hcat(Q_mn_m_n_21, Q_mn_m_n_22))
    )
    get_real_of_nested(x) = x[1] # TODO: is there a better way to separate the first element of a nested array?
    get_imag_of_nested(x) = x[2] # TODO: is there a better way to separate the first element of a nested array?
    Q_r = get_real_of_nested.(Q) # TODO: is there a better way to separate the first element of a nested array?
    Q_i = get_imag_of_nested.(Q) # TODO: is there a better way to separate the first element of a nested array?

    # TODO: is it better to return a 2D square matrix, each element is a 2-vector, representing real and imag parts? This will not work with Zygote
    # TODO: or concatenate 2 2D square arrays, the first and second represent the real and imaginary parts?
    return hcat(Q_r, Q_i)
end

"""Main function that calculates T-matrix

"""
function T_matrix_SeparateRealImag(
        n_max::Int, k1_r::R, k1_i::R, k2_r::R, k2_i::R,
        r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}, ϕ_array::AbstractVecOrMat{R}, n̂_array::Any,
        rotationally_symmetric, symmetric_about_plane_perpendicular_z
    ) where {R <: Real}

    # calculating Regular Q and Irregular Q
    RgQ = Q_matrix_SeparateRealImag(n_max, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, "regular", rotationally_symmetric, symmetric_about_plane_perpendicular_z)
    Q   = Q_matrix_SeparateRealImag(n_max, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, "irregular", rotationally_symmetric, symmetric_about_plane_perpendicular_z)

    #Q = Float64.(Q) # TODO: I need to make the concrete type programmable, in case I need Float128
    #RgQ = Float64.(RgQ) # TODO: I need to make the concrete type programmable, in case I need Float128

    # separating real and imaginary parts of Regular Q
    RgQ_r = RgQ[:, 1:Int(end / 2)]
    RgQ_i = RgQ[:, Int(end / 2) + 1:end]

    # separating real and imaginary parts of irregular Q
    Q_r = Q[:, 1:Int(end / 2)]
    Q_i = Q[:, Int(end / 2) + 1:end]

    # calculating the inverse of irregular Q
    invQ = complex_matrix_inversion(Q_r, Q_i)

    # calculating the T-matrix
    T = -1 .* complex_matrix_multiplication(RgQ_r, RgQ_i, invQ[1], invQ[2]) # equation 5.191
    return hcat(T[1], T[2]) # T[1] and T[2] are the real and imaginary parts of T-matrix, respectively.
end

"""
    allowing BigFloat_precision
"""
function T_matrix_SeparateRealImag(
        n_max::Int, k1_r::R, k1_i::R, k2_r::R, k2_i::R,
        r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}, ϕ_array::AbstractVecOrMat{R}, n̂_array::Any,
        rotationally_symmetric::Bool, symmetric_about_plane_perpendicular_z::Bool, BigFloat_precision
    ) where {R <: Real}
    if BigFloat_precision != nothing
        return setprecision(BigFloat_precision) do
            return T_matrix_SeparateRealImag(
                n_max, big(k1_r), big(k1_i), big(k2_r), big(k2_i),
                big.(r_array), big.(θ_array), big.(ϕ_array), [big.(n) for n in n̂_array],
                rotationally_symmetric, symmetric_about_plane_perpendicular_z
            )
        end

    else
        return T_matrix_SeparateRealImag(
            n_max, k1_r, k1_i, k2_r, k2_i,
            r_array, θ_array, ϕ_array, n̂_array,
            rotationally_symmetric, symmetric_about_plane_perpendicular_z
        )
    end
end

"""
    This can allow input of arbitrary mesh points
"""
function T_matrix_SeparateRealImag_arbitrary_mesh(
        n_max::Int, k1_r::R, k1_i::R, k2_r::R, k2_i::R,
        r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}, ϕ_array::AbstractVecOrMat{R},
        rotationally_symmetric::Bool, symmetric_about_plane_perpendicular_z::Bool, BigFloat_precision
    ) where {R <: Real}
    if rotationally_symmetric
        r_theta_n̂ = get_r_θ_n̂_arrays_from_r_θ_arrays_axisymmetrix(r_array, θ_array)
        r_array = r_theta_n̂[:,1]
        θ_array = r_theta_n̂[:,2]
        ϕ_array = zeros(size(θ_array))
        n̂_r_comp = r_theta_n̂[:,3]
        n̂_θ_comp = r_theta_n̂[:,4]
        n̂_ϕ_comp = r_theta_n̂[:,5]
        n̂_array = reshape([Vector([n̂_r_comp[id], n̂_θ_comp[id], n̂_ϕ_comp[id]]) for id in eachindex(n̂_r_comp)], size(n̂_r_comp))
        # println("rotationally symmetric T-matrix calculation for arbitrary mesh is starting ...")
    else
        println("I need to code it!")
    end

    return T_matrix_SeparateRealImag(
        n_max, k1_r, k1_i, k2_r, k2_i,
        r_array, θ_array, ϕ_array, n̂_array,
        rotationally_symmetric, symmetric_about_plane_perpendicular_z, BigFloat_precision
    )
end

"""
    allowing wavelength/frequency and material properties, rather than wavevectors
"""
function T_matrix_SeparateRealImag_arbitrary_mesh(
    n_max::Int, wl_or_freq_input::R, input_unit::String, Eps_r_r_1::R, Eps_r_i_1::R, Mu_r_r_1::R, Mu_r_i_1::R, Eps_r_r_2::R, Eps_r_i_2::R, Mu_r_r_2::R, Mu_r_i_2::R,
    r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}, ϕ_array::AbstractVecOrMat{R},
    rotationally_symmetric::Bool, symmetric_about_plane_perpendicular_z::Bool, BigFloat_precision
) where {R <: Real}
    k1_complex = get_WaveVector(
            wl_or_freq_input;
            input_unit=input_unit,
            Eps_r = Complex(Eps_r_r_1, Eps_r_i_1),
            Mu_r = Complex(Mu_r_r_1, Mu_r_i_1)
    )
    k2_complex = get_WaveVector(
        wl_or_freq_input;
        input_unit=input_unit,
        Eps_r = Complex(Eps_r_r_2, Eps_r_i_2),
        Mu_r = Complex(Mu_r_r_2, Mu_r_i_2)
    )
    return T_matrix_SeparateRealImag_arbitrary_mesh(
        n_max, real(k1_complex), imag(k1_complex), real(k2_complex), imag(k2_complex),
        r_array, θ_array, ϕ_array,
        rotationally_symmetric, symmetric_about_plane_perpendicular_z, BigFloat_precision
    )
end

"""
    allowing wavelength/frequency and material properties, rather than wavevectors
"""
function T_matrix_SeparateRealImag(
        n_max::Int, wl_or_freq_input::R, input_unit::String, Eps_r_r_1::R, Eps_r_i_1::R, Mu_r_r_1::R, Mu_r_i_1::R, Eps_r_r_2::R, Eps_r_i_2::R, Mu_r_r_2::R, Mu_r_i_2::R,
        r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}, ϕ_array::AbstractVecOrMat{R}, n̂_array::Any,
        rotationally_symmetric, symmetric_about_plane_perpendicular_z, BigFloat_precision
    ) where {R <: Real}
    k1_complex = get_WaveVector(
        wl_or_freq_input;
        input_unit=input_unit,
        Eps_r = Complex(Eps_r_r_1, Eps_r_i_1),
        Mu_r = Complex(Mu_r_r_1, Mu_r_i_1)
    )
    k2_complex = get_WaveVector(
        wl_or_freq_input;
        input_unit=input_unit,
        Eps_r = Complex(Eps_r_r_2, Eps_r_i_2),
        Mu_r = Complex(Mu_r_r_2, Mu_r_i_2)
    )

    return T_matrix_SeparateRealImag(
        n_max, real(k1_complex), imag(k1_complex), real(k2_complex), imag(k2_complex),
        r_array, θ_array, ϕ_array, n̂_array,
        rotationally_symmetric, symmetric_about_plane_perpendicular_z, BigFloat_precision
    )
end

"""
HDF5_filename : string, the path to HDF5 file where Tmatrix is going to be saved.
"""
function calculate_Tmatrix_for_spheroid_SeparateRealImag(
        rx::R, rz::R, n_max::Int,
        k1_r::R, k1_i::R, k2_r::R, k2_i::R;
        n_θ_points=10, n_ϕ_points=20, HDF5_filename=nothing,
        rotationally_symmetric=false, symmetric_about_plane_perpendicular_z=false,
        BigFloat_precision=nothing
    ) where {R <: Real}

    # create a grid of θ_ϕ
    θ_array, ϕ_array = meshgrid_θ_ϕ(n_θ_points, n_ϕ_points; min_θ=1e-16, min_ϕ=1e-16, rotationally_symmetric=rotationally_symmetric)

    # calculate r and n̂ for the geometry
    r_array, n̂_array = ellipsoid(rx, rz, θ_array);

    # calculate T-matrix
    T = T_matrix_SeparateRealImag(
        n_max, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array,
        rotationally_symmetric, symmetric_about_plane_perpendicular_z, BigFloat_precision
    )
    return T
end


text:
# wavread test

using WAV

Fs = 48000;

function prepare_input(input_size)
	t = (0:input_size-1)/Fs;
	y = 0.7*sin.(2*pi*50*t) + sin.(2*pi*120*t) + 2*randn(size(t));
	filename = "/home/shyman/tmp/audiofile.wav";
	wavwrite(y/maximum(abs.(y)), filename, Fs=Fs);
	return filename
end

function test_function(input_data)
	y, Fs = wavread(input_data);
	return y
end


text:
using CmdStan
src = "
data {
    real<lower=0, upper=1> theta;
    int<lower=0> K;
    int<lower=0> N;
}
model {

}
generated quantities {
    int<lower=0, upper=K> y[N];
    for (n in 1:N)
        y[n] = binomial_rng(K, theta);
}
"
model = Stanmodel(method = Sample(algorithm=fixed_param), model = src)
data = Dict("theta" => 0.5, "K" => 100, "N" => 50)
stan(model, data)

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
using Treesitter

parser = Parser()
set_language(parser, tree_sitter_c())

tree = Treesitter.parse_filewith(
    parser,
    joinpath(@__DIR__, "corpus/easystruct.h"),
)
nstring = parser.buffer |> String
rtcrsr = tree_getroot(tree)

rchildren = node_namedchildren(rtcrsr)

rcstr = node_type.(rchildren)

h = node_text.(Ref(nstring), rchildren)

decl = rcstr .=> h

R1 = node_namedchildren.(rchildren)

𝑅 = [node_type.(r) for r in R1]
𝐻 = [node_text.(Ref(nstring), r) for r in R1]

𝑄 = 𝑅 .=> 𝐻


text:
import Base: *

toarray(xs::AbstractArray, ys::AbstractArray) = ys
toarray(xs::AbstractArray, y) = similar(xs, typeof(y), ()) .= y

unarray(xs) = xs
unarray(xs::AbstractArray{T,0} where T) = xs[]

Base.getindex(xs::TrackedArray, i...) =
  TrackedArray(Call(getindex, xs, i...), toarray(xs.data, xs.data[i...]))

function back(::typeof(getindex), Δ, xs::TrackedArray, i...)
  Δ′ = zeros(xs.data)
  Δ′[i...] = unarray(Δ)
  @back(xs, Δ′)
end

Base.:-(xs::TrackedArray) = TrackedArray(Call(-, xs))

back(::typeof(-), Δ, xs::TrackedArray) = back(xs, -Δ)

Base.transpose(xs::TrackedArray) = TrackedArray(Call(transpose, xs))
Base.ctranspose(xs::TrackedArray) = TrackedArray(Call(ctranspose, xs))

back(::typeof(transpose), Δ, xs) = @back(xs, trim(xs, Δ.'))
back(::typeof(ctranspose), Δ, xs) = @back(xs, trim(xs, Δ'))

Base.repmat(x::TrackedVecOrMat, a::Integer...) = TrackedArray(Call(repmat, x, a...))
Base.repmat(x::TrackedVecOrMat, a::Int64...) = TrackedArray(Call(repmat, x, a...))

Base.vcat(a::TrackedVector, b::TrackedVector)  = TrackedArray(Call(vcat, a, b))
Base.vcat(a::TrackedVector, b::AbstractVector) = TrackedArray(Call(vcat, a, b))
Base.vcat(a::AbstractVector, b::TrackedVector) = TrackedArray(Call(vcat, a, b))

Base.vcat(a::TrackedVecOrMat, b::TrackedVecOrMat)  = TrackedArray(Call(vcat, a, b))
Base.vcat(a::TrackedVecOrMat, b::AbstractVecOrMat) = TrackedArray(Call(vcat, a, b))
Base.vcat(a::AbstractVecOrMat, b::TrackedVecOrMat) = TrackedArray(Call(vcat, a, b))

Base.vcat(a::TrackedMatrix, b::TrackedMatrix)  = TrackedArray(Call(vcat, a, b))
Base.vcat(a::TrackedMatrix, b::AbstractMatrix) = TrackedArray(Call(vcat, a, b))
Base.vcat(a::AbstractMatrix, b::TrackedMatrix) = TrackedArray(Call(vcat, a, b))

function back(::typeof(vcat), Δ, xs, ys)
  i = Base.tail(map(_ -> :, size(Δ)))
  @back(xs, Δ[1:size(xs,1), i...])
  @back(ys, Δ[size(xs,1)+1:end, i...])
end

# Reductions

Base.sum(xs::TrackedArray, dim) = TrackedArray(Call(sum, xs, dim))
Base.sum(xs::TrackedArray) = TrackedArray(Call(sum, xs), toarray(xs.data, sum(xs.data)))
Base.sum(xs::TrackedScalar, dim...) = xs

back(::typeof(sum), Δ, xs::TrackedArray, dim...) = back(xs, similar(xs.data) .= Δ)

Base.maximum(xs::TrackedArray, args...) = maximum(xs.data, args...)
Base.findfirst(xs::TrackedArray, args...) = findfirst(xs.data, args...)

Base.mean(xs::TrackedArray) = TrackedArray(Call(mean, xs), toarray(xs.data, mean(xs.data)))
Base.mean(xs::TrackedArray, region) = TrackedArray(Call(mean, xs, region))

back(::typeof(mean), Δ, xs::TrackedArray) = back(xs, similar(xs.data) .= Δ ./ length(xs.data))
back(::typeof(mean), Δ, xs::TrackedArray, region) =
  back(xs, similar(xs.data) .= Δ ./ prod(size(xs.data, region...)))

# BLAS

a::TrackedMatrix * b::TrackedMatrix  = TrackedArray(Call(*, a, b))
a::TrackedMatrix * b::AbstractMatrix = TrackedArray(Call(*, a, b))
a::AbstractMatrix * b::TrackedMatrix = TrackedArray(Call(*, a, b))

a::TrackedMatrix * b::TrackedVector  = TrackedArray(Call(*, a, b))
a::TrackedMatrix * b::AbstractVector = TrackedArray(Call(*, a, b))
a::AbstractMatrix * b::TrackedVector = TrackedArray(Call(*, a, b))

function back(::typeof(*), Δ, a::AbstractMatrix, b::AbstractVecOrMat)
  @back(a, A_mul_Bt(Δ, data(b)))
  @back(b, At_mul_B(data(a), Δ))
end

# NNlib

import NNlib: softmax, ∇softmax

softmax(xs::TrackedArray) = TrackedArray(Call(softmax, xs))

back(::typeof(softmax), Δ, xs) = @back(xs, ∇softmax(Δ, data(xs)))

# Broadcasting

using ForwardDiff: Dual, partials

struct Broadcasted{T}
  data::T
end

(b::Broadcasted)(xs...) = map(x -> x.value, b.data)

dualify(xs, n) = xs
dualify(xs::TrackedArray, ps) = map(x -> Dual(x, ps), data(xs))

function tracked_broadcast(f, args::Vararg{Any,N}) where N
  dargs = map((x,i) -> dualify(x, ntuple(j -> i==j, Val{N})), args, ntuple(identity, Val{N}))
  # TrackedArray(Call(Broadcasted(broadcast(f, dargs...)), args...))
  # Works around a 0.6 type inference issue
  b = Broadcasted(broadcast(f, dargs...))
  TrackedArray(Call(b, args...), b())
end

trim(x, Δ) = reshape(Δ, ntuple(i -> size(Δ, i), Val{ndims(x)}))

unbroadcast(x, Δ) =
  size(x) == size(Δ) ? Δ :
    trim(x, sum(Δ, filter(n -> size(x, n) == 1, 1:ndims(Δ))))

function getpartial(Δ, x, i)
  @inbounds p = getindex(partials(x), i)
  return Δ * p
end

function back(b::Broadcasted, Δ, args::Vararg{Any,N}) where N
  Δargs = ntuple(i -> getpartial.(Δ, b.data, i), Val{N})
  foreach((x, Δ) -> @back(x, unbroadcast(x, Δ)), args, Δargs)
end

Base.Broadcast._containertype(::Type{<:TrackedArray}) = TrackedArray
Base.Broadcast.promote_containertype(::Type{TrackedArray}, ::Type{TrackedArray}) = TrackedArray
Base.Broadcast.promote_containertype(::Type{Array}, ::Type{TrackedArray}) = TrackedArray
Base.Broadcast.promote_containertype(::Type{TrackedArray}, ::Type{Array}) = TrackedArray
Base.Broadcast.promote_containertype(::Type{TrackedArray}, ct) = TrackedArray
Base.Broadcast.promote_containertype(ct, ::Type{TrackedArray}) = TrackedArray
Base.Broadcast.broadcast_indices(::Type{TrackedArray}, A::Ref) = ()
Base.Broadcast.broadcast_indices(::Type{TrackedArray}, A) = indices(A)

Base.Broadcast.broadcast_c(f, ::Type{TrackedArray}, A, Bs...) = tracked_broadcast(f, A, Bs...)


text:
# Current sheet model.

"""
    getB_CS_harris(B₀, L)

Return the magnetic field at location `r` near a current sheet with magnetic strength `B₀`
and sheet length `L`.
"""
function getB_CS_harris(r, B₀, L)
   x = r[1]
   [0.0, 0.0, B₀ * tanh(x/L)]
end