text:
using TuringPatterns
using Meshes
using GeoStatsBase
using Plots; gr(size=(600,400))
using ReferenceTests, ImageIO
using Test, Random

# workaround GR warnings
ENV["GKSwstype"] = "100"

# environment settings
isCI = "CI" ∈ keys(ENV)
islinux = Sys.islinux()
visualtests = !isCI || (isCI && islinux)
datadir = joinpath(@__DIR__,"data")

@testset "TuringPatterns.jl" begin
  @testset "Basic usage" begin
    # TODO:
  end

  @testset "GeoStats.jl API" begin
    Random.seed!(2019)
    problem = SimulationProblem(CartesianGrid(200,200), :z => Float64, 3)
    solution = solve(problem, TPS())

    if visualtests
      @test_reference "data/GeoStatsAPI.png" plot(solution,size=(900,300))
    end
  end
end


text:
using DifferentialEquations: ODEProblem, DifferentialEquations

@inline function __integrate(u::Function, domain::Tuple{T,T}; kwargs...) where T <: AbstractFloat
	prob = ODEProblem( (f, p, τ) -> u(τ), zero(T), domain)
	return DifferentialEquations.solve(prob; kwargs...)
end

@inline function __getNorm(u::Function, c::Compression; kwargs...)
	σ, μ = c.σ, c.μ
	sol = __integrate(
		τ -> abs(τ) ≈ 1.0 ? 0.0 : σ * abs2( u(σ * c.decomp(τ) + μ) ) * c.inv_jacobian(τ),
		(-1.0, 1.0);
		kwargs...
	)
	return t -> abs(sol(c.comp((t - μ) / σ))), sol(1.0)
end

"""
    generateCouplings(u; compression, kwargs...)

inputs
- u:            mode function
- compression:  which compression of ℝ → (-1, 1) will be used for the differential equation
								Options are the keys of the `Compressions` dictionary
								or any `Compression` Object
- kwargs:				all other kwargs get's passed to the numerical integrator

output:
- (gᵢ, gₒ, norm): the coupling constant calculated from the mode function and its norm function.
"""
# TODO test with (df, f, p, τ) lambda expression

function generateCoupling(u::Function; compression::Union{Symbol, Compression} = :exponential, kwargs...)
	if isa(compression, Symbol)
		@assert compression ∈ keys(Compressions)

		return _generateCoupling(u, Compressions[compression]; kwargs...)
	end

	return _generateCoupling(u, compression; kwargs...)
end

function _generateCoupling(u::Function, c::Compression; kwargs...)
	norm, Norm = __getNorm(u, c; kwargs...)
	gᵢ(t) =   u(t) / √(abs(Norm - norm(t)) + (u(t) == 0.0))
	gₒ(t) = - u(t) / √(norm(t) + (u(t) == 0.0))

	return gᵢ, gₒ, norm
end

# """
#     modeContribution(u, v; compression)
#
# inputs
# - u:            mode function of the cavity
# - v:						mode function falsely contributing to the cavity
# - compression:  which compression of ℝ → (-1, 1) will be used for the differential equation
# 								Options are the keys of the `Compressions` dictionary
# 								or any `Compression` Object
#
# output:
# - α:   the contribution of mode u to the cavity
# """
# # TODO test with (df, f, p, τ) lambda expression
# function modeContribution(u::Mode, v::Mode; compression::Union{Symbol, Compression} = :exponential)
# 	if isa(compression, Symbol)
# 		@assert compression ∈ keys(Compressions)
#
# 		return _modeContribution(u, v, Compressions[compression])
# 	end
#
# 	return _modeContribution(u, v, compression)
# end
#
# function _modeContribution(uM::Mode, vM::Mode, c::Compression)
# 	u, v = uM.modeFunction, vM.modeFunction
# 	dfdτ(f, p, τ) = begin
# 		abs(τ) == 1.0 && return 0.0
#
# 		t = c.decomp(τ)
# 		df = u(t) * v(t) * c.inv_jacobian(τ)
# 	end
#
# 	prob = ODEProblem(dfdτ, 0., (-1.0, 1.0))
#   sol = solve(prob, reltol = 1e-5, abstol = 1e-8)
#
# 	return t -> sol(c.comp(t)) / √(uM.norm(t))
# end
#
# """
#     modifiedMode(u, v; compression)
#
# inputs
# - u:            `Mode` of the first cavity
# - v:						target `Mode` for the second cavity
# - [α]:					predetermined modeContribution of v to the first cavity. Optional.
# 								Useful, if α should be calculated with a different compression than the
# 								modified mode.
# - compression:  which compression of ℝ → (-1, 1) will be used for the differential equation
# 								Options are the keys of the `Compressions` dictionary
# 								or any `Compression` Object
#
# output:
# - ṽ:   modified `Mode` for the second cavity, accounting for the scattering at the first cavity
# """
# # TODO test with (df, f, p, τ) lambda expression
# function modifiedMode(u::Mode, v::Mode, α::Union{Function, Nothing} = nothing;
# 		compression::Union{Symbol, Compression} = :exponential
# 	)
#
# 	if isa(compression, Symbol)
# 		@assert compression ∈ keys(Compressions)
# 		compression = Compressions[compression]
# 	end
#
# 	isnothing(α) && (α = _modeContribution(u, v, compression))
# 	return _modifiedMode(u, v, α, compression)
# end
#
# function _modifiedMode(u::Mode, v::Mode, α::Function, compression::Compression)
# 	vMod(t) = v.modeFunction(t) + u.gₒ(t) * α(t)
#
# 	return Mode(vMod, compression = compression)
# end

ext: jl


text:
# Helper function to get reduced variable info
function _reduced_info(vref::ReducedInfiniteVariableRef)::ReducedInfiniteInfo
    return JuMP.owner_model(vref).reduced_info[JuMP.index(vref)]
end

"""
    infinite_variable_ref(vref::ReducedInfiniteVariableRef)::InfiniteVariableRef

Return the `InfiniteVariableRef` associated with the reduced infinite variable
`vref`.

**Example**
```julia-repl
julia> infinite_variable_ref(vref)
g(t, x)
```
"""
function infinite_variable_ref(vref::ReducedInfiniteVariableRef)::InfiniteVariableRef
    return _reduced_info(vref).infinite_variable_ref
end

"""
    eval_supports(vref::ReducedInfiniteVariableRef)::Dict

Return the evaluation supports associated with the reduced infinite variable
`vref`.

**Example**
```julia-repl
julia> eval_supports(vref)
Dict{Int64,Float64} with 1 entry:
  1 => 0.5
```
"""
function eval_supports(vref::ReducedInfiniteVariableRef)::Dict
    return _reduced_info(vref).eval_supports
end

"""
    parameter_refs(vref::ReducedInfiniteVariableRef)::Tuple

Return the `ParameterRef`(s) associated with the reduced infinite variable
`vref`. This is formatted as a Tuple of containing the parameter references as
they were inputted to define the untracripted infinite variable except, the
evaluated parameters are excluded.

**Example**
```julia-repl
julia> parameter_refs(vref)
(t,   [2]  =  x[2]
  [1]  =  x[1])
```
"""
function parameter_refs(vref::ReducedInfiniteVariableRef)
    orig_prefs = parameter_refs(infinite_variable_ref(vref))
    prefs = Tuple(orig_prefs[i] for i = 1:length(orig_prefs) if !haskey(eval_supports(vref), i))
    return prefs
end

"""
    JuMP.name(vref::ReducedInfiniteVariableRef)::String

Extend `JuMP.name` to return name of reduced infinite variable references. This
is used when displaying measure expansions that contain such variables.

**Exanple**
```julia-repl
julia> name(rvref)
g(1.25, x)
```
"""
function JuMP.name(vref::ReducedInfiniteVariableRef)::String
    root_name = _root_name(infinite_variable_ref(vref))
    prefs = parameter_refs(infinite_variable_ref(vref))
    param_names = [_root_name(first(pref)) for pref in prefs]
    for (k, v) in eval_supports(vref)
        param_names[k] = string(v)
    end
    param_name_tuple = "("
    for i = 1:length(param_names)
        if i != length(param_names)
            param_name_tuple *= string(param_names[i], ", ")
        else
            param_name_tuple *= string(param_names[i])
        end
    end
    param_name_tuple *= ")"
    return string(root_name, param_name_tuple)
end

"""
    JuMP.has_lower_bound(vref::ReducedInfiniteVariableRef)::Bool

Extend [`JuMP.has_lower_bound`](@ref) to return a `Bool` whether the original
infinite variable of `vref` has a lower bound.

**Example**
```julia-repl
julia> has_lower_bound(vref)
true
```
"""
function JuMP.has_lower_bound(vref::ReducedInfiniteVariableRef)::Bool
    return JuMP.has_lower_bound(infinite_variable_ref(vref))
end

"""
    JuMP.lower_bound(vref::ReducedInfiniteVariableRef)::Float64

Extend [`JuMP.lower_bound`](@ref) to return the lower bound of the original
infinite variable of `vref`. Errors if `vref` doesn't have a lower bound.

**Example**
```julia-repl
julia> lower_bound(vref)
0.0
```
"""
function JuMP.lower_bound(vref::ReducedInfiniteVariableRef)::Float64
    if !JuMP.has_lower_bound(vref)
        error("Variable $(vref) does not have a lower bound.")
    end
    return JuMP.lower_bound(infinite_variable_ref(vref))
end

# Extend to return the index of the lower bound constraint associated with the
# original infinite variable of `vref`.
function JuMP._lower_bound_index(vref::ReducedInfiniteVariableRef)::Int
    if !JuMP.has_lower_bound(vref)
        error("Variable $(vref) does not have a lower bound.")
    end
    return JuMP._lower_bound_index(infinite_variable_ref(vref))
end

"""
    JuMP.LowerBoundRef(vref::ReducedInfiniteVariableRef)::GeneralConstraintRef

Extend [`JuMP.LowerBoundRef`](@ref) to extract a constraint reference for the
lower bound of the original infinite variable of `vref`.

**Example**
```julia-repl
julia> cref = LowerBoundRef(vref)
var >= 0.0
```
"""
function JuMP.LowerBoundRef(vref::ReducedInfiniteVariableRef)::GeneralConstraintRef
    return JuMP.LowerBoundRef(infinite_variable_ref(vref))
end

"""
    JuMP.has_upper_bound(vref::ReducedInfiniteVariableRef)::Bool

Extend [`JuMP.has_upper_bound`](@ref) to return a `Bool` whether the original
infinite variable of `vref` has an upper bound.

**Example**
```julia-repl
julia> has_upper_bound(vref)
true
```
"""
function JuMP.has_upper_bound(vref::ReducedInfiniteVariableRef)::Bool
    return JuMP.has_upper_bound(infinite_variable_ref(vref))
end

"""
    JuMP.upper_bound(vref::ReducedInfiniteVariableRef)::Float64

Extend [`JuMP.upper_bound`](@ref) to return the upper bound of the original
infinite variable of `vref`. Errors if `vref` doesn't have a upper bound.

**Example**
```julia-repl
julia> upper_bound(vref)
0.0
```
"""
function JuMP.upper_bound(vref::ReducedInfiniteVariableRef)::Float64
    if !JuMP.has_upper_bound(vref)
        error("Variable $(vref) does not have a upper bound.")
    end
    return JuMP.upper_bound(infinite_variable_ref(vref))
end

# Extend to return the index of the upper bound constraint associated with the
# original infinite variable of `vref`.
function JuMP._upper_bound_index(vref::ReducedInfiniteVariableRef)::Int
    if !JuMP.has_upper_bound(vref)
        error("Variable $(vref) does not have a upper bound.")
    end
    return JuMP._upper_bound_index(infinite_variable_ref(vref))
end

"""
    JuMP.UpperBoundRef(vref::ReducedInfiniteVariableRef)::GeneralConstraintRef

Extend [`JuMP.UpperBoundRef`](@ref) to extract a constraint reference for the
upper bound of the original infinite variable of `vref`.

**Example**
```julia-repl
julia> cref = UpperBoundRef(vref)
var <= 1.0
```
"""
function JuMP.UpperBoundRef(vref::ReducedInfiniteVariableRef)::GeneralConstraintRef
    return JuMP.UpperBoundRef(infinite_variable_ref(vref))
end

"""
    JuMP.is_fixed(vref::ReducedInfiniteVariableRef)::Bool

Extend [`JuMP.is_fixed`](@ref) to return `Bool` whether the original infinite
variable of `vref` is fixed.

**Example**
```julia-repl
julia> is_fixed(vref)
true
```
"""
function JuMP.is_fixed(vref::ReducedInfiniteVariableRef)::Bool
    return JuMP.is_fixed(infinite_variable_ref(vref))
end

"""
    JuMP.fix_value(vref::ReducedInfiniteVariableRef)::Float64

Extend [`JuMP.fix_value`](@ref) to return the fix value of the original infinite
variable of `vref`. Errors if variable is not fixed.

**Example**
```julia-repl
julia> fix_value(vref)
0.0
```
"""
function JuMP.fix_value(vref::ReducedInfiniteVariableRef)::Float64
    if !JuMP.is_fixed(vref)
        error("Variable $(vref) is not fixed.")
    end
    return JuMP.fix_value(infinite_variable_ref(vref))
end

# Extend to return the index of the fix constraint associated with the original
# infinite variable of `vref`.
function JuMP._fix_index(vref::ReducedInfiniteVariableRef)::Int
    if !JuMP.is_fixed(vref)
        error("Variable $(vref) is not fixed.")
    end
    return JuMP._fix_index(infinite_variable_ref(vref))
end

"""
    JuMP.FixRef(vref::ReducedInfiniteVariableRef)::GeneralConstraintRef

Extend [`JuMP.FixRef`](@ref) to return the constraint reference of the fix
constraint associated with the original infinite variable of `vref`. Errors
`vref` is not fixed.

**Examples**
```julia-repl
julia> cref = FixRef(vref)
var == 1.0
```
"""
function JuMP.FixRef(vref::ReducedInfiniteVariableRef)::GeneralConstraintRef
    return JuMP.FixRef(infinite_variable_ref(vref))
end

"""
    JuMP.start_value(vref::ReducedInfiniteVariableRef)::Union{Nothing, Float64}

Extend [`JuMP.start_value`](@ref) to return starting value of the original
infinite variable of `vref` if it has one. Returns `nothing` otherwise.

**Example**
```julia-repl
julia> start_value(vref)
0.0
```
"""
function JuMP.start_value(vref::ReducedInfiniteVariableRef)::Union{Nothing, Float64}
    return JuMP.start_value(infinite_variable_ref(vref))
end

"""
    JuMP.is_binary(vref::ReducedInfiniteVariableRef)::Bool

Extend [`JuMP.is_binary`](@ref) to return `Bool` whether the original infinite
variable of `vref` is binary.

**Example**
```julia-repl
julia> is_binary(vref)
true
```
"""
function JuMP.is_binary(vref::ReducedInfiniteVariableRef)::Bool
    return JuMP.is_binary(infinite_variable_ref(vref))
end

# Extend to return the index of the binary constraint associated with the
# original infinite variable of `vref`.
function JuMP._binary_index(vref::ReducedInfiniteVariableRef)::Int
    if !JuMP.is_binary(vref)
        error("Variable $(vref) is not binary.")
    end
    return JuMP._binary_index(infinite_variable_ref(vref))
end

"""
    JuMP.BinaryRef(vref::ReducedInfiniteVariableRef)::GeneralConstraintRef

Extend [`JuMP.BinaryRef`](@ref) to return a constraint reference to the
constraint constrainting the original infinite variable of `vref` to be binary.
Errors if one does not exist.

**Example**
```julia-repl
julia> cref = BinaryRef(vref)
var binary
```
"""
function JuMP.BinaryRef(vref::ReducedInfiniteVariableRef)::GeneralConstraintRef
    return JuMP.BinaryRef(infinite_variable_ref(vref))
end

"""
    JuMP.is_integer(vref::ReducedInfiniteVariableRef)::Bool

Extend [`JuMP.is_integer`](@ref) to return `Bool` whether the original infinite
variable of `vref` is integer.

**Example**
```julia-repl
julia> is_integer(vref)
true
```
"""
function JuMP.is_integer(vref::ReducedInfiniteVariableRef)::Bool
    return JuMP.is_integer(infinite_variable_ref(vref))
end

# Extend to return the index of the integer constraint associated with the
# original infinite variable of `vref`.
function JuMP._integer_index(vref::ReducedInfiniteVariableRef)::Int
    if !JuMP.is_integer(vref)
        error("Variable $(vref) is not an integer.")
    end
    return JuMP._integer_index(infinite_variable_ref(vref))
end

"""
    JuMP.IntegerRef(vref::ReducedInfiniteVariableRef)::GeneralConstraintRef

Extend [`JuMP.IntegerRef`](@ref) to return a constraint reference to the
constraint constrainting the original infinite variable of `vref` to be integer.
Errors if one does not exist.

**Example**
```julia-repl
julia> cref = IntegerRef(vref)
var integer
```
"""
function JuMP.IntegerRef(vref::ReducedInfiniteVariableRef)::GeneralConstraintRef
    return JuMP.IntegerRef(infinite_variable_ref(vref))
end

"""
    used_by_constraint(vref::ReducedInfiniteVariableRef)::Bool

Return a `Bool` indicating if `vref` is used by a constraint.

**Example**
```julia-repl
julia> used_by_constraint(vref)
false
```
"""
function used_by_constraint(vref::ReducedInfiniteVariableRef)::Bool
    return haskey(JuMP.owner_model(vref).reduced_to_constrs, JuMP.index(vref))
end

"""
    used_by_measure(vref::ReducedInfiniteVariableRef)::Bool

Return a `Bool` indicating if `vref` is used by a measure.

**Example**
```julia-repl
julia> used_by_measure(vref)
true
```
"""
function used_by_measure(vref::ReducedInfiniteVariableRef)::Bool
    return haskey(JuMP.owner_model(vref).reduced_to_meas, JuMP.index(vref))
end

"""
    JuMP.is_valid(model::InfiniteModel, vref::ReducedInfiniteVariableRef)::Bool

Extend [`JuMP.is_valid`](@ref) to accomodate reduced infinite variables.

**Example**
```julia-repl
julia> is_valid(model, vref)
true
```
"""
function JuMP.is_valid(model::InfiniteModel,
                       vref::ReducedInfiniteVariableRef)::Bool
    return (model === JuMP.owner_model(vref) && JuMP.index(vref) in keys(model.reduced_info))
end

"""
    JuMP.delete(model::InfiniteModel, vref::ReducedInfiniteVariableRef)

Extend [`JuMP.delete`](@ref) to delete reduced infinite variables and its
dependencies. Errors if `vref` is invalid, meaning it has already been deleted
or it belongs to another model.

**Example**
```julia
julia> print(model)
Min measure(g(0, t)*t + g(1, t)*t) + z
Subject to
 z >= 0.0
 g(0, t) + g(1, t) == 0
 g(x, t) + z >= 42.0
 g(0.5, 0.5) == 0
 t in [0, 6]
 x in [0, 1]

julia> delete(model, rvref1)

julia> print(model)
Min measure(t + g(1, t)*t) + z
Subject to
 z >= 0.0
 g(1, t) == 0
 g(x, t) + z >= 42.0
 g(0.5, 0.5) == 0
 t in [0, 6]
 x in [0, 1]
```
"""
function JuMP.delete(model::InfiniteModel, vref::ReducedInfiniteVariableRef)
    # check valid reference
    @assert JuMP.is_valid(model, vref) "Invalid variable reference."
    # remove from measures if used
    if used_by_measure(vref)
        for mindex in model.reduced_to_meas[JuMP.index(vref)]
            if isa(model.measures[mindex].func, ReducedInfiniteVariableRef)
                model.measures[mindex] = Measure(zero(JuMP.AffExpr),
                                                 model.measures[mindex].data)
            else
                _remove_variable(model.measures[mindex].func, vref)
            end
            JuMP.set_name(MeasureRef(model, mindex),
                          _make_meas_name(model.measures[mindex]))
        end
        # delete mapping
        delete!(model.reduced_to_meas, JuMP.index(vref))
    end
    # remove from constraints if used
    if used_by_constraint(vref)
        for cindex in model.reduced_to_constrs[JuMP.index(vref)]
            if isa(model.constrs[cindex].func, ReducedInfiniteVariableRef)
                model.constrs[cindex] = JuMP.ScalarConstraint(zero(JuMP.AffExpr),
                                                      model.constrs[cindex].set)
            else
                _remove_variable(model.constrs[cindex].func, vref)
            end
        end
        # delete mapping
        delete!(model.reduced_to_constrs, JuMP.index(vref))
    end
    # remove mapping to infinite variable
    ivref = infinite_variable_ref(vref)
    filter!(e -> e != JuMP.index(vref),
            model.infinite_to_reduced[JuMP.index(ivref)])
    if length(model.infinite_to_reduced[JuMP.index(ivref)]) == 0
        delete!(model.infinite_to_reduced, JuMP.index(ivref))
    end
    # delete the info
    delete!(model.reduced_info, JuMP.index(vref))
    return
end


text:
# n! means n × (n − 1) × ... × 3 × 2 × 1
#
# For example, 10! = 10 × 9 × ... × 3 × 2 × 1 = 3628800,
# and the sum of the digits in the number 10! is 3 + 6 + 2 + 8 + 8 + 0 + 0 = 27.
#
# Find the sum of the digits in the number 100!

using ProjectEulerSolutions

# Trivial with Julia's built-in factorial and digits functions.
function p020solution(n::Integer=2)::Integer
    return sum(digits(factorial(BigInt(n))))
end

p020 = Problems.Problem(p020solution)

Problems.benchmark(p020, 100)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
###########################
# FunctionCell Interface  #
###########################
"""
	ohenc(opts)

Generates a FunctionCell that, when data is piped into it, trains a one-hot encoder (e.g. stores all
distinct values for a subset or all variables). `opts` can be either a string taking the
values `"binary"` or `"integer"` depending on the type of data encoding desired
or a `Dict(idx=>method)` if different methods are to be used for different variables.
In this case,

* `idx` can be an `Int`, `Vector{Int}` or `UnitRange` and specifies the variable indices
* `method` is a string that specifies the method. Available methods: `"binary"` and `"integer"`.
Unrecognized options will be ignored.

#Examples
```
 wu = ohenc("binary")
One-hot encoder (binary), varying I/O dimensions, untrained

julia> wt = ["a","b","c"] |> wu
One-hot encoder (binary), 1->3, trained

julia> (+wt).data
DataStructures.SortedDict{Int64,Tuple{String,Array{String,1}},Base.Order.ForwardOrdering} with 1 entry:
  1 => ("binary", String["a", "b", "c"])

julia> ["a","b","d","c","a","d"] |> wt

3×6 Array{Float64,2}:
 1.0  0.0  NaN  0.0  1.0  NaN
 0.0  1.0  NaN  0.0  0.0  NaN
 0.0  0.0  NaN  1.0  0.0  NaN

julia> a = datacell(["a" "a" "b" "b" "b"; "a" "b" "c" "d" "e"],["a","a","b","b","b"])
DataCell, 5 obs, 2 vars, 1 target(s)/obs, 2 distinct values: "b"(3),"a"(2)

julia> wt = a |> ohenc(Dict(1=>"integer",2=>"binary"))
One-hot encoder (mixed), 2->6, trained

julia> (+wt).data
DataStructures.SortedDict{Int64,Tuple{String,Array{String,1}},Base.Order.ForwardOrdering} with 2 entries:
  1 => ("integer", String["a", "b"])
  2 => ("binary", String["a", "b", "c", "d", "e"])

julia> ["x" "x" "a" "b"; "e" "d" "c" "x"] |> wt
6×4 Array{Float64,2}:
 NaN    NaN    1.0    2.0
   0.0    0.0  0.0  NaN
   0.0    0.0  0.0  NaN
   0.0    0.0  1.0  NaN
   0.0    1.0  0.0  NaN
   1.0    0.0  0.0  NaN
```
"""
ohenc(opts::T where T<:AbstractString) = FunctionCell(ohenc, (opts,), ModelProperties(), "One-hot encoder ("*opts*")")
ohenc(opts::T where T<:Dict) = FunctionCell(ohenc, (opts,), ModelProperties(), "One-hot encoder (mixed)")

############################
# DataCell/Array Interface #
############################
"""
	ohenc(data, opts)

Trains a one-hot encoder (e.g. stores all distinct values for a subset or all variables).
"""
# Training
ohenc(x::T where T<:CellData, opts) = ohenc(strip(x), opts)

ohenc(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractArray, opts) = ohenc(x[1], opts)

ohenc(x::T where T<:AbstractArray, opts) = begin

	# Get dictionary or construct a proper one from original;
	# Inputs: encoding option and total number of variables
	_sdict_(x::T, m) where T<:String = Dict(i=>x for i in 1:m)
	_sdict_(x::T, m) where T<:Dict = begin
		out = Dict{Int, String}()
		for (k,v) in x
			if(k isa Int)
				@assert (k>=1)&&(k<= m) "[ohenc] Index $(k) for a variable out of bounds."
				push!(out, k=>v)
			elseif (k isa UnitRange{Int} || k isa Vector{Int})
				@assert (minimum(k)>=1)&&(maximum(k)<= m) "[ohenc] Index $(k) for a variable out of bounds."
				push!(out, (ki=>v for ki in k)...)
			else
				error("[ohenc] Unsupported variable index syntax.")
			end
		end
		return out
	end
	dopts = _sdict_(opts, nvars(x))

	# Filter new dictionary for unrecognized entries and calculate size_out
	size_out = nvars(x) - length(keys(dopts)) 		# initialize with number of unprocessed variables
	for (k,v) in dopts
		if (v == "integer") size_out += 1
		elseif (v == "binary") size_out += length(unique(_variable_(x,k)))
		else
			delete!(dopts, k)			# non-recognized option, delete processing option from dictionary
			size_out += 1				# variable will still be found in output
		end
	end

	# Determine encoder name (for information only)
	_encodername_(x) = begin
		if isempty(x) return "One-hot encoder ()"
		elseif (length(unique(x)) == 1) return "One-hot encoder ("*collect(x)[1]*")"
		else return "One-hot encoder (mixed)"
		end
	end
	encodername = _encodername_(values(dopts))

	# Get list of variables to process/not to process
	variables_to_process = collect(keys(dopts))
	variables_unprocessed = setdiff(collect(1:nvars(x)), variables_to_process)

	# Create model properties
	modelprops = ModelProperties(nvars(x), size_out, nothing,
			Dict("size_unprocessed" => length(variables_unprocessed),
			"variables_to_process" => variables_to_process, 	# Which variables to process
			"variables_unprocessed" => variables_unprocessed
			)
	)

	# Train
	modeldata = SortedDict(Dict( idx => (v, unique(_variable_(x,idx))) for (idx,v) in dopts ))

	# Return trained cell
	FunctionCell(ohenc, Model(modeldata, modelprops), encodername)
end

# Execution
ohenc(x::T where T<:CellData, model::Model) = datacell(ohenc(strip(x), model))

ohenc(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractArray, model::Model) = (ohenc(x[1], model),x[2])

ohenc(x::T where T<:AbstractVector, model::Model) = ohenc(mat(x, ObsDim.Constant{2}()),model)

ohenc(x::T where T<:AbstractMatrix, model::Model) = begin

	# Read some inputs
	size_in::Int = model.properties.idim
	size_out::Int = model.properties.odim
	size_unprocessed::Int = model.properties.other["size_unprocessed"]
	variables_unprocessed::Vector{Int} = model.properties.other["variables_unprocessed"]

	# Allocate encoding matrix
	B = fill(NaN, size_out - size_unprocessed, nobs(x))

	# Iterate over processing options and write values
	i = 1 #internal row (e.g. new feature) counter
	for (idx, (encoding, v)) in model.data
		xv = view(x,idx,:)
		if (encoding == "integer")
			ohenc_integer!(_variable_(B,i), xv, v)
			i += 1
		elseif (encoding == "binary")
			ohenc_binary!(_variable_(B,i:i+length(v)-1), xv, v)
			i += length(v)
		else
			# Should not arrive here, unknown options are removed from start
			error("[ohenc] Unrecognized encoding option for variable $idx")
		end
	end

	# Build output from processed (B) and unprocessed data
	if !isempty(variables_unprocessed)
		if !isempty(B)
			# there are some variable unprocessed
			return dcat(B, _variable_(x,variables_unprocessed))
		else
			# all variables are unprocessed (B empty)
			return dcat(x)
		end
	else
		# all variables are processed
		return B
	end
end

# Binary encoder
function ohenc_binary(x::T where T<:AbstractVector, v::S where S<:AbstractVector, val::Float64=1.0)::Matrix{Float64}
	B = fill(NaN, length(v), length(x))
	ohenc_binary!(B, x, v, val)
	return B
end

function ohenc_binary!(B::AbstractMatrix{Float64}, x::T where T<:AbstractVector, v::S where S<:AbstractVector, val::Float64=1.0)
	vf::Vector{Int} = findin(x, v) 	# indices of samples in x found training values

	@inbounds for j in vf
		for i in eachindex(v)
			B[i,j] = ifelse(x[j]==v[i],val,0.0)
		end
	end
	return B
end

# Integer encoder
function ohenc_integer(x::T where T<:AbstractVector, v::S where S<:AbstractVector)
	n::Int64 = length(x) 		# number of objects
	B = fill(NaN, n)
	ohenc_integer!(B, x, v)
	return B
end

function ohenc_integer!(B::AbstractVector{Float64}, x::T where T<:AbstractVector, v::S where S<:AbstractVector)

	m::Int64 = length(v)		# number of distinct values to search in
	vf::Vector{Int} = findin(x, v) 	# indices of samples in x found training values

	l=1 # keeps track of the last location of x found in v

	# Function that searches fast a value through a vector of unique values
	function _fastsearch_(v::T,uv::Vector{T})::Int where T
		@inbounds for (i,vi) in enumerate(uv)
			if isequal(vi,v)
				return i
			end
		end
		return 0
	end

	for i in 1:length(vf)
		if (l < vf[i])
			B[l:vf[i]-1]=NaN
		end
		@inbounds B[vf[i]] = (1.0:m)[_fastsearch_(x[vf[i]],v)]
		l = vf[i]+1
	end
	return B
end


text:
function add_default_utctai_offsets_loaders(::Type{TimeScalesFactory})
    return jcall(TimeScalesFactory, "addDefaultUTCTAIOffsetsLoaders", void, ())
end

function add_utctai_offsets_loader(::Type{TimeScalesFactory}, arg0::UTCTAIOffsetsLoader)
    return jcall(TimeScalesFactory, "addUTCTAIOffsetsLoader", void, (UTCTAIOffsetsLoader,), arg0)
end

function clear_utctai_offsets_loaders(::Type{TimeScalesFactory})
    return jcall(TimeScalesFactory, "clearUTCTAIOffsetsLoaders", void, ())
end

function get_bdt(::Type{TimeScalesFactory})
    return jcall(TimeScalesFactory, "getBDT", BDTScale, ())
end

function get_glonass(::Type{TimeScalesFactory})
    return jcall(TimeScalesFactory, "getGLONASS", GLONASSScale, ())
end

function get_gmst(::Type{TimeScalesFactory}, arg0::IERSConventions, arg1::jboolean)
    return jcall(TimeScalesFactory, "getGMST", GMSTScale, (IERSConventions, jboolean), arg0, arg1)
end

function get_gps(::Type{TimeScalesFactory})
    return jcall(TimeScalesFactory, "getGPS", GPSScale, ())
end

function get_gst(::Type{TimeScalesFactory})
    return jcall(TimeScalesFactory, "getGST", GalileoScale, ())
end

function get_irnss(::Type{TimeScalesFactory})
    return jcall(TimeScalesFactory, "getIRNSS", IRNSSScale, ())
end

function get_qzss(::Type{TimeScalesFactory})
    return jcall(TimeScalesFactory, "getQZSS", QZSSScale, ())
end

function get_tai(::Type{TimeScalesFactory})
    return jcall(TimeScalesFactory, "getTAI", TAIScale, ())
end

function get_tcb(::Type{TimeScalesFactory})
    return jcall(TimeScalesFactory, "getTCB", TCBScale, ())
end

function get_tcg(::Type{TimeScalesFactory})
    return jcall(TimeScalesFactory, "getTCG", TCGScale, ())
end

function get_tdb(::Type{TimeScalesFactory})
    return jcall(TimeScalesFactory, "getTDB", TDBScale, ())
end

function get_time_scales(::Type{TimeScalesFactory})
    return jcall(TimeScalesFactory, "getTimeScales", LazyLoadedTimeScales, ())
end

function get_tt(::Type{TimeScalesFactory})
    return jcall(TimeScalesFactory, "getTT", TTScale, ())
end

function get_ut1(::Type{TimeScalesFactory}, arg0::EOPHistory)
    return jcall(TimeScalesFactory, "getUT1", UT1Scale, (EOPHistory,), arg0)
end

function get_ut1(::Type{TimeScalesFactory}, arg0::IERSConventions, arg1::jboolean)
    return jcall(TimeScalesFactory, "getUT1", UT1Scale, (IERSConventions, jboolean), arg0, arg1)
end

function get_utc(::Type{TimeScalesFactory})
    return jcall(TimeScalesFactory, "getUTC", UTCScale, ())
end


text:
function value(obj::StateFunction, arg0::SpacecraftState)
    return jcall(obj, "value", Vector{jdouble}, (SpacecraftState,), arg0)
end


text:
#module IdealGasBackend

# One important thing to do here is to make sure that we can handle cases
# with missing data. For Ideal Gas, the possible missing pieces of data are
# Heat Capacity (Cp) data and transport properties.
# The constructors for state structs that depend on heat capacity data (such
# as IdealGasPotentialsState) are responsible to throw an error if heat capacity
# data is missing.

struct IdealGasData{T_} <: AbstractEquationOfState where T <: AbstractFloat
    CompoundName::String
    MW::T_
    CpCoefficients::Union{Vector{T_}, Missing}
    CpEquationNumber::Union{Int, Missing}
end

struct IdealGasMinimalState{T_} where T_ <: AbstractFloat
    P::T_
    V::T_
    ρ::T_
    T::T_
    Cp::Union{T_, Missing}
end

struct IdealGasPotentialsState{T_} where T_ <: AbstractFloat
    S::T_
    U::T_
    H::T_
    A::T_
    G::T_
end

struct IdealGasChemicalState{T_} where T_ <: AbstractFloat
    xᵢ::T_
    f::T_
    ϕ::T_
end

struct IdealGasTransportState{T_} where T_ <: AbstractFloat
    λ::Union{T_, Missing}
    μ::Union{T_, Missing}
    ν::Union{T_, Missing}
    α::Union{T_, Missing}
    Pr::Union{T_, Missing}
end

struct IdealGasFullState{T_} where T_ <: AbstractFloat
    data::IdealGasData{T_}
    minstate::IdealGasMinimalState{T_}
    potentials::IdealGasPotentialsState{T_}
    chemstate::IdealGasChemicalState{T_}
    transtate::IdealGasTransportState{T_}
end

#end # module  IdealGasBackend


text:
export fixed_field, subfields

# Compute basis for the subfield of K that is generated by the elements of as.
function _subfield_basis(K::S, as::Vector{T}) where {
    S <: Union{AnticNumberField, Hecke.NfRel},
    T <: Union{nf_elem, Hecke.NfRelElem}
   }
  if isempty(as)
    return [gen(K)]
  end

  # Notation: k base field, K the ambient field, F the field generated by as

  k = base_field(K)

  d = degree(K)
  Kvs = VectorSpace(k,d)
  # We transition the coefficients of a in reverse order, so that the
  # first vector in the row reduced echelon form yields the highest
  # degree among all elements of Fas.
  (Fvs,phivs) = sub(Kvs, [Kvs([coeff(a,n) for n in d-1:-1:0])
                          for a in as])
  dF = length(Fvs.gens) # dim(Fvs)
  bs = as
  while !isempty(bs)
    nbs = elem_type(K)[]
    for b in bs
      abs = [a*b for a in as]
      abvs,_ = sub(Kvs, [Kvs([coeff(ab,n) for n in d-1:-1:0])
                         for ab in abs])
      (Fvs,phivs) = sub(Kvs, typeof(Fvs)[Fvs, abvs])
      if dF != length(Fvs.gens) # dim(Fvs)
        dF = length(Fvs.gens) # dim(Fvs)
        append!(nbs, abs)
      end
    end
    bs = nbs
  end

  kx = parent(K.pol)
  return [let Kv = phivs(v)
            K(kx([Kv[n] for n in d:-1:1]))
          end
          for v in gens(Fvs)]
end

function _improve_subfield_basis(K, bas)
  # First compute the maximal order of <bas> by intersecting and saturating
  # Then B_Ok = N * B_LLL_OK
  # Then B' defined as lllN * B_LLL_OK will hopefully be small
  OK = maximal_order(K)
  OKbmatinv = basis_mat_inv(OK, copy = false)
  basinOK = bas * matrix(FlintQQ, OKbmatinv.num) * fmpq(1, OKbmatinv.den)
  deno = fmpz(1)
  for i in 1:nrows(basinOK)
    for j in 1:ncols(basinOK)
      deno = lcm(deno, denominator(basinOK[i, j]))
    end
  end
  S = saturate(matrix(FlintZZ, basinOK * deno))
  SS = S * basis_matrix(OK, copy = false)
  lllOK = lll(OK)
  N = (SS * basis_mat_inv(lllOK)).num
  lllN = lll(N)
  maybesmaller = lllN * basis_matrix(lllOK)
  return maybesmaller
end

# Compute a primitive element given a basis of a subfield
function _subfield_primitive_element_from_basis(K::S, as::Vector{T}) where {
    S <: Union{AnticNumberField, Hecke.NfRel},
    T <: Union{nf_elem, Hecke.NfRelElem}
   }
  if isempty(as)
    return gen(K)
  end

  d = length(as)

  # First check basis elements
  i = findfirst(a -> degree(minpoly(a)) == d, as)
  if i <= d
    return as[i]
  end

  k = base_field(K)

  # Notation: cs the coefficients in a linear combination of the as, ca the dot
  # product of these vectors.
  cs = fmpz[zero(ZZ) for n in 1:d]
  cs[1] = one(ZZ)
  while true
    ca = sum(c*a for (c,a) in zip(cs,as))
    if degree(minpoly(ca)) == d
      return ca
    end

    # increment the components of cs
    cs[1] += 1
    let i = 2
      while i <= d && cs[i-1] > cs[i]+1
        cs[i-1] = zero(ZZ)
        cs[i] += 1
        i += 1
      end
    end
  end
end

################################################################################
#
#  Subfield
#
################################################################################

@doc Markdown.doc"""
    subfield(L::NumField, elt::Vector{<: NumFieldelem};
                          isbasis::Bool = false) -> NumField, Map

The simple number field $k$ generated by the elements of `elt` over the base
field $K$ of $L$ together with the embedding $k \to L$.

If isbasis is true, it is assumed that `elt` holds a $K$-basis of $k$.
"""
function subfield(K::NumField, elt::Vector{<:NumFieldElem}; isbasis::Bool = false)
  if length(elt) == 1
    return _subfield_from_primitive_element(K, elt[1])
  end

  if isbasis
    s = _subfield_primitive_element_from_basis(K, elt)
  else
    bas = _subfield_basis(K, elt)
    s = _subfield_primitive_element_from_basis(K, bas)
  end

  return _subfield_from_primitive_element(K, s)
end

function _subfield_from_primitive_element(K, s)
  f = minpoly(s)
  L,_ = NumberField(f, cached = false)
  return L, hom(L, K, s)
end

################################################################################
#
#  Fixed field
#
################################################################################

@doc Markdown.doc"""
    fixed_field(K::SimpleNumField,
                sigma::Map;
                simplify::Bool = true) -> NumberField, NfToNfMor

Given a number field $K$ and an automorphisms $\sigma$ of $K$, this function
returns the fixed field of $\sigma$ as a pair $(L, i)$ consisting of a number
field $L$ and an embedding of $L$ into $K$.

By default, the function tries to find a small defining polynomial of $L$. This
can be disabled by setting `simplify = false`.
"""
function fixed_field(K::SimpleNumField, sigma::T; simplify::Bool = true) where {T <: Union{NfToNfMor, NfRelToNfRelMor}}
  return fixed_field(K, T[sigma], simplify = simplify)
end

@doc Markdown.doc"""
    fixed_field(K::SimpleNumField, A::Vector{NfToNfMor}) -> NumberField, NfToNfMor

Given a number field $K$ and a set $A$ of automorphisms of $K$, this function
returns the fixed field of $A$ as a pair $(L, i)$ consisting of a number field
$L$ and an embedding of $L$ into $K$.

By default, the function tries to find a small defining polynomial of $L$. This
can be disabled by setting `simplify = false`.
"""
function fixed_field(K::SimpleNumField, A::Vector{T}; simplify::Bool = true) where {T <: Union{NfToNfMor, NfRelToNfRelMor}}
  autos = A

  # Everything is fixed by nothing :)
  if length(autos) == 0
    return K, id_hom(K)
  end

  F = base_field(K)
  a = gen(K)
  n = degree(K)
  ar_mat = Vector{dense_matrix_type(elem_type(F))}()
  v = Vector{elem_type(K)}(undef, n)
  for i in 1:length(autos)
    domain(autos[i]) !== codomain(autos[i]) && throw(error("Maps must be automorphisms"))
    domain(autos[i]) !== K && throw(error("Maps must be automorphisms of K"))
    o = one(K)
    # Compute the image of the basis 1,a,...,a^(n - 1) under autos[i] and write
    # the coordinates in a matrix. This is the matrix of autos[i] with respect
    # to 1,a,...a^(n - 1).
    as = autos[i](a)
    if a == as
      continue
    end
    v[1] = o
    for j in 2:n
      o = o * as
      v[j] = o
    end

    if isabsolute(K)
      bm = basis_matrix(v, FakeFmpqMat)
      # We have to be a bit careful (clever) since in the absolute case the
      # basis matrix is a FakeFmpqMat

      m = matrix(FlintQQ, bm.num)
      for j in 1:n
        m[j, j] = m[j, j] - bm.den # This is autos[i] - identity
      end
    else
      bm = basis_matrix(v)
      # In the generic case just subtract the identity
      m = bm - identity_matrix(F, degree(K))
    end

    push!(ar_mat, m)
  end

  if length(ar_mat) == 0
    return K, id_hom(K)
  else
    bigmatrix = hcat(ar_mat)
    k, Ker = kernel(bigmatrix, side = :left)
    bas = Vector{elem_type(K)}(undef, k)
    if isabsolute(K)
      # Try to find a small basis for absolute simple number fields
      if simplify
        KasFMat = _improve_subfield_basis(K, Ker)
        for i in 1:k
          bas[i] = elem_from_mat_row(K, KasFMat.num, i, KasFMat.den)
        end
      else
        KasFMat = FakeFmpqMat(Ker)
        Ksat = saturate(KasFMat.num)
        Ksat = lll(Ksat)
        onee = one(fmpz)
        for i in 1:k
          bas[i] = elem_from_mat_row(K, Ksat, i, onee)
        end
      end
    else
      for i in 1:k
        bas[i] = elem_from_mat_row(K, Ker, i)
      end
    end

    return subfield(K, bas, isbasis = true)
  end
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations3D(1.4)

initial_condition = initial_condition_convergence_test

# you can either use a single function to impose the BCs weakly in all
# 2*ndims == 4 directions or you can pass a tuple containing BCs for each direction
boundary_condition = BoundaryConditionDirichlet(initial_condition)
boundary_conditions = (x_neg=boundary_condition,
                       x_pos=boundary_condition,
                       y_neg=boundary_condition,
                       y_pos=boundary_condition,
                       z_neg=boundary_condition,
                       z_pos=boundary_condition,)

solver = DGSEM(polydeg=4, surface_flux=flux_hll,
               volume_integral=VolumeIntegralWeakForm())

# Mapping as described in https://arxiv.org/abs/2012.12040 but with less warping.
function mapping(xi, eta, zeta)
  # Don't transform input variables between -1 and 1 onto [0,3] to obtain curved boundaries
  # xi = 1.5 * xi_ + 1.5
  # eta = 1.5 * eta_ + 1.5
  # zeta = 1.5 * zeta_ + 1.5

  y = eta + 1/6 * (cos(1.5 * pi * (2 * xi - 3)/3) *
                   cos(0.5 * pi * (2 * eta - 3)/3) *
                   cos(0.5 * pi * (2 * zeta - 3)/3))

  x = xi + 1/6 * (cos(0.5 * pi * (2 * xi - 3)/3) *
                  cos(2 * pi * (2 * y - 3)/3) *
                  cos(0.5 * pi * (2 * zeta - 3)/3))

  z = zeta + 1/6 * (cos(0.5 * pi * (2 * x - 3)/3) *
                    cos(pi * (2 * y - 3)/3) *
                    cos(0.5 * pi * (2 * zeta - 3)/3))

  # Transform the weird deformed cube to be approximately the cube [0,2]^3
  return SVector(x + 1, y + 1, z + 1)
end

cells_per_dimension = (4, 4, 4)
mesh = StructuredMesh(cells_per_dimension, mapping, periodicity=false)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    source_terms=source_terms_convergence_test,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=0.6)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
"""

```
contrib(s::Selection, atom_contributions::Array{Float64}, selection)
```

Extract the contribution of a given atom type selection from the
solute or solvent atomic contributions to the MDDF.

`s` here is the solute or solvent selection (type `ComplexMixtures.Selection`)
`atom_contributions` is the `R.solute_atom` or `R.solvent_atom` arrays of the `Result` structure,
and the last argument is the selection of atoms from the solute to be considered, given as a list
of indexes, list of atom names, vector of `PDBTools.Atom`s, or a `PDBTools.Residue`.

"""
function contrib(s::Selection, atom_contributions::Array{Float64}, indexes::Vector{Int})
  nbins = size(atom_contributions,1)
  c = zeros(nbins)
  for it in indexes
    if it > s.natomspermol
      error("The index list contains atoms with indexes greater than the number of atoms of the molecule.")
    end
    c += atom_contributions[:,it]
  end
  return c
end

#
# If a list of atom names is provided
#
function contrib(s::Selection, atom_contributions::Array{Float64}, names::Vector{String})
  indexes = Vector{Int}(undef,0)
  for name in names
    index = findall(isequal(name), s.names)
    if length(index) == 0
      error(" Atom in input list is not part of solvent (or solute): $name")
    end
    append!(indexes,index)
  end
  return contrib(s, atom_contributions, indexes)
end

#
# If a list of atoms of PDBTools.Atom is provided
#
function contrib(s::Selection,
                 atom_contributions::Array{Float64},
                 atoms::Vector{PDBTools.Atom};
                 warning=true)
  (warning && s.nmols > 1) && warning_nmols_types()
  indexes = [ atom.index for atom in atoms ]
  # Check which types of atoms belong to this selection
  selected_types = which_types(s, indexes, warning=warning)
  return contrib(s, atom_contributions, selected_types)
end

#
# If a residue of type PDBTools.Residue is provided
#
function contrib(s::Selection,
                 atom_contributions::Array{Float64},
                 residue::Residue;
                 warning=true)
  (warning && s.nmols > 1) && warning_nmols_types()
  indexes = collect(residue.range)
  # Check which types of atoms belong to this selection
  selected_types = which_types(s, indexes, warning=warning)
  return contrib(s, atom_contributions, selected_types)
end

function warning_nmols_types()
  println("WARNING: there is more than one molecule in this selection. "*
          "Contributions are summed over all atoms of the same type.")
end


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

#=
Implementation of submodule manager. Submodules are optional modules with heavy
dependencies. See documentation for usage and implementation details.
=#

using Pkg
using Scratch
using Suppressor
using TOML

ENV_DIR = ""                                                # directory of submodule environment
const DEV_DIR = mkpath(joinpath(first(DEPOT_PATH), "dev"))  # directory of packages added via `dev`
const PKG_NAME = "$(@__MODULE__)"                           # top-level package name (AdaStress)
const PKG_PATH = dirname(dirname(pathof(@__MODULE__)))      # top-level package directory
const SUBMODULES = Dict{String, String}()                   # submodule table

"""
Suppresses output if verbosity condition is not met.
"""
macro verboseif(v, expr)
    :($(esc(v)) ? $(esc(expr)) : @suppress $(esc(expr)))
end

"""
    exclude(dir::String)

Associate package directory to main package via submodule table instead of direct code
loading. Equivalent of `include` for submodules / optional dependencies.
"""
function exclude(dir::String)
    fprev = stacktrace()[2].file
    path = joinpath(dirname(String(fprev)), dir)
    name = basename(dir) # directory name must match corresponding module name
    SUBMODULES[name] = path
    return
end

"""
    submodules()

List all available submodules.
"""
submodules() = collect(keys(SUBMODULES))

"""
    enabled()

List enabled submodules.
"""
function enabled(; verbose::Int64=1)
    curr = Pkg.project().path
    @verboseif (verbose >= 2) Pkg.activate(ENV_DIR)
    deps = filter(p -> p[2].is_direct_dep, Pkg.dependencies())
    names = (d -> d.name).(values(deps))
    @verboseif (verbose >= 2) Pkg.activate(curr)
    return filter(d -> d != PKG_NAME && d in keys(SUBMODULES), names)
end

"""
List unregistered dependencies of submodule.
"""
function unregistered_deps(submodule::String)
    path = SUBMODULES[submodule]
    toml = TOML.parsefile(joinpath(path, "Project.toml"))
    return "unreg" in keys(toml) ? toml["unreg"] : Dict{String, Any}()
end

"""
    enable(submodule)

Enable submodule(s). Accepts string or vector of strings. With zero arguments defaults to
all associated submodules. Takes effect immediately.
"""
function enable(submodules::Vector{String}; verbose::Int64=1)
    curr = Pkg.project().path
    pkgs = [PackageSpec(path=PKG_PATH)] # update AdaStress if necessary
    dev_pkgs = readdir(DEV_DIR)
    for submodule in submodules
        for (dep, url) in unregistered_deps(submodule)
            push!(pkgs, dep in dev_pkgs ? PackageSpec(name=dep) : PackageSpec(url=url))
        end
        push!(pkgs, PackageSpec(path=SUBMODULES[submodule]))
    end

    try
        @verboseif (verbose >= 2) begin
            Pkg.activate(ENV_DIR)
            Pkg.develop(pkgs)
            for submodule in submodules
                @eval using $(Symbol(submodule))
            end
        end
        @verboseif (verbose >= 1) foreach(s -> @info("Enabled submodule $s."), submodules)
    catch e
        @verboseif (verbose >= 1) @error "Cannot enable submodule(s)."
        @verboseif (verbose >= 2) foreach(pkg -> try Pkg.rm(pkg) catch end, pkgs[2:end])
        throw(e)
    finally
        @verboseif (verbose >= 2) Pkg.activate(curr)
    end
    return
end

enable(submodule::String; kwargs...) = enable([submodule]; kwargs...)
enable(; kwargs...) = enable(submodules(); kwargs...)

"""
    disable(submodule)

Disable submodule(s). Accepts string or vector of strings. With zero arguments defaults to
all enabled submodules. Takes effect after Julia restart.
"""
function disable(submodules::Vector{String}; verbose::Int64=1)
    curr = Pkg.project().path
    pkgs = String[]
    for submodule in submodules
        push!(pkgs, submodule)
        for (dep, _) in unregistered_deps(submodule)
            # NOTE: will cause problems if multiple enabled submodules use the same
            #       unregistered dependencies, but this is unlikely, given the ideally
            #       rare usage of such dependencies.
            push!(pkgs, dep)
        end
    end

    try
        @verboseif (verbose >= 2) begin
            Pkg.activate(ENV_DIR)
            Pkg.rm(pkgs)
        end
        @verboseif (verbose >= 1) foreach(s -> @info("Disabled submodule $s."), submodules)
    catch e
        @verboseif (verbose >= 1) @error "Cannot disable submodule(s)."
        throw(e)
    finally
        @verboseif (verbose >= 2) Pkg.activate(curr)
    end
    return
end

disable(submodule::String; kwargs...) = disable([submodule]; kwargs...)
disable(; kwargs...) = disable(enabled(; kwargs...); kwargs...)

"""
    load()

Load enabled submodules (necessary after Julia restart). Takes effect immediately.
"""
function load(; verbose::Int64=1)
    curr = Pkg.project().path
    try
        @verboseif (verbose >= 2) Pkg.activate(ENV_DIR)
        for submodule in enabled()
            @suppress @eval using $(Symbol(submodule))
            @verboseif (verbose >= 1) @info "Loaded submodule $submodule."
        end
    catch e
        @verboseif (verbose >= 1) @error "Unable to load enabled submodules."
        throw(e)
    finally
        @verboseif (verbose >= 2) Pkg.activate(curr)
    end
    return
end

"""
    clean()

Forcibly remove temporary environment, purging all enabled submodules. Only necessary if
submodule manager is corrupted and `disable` cannot restore functionality. Takes effect
after Julia restart.
"""
function clean(; verbose::Int64=1)
    rm.(joinpath.(ENV_DIR, readdir(ENV_DIR)))
    @verboseif (verbose >= 1) @info "Cleaned submodule environment."
end

"""
Invoke any submodule-related initializations. Called by top-level `__init__` function.
"""
function init_submodules()
    # scratchspace for intermediate projects
    global ENV_DIR = @get_scratch!("env")
end


text:
using PtFEM
using Test

ProjDir = dirname(@__FILE__)

#=
Compare formulas at:
http://www.awc.org/pdf/codes-standards/publications/design-aids/AWC-DA6-BeamFormulas-0710.pdf
=#

data = Dict(
  # Frame(nels, nn, ndim, nst, nip, finite_element(nod, nodof))
  :struc_el => Frame(200, 201, 3, 1, 1, Line(2, 3)),
  :properties => [1.0e6 1.0e6 1.0e6 3.0e5;],
  :x_coords => collect(range(0, stop=4, length=201)),
  :y_coords => zeros(201),
  :z_coords => zeros(201),
  :g_num => [
    collect(1:200)';
    collect(2:201)'],
  :support => [
    (1, [0 0 0 0 0 0]),
    (201, [0 0 0 0 0 0]),
    ],
  :loaded_nodes => [
    (101, [0.0 -10000.0 0.0 0.0 0.0 0.0])]
)

data |> display
println()

m, dis_df, fm_df = p44(data)

println("Displacements:")
m.displacements |> display
println()

println("Actions:")
m.actions |> display
println()

println("y displacements:")
m.displacements[2,:] |> display
println()

println("y moment actions:")
m.actions[12,:] |> display
println()

using Plots
gr(size=(400,600))

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 3)
p[1] = plot(m.x_coords, m.displacements[2,:], ylim=(-0.004, 0.001), lab="Displacement",
 xlabel="x [m]", ylabel="deflection [m]", color=:red)
p[2] = plot(m.actions[2,:], lab="Shear force", ylim=(-6000, 6000), xlabel="element",
  ylabel="shear force [N]", palette=:greens,fill=(0,:auto),α=0.6)
p[3] = plot(m.actions[12,:], lab="Moment", ylim=(-6000, 6000), xlabel="element",
  ylabel="moment [Nm]", palette=:grays,fill=(0,:auto),α=0.6)

plot(p..., layout=(3, 1))
#savefig(ProjDir*"/figure-24.pdf")
gui()

# See figure 24 in above reference (Δmax):
@test round.(m.displacements[2,101], digits=11) ≈ -0.00333333333 atol = 0.001

# See figure 24 in above reference (Mmax):
@test m.actions[12,101] ≈ 4899.999997801155

max_stars_repo_path: examples/Nodal Formulas/figure-24.jl

max_issues_repo_path: examples/Nodal Formulas/figure-24.jl

max_forks_repo_path: examples/Nodal Formulas/figure-24.jl


text:
module Constants

    export cee,
           emgm,
           pmgm,
           kboltz,
           kboltz_kev2erg,
           gr_to_kev,
           me_kev,
           emerg,
           charg,
           sigtom,
           herg,
           hkev,
           mjy,
           re0,
           gconst,
           sbconst,
           aconst,
           kpc,
           msun

    using GSL
    const GSL_CONST = GSL.Const

    const cee  =           GSL_CONST.CGSM_SPEED_OF_LIGHT
    const emgm=            GSL_CONST.CGSM_MASS_ELECTRON
    const pmgm  =          GSL_CONST.CGSM_MASS_PROTON
    const kboltz    =      GSL_CONST.CGSM_BOLTZMANN
    const kboltz_kev2erg = 1.6022e-9
    const gr_to_kev   =    5.6095883571872e+29
    const me_kev    =      511.
    const emerg   =        (GSL_CONST.CGSM_MASS_ELECTRON*GSL_CONST.CGSM_SPEED_OF_LIGHT^2.0)
    const charg   =        4.8e-10
    const sigtom =         GSL_CONST.CGSM_THOMSON_CROSS_SECTION
    const herg   =         GSL_CONST.CGSM_PLANCKS_CONSTANT_H
    const hkev   =         (GSL_CONST.CGSM_PLANCKS_CONSTANT_H*6.2415e8)
    const mjy   =          1.e-26
    const re0    =         2.81794e-13
    const gconst  =        GSL_CONST.CGSM_GRAVITATIONAL_CONSTANT
    const sbconst  =       GSL_CONST.CGSM_STEFAN_BOLTZMANN_CONSTANT
    const aconst   =       7.56e-15
    const kpc =             1e3*GSL_CONST.CGSM_PARSEC
    const msun =            GSL_CONST.CGSM_SOLAR_MASS

end


text:
function OutlineExtractor(arg0::Vector3D, arg1::Vector3D)
    return OutlineExtractor((Vector3D, Vector3D), arg0, arg1)
end

function get_outline(obj::OutlineExtractor, arg0::PolyhedronsSet)
    return jcall(obj, "getOutline", Vector{Vector{Vector2D}}, (PolyhedronsSet,), arg0)
end


text:
using FastRef
using Test

@testset "FastRef.jl" begin
    # Write your tests here.
end


text:
# hyperbolic_expansion.jl

module hypergeometric_expansion
export hypergeom_exp

include("arb_hypgeom.jl")
import .arb_hypgeom
using ArbNumerics

function arb_binomial(x, y; prec = 64)
    setextrabits(0)
    setprecision(ArbReal, prec)
    x = ArbReal("$x")
    y = ArbReal("$y")
    one = ArbReal("1")
    ArbNumerics.gamma(x + one) / (ArbNumerics.gamma(y + one) * ArbNumerics.gamma(x - y + one))
end

"""
This function implements the binomial expansion of the denominator of the integrals:

    ∫_0^1 cosh(zt) / (4a^2/β^2 - t^2)^(n+3/2) dt
    ∫_0^1 sinh(zt) / (4a^2/β^2 - t^2)^(n+3/2) dt

which is:

    ∑_{m=0}^∞ C(-n-3/2, m) (-1)^m (β/2a)^{2n+2m+3} ∫_0^1 cosh(zt) t^{2m} dt
    ∑_{m=0}^∞ C(-n-3/2, m) (-1)^m (β/2a)^{2n+2m+3} ∫_0^1 sinh(zt) t^{2m} dt.

and also substitutes the integrals for:

    ∫_0^1 cosh(zt) t^{2m} dt = ∑_{t=0}^∞ z^{2t} / ((2t+2m+1)⋅(2t)!)
    ∫_0^1 sinh(zt) t^{2m} dt = ∑_{t=0}^∞ z^{2t+1} / ((2t+2m+2)⋅(2t+1)!)

by replacing those integrals with one_F_two_fast() from arb_hypgeom.jl. h=0 gives cosh version, h=1 gives sinh version.
"""
function hypergeom_exp(z, n, β, a, h; prec = 64)
# h = 0 gives cosh, h = 1 gives sinh.

    # Initialise precision of ArbReal to prec.
    p = prec
    setextrabits(0)
    setprecision(ArbReal, p)

    z = ArbReal("$z")
    n = ArbReal("$n")
    β = ArbReal("$β")
    a = ArbReal("$a")

    m = ArbReal("0")
    result = ArbReal("0.0")
    term = ArbReal("1.0")
    err = eps(result)  # Machine accuracy of specified precision prec.

    while abs(midpoint(term)) > err * abs(midpoint(result))

        term = ArbReal(β * arb_binomial(-n - 3/2, m; prec = p) * (-1)^m * (β / (2 * a))^(2 * m) * arb_hypgeom.one_f_two_fast(z, m, h; prec = p) / a^(n + 2))
        result += term
        # println("term: m = ", m, "\nterm value: ", ball(ArbReal(term, bits = prec)), "\ncumulant result: ", ball(ArbReal(result, bits = prec)), "\n")
        m += ArbReal("1")

        # Double precision if rounding error in result exceeds accuracy specified by prec.
        if radius(result) > err * abs(midpoint(result))
            p *= 2
            setprecision(ArbReal, p)
            # println("Not precise enough. Error = ", abs(radius(result)/midpoint(result)), " > ", err, ". Increasing precision to ", p, " bits.\n")

            n = ArbReal("$n")
            z = ArbReal("$z")
            β = ArbReal("$β")
            a = ArbReal("$a")

            m = ArbReal("0")
            result = ArbReal("0.0")
            term = ArbReal("1.0")
        end
    end
    # println("z: ", ArbReal(z, bits = prec), ". Final result: ", ArbReal(result, bits = prec))
    ArbReal(result, bits = prec)
end

# β = ArbReal("2.0")
# α = ArbReal("7.0")
# v = ArbReal("5.8")
# w = ArbReal("1.6")
# R = ArbReal((v^2 - w^2) / (w^2 * v))
# a = ArbReal(sqrt(β^2 / 4 + R * β * coth(β * v / 2)))
# z = ArbReal("90.61")
# n = ArbReal("12.0")
# @time c = hypergeom_exp(z, n, β, a, 0; prec = 64)
# @show(c)
end # end module


text:
using Docile, Lexicon, ML

# config = Config()
config = Config(md_subheader=:skip, mathjax = true)

index = save("ML.md", ML, config)
save("index.md", Index([index]), config)

# run(`../mkdocs build`)
# mkdocs gh-deploy --clean


text:
import FunctionWrappers.FunctionWrapper;
using BenchmarkTools;

const RheoFloat = Float64

#
# Comment / uncomment the paramteric types to test their effect
#
struct RheoModel#{GT, GaT}

    params::Tuple

    _G#::GT
    _Ga#::GaT
end

#
#  Helper functions to shape expressions
#
expr_replace!(ex, s, v) = ex
expr_replace!(ex::Symbol, s, v) = s == ex ? v : ex

function expr_replace!(ex::Expr, s, v)
    for i=1:length(ex.args)
        ex.args[i] = expr_replace!(ex.args[i], s, v)
    end
    return ex
end

expr_replace(ex, s, v) = expr_replace!(copy(ex), s, v)

function expr_replace(ex, nt)
    e = copy(ex)
    k=keys(nt)
    v=values(nt)
    for i in 1:length(nt)
        expr_replace!(e, k[i], v[i])
    end
    return e
end

#
#  Create the scalar and vector function wrappers for the moduli.
#
function buildmoduli(G::Expr, paramsymbs::Tuple, Ga=nothing)
	#Ga=quote k.*exp.((-k./η).*ta) end

	if Ga==nothing
		Ge = expr_replace(G,NamedTuple{paramsymbs}((:(p_arr[1]), :(p_arr[2]))) )
		@eval return( (         ((t,p_arr) -> $Ge)                 |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}},
		        ((ta,p_arr) -> [$Ge for t in ta] )  |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}} ) )
	else
		Gae = expr_replace(Ga,NamedTuple{paramsymbs}((:(p_arr[1]), :(p_arr[2]))) )
		@eval return( (         ((t,p_arr) -> $Gae)                 |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}},
				((t,p_arr) -> $Gae )  |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}} ) )
	end
end

function buildmoduli_fixed(G::Expr, pvals::NamedTuple, Ga=nothing)

	if Ga==nothing
		Ge = expr_replace(G,pvals )
		println(Ge)
		@eval return( (         (t -> $Ge)                 	|> FunctionWrapper{RheoFloat,Tuple{RheoFloat}},
						(ta -> [$Ge for t in ta])  			|> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}} ) )
	else
		println("here")
		Gae = expr_replace(Ga,pvals )
		@eval return( (         (t -> $Gae)                 |> FunctionWrapper{RheoFloat,Tuple{RheoFloat}},
						(t -> $Gae)  						|> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}} ) )
	end
end
# julia> Meta.parse(":(p_arr[$i])").args[1]
# :(p_arr[1])

const nanexp = quote NaN end

function RheoModel1(;p::Tuple,         G::Expr = nanexp, Ga = nothing       )
    return(RheoModel(p,buildmoduli(G,p,Ga)... ))
end

function RheoModel1_fixed(;p::Tuple,         G::Expr = nanexp , pvals, Ga = nothing          )
    return(RheoModel(p,buildmoduli_fixed(G,pvals,Ga)... ))
end

function RheoModel2(;p::Tuple,   G = nothing, Ga = nothing )
	if Ga==nothing
    	return(RheoModel(p, (t::RheoFloat,p_arr) -> G(t; NamedTuple{p}(p_arr)...), (t_arr::Vector{RheoFloat},p_arr) -> begin nt=NamedTuple{p}(p_arr); [G(t; nt...) for t in t_arr]; end ))
	else
		return(RheoModel(p, (t::RheoFloat,p_arr) -> Ga(t; NamedTuple{p}(p_arr)...), (t::Vector{RheoFloat},p_arr) -> Ga(t; NamedTuple{p}(p_arr)...)  ))
	end
end

function RheoModel2_fixed(;p::Tuple,   G = nothing, pvals , Ga = nothing   )
	if Ga==nothing
    	return(RheoModel(p, t::RheoFloat -> G(t; pvals...), t_arr::Vector{RheoFloat} -> [G(t; pvals...) for t in t_arr] ))
	else
		return(RheoModel(p, t::RheoFloat -> Ga(t; pvals...), t::Vector{RheoFloat}-> Ga(t; pvals...) ))
	end
end

function RheoModel3(;p::Tuple,   G::Function  )
    return(RheoModel(p, G, (t_arr,p_arr) -> [G(t, p_arr) for t in t_arr] ))
end

function RheoModel3_fixed(;p::Tuple,   G::Function, pvals  )
    return(RheoModel(p, t->G(t,pvals), t_arr -> [G(t, pvals) for t in t_arr] ))
end

function RheoModel4(;p::Tuple,   G::Function  )
    return(RheoModel(p, (t,p_arr) -> G([t], p_arr)[1], G ))
end

function RheoModel4_fixed(;p::Tuple,   G::Function  , pvals)
    return(RheoModel(p, t -> G([t], pvals)[1], t_arr->G(t_arr,pvals) ))
end

#
#  Choose between the scalar/vectorized for of the modulus
#

#m1 = RheoModel1(p=(:k, :η), G =quote k*exp((-k/η)*t) end)
m1 = RheoModel1(p=(:k, :η), Ga=quote k.*exp.((-k./η).*t) end)

#m1_fixed = RheoModel1_fixed(p=(:k, :η), G=quote k*exp((-k/η)*t) end, pvals=(k=1.0, η=2.0))
m1_fixed = RheoModel1_fixed(p=(:k, :η), Ga = quote k.*exp.((-k./η).*t) end, pvals=(k=1.0, η=2.0))

function G_Maxwell2(t::RheoFloat; k, η)::RheoFloat
    return( k * exp( (-k/η)*t) )
end

function Ga_Maxwell2(t::Union{RheoFloat,Vector{RheoFloat}}; k, η)::Union{RheoFloat,Vector{RheoFloat}}
    return( k .* exp.( (-k./η).*t) )
end

#
#  Choose between the scalar/vectorized for of the modulus
#
#m2 = RheoModel2( p=(:k, :η), G=G_Maxwell2 )
m2 = RheoModel2( p=(:k, :η), Ga=Ga_Maxwell2 )

#m2_fixed = RheoModel2_fixed( p=(:k, :η), G=G_Maxwell2, pvals=(k=1.0, η=2.0) )
m2_fixed = RheoModel2_fixed( p=(:k, :η), Ga=Ga_Maxwell2, pvals=(k=1.0, η=2.0) )

function G_Maxwell3(t, params)
    return(params[1]*exp((-params[1]/params[2])*t))
end
m3 = RheoModel3( p=(:k, :η), G=G_Maxwell3 )
m3_fixed = RheoModel3_fixed( p=(:k, :η), G=G_Maxwell3, pvals=[1.0,2.0] )

function G_Maxwell4(t_arr, params)
    return(params[1].*exp.((-params[1]/params[2]).*t_arr))
end

m4 = RheoModel4( p=(:k, :η), G=G_Maxwell4 )
m4_fixed = RheoModel4_fixed( p=(:k, :η), G=G_Maxwell4, pvals=[1.0,2.0]  )

t=RheoFloat.(0.0:0.001:1.0);
t0=1.0; k=1.0; η=2.0

function dosomething(m,t::RheoFloat)
	return((m._G(t)))
end

function dosomething(m,t::Vector{RheoFloat})
	return((m._Ga(t)))
end

function dosomething(m,t::Vector{RheoFloat},p)
	return((m._Ga(t,p)))
end

function dosomething(m,t::RheoFloat,p)
	return((m._G(t,p)))
end

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

println("")
println("Tests...")
println("")
println("")

println("Direct function calls")
#@btime k*exp((-k/η)*t0)
#@btime 1.0*exp((-1.0/2.0)*t0)
println("Scalar function with keyword params with set parameters")
@btime G_Maxwell2(t0, k=1.0, η=1.0)
println("Scalar function with keyword params with passed parameters")
@btime G_Maxwell2(t0, k=k, η=η)
println("Scalar function with array params")
@btime G_Maxwell3(t0, [1.,2.])
println("Vector function with array params")
@btime G_Maxwell4([t0], [1.,2.])[1]
println("")

println("")
println("RheoModel function on scalar time values with free params")
println("")

println("Function wrappers")
print("   Direct call: "); @btime m1._G(1.,[1.,2.])
print("   In function: "); @btime dosomething(m1,2.0,[1.0,2.0]);
println("")

println("Scalar function with keyword params")
print("   Direct call: "); @btime m2._G(1.,[1.,2.])
print("   In function: "); @btime dosomething(m2,2.0,[1.0,2.0]);
println("")

println("Scalar function with array params")
print("   Direct call: "); @btime m3._G(1.,[1.,2.])
print("   In function: "); @btime dosomething(m3,2.0,[1.0,2.0]);
println("")

println("Vector function with array params")
print("   Direct call: "); @btime m4._G(1.,[1.,2.])
print("   In function: "); @btime dosomething(m4,2.0,[1.0,2.0]);
println("")

println("")
println("RheoModel function on array time values with free params")
println("")

println("Function wrappers")
print("   Direct call: "); @btime m1._Ga(t,[1.,2.]);
print("   In function: "); @btime dosomething(m1,t,[1.0,2.0]);
println("")

println("Scalar function with keyword params")
print("   Direct call: "); @btime m2._Ga(t,[1.,2.]);
print("   In function: "); @btime dosomething(m2,t,[1.0,2.0]);
println("")

println("Scalar function with array params")
print("   Direct call: "); @btime m3._Ga(t,[1.,2.]);
print("   In function: "); @btime dosomething(m3,t,[1.0,2.0]);
println("")

println("Vector function with array params")
print("   Direct call: "); @btime m4._Ga(t,[1.,2.]);
print("   In function: "); @btime dosomething(m4,t,[1.0,2.0]);
println("")

println("")
println("RheoModel function on scalar time values with fixed params")
println("")

println("Function wrappers")
print("   Direct call: "); @btime m1_fixed._G(1.0);
print("   In function: "); @btime dosomething(m1_fixed,2.0);
println("")
println("Scalar function with keyword params")
print("   Direct call: "); @btime m2_fixed._G(1.0);
print("   In function: "); @btime dosomething(m2_fixed,2.0);
println("")
println("Scalar function with array params")
print("   Direct call: "); @btime m3_fixed._G(1.0);
print("   In function: "); @btime dosomething(m3_fixed,2.0);
println("")
println("Vector function with array params")
print("   Direct call: "); @btime m4_fixed._G(1.0);
print("   In function: "); @btime dosomething(m4_fixed,2.0);
println("")

println("")
println("RheoModel function on array time values with fixed params")
println("")

println("Function wrappers")
print("   Direct call: "); @btime m1_fixed._Ga(t);
print("   In function: "); @btime dosomething(m1_fixed,t);
#println(m1_fixed._Ga(t)[end])
println()

println("Scalar function with keyword params")
print("   Direct call: "); @btime m2_fixed._Ga(t);
print("   In function: "); @btime dosomething(m2_fixed,t);
#println(m2_fixed._Ga(t)[end])
println()

println("Scalar function with array params")
print("   Direct call: "); @btime m3_fixed._Ga(t);
print("   In function: "); @btime dosomething(m3_fixed,t);
#println(m3_fixed._Ga(t)[end])
println()
println("Vector function with array params")
print("   Direct call: "); @btime m4_fixed._Ga(t);
print("   In function: "); @btime dosomething(m4_fixed,t);
#println(m4_fixed._Ga(t)[end])
println()

nothing

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
# -----------------------------------------------------------
#
# An example of interpolating a given function in 1D and 2D
#
# -----------------------------------------------------------

using GalerkinSparseGrids
using LinearAlgebra

D = 1
# Random function to interpolate:
truesoln = x -> 1.152 * norm(x)^2 * sin(dot([1.214], x))^3
k_max = 5
n_max = 5

println("Interpolation in ", D, "D.")
println("Going to k_max = ", k_max, ", n_max = ", n_max, ":")

for k_used in 1:k_max
    for n_used in 1:n_max
        # By default, sparse grids are used. Change scheme to "full"
        # to compare the different accuraces and costs of the
        # interpolations
        dict = coeffs_DG(D, k_used, n_used, truesoln)
        err = mcerr(x->reconstruct_DG(dict, [x...]), truesoln, D)
        println("(k = ", k_used, ", n = ", n_used, ") : err = ", err)
    end
end

D = 2
# Random function to interpolate:
truesoln = x -> 1.152 * norm(x)^2 * sin(dot([1.214, -0.5], x))^3
k_max = 5
n_max = 5

println("Interpolation in ", D, "D.")
println("Going to k_max = ", k_max, ", n_max = ", n_max, ":")

for k_used in 1:k_max
    for n_used in 1:n_max
        # By default, sparse grids are used. Change scheme to "full"
        # to compare the different accuraces and costs of the
        # interpolations
        dict = coeffs_DG(D, k_used, n_used, truesoln)
        err = mcerr(x->reconstruct_DG(dict, [x...]), truesoln, D)
        println("(k = ", k_used, ", n = ", n_used, ") : err = ", err)
    end
end


text:
using Compat
# Do not delete this, this will let you load other codes locally
push!(LOAD_PATH, joinpath(@__DIR__, "modules"))
push!(LOAD_PATH, joinpath(@__DIR__, "modules/dataset"))
push!(LOAD_PATH, joinpath(@__DIR__, "."))

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
using Diana

struct QueryOptions
    name::String
    query::String
    variables::Dict{String,Any}
end

struct MutationOptions
    name::String
    mutation::String
    variables::Dict{String,Any}
end

struct NavAbilityClient
    query::Function
    mutate::Function
end

function NavAbilityWebsocketClient( apiUrl::String="wss://api.navability.io/graphql")::NavAbilityClient
    throw("Not implemented")
end

function NavAbilityHttpsClient(
        apiUrl::String="https://api.navability.io";
        authorize::Bool=false
    )::NavAbilityClient
    #
    dianaClient = GraphQLClient(apiUrl)

    # auth
    if authorize
        # FIXME, use Base.getpass instead of readline once VSCode supports getpass.
            # st = Base.getpass("Copy-paste auth token")
            # seekstart(st)
            # tok = read(st, String)
            # Base.shred!(st)
        println("  > VSCode ONLY WORKAROUND, input issue, see https://github.com/julia-vscode/julia-vscode/issues/785")
        println("  >  Workaround: first press 0 then enter, and then paste the token and hit enter a second time.")
        println("Copy-paste auth token: ")
        tok = readline(stdin)
        dianaClient.serverAuth("Bearer "*tok)
    end

    function query(options::QueryOptions)
        # NOTE, the query client library used is synchronous, locally converted to async for package consistency
        @async dianaClient.Query(options.query, operationName=options.name, vars=options.variables)
    end
    function mutate(options::MutationOptions)
        # NOTE, the query client library used is synchronous, locally converted to async for package consistency
        @async dianaClient.Query(options.mutation, operationName=options.name, vars=options.variables)
    end
    return NavAbilityClient(query, mutate)
end


text:
using .Flux
using .Flux.Tracker: TrackedArray, TrackedReal

if !usetyped
  unwrap(x::Union{TrackedArray,TrackedReal}) = Flux.data(x)
end

forward(f, ps::Flux.Tracker.Params) = forward(f, Params(ps))


text:
using sb, LinearAlgebra, Statistics, StaticArrays

struct storages
    time :: Vector{Float64}
    energy :: Vector{Float64}
    function storages()
        new(Vector{Float64}(),Vector{Float64}())
    end
end

function store!(storage, t, ops, sbm)
    append!(storage.time,t)
    append!(storage.energy  ,sb.ECEIDenergy(ops,sbm))
end

function my_run(thermostat)
    storage = storages()
    # Create a model
    sbm = SBModel(OhmicSBParams(ϵ = 0.0, Δ = 1.0, ωc = 7.5, α = 0.1, No = 60, β = 5.0))
    # Create containers starting from a localized state
    ops = ECEIDOps(SA[1.0+0.0im 0.0im; 0.0im 0.0im],sbm)
    # run dynamics
    RunECEID!(sbm,ops,2000,0.001,store!, storage, thermostat = thermostat)
    return storage.time,storage.energy
end

ts, energies = my_run(false)
@test mean(energies .- energies[1]) < 1e-13

ts, energies = my_run(true)
@test (energies[2000] - energies[1]) < 0.0


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, "..", ".."))
using ZigZagBoomerang

Pkg.activate(@__DIR__)
cd(@__DIR__)
#using Revise
using DataStructures
using LinearAlgebra
using Random
using SparseArrays
using Test
using FileIO
using Statistics
using Makie, AbstractPlotting
using ForwardDiff
using StructArrays
const ρ0 = 0.0
Σ = [1 ρ0; ρ0 1]
const d = 2
const dist = 1.5

phi(x, y, rho) =  1/(2*pi*sqrt(1-rho^2))*exp(-0.5*(x^2 + y^2 - 2x*y*rho)/(1-rho^2))
logdensity(x, y) = log(phi(x - dist, y - dist, 0.0) + phi(x + dist, y + dist, 0.0))
#r = -4.5:0.01:4.5

#surf = surface(r, r, [10*exp(logdensity(x1, x2)) for x1 in r, x2 in r])
#=
function partiali()
    ith = zeros(d)
    function (x,i)
        ith[i] = 1
        sa = StructArray{ForwardDiff.Dual{}}((x, ith))
        δ = density(sa).partials[]
        ith[i] = 0
        return δ
    end
end
=#

function gradϕ(x,i)
    -ForwardDiff.partials(logdensity(ForwardDiff.Dual{}(x[1], 1.0*(i==1)), ForwardDiff.Dual{}(x[2], 1.0*(i==2))))[]
end

function gradϕ!(y, x)
    y[1] = gradϕ(x,1)
    y[2] = gradϕ(x,2)
    y
end

κ = [0.1, 0.1]
t0 = 0.0
x0 = randn(d)
θ0 = [1.0, 1.0]
c = 20.0
T = 1000.0

#sspdmp(∇ϕi, t0, x0, θ0, T, c*ones(d), ZigZag(sparse(I(d)), 0*x0), κ; adapt=false)
trace, (tT, xT, θT), (acc, num) = sspdmp(gradϕ, t0, x0, θ0, T, c*ones(d), ZigZag(sparse(I(d)), 0*x0), κ; adapt=false)
#ts, xs = ZigZagBoomerang.sep(collect(discretize(trace, 0.05)))
ts, xs = ZigZagBoomerang.sep(collect(trace))

#p1 = scatter(first.(xs), last.(xs), color=ts, markersize = 3, strokewidth = 0.1, alpha = 0.5)
p1 = lines(first.(xs), last.(xs), color=ts)

scene, layout = layoutscene(resolution = (1200, 900))
layout[1, 1] = ax1 = Axis(scene)
layout[2, 1] = ax2 = Axis(scene)

linkyaxes!(ax1, ax2)
linkxaxes!(ax1, ax2)

lines!(ax1, ts, first.(xs))
lines!(ax2, ts, last.(xs))
p1a = scene
trace, (tT, xT, θT), (acc, num) = sspdmp(gradϕ!, t0, x0, θ0, T, c, BouncyParticle(sparse(I(d)), 0*x0, 0.1), κ; adapt=false)

#ts, xs = ZigZagBoomerang.sep(collect(discretize(trace, 0.05)))
ts, xs = ZigZagBoomerang.sep(collect(trace))

p2 = lines(first.(xs), last.(xs), color=ts)

scene, layout = layoutscene(resolution = (1200, 900))
layout[1, 1] = ax1 = Axis(scene)
layout[2, 1] = ax2 = Axis(scene)

linkyaxes!(ax1, ax2)
linkxaxes!(ax1, ax2)

lines!(ax1, ts, first.(xs))
lines!(ax2, ts, last.(xs))
p2a = scene

Random.seed!(4)
BB = Boomerang(sparse(I(d)), 0*x0, 1.0)
trace, (tT, xT, θT), (acc, num) = sspdmp(gradϕ!, t0, x0, θ0, T, c, BB, κ; adapt=false)

tsxs = sdiscretize(trace.events, BB, 0.05)

ts, xs= tsxs.t, tsxs.x

p3 = scatter(first.(xs), last.(xs), color=ts, markersize = 3, strokewidth = 0.1, alpha = 0.5)

scene, layout = layoutscene(resolution = (1200, 900))
    layout[1, 1] = ax1 = Axis(scene)
    layout[2, 1] = ax2 = Axis(scene)

    linkyaxes!(ax1, ax2)
    linkxaxes!(ax1, ax2)

    lines!(ax1, ts, first.(xs))
    lines!(ax2, ts, last.(xs))
    p3a = scene

save("zigzagphase.png", p1)
save("zigzagtrace.png", p1a)

save("bouncyphase.png", p2)
save("bouncytrace.png", p2a)

save("boomyphase.png", p3)
save("boomytrace.png", p3a)


text:
using MLJ
import RDatasets: dataset
using PrettyPrinting
using Random

Random.seed!(3203)
X = randn(20, 2)
y = vcat(-ones(10), ones(10))

using PyPlot
figure(figsize=(8,6))

ym1 = y .== -1
ym2 = .!ym1
plot(X[ym1, 1], X[ym1, 2], ls="none", marker="o")
plot(X[ym2, 1], X[ym2, 2], ls="none", marker="x")

X = MLJ.table(X)
y = categorical(y);

SVC = @load SVC pkg=LIBSVM

svc_mdl = SVC()
svc = machine(svc_mdl, X, y)

fit!(svc);

ypred = MLJ.predict(svc, X)
misclassification_rate(ypred, y)

rc = range(svc_mdl, :cost, lower=0.1, upper=5)
tm = TunedModel(model=svc_mdl, ranges=[rc], tuning=Grid(resolution=10),
                resampling=CV(nfolds=3, rng=33), measure=misclassification_rate)
mtm = machine(tm, X, y)

fit!(mtm)

ypred = MLJ.predict(mtm, X)
misclassification_rate(ypred, y)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
@testset "Test symmetric version of pressure for point vortices" begin
    atol = 1000*eps()
    Nv = 10
    zv = randn(Nv) + im *(0.05 .+ rand(Nv))
    Sv = randn(Nv)

    points₊ = Vortex.Point.(zv, Sv)
    points₋ = Vortex.Point.(conj(zv), -Sv)

    sys = deepcopy(vcat(points₊, points₋))

    xsensors = collect(-2.0:0.1:2.0)
    sensors = complex(xsensors)
    Ny = size(xsensors, 1)

    press_sym = symmetric_pressure(xsensors, points₊, 0.0)

    press = pressure_AD(sensors, deepcopy(sys), 0.0)
    @test isapprox(press_sym, press, atol = atol)
end

@testset "Test symmetric version of pressure for point vortices with real freestream" begin
    atol = 1000*eps()

    U = randn(ComplexF64)
    freestream = Freestream(real(U))
    Nv = 10
    zv = randn(Nv) + im *(0.05 .+ rand(Nv))
    Sv = randn(Nv)

    points₊ = Vortex.Point.(zv, Sv)
    points₋ = Vortex.Point.(conj(zv), -Sv)

    sys = deepcopy(vcat(points₊, points₋))

    xsensors = collect(-2.0:0.1:2.0)
    sensors = complex(xsensors)
    Ny = size(xsensors, 1)

    press_sym = symmetric_pressure(xsensors, points₊, freestream, 0.0)

    press_AD = pressure_AD(sensors, deepcopy(sys), freestream, 0.0)
    press = pressure(sensors, deepcopy(sys), freestream, 0.0)

    @test isapprox(press_AD, press, atol = atol)
    @test isapprox(press_sym, press, atol = atol)

end

@testset "Test symmetric version of pressure for regularized point vortices" begin
    atol = 1000*eps()
    Nv = 10
    zv = randn(Nv) + im *rand(Nv)
    Sv = randn(Nv)
    δ = 0.5

    blobs₊ = Vortex.Blob.(zv, Sv, δ*ones(Nv))
    blobs₋ = Vortex.Blob.(conj(zv), -Sv, δ*ones(Nv))

    sys = deepcopy(vcat(blobs₊, blobs₋))

    xsensors = collect(-2.0:0.1:2.0)
    sensors = complex(xsensors)
    Ny = size(xsensors, 1)

    press_sym = symmetric_pressure(xsensors, blobs₊, 0.0)

    press = pressure_AD(sensors, deepcopy(sys), 0.0)
    @test isapprox(press_sym, press, atol = atol)
end

@testset "Test symmetric version of pressure for regularized point vortices with real freestream" begin
    atol = 1000*eps()

    U = randn(ComplexF64)
    freestream = Freestream(real(U))
    Nv = 10
    zv = randn(Nv) + im *rand(Nv)
    Sv = randn(Nv)
    δ = 0.5

    blobs₊ = Vortex.Blob.(zv, Sv, δ*ones(Nv))
    blobs₋ = Vortex.Blob.(conj(zv), -Sv, δ*ones(Nv))

    sys = deepcopy(vcat(blobs₊, blobs₋))

    xsensors = collect(-2.0:0.1:2.0)
    sensors = complex(xsensors)
    Ny = size(xsensors, 1)

    press_sym = symmetric_pressure(xsensors, blobs₊, freestream, 0.0)

    press = pressure(sensors, deepcopy(sys), freestream, 0.0)
    press_AD = pressure_AD(sensors, deepcopy(sys), freestream, 0.0)

    @test isapprox(press_sym, press_AD, atol = atol)
    @test isapprox(press_sym, press, atol = atol)
end

@testset "Test measure_state_symmetric and measure_state" begin
    Nv = 15
    Nx = 3*Nv
    x = rand(Nx)

    atol = 1000*eps()

    xsensors = collect(-2.0:0.5:10)
    sensors = complex(xsensors)
    U = complex(0.0)

    config = let Nv = Nv,
         U = U,
         ss = sensors, Δt = 1e-3, δ = 5e-2,
         ϵX = 1e-3, ϵΓ = 1e-4,
         β = 1.0,
         ϵY = 1e-16
         VortexConfig(Nv, U, ss, Δt, δ, ϵX, ϵΓ, β, ϵY)
    end

    sys = state_to_lagrange(x, config)
    press_truth = pressure(sensors, sys, 0.0)

    press = measure_state(x, 0.0, config)

    press_symmetric = measure_state_symmetric(x, 0.0, config)

    @test isapprox(press, press_truth, atol = atol)
    @test isapprox(press_symmetric, press_truth, atol = atol)
end

@testset "Test measure_state_symmetric with real freestream" begin
    Nv = 15
    Nx = 3*Nv
    x = rand(Nx)

    atol = 1000*eps()

    xsensors = collect(-2.0:0.5:10)
    sensors = complex(xsensors)

    U = randn() + 0.0*im
    freestream = Freestream(U)

    config = let Nv = Nv,
         U = U,
         ss = sensors, Δt = 1e-3, δ = 5e-2,
         ϵX = 1e-3, ϵΓ = 1e-4,
         β = 1.0,
         ϵY = 1e-16
         VortexConfig(Nv, U, ss, Δt, δ, ϵX, ϵΓ, β, ϵY)
    end

    sys = state_to_lagrange(x, config)

    press_symmetric = measure_state_symmetric(x, 0.0, config; withfreestream = true)
    press_sym = symmetric_pressure(real.(sensors), sys[1], freestream, 0.0)
    press_truth = pressure(sensors, sys, freestream, 0.0)

    @test isapprox(press_symmetric, press_truth, atol = atol)
    @test isapprox(press_sym, press_truth, atol = atol)
end


text:
len(m::single_part_motif) = m.len;
# view(m::single_part_motif) = @pt :header = DNA_char_cap view(m.P);
sample_m(m::single_part_motif)::String = join([DNA_char_cap[i] for i in rand(m.P)]);
sample_whole(m::single_part_motif) = (motif=sample_m(m), mode=1);

function sample_contiguous_subset(m::single_part_motif, r::UnitRange{Int})::String
    if maximum(r) > m.len || minimum(r) < 1
        error("unit range specified is not covered by 1:$(m.len)");
    end
    sample_m(m)[r] # make it Vector{String} for concatenation purposes later (sample fake data)
end

view(m::k_parts_motif) = for v in m.P view(v) end;
sample_m(m::k_parts_motif)::String = join([sample_m(v) for v in m.P]);
sample_vec(m::k_parts_motif) = [sample_m(v) for v in m.P];
sample_whole(m::k_parts_motif) = (motif=sample_m(m), mode=1);

# distribution on the number of gap alphbet for now is assumed to be discrete uniform (for now)
function sample_m(m::gapped_k_parts_motif)::String
    s = [];
    parts = sample_vec(m.P_motifs);
    for i = 1:(m.K-1)
        push!(s, parts[i]);
        num_alphabet = rand(DiscreteUniform(0, m.gap_len[i]));
        push!(s, join([DNA_char_nocap[rand(m.P_gap)] for j = 1:num_alphabet]));
    end
    push!(s,parts[m.K]);
    return join(s);
end

sample_whole(m::gapped_k_parts_motif) = (motif=sample_m(m), mode=1);

"""
sample a (contiguous) subset of gapped_p_fam_k_motif_PC
according to a unitrange r
"""

function sample_contiguous_subset(m::gapped_k_parts_motif, r::UnitRange{Int})::String
    r_max = maximum(r); r_min = minimum(r);
    if r_max > m.K || r_min < 1
        error("unit range specified is not covered by 1:$(m.len)");
    end
    s = [];
    parts = sample_vec(m.P_motifs);
    for i = r_min:(r_max-1)
        push!(s, parts[i]);
        num_alphabet = rand(DiscreteUniform(0,m.gap_len[i]));
        push!(s, join([DNA_char_nocap[rand(m.P_gap)] for j = 1:num_alphabet]));
    end
    push!(s, parts[r_max]);
    return join(s);
end

view(m::gapped_k_parts_motif) = view(m.P_motifs);

# example ----------------------------------
# a = gapped_k_parts_motif([3,3,3],[2,2]);
# sample_m(a)
# view(a)
# ------------------------------------------

# Sample the ith mode of the motif, a gapped_p_fam_k_motif

function sample_m(m::mixture_gapped_k_parts_motifs)::String
    i = rand(m.mixture_weights);
    sample_contiguous_subset(m.motif, m.modes[i]);
end

function sample_whole(m::mixture_gapped_k_parts_motifs)
    i = rand(m.mixture_weights);
    str = sample_contiguous_subset(m.motif, m.modes[i]);
    (motif=str, mode=i);
end

# example ----------------------------------
# a=mixture_gapped_k_parts_motifs([7,7,7],[4,4],[1:2,2:3]);
# a.modes
# a.num_modes
# a.mixture_weights.p
# rand(a.mixture_weights)
# sample_m(a)

# b = mixture_gapped_k_parts_motifs([1:9, 6:12, 8:15], "uniform");
# b.modes
# b.num_modes
# b.mixture_weights.p
# rand(b.mixture_weights)
# sample_m(b)
# ------------------------------------------

"""
Sample the ith mode of the motif, a mixture_p_fam_motifs
"""
function sample_m(m::mixture_k_parts_motifs)::String
    i = rand(m.mixture_weights);
    sample_contiguous_subset(m.motif, m.modes[i]);
end

function sample_whole(m::mixture_k_parts_motifs)
    i = rand(m.mixture_weights);
    str = sample_contiguous_subset(m.motif, m.modes[i]);
    # println("mode $i");
    (motif=str, mode=i);
end

# example ----------------------------------
# a = mixture_k_parts_motifs([1:9, 3:8, 6:15]);
# a.modes
# a.num_modes
# a.mixture_weights.p
# rand(a.mixture_weights)
# sample_m(a)
# sample_whole(a)

# b = mixture_gapped_p_fam_k_motifs([1:7, 6:8, 6:10], "uniform");
# b.modes
# b.num_modes
# b.mixture_weights.p
# rand(b.mixture_weights)
# sample_m(b)
# ------------------------------------------

# Sample a background-string of length len (default: len = data_pt_len)
sample_background_single(len::Int = data_pt_len)::String =
        join([DNA_char_nocap[rand(background)] for _ in 1:len]);

#=
sample a background string of length len (default: len = data_pt_len)
with a motif embedded in it;

returns a named-tuple:
(str=dna_string, motif_where=motif_start:motif_end, mode=binding_mode),
which has a motif-embedded background string , an UnitRange that
indicates where motif is embedded, and the mode of binding
=#

function sample_backgound_with_motif_single(motif::motif_type,
                        len::Int=data_pt_len
                        )::sim_dna_str_w_motif

    motif_realization = sample_whole(motif);

    ########### either make it complement or not ###########
    motif_str = motif_realization.motif; motif_mode = motif_realization.mode;
    is_complement = rand() > percent_complement;
    motif_realization = is_complement ? (motif=reverse_complement(motif_str),mode=motif_mode) : motif_realization;
    ########################################################

    motif_len = length(motif_realization.motif);
    if motif_len > data_pt_len
        error("Please use a motif that has width shorter than $data_pt_len");
    end

    avail_positions::Int = (len-motif_len+1)-sample_near_center;
    at::Int = rand(DiscreteUniform(1+sample_near_center, avail_positions));
    before_at::Int = (at-1);
    dna_str = join([sample_background_single(before_at),
        motif_realization.motif::String,
        sample_background_single((data_pt_len-(before_at+motif_len)))]);

    return  (str=dna_str,
             motif_where=at:(at+motif_len-1),
             mode=motif_realization.mode,
             complement=is_complement);
end

# example ----------------------------------
# a = gapped_k_parts_motif([3,3,3],[2,2]);
# sample_a = sample_backgound_with_motif_single(a);
# sample_a
# sample_a.mode
# sample_a[1][sample_a.motif_where]

# b = mixture_gapped_k_parts_motifs([3,3,3],[2,2],[1:2,2:3]);
# sample_b = sample_backgound_with_motif_single(b);
# sample_b
# sample_b.mode
# sample_b[1][sample_b.motif_where]
# ------------------------------------------

#=
sample multiple background string of length len (default: len = data_pt_len)
with a motif embedded in it; see sample_backgound_with_motif_single.
Instead of returning an instance of sim_dna_str_w_motif, it returns an
instance of Vector{sim_dna_str_w_motif}.
=#

function sample_backgound_with_motif_multiple(
                            motif::motif_type,
                            how_many::Int,
                            len::Int=data_pt_len
                            )::Vector{sim_dna_str_w_motif}
    return [sample_backgound_with_motif_single(motif, len)
        for _ in 1:how_many];
end

# example ----------------------------------
# b = mixture_gapped_k_parts_motifs([3,3,3],[2,2],[1:2,2:3]);
# sample_b = sample_backgound_with_motif_multiple(b, 5);
# sample_b
# ------------------------------------------


text:
using QInterpolations, Test

@testset "QInterpoations" begin
x = range(0,stop=10,length=100)
y1 = Array(x) + 1.0im*Array(x)
y2 = (10.0+10.0im) .- Array(x)
inter_y1 = construct_interpolations(x,y1)
test_x = range(0,stop=10,length=50)
exp_y1 = Array(test_x) + 1.0im*Array(test_x)
exp_y2 = (10.0+10.0im) .- Array(test_x)
res_y1 = [inter_y1(i) for i in test_x]
@test isapprox(exp_y1, res_y1)
# Test for complex vector interpolation
y_vec_test = [[y1[i],y2[i]] for i in eachindex(y1)]
inter_yvec = construct_interpolations(x, y_vec_test)
res_y_vec = [inter_yvec(i) for i in test_x]
exp_y_vec = [[exp_y1[i],exp_y2[i]] for i in eachindex(exp_y1)]
@test isapprox(res_y_vec, exp_y_vec)
# Test for complex matrix interolation
y_mat_test = [[y1[i] y2[i]; y2[i] y1[i]] for i in eachindex(y1)]
inter_ymat = construct_interpolations(x, y_mat_test)
res_y_mat = [inter_ymat(i) for i in test_x]
exp_y_mat = [[exp_y1[i] exp_y2[i]; exp_y2[i] exp_y1[i]] for i in eachindex(exp_y1)]
@test isapprox(res_y_mat, exp_y_mat)
# Test for complex vector gradient
res_g_yvec = gradient(inter_yvec, 1.0)
@test isapprox(res_g_yvec, [1.0+1.0im,-1], atol=1e-6)
# Test for real number extrapolation
x = range(1.0,stop=10.0)
eitp = construct_interpolations(x, Array(x) , extrapolation="line")
@test isapprox(eitp(0.0),0.0, atol=1e-8)
end

max_stars_repo_name: neversakura/QInterpolate

max_issues_repo_name: neversakura/QInterpolate

max_forks_repo_name: neversakura/QInterpolate


text:
"""

    XYZ = L63(Float32)

Integrates the Lorenz 1963 model in any given DataType. Default parameters

    L63(::Type{T}=Float64;                      # number format
    N::Int=10_000,                              # number of time steps
    xyz::Array{Float64,1}=[5.0,5.0,20.0],       # initial conditions
    σ::Float64=10.0,                            # σ parameter
    ρ::Float64=28.0,                            # ρ parameter
    β::Float64=8/3,                             # β parameter
    s::Float64=1.0,                             # scaling
    Δt::Float64=0.002,                          # time step
    scheme::String="RK4"                        # time integration scheme

# Examples
```jldoc
julia> X1 = L63(Float64,N=100_000);
julia> X2 = L63(Float32,Δt=0.05);
```
"""
function L63(::Type{T}=Float64;
            N::Int=10_000,
            xyz::Array{Float64,1}=[5.0,5.0,20.0],
            σ::Float64=10.0,
            ρ::Float64=28.0,
            β::Float64=8/3,
            s::Float64=1.0,
            Δt::Float64=0.002,
            scheme::String="RK4") where {T<:AbstractFloat}

            if scheme == "RK4"
                return RK4(T,N,xyz,σ,ρ,β,s,Δt)
            else
                throw(error("Other schemes than RK4 not implemented yet."))
            end
end


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

mutable struct RegStruct

    boxlo::Array{Float64,2};
    boxhi::Array{Float64,2};
    boxtilt::Array{Float64,2};
    triclinic::Int64;

    RegStruct() = new();
end

function setregion(boxhi, boxtilt=[0.0 0.0 0.0])

reg = RegStruct();

reg.boxlo = [0.0 0.0 0.0];
reg.boxhi = boxhi;
reg.boxtilt = boxtilt;

if (reg.boxtilt[1] == 0 && reg.boxtilt[2] == 0 && reg.boxtilt[3] == 0)
    reg.triclinic = 0;
else
    reg.triclinic = 1;
end

return reg

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
module MKLsparseinspector

using Base.Enums
using LinearAlgebra
using MKL_jll
using SparseArrays

using LinearAlgebra: BlasInt

function __init__()
    ccall((:MKL_Set_Interface_Layer, libmkl_rt), Cint, (Cint,), Base.USE_BLAS64 ? 1 : 0)
end

include("enums.jl")
include("types.jl")
include("toplevel.jl")

export
    MKLcsr,
    MKLFloats,

    csrptr,
    syrkd!

end


text:
using ExperimentalDesign, DataStructures

ranges = [k:k for k = 20:1:24]

factors = OrderedDict{Symbol, Any}()

for i = 1:10
    factors[Symbol(:f, i)] = [:a, :b, :c]
end

factor_list = [factors for i in ranges]

designs = 4000

sampled_subsets = sample_subsets(factor_list, ranges, designs, check_bounds = true)


text:
using PureSeq
using Base.Test

## BamReader

# test forward reads
reader = BamReader("data/small.bam", :forward, ReferenceContigs_hg38)
@test position(reader) == 10544
advance!(reader)
@test position(reader) == 11247
close(reader)

# test reverse reads
reader = BamReader("data/small.bam", :reverse, ReferenceContigs_hg38)
@test position(reader) == 10056
advance!(reader)
@test position(reader) == 10102
close(reader)

# test reading through a whole file
reader = BamReader("data/small.bam", :any, ReferenceContigs_hg38)
lastPos = -1
while !eof(reader)
	lastPos = position(reader)
	advance!(reader)
end
close(reader)

# test reading through a whole file with an iterator
reader = BamReader("data/small.bam", :reverse, ReferenceContigs_hg38)
for pos in eachposition(reader)
end
close(reader)

## FeatureMap

# test forward reads
reader = BamReader("data/small.bam", :forward, ReferenceContigs_hg38)
fm = FeatureMap(reader, ones(2001))
@test position(fm) == 9544 # where data first enters the window
@test value(fm) == 1.0
while !eof(fm)
	if position(fm) == 10544
		@test value(fm) == 3.0
	end
	advance!(fm)
end
close(fm)

## DenseBlockIterator

# test forward reads
fm1 = FeatureMap(BamReader("data/small.bam", :forward, ReferenceContigs_hg38), ones(21)*2)
fm2 = FeatureMap(BamReader("data/small.bam", :forward, ReferenceContigs_hg38), ones(21))
count = 0
for block in denseblocks([fm1, fm2], 100)
	if count == 105
		@test block[33,1] == 0.0
		@test block[33,2] == 0.0
		@test block[34,1] == 2.0
		@test block[34,2] == 1.0
		@test block[54,1] == 2.0
		@test block[54,2] == 1.0
		@test block[55,1] == 0.0
		@test block[55,2] == 0.0
	end
	count += 1
	if count > 120 break end
end
close(fm1)
close(fm2)

# test infinite looping of blocks
fm1 = FeatureMap(BamReader("data/small.bam", :forward, ReferenceContigs_hg38), ones(21)*2)
fm2 = FeatureMap(BamReader("data/small.bam", :forward, ReferenceContigs_hg38), ones(21))
count = 0
for block in denseblocks([fm1, fm2], 1000, loop=true)
	count += 1
	if count > 12000 break end
end
@test count == 12001
close(fm1)
close(fm2)

## SamWriter

f = open("data/tmp.sam", "w")
sw = SamWriter(f, ReferenceContigs_hg38)
writeRead(sw, 10, 16)
close(f)

## write_binned and BinnedReader

# test forward reads
write_binned("data/small.bam", 1000, :reverse, skipDup=false)
reader = BinnedReader("data/small.bam.rbin1000")
@test position(reader) == 11 # where data first enters the window
@test value(reader) == 3.0

while !eof(reader)
	if position(reader) == 12
		@test value(reader) == 5.0
	end
	advance!(reader)
end

# now skip duplicate reads
write_binned("data/small.bam", 1000, :reverse, skipDup=true)
reader = BinnedReader("data/small.bam.rbin1000")
@test position(reader) == 11 # where data first enters the window
@test value(reader) == 3.0
while !eof(reader)
	if position(reader) == 12
		@test value(reader) == 4.0
	end
	advance!(reader)
end
close(reader)

## ContextMap

# test forward reads
reader = BamReader("data/small.bam", :forward, ReferenceContigs_hg38)
cm = ContextMap(reader, 1000, 1000)
@test position(cm) == 9544 # where data first enters the window
@test value(cm) == 1.0
while !eof(cm)
	if position(cm) == 10544
		@test value(cm) == 3.0
	end
	advance!(cm)
end
close(cm)

# test bins
reader = BinnedReader("data/small.bam.fbin1000")
cm = ContextMap(reader, 2, 2)
@test position(cm) == 9 # where data first enters the window
@test value(cm) == 1.0
while !eof(cm)
	if position(cm) == 12
		@test value(cm) == 13.0
	end
	advance!(cm)
end
close(cm)