text:
# Medical system utility
# There are three kinds of diseases, each comes in two types
#
#  *  Heart disease
#    *  h1
#    *  h2
#  *  Lung disease 1
#    *  l1
#    *  l2
#  *  Lung disease 2
#    *  l3
#    *  l4
#
#
# For each specific type, there is a specific treatment that works best (=highest utility)
# Applying the specific treatment to the wrong type of the same disease is less effective
# Applying the specific treatment to the wrong disease is bad (low utility)
#
# There are also general treatments for each disease: ``treat h``, ``treat l12``, ``treat l34``
# The general treatement works for both types of the corresponding disease equally well,
# but has no effect when applied to the wrong disease.
# The general treatement for the heart is less effective (lower utility) than the general treatments
# for lung disease 1 or lung disease 2.
#
# There is one general lung treatement ``treat l`` that has an effect when applied to ``l1, l2, l3`` or ``l4``
# However, it is less efficient than the general treatments for the corresponding lung disease ``treat l12`` or ``treat l34``
# The general lung treatment has no effect when applied to a heart disease
#
#
#
# The function ``setup_medical_example`` has an optional argument that allows to choose between
#   -) a uniform environment where all diseases appear with equal probability and
#   -) a nonuniform environment where both heart disease types ``h1, h2`` appear with inceased chance.
#
#
function setup_medical_example(;uniform_w=true)

    w_strings = ["h1","h2", "l1","l2", "l3", "l4"]
    w_values = collect(1:length(w_strings))

    numw = length(w_values)

    if uniform_w
        #uniform distribution over w's
        pw = ones(numw)/numw
    else
        #non-uniform p(w)
        pw = medical_nonuniform_pw(numw)
    end

    a_strings = [ map((x)->("treat "*w_strings[x]), w_values);
                  "treat l12"; "treat l34"; "treat h"; "treat l" ]
    a_values = collect(1:length(a_strings))

    function U(a,w)
        #correct treatment
        correct_utility = 3

        #wrong treatment for cause: heard-disease
        wrong_utility_heart = correct_utility * 0.3

        #general heart-treatment
        utility_heart_general = 1.5

        #wrong treatment for cause: lung-disease
        wrong_utility_lung = correct_utility * 0.5

        #wrong treatment for cause: lung-disease
        wrong_utility_lung_2 = correct_utility * 0

        #general lung treatment
        utility_lung_general = 1.5

        #general lung treatment l12
        utiliy_lung_general_12 = 2.5

        #general lung treatment l12
        utiliy_lung_general_34 = 2.5

        #correct treatment
        if a==w
            return correct_utility
        end

        #heart-disease, heart treatment (but not correct one)
        if w<3 && a<3
            return wrong_utility_heart
        end

        #lung-disease, lung treatment (but not correct one)
        if w>2 && w<5 && a>2 && a<5
            return wrong_utility_lung
        end

        #lung-disease, lung treatment (but not correct one)
        if w>2 && w<5 && a>4 && a<7
            return wrong_utility_lung_2
        end

        #lung-disease, lung treatment (but not correct one)
        if w>4 && w<7 && a>4 && a<7
            return wrong_utility_lung
        end

        #lung-disease, lung treatment (but not correct one)
        if w>4 && w<7 && a>2 && a<5
            return wrong_utility_lung_2
        end

        #general heart treatment
        if w<3 && a==9
            return utility_heart_general
        end

        #general lung treatments
        if w>2 && w<7
            if a==7 && w<5
                return utiliy_lung_general_12
            end

            if a==8 && w>4
                return utiliy_lung_general_34
            end

            if a==10
                return utility_lung_general
            end
        end

        #wrong treatment for wrong cause
        return 0

    end

    return w_values, w_strings, a_values, a_strings, pw, U

end

function medical_nonuniform_pw(num_w_values)
        #increased probability of one heart-disease and pancreatic disease
        pw = ones(num_w_values)
        pw[1:2] = 3
        #pw[3:6] = 3
        #pw[6:7] = 3
        #pw[8:9] = 1
        pw /= sum(pw)  #re-normalize
    return pw
end


text:
module PhotometricFilters

using DataDeps
using DelimitedFiles

export PhotometricFilter,
       wave,
       throughput,
       apply,
       apply!

include("core.jl") # types and base utilities (like size, length)
include("library.jl") # the library of filters
include("plots.jl") # plotting recipes
include("vega.jl") # vega calibration spectrum

function __init__()
    register(PYPHOT_DATADEP)
    register(VEGA_DATADEP)
end

end

ext: jl


text:
#
# The type of the query input.
#

# InputSlots with their types.

typealias InputSlot Pair{Symbol, Output}
typealias InputSlots Vector{InputSlot}

const NO_SLOTS = InputSlot[]

# Input context.

immutable InputMode
    # Depends on the preceeding and subsequent input values.
    relative::Bool
    # List of query slots.
    slots::Vector{InputSlot}
end

InputMode() = InputMode(false, NO_SLOTS)

setrelative(imode::InputMode, relative::Bool=true) =
    InputMode(relative, imode.slots)

setslots(imode::InputMode, slots::InputSlots=NO_SLOTS) =
    InputMode(imode.relative, slots)

setrelative(relative::Bool=true) =
    imode -> setrelative(imode, relative)

setslots(slots::InputSlots=NO_SLOTS) =
    imode -> setslots(imode, slots)

# Predicates and properties.

isfree(imode::InputMode) = !imode.relative && isempty(imode.slots)
isrelative(imode::InputMode) = imode.relative
slots(imode::InputMode) = imode.slots

function slot(imode::InputMode, tag::Symbol)
    for param in imode.slots
        if param.first == tag
            return param.second
        end
    end
    return Output(None)
end

# The type and the context of the input.

immutable Input
    dom::Domain
    mode::InputMode
end

Input(dom) = Input(convert(Domain, dom), InputMode())

decorate(ity::Input, d::Decoration) =
    Input(decorate(ity.dom, d), ity.mode)

setrelative(ity::Input, relative::Bool=true) =
    Input(ity.dom, setrelative(ity.mode, relative))

setslots(ity::Input, slots::InputSlots=NO_SLOTS) =
    Input(ity.dom, setslots(ity.mode, slots))

convert(::Type{Input}, dom::Union{Type, Symbol, Tuple, Vector, Domain}) =
    Input(convert(Domain, dom))

function show(io::IO, ity::Input)
    mode = ity.mode
    if !isempty(mode.slots)
        print(io, "{")
    end
    if mode.relative
        print(io, "(")
    end
    print(io, ity.dom)
    if mode.relative
        print(io, "...)")
    end
    for (n, o) in mode.slots
        print(io, ", $n => $o")
    end
    if !isempty(mode.slots)
        print(io, "}")
    end
end

# Predicates and properties.

domain(ity::Input) = ity.dom
mode(ity::Input) = ity.mode

isdata(ity::Input) = isdata(ity.dom)
isany(ity::Input) = isany(ity.dom)
isvoid(ity::Input) = isvoid(ity.dom)
isnone(ity::Input) = isnone(ity.dom)
isentity(ity::Input) = isentity(ity.dom)
isrecord(ity::Input) = isentity(ity.dom)

isfree(ity::Input) = isfree(ity.mode)
isrelative(ity::Input) = isrelative(ity.mode)
slots(ity::Input) = slots(ity.mode)
slot(ity::Input, tag::Symbol) = slot(ity, tag)

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
module ParameterDistributions

## Usings
using Distributions
using Statistics
using Random

#import (to add definitions)
import StatsBase: mean, var, cov, sample
import Base: size, length, ndims
## Exports

#types
export ParameterDistributionType
export ConstraintType

#objects
export Parameterized, Samples, VectorOfParameterized
export ParameterDistribution
export Constraint

#functions
export get_name, get_distribution, ndims, get_dimensions, get_all_constraints, get_n_samples
export sample
export no_constraint, bounded_below, bounded_above, bounded
export transform_constrained_to_unconstrained, transform_unconstrained_to_constrained
export get_logpdf, batch

export combine_distributions
## Objects
# for the Distribution
abstract type ParameterDistributionType end

"""
    Parameterized <: ParameterDistributionType

A distribution constructed from a parametrized formula (e.g Julia Distributions.jl)
"""
struct Parameterized <: ParameterDistributionType
    distribution::Distribution
end

"""
    Samples{FT <: Real} <: ParameterDistributionType

A distribution comprised of only samples, stored as columns of parameters.
"""
struct Samples{FT <: Real} <: ParameterDistributionType
    distribution_samples::AbstractMatrix{FT} #parameters are columns
    Samples(distribution_samples::AbstractMatrix{FT}; params_are_columns = true) where {FT <: Real} =
        params_are_columns ? new{FT}(distribution_samples) : new{FT}(permutedims(distribution_samples, (2, 1)))
    #Distinguish 1 sample of an ND parameter or N samples of 1D parameter, and store as 2D array
    Samples(distribution_samples::AbstractVector{FT}; params_are_columns = true) where {FT <: Real} =
        params_are_columns ? new{FT}(reshape(distribution_samples, 1, :)) : new{FT}(reshape(distribution_samples, :, 1))
end

"""
    VectorOfParameterized <: ParameterDistributionType

A distribution built from an array of Parametrized distributions.
A utility to help stacking of distributions where a multivariate equivalent doesn't exist
"""
struct VectorOfParameterized{DT <: Distribution} <: ParameterDistributionType
    distribution::AbstractVector{DT}
end

# For the transforms
abstract type ConstraintType end

"""
    Constraint <: ConstraintType

Contains two functions to map between constrained and unconstrained spaces.
"""
struct Constraint <: ConstraintType
    constrained_to_unconstrained::Function
    unconstrained_to_constrained::Function
end

"""
    no_constraint()

Constructs a Constraint with no constraints, enforced by maps x -> x and x -> x.
"""
function no_constraint()
    c_to_u = (x -> x)
    u_to_c = (x -> x)
    return Constraint(c_to_u, u_to_c)
end

"""
    bounded_below(lower_bound::FT) where {FT <: Real}

Constructs a Constraint with provided lower bound, enforced by maps `x -> log(x - lower_bound)`
and `x -> exp(x) + lower_bound`.
"""
function bounded_below(lower_bound::FT) where {FT <: Real}
    c_to_u = (x -> log(x - lower_bound))
    u_to_c = (x -> exp(x) + lower_bound)
    return Constraint(c_to_u, u_to_c)
end

"""
    bounded_above(upper_bound::FT) where {FT <: Real}

Constructs a Constraint with provided upper bound, enforced by maps `x -> log(upper_bound - x)`
and `x -> upper_bound - exp(x)`.
"""
function bounded_above(upper_bound::FT) where {FT <: Real}
    c_to_u = (x -> log(upper_bound - x))
    u_to_c = (x -> upper_bound - exp(x))
    return Constraint(c_to_u, u_to_c)
end

"""
    bounded(lower_bound::FT, upper_bound::FT) where {FT <: Real}

Constructs a Constraint with provided upper and lower bounds, enforced by maps
`x -> log((x - lower_bound) / (upper_bound - x))`
and `x -> (upper_bound * exp(x) + lower_bound) / (exp(x) + 1)`.

"""
function bounded(lower_bound::FT, upper_bound::FT) where {FT <: Real}
    if (upper_bound <= lower_bound)
        throw(DomainError("upper bound must be greater than lower bound"))
    end
    c_to_u = (x -> log((x - lower_bound) / (upper_bound - x)))
    u_to_c = (x -> (upper_bound * exp(x) + lower_bound) / (exp(x) + 1))
    return Constraint(c_to_u, u_to_c)
end

#extending Base.length
"""
    length(c<:ConstraintType)

A constraint has length 1.
"""
length(c::CType) where {CType <: ConstraintType} = length([c])

#extending Base.size
"""
    size(c<:ConstraintType)

A constraint has size 1.
"""
size(c::CType) where {CType <: ConstraintType} = size([c])

"""
    ndims(d<:ParametrizedDistributionType)

The number of dimensions of the parameter space
"""
ndims(d::Parameterized) = length(d.distribution)

ndims(d::Samples) = size(d.distribution_samples, 1)

ndims(d::VectorOfParameterized) = sum(length.(d.distribution))
"""
    n_samples(d<:Samples)

The number of samples in the array.
"""
n_samples(d::Samples) = size(d.distribution_samples)[2]

n_samples(d::Parameterized) = "Distribution stored in Parameterized form, draw samples using `sample` function"

n_samples(d::VectorOfParameterized) = "Distribution stored in Parameterized form, draw samples using `sample` function"

"""
    ParameterDistribution

Structure to hold a parameter distribution, always stored as an array internally.
"""
struct ParameterDistribution{PDType <: ParameterDistributionType, CType <: ConstraintType, ST <: AbstractString}
    distribution::AbstractVector{PDType}
    constraint::AbstractVector{CType}
    name::AbstractVector{ST}
end

"""
    ParameterDistribution(param_dist_dict::Union{Dict,AbstractVector})

Constructor taking in a Dict or array of Dicts. Each dict must contain the key-val pairs:
- `"distribution"` - a distribution of `ParameterDistributionType`
- `"constraint"` - constraint(s) given as a `ConstraintType` or array of `ConstraintType`s with length equal to the dims of the distribution
- `"name"` - a name of the distribution as a String.
"""
function ParameterDistribution(param_dist_dict::Union{Dict, AbstractVector})

    #check type
    if !(isa(param_dist_dict, Dict) || eltype(param_dist_dict) <: Dict)
        throw(ArgumentError("input argument must be a Dict, or <:AbstractVector{Dict}"))
    end

    # make copy as array
    param_dist_dict_array = !isa(param_dist_dict, AbstractVector) ? [param_dist_dict] : param_dist_dict
    # perform checks on the individual distributions
    for pdd in param_dist_dict_array
        # check all keys are present
        if !all(["distribution", "name", "constraint"] .∈ [collect(keys(pdd))])
            throw(
                ArgumentError("input dictionaries must contain the keys: \"distribution\", \"name\", \"constraint\" "),
            )
        end

        distribution = pdd["distribution"]
        name = pdd["name"]
        constraint = pdd["constraint"]

        # check key types
        if !isa(distribution, ParameterDistributionType)
            throw(
                ArgumentError(
                    "Value of \"distribution\" must be a valid ParameterDistribution object: Parameterized or Samples",
                ),
            )
        end
        if !isa(constraint, ConstraintType)
            if !isa(constraint, AbstractVector) #it's not a vector either
                throw(
                    ArgumentError(
                        "Value of \"constraint\" must be a ConstraintType, or <:AbstractVector(ConstraintType)",
                    ),
                )
            elseif !(eltype(constraint) <: ConstraintType) #it is a vector, but not of constraint
                throw(
                    ArgumentError(
                        "Value of \"constraint\" must be a ConstraintType, or <:AbstractVector(ConstraintType)",
                    ),
                )
            end
        end
        if !isa(name, String)
            throw(ArgumentError("Value of \"name\" must be a String"))
        end

        # 1 constraint per dimension check
        constraint_array = isa(constraint, ConstraintType) ? [constraint] : constraint
        n_parameters = ndims(distribution)

        if !(n_parameters == length(constraint_array))
            throw(
                DimensionMismatch(
                    "There must be one constraint dimension in a parameter distribution, use no_constraint() type if no constraint is required",
                ),
            )
        end
    end

    # flatten the structure
    distribution = getindex.(param_dist_dict_array, "distribution")
    flat_constraint = cat(getindex.(param_dist_dict_array, "constraint")..., dims = 1)
    name = getindex.(param_dist_dict_array, "name")

    # build the object
    return ParameterDistribution(distribution, flat_constraint, name)

end

"""
    function ParameterDistribution(distribution::ParameterDistributionType,
                                   constraint::Union{ConstraintType,AbstractVector{ConstraintType}},
                                   name::AbstractString)

constructor of a ParameterDistribution from a single `distribution`, (array of) `constraint`, `name`.
these can used to build another ParameterDistribution
"""
function ParameterDistribution(
    distribution::ParameterDistributionType,
    constraint::Union{ConstraintType, AbstractVector},
    name::AbstractString,
)

    if !(typeof(constraint) <: ConstraintType || eltype(constraint) <: ConstraintType) # if it is a vector, but not of constraint
        throw(ArgumentError("`constraint` must be a ConstraintType, or Vector of ConstraintType's "))
    end
    # 1 constraint per dimension check
    constraint_vec = isa(constraint, ConstraintType) ? [constraint] : constraint
    n_parameters = ndims(distribution)

    if !(n_parameters == length(constraint_vec))
        throw(
            DimensionMismatch(
                "There must be one constraint dimension in a parameter distribution, use no_constraint() type if no constraint is required",
            ),
        )
    end

    # flatten the structure
    distribution_vec = [distribution]
    flat_constraint_vec = cat(constraint_vec..., dims = 1)
    name_vec = [name]

    # build the object
    return ParameterDistribution(distribution_vec, flat_constraint_vec, name_vec)

end

"""
    combine_distributions(pd_vec::AbstractVector{ParameterDistribution}...; dims = 1)

concatenate from vector of single ParameterDistributions, dims keyword is unused
"""
function combine_distributions(pd_vec::AbstractVector{PD}; dims = 1) where {PD <: ParameterDistribution}
    # flatten the structure
    distribution = cat([d.distribution for d in pd_vec]..., dims = 1)
    constraint = cat([d.constraint for d in pd_vec]..., dims = 1)
    name = cat([d.name for d in pd_vec]..., dims = 1)
    return ParameterDistribution(distribution, constraint, name)
end

## Functions

"""
    get_name(pd::ParameterDistribution)

Returns a list of ParameterDistribution names.
"""
get_name(pd::ParameterDistribution) = pd.name

"""
    get_dimensions(pd::ParameterDistribution)

The number of dimensions of the parameter space.
"""
function get_dimensions(pd::ParameterDistribution)
    return [ndims(d) for d in pd.distribution]
end
function get_dimensions(d::VectorOfParameterized)
    return [length(dd) for dd in d.distribution]
end

function ndims(pd::ParameterDistribution)
    return sum(get_dimensions(pd))
end

"""
    get_n_samples(pd::ParameterDistribution)

The number of samples in a Samples distribution
"""
function get_n_samples(pd::ParameterDistribution)
    return Dict{String, Any}(pd.name[i] => n_samples(d) for (i, d) in enumerate(pd.distribution))
end

"""
    get_all_constraints(pd::ParameterDistribution)

Returns the (flattened) array of constraints of the parameter distribution.
"""
get_all_constraints(pd::ParameterDistribution) = pd.constraint

"""
    batch(pd::ParameterDistribution)

Returns a list of contiguous `[collect(1:i), collect(i+1:j),... ]` used to split parameter arrays by distribution dimensions.
"""
function batch(pd::Union{ParameterDistribution, VectorOfParameterized})
    #chunk xarray to give to the different distributions.
    d_dim = get_dimensions(pd) #e.g [4,1,2]
    d_dim_tmp = Array{Int64}(undef, size(d_dim)[1] + 1)
    d_dim_tmp[1] = 0
    for i in 2:(size(d_dim)[1] + 1)
        d_dim_tmp[i] = sum(d_dim[1:(i - 1)]) # e.g [0,4,5,7]
    end

    return [collect((d_dim_tmp[i] + 1):d_dim_tmp[i + 1]) for i in 1:size(d_dim)[1]] # e.g [1:4, 5:5, 6:7]
end

"""
    get_distribution(pd::ParameterDistribution)

Returns a `Dict` of `ParameterDistribution` distributions, with the parameter names
as dictionary keys. For parameters represented by `Samples`, the samples are returned
as a 2D (`parameter_dimension x n_samples`) array.
"""
function get_distribution(pd::ParameterDistribution)
    return Dict{String, Any}(pd.name[i] => get_distribution(d) for (i, d) in enumerate(pd.distribution))
end

get_distribution(d::Samples) = d.distribution_samples

get_distribution(d::Parameterized) = d.distribution

get_distribution(d::VectorOfParameterized) = d.distribution

"""
    sample([rng], pd::ParameterDistribution, [n_draws])

Draws `n_draws` samples from the parameter distributions `pd`. Returns an array, with
parameters as columns. `rng` is optional and defaults to `Random.GLOBAL_RNG`. `n_draws` is
optional and defaults to 1.
"""
function sample(rng::AbstractRNG, pd::ParameterDistribution, n_draws::IT) where {IT <: Integer}
    return cat([sample(rng, d, n_draws) for d in pd.distribution]..., dims = 1)
end

# define methods that dispatch to the above with Random.GLOBAL_RNG as a default value for rng
sample(pd::ParameterDistribution, n_draws::IT) where {IT <: Integer} = sample(Random.GLOBAL_RNG, pd, n_draws)
sample(rng::AbstractRNG, pd::ParameterDistribution) = sample(rng, pd, 1)
sample(pd::ParameterDistribution) = sample(Random.GLOBAL_RNG, pd, 1)

"""
    sample([rng], d::Samples, [n_draws])

Draws `n_draws` samples from the parameter distributions `d`. Returns an array, with
parameters as columns. `rng` is optional and defaults to `Random.GLOBAL_RNG`. `n_draws` is
optional and defaults to 1.
"""
function sample(rng::AbstractRNG, d::Samples, n_draws::IT) where {IT <: Integer}
    n_stored_samples = n_samples(d)
    samples_idx = sample(rng, collect(1:n_stored_samples), n_draws)
    if ndims(d) == 1
        return reshape(d.distribution_samples[:, samples_idx], :, n_draws) #columns are parameters
    else
        return d.distribution_samples[:, samples_idx]
    end
end

# define methods that dispatch to the above with Random.GLOBAL_RNG as a default value for rng
sample(d::Samples, n_draws::IT) where {IT <: Integer} = sample(Random.GLOBAL_RNG, d, n_draws)
sample(rng::AbstractRNG, d::Samples) = sample(rng, d, 1)
sample(d::Samples) = sample(Random.GLOBAL_RNG, d, 1)

"""
    sample([rng], d::Parameterized, [n_draws])

Draws `n_draws` samples from the parameter distributions `d`. Returns an array, with
parameters as columns. `rng` is optional and defaults to `Random.GLOBAL_RNG`. `n_draws` is
optional and defaults to 1.
"""
function sample(rng::AbstractRNG, d::Parameterized, n_draws::IT) where {IT <: Integer}
    if ndims(d) == 1
        return reshape(rand(rng, d.distribution, n_draws), :, n_draws) #columns are parameters
    else
        return rand(rng, d.distribution, n_draws)
    end
end

# define methods that dispatch to the above with Random.GLOBAL_RNG as a default value for rng
sample(d::Parameterized, n_draws::IT) where {IT <: Integer} = sample(Random.GLOBAL_RNG, d, n_draws)
sample(rng::AbstractRNG, d::Parameterized) = sample(rng, d, 1)
sample(d::Parameterized) = sample(Random.GLOBAL_RNG, d, 1)

"""
    sample([rng], d::VectorOfParameterized, [n_draws])

Draws `n_draws` samples from the parameter distributions `d`. Returns an array, with
parameters as columns. `rng` is optional and defaults to `Random.GLOBAL_RNG`. `n_draws` is
optional and defaults to 1.
"""
function sample(rng::AbstractRNG, d::VectorOfParameterized, n_draws::IT) where {IT <: Integer}
    samples = zeros(ndims(d), n_draws)
    batches = batch(d)
    dimensions = get_dimensions(d)
    for (i, dd) in enumerate(d.distribution)
        samples[batches[i], :] = rand(rng, dd, n_draws) #columns are parameters
    end
    return samples
end

sample(d::VectorOfParameterized, n_draws::IT) where {IT <: Integer} = sample(Random.GLOBAL_RNG, d, n_draws)
sample(rng::AbstractRNG, d::VectorOfParameterized) = sample(rng, d, 1)
sample(d::VectorOfParameterized) = sample(Random.GLOBAL_RNG, d, 1)

"""
    logpdf(pd::ParameterDistribution, xarray::Array{<:Real,1})

Obtains the independent logpdfs of the parameter distributions at `xarray`
(non-Samples Distributions only), and returns their sum.
"""
get_logpdf(d::Parameterized, xarray::AbstractVector{FT}) where {FT <: Real} = logpdf.(d.distribution, xarray)

function get_logpdf(d::VectorOfParameterized, xarray::AbstractVector{FT}) where {FT <: Real}
    # get the index of xarray chunks to give to the different distributions.
    batches = batch(d)
    dimensions = get_dimensions(d)
    lpdfsum = 0.0
    # perform the logpdf of each of the distributions, and returns their sum
    for (i, dd) in enumerate(d.distribution)
        if dimensions[i] == 1
            lpdfsum += logpdf.(dd, xarray[batches[i]])[1]
        else
            lpdfsum += logpdf(dd, xarray[batches[i]])
        end
    end
    return lpdfsum
end

function get_logpdf(pd::ParameterDistribution, xarray::AbstractVector{FT}) where {FT <: Real}
    #first check we don't have sampled distribution
    for d in pd.distribution
        if typeof(d) <: Samples
            throw(
                ErrorException(
                    "Cannot compute get_logpdf of Samples distribution. Consider using a Parameterized type for your prior.",
                ),
            )
        end
    end
    #assert xarray correct dim/length
    if size(xarray)[1] != ndims(pd)
        throw(DimensionMismatch("xarray must have dimension equal to the parameter space"))
    end

    # get the index of xarray chunks to give to the different distributions.
    batches = batch(pd)

    # perform the logpdf of each of the distributions, and returns their sum
    return sum(cat([get_logpdf(d, xarray[batches[i]]) for (i, d) in enumerate(pd.distribution)]..., dims = 1))
end

#extending StatsBase cov,var
"""
    var(pd::ParameterDistribution)
Returns a flattened variance of the distributions
"""
var(d::Parameterized) = var(d.distribution)
var(d::Samples) = var(d.distribution_samples, dims = 2)
function var(d::VectorOfParameterized)
    d_dims = get_dimensions(d)
    block_var = Array{Any}(undef, size(d_dims)[1])

    for (i, dimension) in enumerate(d_dims)
        block_var[i] = var(d.distribution[i])
    end
    return cat(block_var..., dims = 1)
end

function var(pd::ParameterDistribution)
    d_dims = get_dimensions(pd)
    block_var = Array{Any}(undef, size(d_dims)[1])

    for (i, dimension) in enumerate(d_dims)
        block_var[i] = var(pd.distribution[i])
    end
    return cat(block_var..., dims = 1) #build the flattened vector
end

"""
    cov(pd::ParameterDistribution)

Returns a dense blocked (co)variance of the distributions.
"""
cov(d::Parameterized) = cov(d.distribution)
cov(d::Samples) = cov(d.distribution_samples, dims = 2) #parameters are columns
function cov(d::VectorOfParameterized)
    d_dims = get_dimensions(d)

    # create each block (co)variance
    block_cov = Array{Any}(undef, size(d_dims)[1])
    for (i, dimension) in enumerate(d_dims)
        if dimension == 1
            block_cov[i] = var(d.distribution[i])
        else
            block_cov[i] = cov(d.distribution[i])
        end
    end

    return cat(block_cov..., dims = (1, 2)) #build the block diagonal (dense) matrix

end

function cov(pd::ParameterDistribution)
    d_dims = get_dimensions(pd)

    # create each block (co)variance
    block_cov = Array{Any}(undef, size(d_dims)[1])
    for (i, dimension) in enumerate(d_dims)
        if dimension == 1
            block_cov[i] = var(pd.distribution[i])
        else
            block_cov[i] = cov(pd.distribution[i])
        end
    end

    return cat(block_cov..., dims = (1, 2)) #build the block diagonal (dense) matrix

end

#extending mean
"""
    mean(pd::ParameterDistribution)

Returns a concatenated mean of the parameter distributions.
"""
mean(d::Parameterized) = mean(d.distribution)
mean(d::Samples) = mean(d.distribution_samples, dims = 2)
function mean(d::VectorOfParameterized)
    return cat([mean(dd) for dd in d.distribution]..., dims = 1)
end
function mean(pd::ParameterDistribution)
    return cat([mean(d) for d in pd.distribution]..., dims = 1)
end

#apply transforms

"""
    transform_constrained_to_unconstrained(pd::ParameterDistribution, xarray::Array{<:Real,1})

Apply the transformation to map (possibly constrained) parameters `xarray` into the unconstrained space.
"""
function transform_constrained_to_unconstrained(
    pd::ParameterDistribution,
    xarray::AbstractVector{FT},
) where {FT <: Real}
    return cat([c.constrained_to_unconstrained(xarray[i]) for (i, c) in enumerate(pd.constraint)]..., dims = 1)
end

"""
    transform_constrained_to_unconstrained(pd::ParameterDistribution, xarray::Array{<:Real,2})

Apply the transformation to map (possibly constrained) parameter samples `xarray` into the unconstrained space.
Here, `xarray` contains parameters as columns and samples as rows.
"""
function transform_constrained_to_unconstrained(
    pd::ParameterDistribution,
    xarray::AbstractMatrix{FT},
) where {FT <: Real}
    return Array(hcat([c.constrained_to_unconstrained.(xarray[i, :]) for (i, c) in enumerate(pd.constraint)]...)')
end

"""
    transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{<:Real,1})

Apply the transformation to map parameters `xarray` from the unconstrained space into (possibly constrained) space.
"""
function transform_unconstrained_to_constrained(
    pd::ParameterDistribution,
    xarray::AbstractVector{FT},
) where {FT <: Real}
    return cat([c.unconstrained_to_constrained(xarray[i]) for (i, c) in enumerate(pd.constraint)]..., dims = 1)
end

"""
    transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{<:Real,2})

Apply the transformation to map parameter samples `xarray` from the unconstrained space into (possibly constrained) space.
Here, `xarray` contains parameters as columns and samples as rows.
"""
function transform_unconstrained_to_constrained(
    pd::ParameterDistribution,
    xarray::AbstractMatrix{FT},
) where {FT <: Real}
    return Array(hcat([c.unconstrained_to_constrained.(xarray[i, :]) for (i, c) in enumerate(pd.constraint)]...)')
end

"""
    transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{Array{<:Real,2},1})

Apply the transformation to map parameter sample ensembles `xarray` from the unconstrained space into (possibly constrained) space.
Here, `xarray` is an iterable of parameters sample ensembles for different EKP iterations.
"""
function transform_unconstrained_to_constrained(
    pd::ParameterDistribution,
    xarray, # ::Iterable{AbstractMatrix{FT}},
) where {FT <: Real}
    transf_xarray = []
    for elem in xarray
        push!(transf_xarray, transform_unconstrained_to_constrained(pd, elem))
    end
    return transf_xarray
end

end # of module


text:
println("hello world")

function quadratic2(a::Float64, b::Float64, c::Float64)
    # unlike other languages 2a is equivalent to 2*a
    # a^2 is used instead of a**2 or pow(a,2)
    sqr_term = sqrt(b^2-4a*c)
    r1 = quadratic(a, sqr_term, b)
    r2 = quadratic(a, -sqr_term, b)
    # multiple values can be returned from a function using tuples
    # if the return keyword is omitted, the last term is returned
    r1, r2
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to conversion from direction cosine matrices to quaternion.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/conversions/dcm_to_quaternion.jl
# ============================================

# Functions: dcm_to_quat
# ----------------------

@testset "DCM => Quaternion (Float64)" begin
    # The conversion between DCM and quaternion is tested by observing if the
    # rotation of a vector is consistent in both representations. In the
    # testset, we add all the special cases in the conversion.
    T = Float64
    testset = [
        (_rand_ang(), _rand_ang2(), _rand_ang(), :Z, :Y, :X)
        (-0.3,        +0.5,         +π,          :Z, :Y, :X)
        (+0.3,        +0.5,         +π,          :Z, :Y, :X)
        (+0.5,        +π,           -0.3,        :Z, :Y, :X)
        (+0.5,        +π,           +0.3,        :Z, :Y, :X)
        (+π,          +0.5,         -0.3,        :Z, :Y, :X)
        (+π,          +0.5,         +0.3,        :Z, :Y, :X)
    ]

    for test in testset
        # Unpack values in tuple.
        a₁, a₂, a₃, r₁, r₂, r₃ = test

        # Create the DCM.
        D = angle_to_dcm(a₃, r₃) * angle_to_dcm(a₂, r₂) * angle_to_dcm(a₁, r₁)

        # Convert it to a quaternion.
        q = dcm_to_quat(D)

        # Test quaternion type.
        @test eltype(q) === T

        # The real part must always be positive.
        @test q.q0 > 0

        # The quaternion must be unitary.
        @test √(q.q0^2 + q.q1^2 + q.q2^2 + q.q3^2) ≈ 1

        # Compare the rotations between the DCM and quaternion.
        v = @SVector randn(3)

        vrd = D * v
        vrq = vect(q \ v * q)

        @test vrd ≈ vrq
        @test eltype(vrd) === eltype(vrq) === T
    end
end

@testset "DCM => Quaternion (Float32)" begin
    # The conversion between DCM and quaternion is tested by observing if the
    # rotation of a vector is consistent in both representations. In the
    # testset, we add all the special cases in the conversion.
    T = Float32
    testset = [
        (_rand_ang(T), _rand_ang2(T), _rand_ang(T), :Z, :Y, :X)
        (T(-0.3),      T(+0.5),       T(+π),        :Z, :Y, :X)
        (T(+0.3),      T(+0.5),       T(+π),        :Z, :Y, :X)
        (T(+0.5),      T(+π),         T(-0.3),      :Z, :Y, :X)
        (T(+0.5),      T(+π),         T(+0.3),      :Z, :Y, :X)
        (T(+π),        T(+0.5),       T(-0.3),      :Z, :Y, :X)
        (T(+π),        T(+0.5),       T(+0.3),      :Z, :Y, :X)
    ]

    for test in testset
        # Unpack values in tuple.
        a₁, a₂, a₃, r₁, r₂, r₃ = test

        # Create the DCM.
        D = angle_to_dcm(a₃, r₃) * angle_to_dcm(a₂, r₂) * angle_to_dcm(a₁, r₁)

        # Convert it to a quaternion.
        q = dcm_to_quat(D)

        # Test quaternion type.
        @test eltype(q) === T

        # The real part must always be positive.
        @test q.q0 > 0

        # The quaternion must be unitary.
        @test √(q.q0^2 + q.q1^2 + q.q2^2 + q.q3^2) ≈ 1

        # Compare the rotations between the DCM and quaternion.
        v = @SVector randn(T, 3)

        vrd = D * v
        vrq = vect(q \ v * q)

        @test vrd ≈ vrq
        @test eltype(vrd) === eltype(vrq) === T
    end
end


text:
import Base: show

using Parquet
using DataFrames

mutable struct VersionedFile
    # Versioned Model Data
    versionstoredf::DataFrame

    # Versioning Detail Fields
    # Storage Information
    version_store_file::AbstractString
end

function VersionedFile(version_store_file::AbstractString)
    if isfile(version_store_file)
        versionstoredf = DataFrame(read_parquet(version_store_file))
    else
        versionstoredf = DataFrame(Name = Any[], Model=Any[], Version=Int[], Hash=Any[])
    end
    VersionedFile(versionstoredf, version_store_file)
end

function saveFile(versionfile::VersionedFile)
    # Function to write the data to the file
    write_parquet(versionfile.version_store_file, versionfile.versionstoredf)
end

function addModelVersion!(versionfile::VersionedFile, versionedmodel::VersionedModel)
    modelfromfile = filter(:name => modelnames -> modelnames == name, versionfile.versionstoredf)
    if length(modelfromfile) == 0
        df = DataFrame(Name = [versionedmodel.name], Model=[versionedmodel.stan_model], Version=[versionedmodel.version], Hash=[versionedmodel.hash])
        append!(versionfile.versionstoredf, df)
    else
        duplicatemodelfromfile = filter(:hash => modelhash -> modelhash == hash, modelfromfile)
        if length(duplicatemodelfromfile) > 0
            versionedmodel.version = duplicatemodelfromfile[:hash, 1]
        else
            versionedmodel.version = max(duplicatemodelfromfile[:Version]) + 1
            df = DataFrame(Name = [name], Model=[model], Version=[version_number], Hash=[hash])
            append!(versionstoredf, df)
        end
    end
end


text:
#export top

#### preciso editar essa funcação daqui

# one IL + Protein + water
function top(dict,top_dir, prot_dir, data, nions, nwater)
  anion  = data.anion
  cation = data.cation
  file = open("$prot_dir/topol.top","r")
  topol  = open("./topol.top","w") # creation of the topol.top which will be used in the simulations
  n = 0
  for line in eachline(file)
     n = n + 1
     if n==19
       println(topol,"#define _ff_opls")
       println(topol,"#define _ff_oplsaa")
       println(topol,"                 ")
     elseif n==22
       println(topol,"      ")
       println(topol,"; load atom types and molecules")
       println(topol,"#include \"$top_dir/$(data.an_atomtypes)\" ")
       println(topol,"#include \"$top_dir/$(data.cat_atomtypes)\" ")
       println(topol,"#include \"$top_dir/$(data.an_par)\" ")
       println(topol,"#include \"$top_dir/$(data.cat_par)\" ")
       println(topol,"             ")
     elseif n>19 && n != 22 && occursin("SOL",line)==false
       println(topol, line)
     elseif occursin("SOL",line)==true
       println(topol,line)
       println(topol,"SOL                $nwater")
       println(topol,"$(dict[cation])                $nions")
       println(topol,"$(dict[anion])                $nions")
     end
  end
  close(file)
  close(topol)
end

# Two ILs +  water
function topelec(dict,top_dir, ionic_dir, data, ncat, nan, nwater)
  anion1  = data.anion1
  anion2  = data.anion2
  cation  = data.cation
  file = open("$ionic_dir/topol.top","r")
  topol  = open("./topol.top","w") # creation of the topol.top which will be used in the simulations
  n = 0
  for line in eachline(file)
     n = n + 1
     if n==7
       println(topol,"      ")
       println(topol,"; load atom types and molecules")
       println(topol,"#include \"$top_dir/$(data.an1_atomtypes)\" ")
       println(topol,"#include \"$top_dir/$(data.an2_atomtypes)\" ")
       println(topol,"#include \"$top_dir/$(data.cat_atomtypes)\" ")
       println(topol,"#include \"$top_dir/$(data.an1_par)\" ")
       println(topol,"#include \"$top_dir/$(data.an2_par)\" ")
       println(topol,"#include \"$top_dir/$(data.cat_par)\" ")
       println(topol,"             ")
     elseif  n != 7 && occursin("Compound",line)==false
       println(topol, line)
     elseif occursin("Compound",line)==true
       println(topol,line)
       println(topol,"SOL                $nwater")
       println(topol,"$(dict[cation])                $ncat")
       println(topol,"$(dict[anion1])                $nan")
       println(topol,"$(dict[anion2])                $nan")
     end
  end
  close(file)
  close(topol)
end

# one IL +  water
function topelec(dict,top_dir, ionic_dir, data, nions, nwater)
  anion   = data.anion
  cation  = data.cation
  file = open("$ionic_dir/topol.top","r")
  topol  = open("./topol.top","w") # creation of the topol.top which will be used in the simulations
  n = 0
  for line in eachline(file)
     n = n + 1
     if n==7
       println(topol,"      ")
       println(topol,"; load atom types and molecules")
       println(topol,"#include \"$top_dir/$(data.an_atomtypes)\" ")
       println(topol,"#include \"$top_dir/$(data.cat_atomtypes)\" ")
       println(topol,"#include \"$top_dir/$(data.an_par)\" ")
       println(topol,"#include \"$top_dir/$(data.cat_par)\" ")
       println(topol,"             ")
     elseif  n != 7 && occursin("Compound",line)==false
       println(topol, line)
     elseif occursin("Compound",line)==true
       println(topol,line)
       println(topol,"SOL                $nwater")
       println(topol,"$(dict[cation])                $nions")
       println(topol,"$(dict[anion])                $nions")
     end
  end
  close(file)
  close(topol)
end

# Two IL + Protein + water
function top(dict,top_dir, prot_dir, data, ncat, nan, nwater)
  anion1 = data.anion1
  anion2 = data.anion2
  cation = data.cation
  file = open("$prot_dir/topol.top","r")
  topol  = open("./topol.top","w") # creation of the topol.top which will be used in the simulations
  n = 0
  for line in eachline(file)
     n = n + 1
     if n==19
       println(topol,"#define _ff_opls")
       println(topol,"#define _ff_oplsaa")
       println(topol,"                 ")
     elseif n==22
       println(topol,"      ")
       println(topol,"; load atom types and molecules")
       println(topol,"#include \"$top_dir/$(data.an1_atomtypes)\" ")
       println(topol,"#include \"$top_dir/$(data.an2_atomtypes)\" ")
       println(topol,"#include \"$top_dir/$(data.cat_atomtypes)\" ")
       println(topol,"#include \"$top_dir/$(data.an1_par)\" ")
       println(topol,"#include \"$top_dir/$(data.an2_par)\" ")
       println(topol,"#include \"$top_dir/$(data.cat_par)\" ")
       println(topol,"             ")
     elseif n>19 && n != 22 && occursin("SOL",line)==false
       println(topol, line)
     elseif occursin("SOL",line)==true
       println(topol,line)
       println(topol,"SOL                $nwater")
       println(topol,"$(dict[cation])                $ncat")
       println(topol,"$(dict[anion1])                $nan")
       println(topol,"$(dict[anion2])                $nan")
     end
  end
  close(file)
  close(topol)
end


text:
using Printf

include("gen_rmesh_exp.jl")
include("lagrange_interp.jl")
include("radial_interp.jl")
include("rsch_integ_rk4.jl")
include("sch_rk4_step.jl")
include("get_min_idx.jl")
include("get_n_nodes.jl")

function main()
    # Build radial function
    r_min = 1e-8
    r_max = 50.0
    a = 1e6
    Nr = 3001
    rmesh = gen_rmesh_exp(r_min, r_max, a, Nr)

    Z = 3
    n = 1
    l = 0

    P = zeros(Float64,Nr)
    Q = zeros(Float64,Nr)
    V = -Z ./ rmesh

    Emin = -5000.0
    Emax = 0.0
    #E = 0.5*(Emin + Emax)
    E = -1000.0

    for iterShoot in 1:30

        println("\niterShoot = ", iterShoot)
        println("E = ", E)
        dE = abs(Emax - Emin)
        if dE < 1e-10
            println("Converged by dE")
            println("dE = ", dE)
            break
        end

        ctp = Nr # no perturbation correction

        imax = rsch_integ_rk4!( E, Z, l, rmesh, V, P, Q )
        println("imax = ", imax)

        minidx = get_min_idx(imax, P)
        println("minidx = ", minidx)

        #@views P[minidx:Nr] .= 0.0
        #@views Q[minidx:Nr] .= 0.0

        Nnodes = get_n_nodes( imax, P )
        println("Nnodes = ", Nnodes)
        println("minidx = ", minidx)

        cond1 = Nnodes != n - l - 1
        cond2 = ctp == Nr
        cond3 = imax < ctp

        if cond1 || cond2 || cond3
            is_big = Nnodes > (n - l - 1)
            println("n = ", n, " l = ", l)
            println("is_big = ", is_big)
            @printf("before: %18.10f %18.10f\n", Emin, Emax)
            if is_big
                Emax = E
            else
                Emin = E
            end
            @printf("After: %18.10f %18.10f\n", Emin, Emax)
            #
            E = 0.5*(Emin + Emax)
            println("E = ", E)
        end

    end

end

@time main()


text:
# Univariate Options
function optimize(f,
     lower::T,
     upper::T;
     method = Brent(),
     rel_tol::Real = sqrt(eps(float(T))),
     abs_tol::Real = eps(float(T)),
     iterations::Integer = 1_000,
     store_trace::Bool = false,
     show_trace::Bool = false,
     callback = nothing,
     show_every = 1,
     extended_trace::Bool = false) where T <: Real
    show_every = show_every > 0 ? show_every : 1
    if extended_trace && callback === nothing
        show_trace = true
    end

    show_trace && print_header(method)
    Tf = float(T)
    optimize(f, Tf(lower), Tf(upper), method;
             rel_tol = Tf(rel_tol),
             abs_tol = Tf(abs_tol),
             iterations = iterations,
             store_trace = store_trace,
             show_trace = show_trace,
             show_every = show_every,
             callback = callback,
             extended_trace = extended_trace)
end

function optimize(f,
    lower::Union{Integer, Real},
    upper::Union{Integer, Real};
    kwargs...)

    T = promote_type(typeof(lower/1), typeof(upper/1))
    optimize(f,
             T(lower),
             T(upper);
             kwargs...)
end

function optimize(f,
    lower::Union{Integer, Real},
    upper::Union{Integer, Real},
    method::Union{Brent, GoldenSection};
    kwargs...)

    T = promote_type(typeof(lower/1), typeof(upper/1))
    optimize(f,
             T(lower),
             T(upper),
             method;
             kwargs...)
end


text:
using SpecialFunctions: gamma

"
Evaluate Jacobi polynomial of type
(alpha, beta) > -1, (alpha+beta != -1)
at points x for order N. Returns array of values, same length as x

See: https://en.wikipedia.org/wiki/Jacobi_polynomials

NOTE: This calculates the normalised Jacobi polynomial

"
function JacobiP(x, alpha, beta, N::Int)

    nvalues = length(x)

    PL = zeros(Float64, nvalues, N+1)

    gamma0 = 2^(alpha+beta+1)/(alpha+beta+1) * gamma(alpha+1) * gamma(beta+1)/gamma(alpha+beta+1)
    PL[:,1] .= 1.0 / sqrt(gamma0)

    if N == 0
        return PL[:,1]
    end

    gamma1 = (alpha+1)*(beta+1)/(alpha+beta+3)*gamma0
    PL[:,2] = ((alpha+beta+2)*x/2 .+ (alpha-beta)/2)/sqrt(gamma1)

    if N == 1
        return PL[:,2]
    end

    aold = 2/(2+alpha+beta)*sqrt((alpha+1)*(beta+1)/(alpha+beta+3))

    # Forward recurrence
    for i = 1:(N-1)
        h1 = 2*i + alpha + beta
        anew = 2/(h1+2)*sqrt( (i+1)*(i+1+alpha+beta)*(i+1+alpha)*(i+1+beta)/(h1+1)/(h1+3))
        bnew = - (alpha^2 - beta^2)/h1/(h1+2)
        PL[:,i+2] = (1/anew)*(-aold * PL[:,i] + (x .- bnew) .* PL[:,i+1])
        aold = anew
    end
    return PL[:,end]
end

"Evaluate the derivative of the Jacobi polynomial
of type (alpha,beta) > -1 at points r for order N

## Returns
array of the same length as r
"
function GradJacobiP(r, alpha, beta, N::Int)
    if N == 0
        zeros(length(r))
    else
        sqrt(N * (N+alpha+beta+1))*JacobiP(r, alpha+1, beta+1, N-1)
    end
end

struct JacobiBasis1D
    order::Int  # Number of polynomials
end

function basisFunction(basis::JacobiBasis1D, locations, index::Int)
    JacobiP(locations, 0, 0, index-1)
end

function basisGradient(basis::JacobiBasis1D, locations, index::Int)
    GradJacobiP(locations,0,0,index-1)
end


text:
using Luxor, Colors

function main(fname)
    Drawing(1280, 640, fname)
    origin()

    bbx= BoundingBox()

    # background
    setblend(blend(
    boxtopcenter(bbx),
    boxbottomcenter(bbx),
    RGB(0.18, 0, .24),
    RGB(0.02, 0, 0.14)))
    paint()

    setblend(blend(
    boxtopcenter(bbx),
    boxbottomcenter(bbx),
    RGB(0.02, 0, 0.14),
    RGB(0.18, 0, .24)))
    poly([boxtopleft(bbx),
    between(boxtopleft(bbx), boxtopright(bbx), 360/1280),
    between(boxbottomleft(bbx), boxbottomright(bbx), 590/1280),
    boxbottomleft(bbx)], :fill)

    # show croppable zones
    #@layer begin
    # sethue("red")
    # setline(0.5)
    # rule(boxtopcenter(bbx) + (0, 80))
    # rule(boxbottomcenter(bbx) + (0, -80))
    # rule(boxmiddleleft(bbx) + (80, 0), π/2)
    # rule(boxmiddleright(bbx) + (-80, 0), π/2)
    #end

    # logo
    try
        logo = readsvg("symbolics/symbolics-logo-black-white.svg")
        @layer begin
            translate(boxmiddleleft(bbx) + (250, 0))
            scale(0.4)
            placeimage(logo, centered=true)
        end
    catch e
        @error "Can't load logo"
    end

    # text
    sethue("white")

    # top
    fontsize(90)
    fontface("WorkSans-ExtraBold")
    text("JuliaSymbolics", O + (-130, -170), halign=:left)

    fontsize(29)
    fontface("WorkSans-ExtraBoldItalic")
    text("A fast and modern Computer Algebra System ...", O + (-130, -120), halign=:left)

    # middle
    fontsize(19)
    fontface("WorkSans-Regular")
    textwrap("Symbolics.jl is a pure Julia CAS which uses the Julia core library to its fullest. It is built from the ground up with performance in mind. We use specialized structures for automatic simplification to match the performance of the most fully optimized C++ libraries. It exploits parallelism at every level; our symbolic simplification takes advantage of Julia's task-based multithreading to transform symbolic equations into parallelized Julia code.", 650, O + (-60, -80), leading=30, rightgutter=30)

    # bottom
    fontsize(29)
    fontface("WorkSans-ExtraBoldItalic")
    text("... for a fast and modern language", boxbottomright(bbx) + (-80, -90), halign=:right)

    finish()
    preview()
end

main("/tmp/social-media-preview.png")

max_stars_repo_name: JuliaSymbolics/Resources

max_issues_repo_name: JuliaSymbolics/Resources

max_forks_repo_name: JuliaSymbolics/Resources


text:
include(joinpath("deps", "ReadMolecule.jl"))
include(joinpath("deps", "Shapes.jl"))
import .ReadMoelcule
import .Shapes
import SimpleMolecule: Molecule, obtain
import MolecularPDB

function main(argd)
    systems = ReadMolecule.read(convert(Array{String,1}, argd["inputs"]))

    fn_assemble = Shapes.choose(Symbol(argd["shape"]))

    @time molecules = fn_assemble(systems, argd["params"])
    new_system = Molecule(vcat([obtain(mol, :atoms) for mol in molecules]...))

    for file_output in argd["outputs"]
        @time MolecularPDB.save(
            MolecularPDB.VMD,
            file_output,
            new_system
        )
    end
end

max_stars_repo_name: yuhangwang/SimShape

max_issues_repo_name: yuhangwang/SimShape

max_forks_repo_name: yuhangwang/SimShape


text:
function calc_W_from_B(Tk, B, P)

    pws = Pws(B)
    efac = efactor(B,P)
    xsv = efac * pws / P
    w2 = Mv/Ma * xsv/(1-xsv)

    w = (enthalpyair(B,P) - enthalpyair(Tk,P) - w2*(enthalpywi(B) - enthalpyvapor(B))) /
        (enthalpyvapor(Tk) - enthalpywi(B))

    f0(w0) = aux_WB(w0, Tk, B, P, pws, efac)
    return Roots.find_zero(f0,w)
end

function aux_WB(w, Tk, B, P, pws, efac)
    xv1 = w / (Mv/Ma+w)
    xv2 = efac * pws / P
    w2 = Mv / Ma * xv2 / (1.0-xv2)
    #(1.0+w)*enthalpymoist(Tk,P,xv1) + (w2-w)*enthalpywi(B) - (1.0+w2)*enthalpymoist(B,P,xv2)
    enthalpymoist(Tk,P,xv1) + (w2-w)*enthalpywi(B) - enthalpymoist(B,P,xv2)
end

function aux_WB(w, Tk, B, P)
    pws = Pws(B)
    efac = efactor(B,P)
    xv1 = w / (Mv/Ma+w)
    xv2 = efac * pws / P
    w2 = Mv / Ma * xv2 / (1.0-xv2)
    #(1.0+w)*enthalpymoist(Tk,P,xv1) + (w2-w)*enthalpywi(B) - (1.0+w2)*enthalpymoist(B,P,xv2)
    enthalpymoist(Tk,P,xv1) + (w2-w)*enthalpywi(B) - enthalpymoist(B,P,xv2)
end

function calcwetbulb(Tk, P, xv)
    w = humrat(xv)
    B = Tk - 1.0 # Initial guess
    h = 1e-7
    dB = 0.0
    f0(_B) = aux_WB(w, Tk,_B, P)
    return Roots.find_zero(f0,B)
end

humrat(xv) = xv / (1-xv) * (Mv/Ma)


text:
import Base.pop!, Base.push!, Base.length, Base.*
using LinearAlgebra, Random

struct ProjectorStack{P <: AbstractVector}
    vs::Vector{P}
end
Base.push!(p::ProjectorStack, v) = (push!(p.vs, v); p)
Base.pop!(p::ProjectorStack) = (pop!(p.vs); p)
Base.:*(p::ProjectorStack, v) = (s = zero(v); for vi in p.vs s += vi * dot(vi, v) end; s)

"""
    GärtnerBoundary

See the passage regarding M_B in Section 4 of Gärtner's paper.
"""
mutable struct GärtnerBoundary{P<:AbstractVector, F<:AbstractFloat}
    centers::Vector{P}
    square_radii::Vector{F}
    # projection onto of affine space spanned by points
    # shifted such that first point becomes origin
    projector::ProjectorStack{P}
    empty_center::P # center of nsphere spanned by empty boundary
end

function GärtnerBoundary(pts)
    P = eltype(pts)
    F = eltype(P)
    projector, centers, square_radii = ProjectorStack(P[]), P[], F[]
    empty_center = F(NaN) * first(pts)
    GärtnerBoundary(centers, square_radii, projector, empty_center)
end

function push_if_stable!(b::GärtnerBoundary, pt)
    if isempty(b)
        push!(b.square_radii, zero(eltype(pt)))
        push!(b.centers, pt)
        dim = length(pt)
        return true
    end
    q0, center = first(b.centers), b.centers[end]
    C, r2  = center - q0, b.square_radii[end]
    Qm, M = pt - q0, b.projector
    Qm_bar = M*Qm
    residue, e = Qm - Qm_bar, sqdist(Qm, C) - r2
    z, tol = 2*sqnorm(residue), eps(eltype(pt)) * max(r2, one(r2))
    isstable = abs(z) > tol
    if isstable
        center_new  = center + (e/z) * residue
        r2new = r2 + (e^2)/(2z)
        push!(b.projector, residue / norm(residue))
        push!(b.centers, center_new)
        push!(b.square_radii, r2new)
    end
    return isstable
end

function Base.pop!(b::GärtnerBoundary)
    n = length(b)
    pop!(b.centers)
    pop!(b.square_radii)
    if n >= 2
        pop!(b.projector)
    end
    return b
end

struct NSphere{P,F}
    center::P
    sqradius::F
end

function isinside(pt, nsphere::NSphere; atol=0, rtol=0)
    r2, R2 = sqdist(pt, center(nsphere)), sqradius(nsphere)
    return r2 <= R2 || isapprox(r2, R2;atol=atol^2,rtol=rtol^2)
end
allinside(pts, nsphere; kw...) = all(pt -> isinside(pt, nsphere; kw...), pts)

function move_to_front!(pts, i)
    pt = pts[i]
    for j in eachindex(pts)
        pts[j], pt = pt, pts[j]
        j == i && break
    end
    pts
end

prefix(pts, i) = view(pts, 1:i)
Base.length(b::GärtnerBoundary) = length(b.centers)
Base.isempty(b::GärtnerBoundary) = length(b) == 0
center(b::NSphere) = b.center
radius(b::NSphere) = sqrt(b.sqradius)
sqradius(b::NSphere) = b.sqradius
dist(p1,p2) = norm(p1-p2)
sqdist(p1::AbstractVector, p2::AbstractVector) = sqnorm(p1-p2)
sqnorm(p) = sum(abs2,p)
ismaxlength(b::GärtnerBoundary) = length(b) == length(b.empty_center) + 1

function NSphere(b::GärtnerBoundary)
    return isempty(b) ? NSphere(b.empty_center, 0.0) :
        NSphere(b.centers[end], b.square_radii[end])
end

function welzl!(pts, bdry)
    bdry_len, support_count, nsphere = length(bdry), 0, NSphere(bdry)
    ismaxlength(bdry) && return nsphere, 0
    for i in eachindex(pts)
        pt = pts[i]
        if !isinside(pt, nsphere)
            pts_i = prefix(pts, i-1)
            isstable = push_if_stable!(bdry, pt)
            if isstable
                nsphere, s = welzl!(pts_i, bdry)
                pop!(bdry)
                move_to_front!(pts, i)
                support_count = s + 1
            end
        end
    end
    return nsphere, support_count
end

function find_max_excess(nsphere, pts, k1)
    T = eltype(first(pts))
    e_max, k_max = T(-Inf), k1 -1
    for k in k1:length(pts)
        pt = pts[k]
        e = sqdist(pt, center(nsphere)) - sqradius(nsphere)
        if  e > e_max
            e_max = e
            k_max = k
        end
    end
    return e_max, k_max
end

function welzl(points, maxiterations=2000)
    pts = deepcopy(points)
    bdry, t = GärtnerBoundary(pts), 1
    nsphere, s = welzl!(prefix(pts, t), bdry)
    for i in 1:maxiterations
        e, k = find_max_excess(nsphere, pts, t + 1)
        P = eltype(pts)
        F = eltype(P)
        e <= eps(F) && break
        pt = pts[k]
        push_if_stable!(bdry, pt)
        nsphere_new, s_new = welzl!(prefix(pts, s), bdry)
        pop!(bdry)
        move_to_front!(pts, k)
        nsphere = nsphere_new
        t, s = s + 1, s_new + 1
    end
    return nsphere
end

function testwelzl()
    testdata =[
        [[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]],
        [[5.0, -2.0], [-3.0, -2.0], [-2.0, 5.0], [1.0, 6.0], [0.0, 2.0]],
        [[2.0, 4.0, -1.0], [1.0, 5.0, -3.0], [8.0, -4.0, 1.0], [3.0, 9.0, -5.0]],
        [randn(5) for _ in 1:8]
    ]
    for test in testdata
        nsphere = welzl(test)
        println("For points: ", test)
        println("    Center is at: ", nsphere.center)
        println("    Radius is: ", radius(nsphere), "\n")
    end
end

testwelzl()

max_stars_repo_name: ashwani-rathee/Sigil

max_issues_repo_name: ashwani-rathee/Sigil

max_forks_repo_name: ashwani-rathee/Sigil


text:
@testset "ch07_kld_mcl" begin
    dt = 0.1
    # environment
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    # id of landmark must start from 0 with 1 step
    landmarks = [Landmark([2.0, -3.0], 0), Landmark([3.0, 3.0], 1)]
    envmap = Map()
    push!(envmap, landmarks)
    world = World(xlim, ylim)
    push!(world, envmap)
    # robot side
    initial_pose = [0.0, 0.0, 0.0]
    estimator = KldMcl(initial_pose, 1000)
    circling_agent = EstimatorAgent(0.2, 10.0 * pi / 180, dt, estimator)
    robot = RealRobot(initial_pose, circling_agent, RealCamera(landmarks); color = "red")
    push!(world, robot)
    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        z = observations(robot.sensor_, robot.pose_; noise = true, bias = true)
        p = draw(world, annota)
        v, ω = decision(circling_agent, z, envmap; resample = true)
        state_transition(robot, v, ω, dt; move_noise = true, vel_bias_noise = true)
    end
    if GUI
        gif(anim, "ch07_kld_mcl.gif", fps = 20)
    end
end

@testset "ch07_mcl_global" begin
    dt = 0.1
    # environment
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    # id of landmark must start from 0 with 1 step
    landmarks =
        [Landmark([2.0, -3.0], 0), Landmark([3.0, 3.0], 1), Landmark([-4.0, 2.0], 2)]
    envmap = Map()
    push!(envmap, landmarks)
    world = World(xlim, ylim)
    push!(world, envmap)
    # robot side
    initial_pose = [0.0, 0.0, 0.0]
    estimator = Mcl(initial_pose, 100; glob = true, xlim = xlim, ylim = ylim)
    circling_agent = EstimatorAgent(0.2, 10.0 * pi / 180, dt, estimator)
    robot = RealRobot(initial_pose, circling_agent, RealCamera(landmarks); color = "red")
    push!(world, robot)
    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        z = observations(robot.sensor_, robot.pose_; noise = true, bias = true)
        p = draw(world, annota)
        v, ω = decision(circling_agent, z, envmap; resample = true)
        state_transition(robot, v, ω, dt; move_noise = true, vel_bias_noise = true)
    end
    if GUI
        gif(anim, "ch07_mcl_global.gif", fps = 20)
    end
end

@testset "ch07_kf_global" begin
    dt = 0.1
    # environment
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    # id of landmark must start from 0 with 1 step
    landmarks =
        [Landmark([-4.0, 2.0], 0), Landmark([2.0, -3.0], 1), Landmark([3.0, 3.0], 2)]
    envmap = Map()
    push!(envmap, landmarks)
    world = World(xlim, ylim)
    push!(world, envmap)
    # robot side
    initial_pose = [0.0, 0.0, 0.0]
    estimator = KalmanFilter(envmap, initial_pose; glob = true, xlim = xlim, ylim = ylim)
    agent = EstimatorAgent(0.2, 10.0 * pi / 180, dt, estimator)
    robot = RealRobot(initial_pose, agent, RealCamera(landmarks); color = "red")
    push!(world, robot)

    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        z = observations(robot.sensor_, robot.pose_; noise = true, bias = true)
        p = draw(world, annota)
        v, ω = decision(agent, z, envmap)
        state_transition(robot, v, ω, dt; move_noise = true, vel_bias_noise = true)
    end
    if GUI
        gif(anim, "ch07_kf_global.gif", fps = 20)
    end
end

@testset "ch07_kf_kidnap" begin
    dt = 0.1
    # environment
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    # id of landmark must start from 0 with 1 step
    landmarks =
        [Landmark([-4.0, 2.0], 0), Landmark([2.0, -3.0], 1), Landmark([3.0, 3.0], 2)]
    envmap = Map()
    push!(envmap, landmarks)
    world = World(xlim, ylim)
    push!(world, envmap)
    # robot side
    initial_pose = [0.0, 0.0, 0.0]
    estimator = KalmanFilter(envmap, initial_pose)
    agent = EstimatorAgent(0.2, 10.0 * pi / 180, dt, estimator)
    robot = RealRobot(
        initial_pose,
        agent,
        RealCamera(landmarks);
        color = "red",
        expected_kidnap_time = 10,
    )
    push!(world, robot)

    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        z = observations(robot.sensor_, robot.pose_; noise = true, bias = true)
        p = draw(world, annota)
        v, ω = decision(agent, z, envmap)
        state_transition(
            robot,
            v,
            ω,
            dt;
            move_noise = true,
            vel_bias_noise = true,
            kidnap = true,
        )
    end
    if GUI
        gif(anim, "ch07_kf_global.gif", fps = 20)
    end
end

@testset "ch07_mcl_kidnap" begin
    dt = 0.1
    # environment
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    # id of landmark must start from 0 with 1 step
    landmarks =
        [Landmark([2.0, -3.0], 0), Landmark([3.0, 3.0], 1), Landmark([-4.0, 2.0], 2)]
    envmap = Map()
    push!(envmap, landmarks)
    world = World(xlim, ylim)
    push!(world, envmap)
    # robot side
    initial_pose = [0.0, 0.0, 0.0]
    estimator = Mcl(initial_pose, 100)
    circling_agent = EstimatorAgent(0.2, 10.0 * pi / 180, dt, estimator)
    robot = RealRobot(
        initial_pose,
        circling_agent,
        RealCamera(landmarks);
        color = "red",
        expected_kidnap_time = 20,
    )
    push!(world, robot)
    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        z = observations(robot.sensor_, robot.pose_; noise = true, bias = true)
        p = draw(world, annota)
        v, ω = decision(circling_agent, z, envmap; resample = true)
        state_transition(
            robot,
            v,
            ω,
            dt;
            move_noise = true,
            vel_bias_noise = true,
            kidnap = true,
        )
    end
    if GUI
        gif(anim, "ch07_mcl_kidnap.gif", fps = 20)
    end
end

@testset "ch07_reset_mcl" begin
    dt = 0.1
    # environment
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    # id of landmark must start from 0 with 1 step
    landmarks =
        [Landmark([2.0, -3.0], 0), Landmark([3.0, 3.0], 1), Landmark([-4.0, 2.0], 2)]
    envmap = Map()
    push!(envmap, landmarks)
    world = World(xlim, ylim)
    push!(world, envmap)
    # robot side
    initial_pose = [0.0, 0.0, 0.0]
    estimator = ResetMcl(initial_pose, 100)
    circling_agent = EstimatorAgent(0.2, 10.0 * pi / 180, dt, estimator)
    robot = RealRobot(initial_pose, circling_agent, RealCamera(landmarks); color = "red")

    push!(world, robot)
    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        z = observations(robot.sensor_, robot.pose_; noise = true, bias = true)
        p = draw(world, annota)
        v, ω = decision(circling_agent, z, envmap; resample = true, sensor_reset = true)
        state_transition(robot, v, ω, dt; move_noise = true, vel_bias_noise = true)
    end
    if GUI
        gif(anim, "ch07_reset_mcl.gif", fps = 20)
    end
end

@testset "ch07_sensor_reset" begin
    dt = 0.1
    # environment
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    # id of landmark must start from 0 with 1 step
    landmarks = [
        Landmark([2.0, -3.0], 0),
        Landmark([3.0, 3.0], 1),
        Landmark([-4.0, 2.0], 2),
        Landmark([0.0, 0.0], 3),
        Landmark([-4.0, -4.0], 4),
    ]
    envmap = Map()
    push!(envmap, landmarks)
    world = World(xlim, ylim)
    push!(world, envmap)
    # robot side
    initial_pose = uniform(PoseUniform(xlim, ylim))
    estimator = ResetMcl(initial_pose, 100; xlim = xlim, ylim = ylim, α_threshold = 0.005)
    circling_agent = EstimatorAgent(0.2, 10.0 * pi / 180, dt, estimator)
    robot = RealRobot(
        initial_pose,
        circling_agent,
        RealCamera(landmarks);
        color = "red",
        expected_kidnap_time = 10,
    )
    push!(world, robot)
    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        z = observations(robot.sensor_, robot.pose_; noise = true, bias = true)
        p = draw(world, annota)
        v, ω = decision(circling_agent, z, envmap; resample = true, sensor_reset = true)
        state_transition(
            robot,
            v,
            ω,
            dt;
            move_noise = true,
            vel_bias_noise = true,
            kidnap = true,
        )
    end
    if GUI
        gif(anim, "ch07_sensor_reset2.gif", fps = 20)
    end
end

@testset "ch07_reset_mcl2" begin
    dt = 0.1
    # environment
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    # id of landmark must start from 0 with 1 step
    landmarks =
        [Landmark([2.0, -3.0], 0), Landmark([3.0, 3.0], 1), Landmark([-4.0, 2.0], 2)]
    envmap = Map()
    push!(envmap, landmarks)
    world = World(xlim, ylim)
    push!(world, envmap)
    # robot side
    initial_pose = [0.0, 0.0, 0.0]
    estimator = ResetMcl(initial_pose, 100; xlim = xlim, ylim = ylim, α_threshold = 0.005)
    circling_agent = EstimatorAgent(0.2, 10.0 * pi / 180, dt, estimator)
    robot = RealRobot(
        initial_pose,
        circling_agent,
        RealCamera(landmarks);
        color = "red",
        expected_kidnap_time = 30,
    )
    push!(world, robot)
    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        z = observations(robot.sensor_, robot.pose_; noise = true, bias = true)
        p = draw(world, annota)
        v, ω = decision(circling_agent, z, envmap; resample = true, sensor_reset = true)
        state_transition(
            robot,
            v,
            ω,
            dt;
            move_noise = true,
            vel_bias_noise = true,
            kidnap = true,
        )
    end
    if GUI
        gif(anim, "ch07_reset_mcl2.gif", fps = 20)
    end
end

@testset "ch07_adaptive_mcl" begin
    dt = 0.1
    # environment
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    # id of landmark must start from 0 with 1 step
    landmarks = [Landmark([2.0, -3.0], 0), Landmark([3.0, 3.0], 1)]
    envmap = Map()
    push!(envmap, landmarks)
    world = World(xlim, ylim)
    push!(world, envmap)
    # robot side
    initial_pose = [0.0, 0.0, 0.0]
    estimator = AMcl(initial_pose, 100; xlim = xlim, ylim = ylim)
    circling_agent = EstimatorAgent(0.2, 10.0 * pi / 180, dt, estimator)
    robot = RealRobot(
        initial_pose,
        circling_agent,
        RealCamera(
            landmarks,
            phantom_prob = 0.1,
            phantom_range_x = xlim,
            phantom_range_y = ylim,
        );
        color = "red",
    )
    push!(world, robot)
    anim = @animate for i = 1:300
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        z = observations(
            robot.sensor_,
            robot.pose_;
            noise = true,
            bias = true,
            phantom = true,
        )
        p = draw(world, annota)
        v, ω = decision(circling_agent, z, envmap)
        state_transition(
            robot,
            v,
            ω,
            dt;
            move_noise = true,
            vel_bias_noise = true,
            kidnap = false,
        )
    end
    if GUI
        gif(anim, "ch07_adaptive_mcl.gif", fps = 20)
    end
end

max_stars_repo_name: soblin/JuliaProbo

max_issues_repo_name: soblin/JuliaProbo

max_forks_repo_name: soblin/JuliaProbo


text:
# Author: AJ Tyler
# Date: 24/01/22
module IO # This module contains functions that read in and write out files.
include("CNA.jl") # This file contains the CommonNeighbourAnalysis module.
include("atomicShells.jl") # This file contains the AtomicShells module.
using .CommonNeighbourAnalysis # This module computes the total and atomic CNA profile of a configuration.
using .AtomicShells # This module computes the atomic shell number of a cluster.
using Graphs # This module contains graph theory functions which allow for a graphical representation of configurations.
# These functions either read in or write out data.
export processFileName, processFile, writeClassification, writeComparison, writeProfiles, writeUnique

"""
This function extracts the information embedded in the .xyz filenames.

Inputs:
configDir: (String) The path to the directory containing the .xyz files.
fileName: (String) A filename.

Outputs:
isXYZ: (Boolean) Whether or not the file is a .xyz file
prefix: (String) The unique prefix of the fileName.
fp: (IOStream) The file pointer for the file.
L: (Int64) The number of configurations in the file.
N: (Int64) The number of atoms in each configuration.
B :(Bool) Whether configurations were found under a magnetic field of 0.3 au or not.
"""
function processFileName(configDir,fileName)
	parts = split(fileName,".xy") # Process filename to be able to tell if .xyz file
	if (last(parts) == "z") # If is an .xyz file
		fp = open("$configDir/$fileName") # Open the file
		parts = split(fileName,'_') # Process filename
		prefix = parts[1] # Get filename prefix
		N,L = (parse.(Int64,parts[2:3]))[:] # Number of configurations in the file (part of file name)
		B = false
		if parts[4] == "BField"
			B = true
		end

		return true,prefix,fp,L,N,B
	else
		return false,nothing,nothing,nothing,nothing,nothing
	end
end

"""
This function reads in all the .xyz files in a directory and performs CNA on each of the configurations for a range of rCut values.
The atomic shell numbers for configurations that are to be classified.

Inputs:
fp: (String) The file path to the .xyz file to process.
L: (Int64) The number of configurations in the file.
N: (Int64) The number of atoms in each configuration.
B: (Bool) Whether configurations were found under a magnetic field of 0.3 au or not.
compare: (Bool) Whether the configurations are to be compared or classified.
rCutRange: (Vector{Float64}) A range of cut-off radii to test, in units of equilibriumBondLength.
M: (Int64) The number of rCut values to test. Default of 1 for classification.

Outputs:
energies: (Vector{Float64}) The energy of each configuration in each the file.
totalProfiles: Array{Dict{String,Int}} 2D array containg the totalCNA profiles of each configuration, for each rCut value.
atomicProfiles: Array{Dict{String,Int}} 3D array containing the CNA profile of each configuration, for each rCut value, for each atom.
configurations: (Array{Float64}) 3D array containing the atom coordinates for each configuration in the file. Only used if classifying.
shells: (Vector{Vector{Int64}}) A vector of shell numbers of each atom in each configuration. Empty if comparing.
bondGraphs: Vector{SimpleGraph} A vector of graphical representation of each configuration. Empty if comparing.
"""
function processFile(fp,L,N,B,compare,EBL,rCutRange;M=1)

	energies = Vector{Float64}() # Initialise list of configuration energies
	configurations = Array{Float64}(undef, (L,N,3)) # Initialise 3D array for the atom coordinates for each configuration in the file.
	totalProfiles = Array{Dict{String,Int}}(undef,M,L) # Initialise array to hold total CNA profiles
	atomicProfiles = Array{Dict{String,Int}}(undef,M,L,N) # Initialise array to hold atomic CNA profiles
	shells = Vector{Vector{Int64}}() # Initialise array to hold shell number of each atom in cluster
	bondGraphs = Vector{SimpleGraph}() # Initialise vector of graphical representation of each configuration.
	coordinatesNext = false # Does next line contain coordinates?
	headerNext = true # Does next line contain the header for a new configuration?
	skipLine = false # Does next line contain useless/no information?
	configuration, E, atomNum, configNum::Int = Array{Float64}(undef, (N,3)), 0.0, 0, 0 # Initialise variables
	while (true) # While not at end of file
		if (headerNext) # If next line is the start of a new configuration
			atomNum = 0 # Initialise number of atoms processed in current configuration
			line = readline(fp) # Read first line
			configNum = parse(Int,line) # Get configuration number
			line = readline(fp) # Read blank line
			line = readline(fp) # Skip over number of atoms in the configuration as known from filename
			line = readline(fp) # Read blank line
			line = readline(fp) # Read next line
			E = round(parse(Float64,line); digits = 15) # Get the number of atoms in the configuration
			line = readline(fp) # Read blank line
			headerNext = false # Have completed reading the header
			coordinatesNext = true # Cooridinate information is on the next line
		elseif (coordinatesNext) # If next line is coordinates
			line = readline(fp) # Read the next line
			for (j,x) in enumerate((split(line, r" +"))[2:4]) # For each component of the coordinate
				configuration[atomNum+1,j] = parse(Float64,x) # Store the component
			end
			atomNum += 1 # Increment the number of atoms processed in the current configuration
			if (atomNum==N) # If have processed all the atoms in the current configuration
				coordinatesNext = false # Have finished reading in coordinates
				skipLine = true # Either eof or another configuration
			end
		elseif (skipLine) # If need to determine if at eof or there is another configuration
			line = readline(fp) # Read next line
			for (i,rCut) in enumerate(rCutRange) # For each rCut value
				if !compare # If classifying
					labels,graph = ShellLabelling(configuration,N,B,rCut,EBL) # Compute atomic shell numbers
					push!(shells,labels) # Add shell numbers for current configuration to list.
					push!(bondGraphs,graph) # Add graph for current configuration to list.
				end
				totalProfile,atomicProfile = CNA(configuration,N,rCut,B,EBL) # Compute the CNA profile of the current configuration
				totalProfiles[i,configNum] = totalProfile # Add the total CNA profile to vector
				atomicProfiles[i,configNum,:] = atomicProfile[:] # Store the atomic Profile
			end
			push!(energies,E) # Add the configuration energy to vector
			configurations[configNum,:,:] = configuration
			if (eof(fp)) # If at end of file
				break # end while loop
			elseif (isempty(line)) # If just a blank line -> another configuration follows
				line = readline(fp) # Read a blank line
				headerNext = true # New configuration header on next line
				configuration = Array{Float64}(undef, (N,3)) # Preallocate space for the configuration coordinates
			end
		end
	end
	return energies, totalProfiles, atomicProfiles, configurations, shells, bondGraphs # Return outputs
end

"""
Write out symmetry classifications for all configurations.

Inputs:
configDir: (String) The path to (or name of if in same location as this program) the directory containing the .xyz files.
filename: (String) The prefix of the file.
classifications: Vector{Dict{Vector{Int64},String}} Vector of dictionaries mapping a list of atoms to their symmetry.
capSymmetries: Vector{Dict{Dict{Int64,Vector{Int64}},String}} Vector of dictionaries mapping list of atoms forming a cap to their symmetry.
L: (Int64) The number of configurations in the file.
"""
function writeClassification(configDir,filename,classifications,capSymmetries,L)
	# Write output file contatining classifications of all configurations in the file
	open("$configDir/classification/$(filename)_classification.txt","a") do io # Open the output file corresponding to the input file
		# Write file headers
		println(io,"Config #    classification")
		for j in 1:L # all the configurations in the file
			# Print out the classifcations
			println(io,"$(rpad(j,8," "))    $(classificationToString(classifications[j],capSymmetries[j]))")
		end
	end
end

"""
This function converts the classification and capSymmetry dictionaries to a string.

Inputs:
classification: Dict{Vector{Int64},String} Dictionary mapping a list of atoms to their symmetry.
capSymmetries: Dict{Dict{Int64,Vector{Int64}},String} Dictionary mapping list of atoms forming a cap to their symmetry.

Outputs:
string: (String) Formatted string.
"""
function classificationToString(classification,capSymmetry)
	if isempty(classification) # If no symmetries identified
		string = "No core or symmetries detected"
	else # If symmetries detected
		if isempty(capSymmetry) # If no caps detected
			capSymmetry = "None dectected"
		end
		string = "Symmetries: $classification; Caps: $capSymmetry" # Format string
	end
	return string # Return string
end

"""
Write out total CNA profiles for all configurations for a particular rCut value.

Inputs:
rCutRange: (Vector{Float64}) A range of cut-off radii to test, in units of equilibriumBondLength.
configDir: (String) The path to (or name of if in same location as this program) the directory containing the .xyz files.
filename: (String) The prefix of the file.
energies: (Vector{Float64}) The energy of each configuration in each the file.
totalProfiles: Array{Dict{String,Int}} 2D array containg the totalCNA profiles of each configuration, for each rCut value.
i: (Int64) The index of the current rCut value.
L: (Int64) The number of configurations in the file.
"""
function writeProfiles(rCutRange,configDir,filename,energies,totalProfiles,i,L)
	if !isdir("$configDir/comparison/rCut_$(rCutRange[i])")
		mkdir("$configDir/comparison/rCut_$(rCutRange[i])") # Create the output directory
	end
	# Write output files for each rCut value
	open("$configDir/comparison/rCut_$(rCutRange[i])/profile_$(filename).txt","a") do io # Open the output file corresponding to the input file
		# Write file headers
		println(io,"Config #    Energy                 Total CNA Profile")
		for j in 1:L # all the configurations in the file
			# Print out the results
			println(io,"$(rpad(j,8," "))    $(rpad(energies[j], 19, " "))    $(totalProfiles[i,j])")
		end
	end
end

"""
Write out the configuration comparision results.

	Inputs:
	configDir: (String) The path to (or name of if in same location as this program) the directory containing the .xyz files.
	filename: (String) The prefix of the file.
	sortingArray: Vector{Int64} A vector of configuration numbers, sorted by maximum similarity score.
	sortedEnergies: (Vector{Float64}) The energy of each configuration, sorted by maximum similarity score.
	maxSims: Vector{Float64} A Vector of the maximum similarity score for each configuration when compared to every other configuration
							 over all rCut values, in decreasing order.
	similarConfigs: Vector{Dict{Int64,Vector{Float64}}} A vector of dictionaries which map the most similar configurations to the rCut values
														for which it is maximised, for each configuraiton, sorted by similarity.
	L: (Int64) The number of configurations in the file.
"""
function writeComparison(configDir,filename, sortingArray, sortedEnergies, maxSims, similarConfigs,L) # Write comparison output file
	open("$configDir/comparison/comparison_$(filename).txt","a") do io # Open the output file corresponding to the input file
		# Write file headers
		println(io,"Config #    Energy                 Similarity Score    $("Most similar to (Config #; rCut)")")
		for j in 1:L # all the configurations in the file
			# Print out the results
			println(io,"$(rpad(sortingArray[j],8," "))    $(rpad(sortedEnergies[j], 19, " "))    $(rpad(maxSims[j],17, " "))    $(similarConfigs[j])")
		end
	end
end

"""
Write out the distinguisable configuration set results.

	Inputs:
	configDir: (String) The path to (or name of if in same location as this program) the directory containing the .xyz files.
	filename: (String) The prefix of the file.
	numUnique (Int64) Number of distinguiable sets of configurations.
	uniqueConfigs: Vector{Vector{Int64}} Vector of configuration numbers of the distinguisable sets of configurations.
	uniqueSimilarities: Vector{Float64} Vector of similarity scores of each of the distinguiable sets.
	uniqueEnergyStats: Array{Float64} 2D array of energy statistics. First row is the mean, second row is standard deviation.
"""
function writeUnique(configDir,filename,numUnique,uniqueConfigs,uniqueSimilarities, uniqueEnergyStats)
	# Write unique configurations output file
	open("$configDir/comparison/unique_$(filename).txt","a") do io # Open the output file corresponding to the input file
		# Write file headers
		println(io,"Unique Config    Similarity Score     Energy Mean            Energy Std")
		for j in 1:numUnique # all the configurations in the file
			# Print out the results
			println(io,"$(rpad(uniqueConfigs[j],13," "))    $(rpad(uniqueSimilarities[j],17, " "))    $(rpad(uniqueEnergyStats[1,j], 19, " "))    $(uniqueEnergyStats[2,j])")
		end
	end
end

end


text:
# Autogenerated wrapper script for VMEC_jll for x86_64-linux-gnu-libgfortran4
export libvmec_mkl, libvmec_openblas

using MPICH_jll
using OpenBLAS_jll
using SCALAPACK_jll
using MKL_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("VMEC")
JLLWrappers.@declare_library_product(libvmec_mkl, "libvmec_mkl.so.0")
JLLWrappers.@declare_library_product(libvmec_openblas, "libvmec_openblas.so.0")
function __init__()
    JLLWrappers.@generate_init_header(MPICH_jll, OpenBLAS_jll, SCALAPACK_jll, MKL_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libvmec_mkl,
        "lib/libvmec_mkl.so",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libvmec_openblas,
        "lib/libvmec_openblas.so",
        nothing,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using BlockArrays, Test

@testset "block slice" begin
    A = BlockArray(1:6,1:3)
    b = parentindices(view(A, Block(2)))[1] # A BlockSlice

    @test first(b) == 2
    @test last(b) == 3
    @test length(b) == 2
    @test step(b) == 1
    @test Base.unsafe_length(b) == 2
    @test axes(b) == (Base.OneTo(2),)
    @test Base.axes1(b) == Base.OneTo(2)
    @test Base.unsafe_indices(b) == (Base.OneTo(2),)
    @test size(b) == (2,)
    @test collect(b) == [2,3]
    @test b[1] == 2
    @test b[1:2] == 2:3
end

@testset "block view" begin
    A = BlockArray(collect(1:6), 1:3)
    @test view(A, Block(2)) == [2,3]
    view(A, Block(2))[2] = -1
    @test A[3] == -1

    # backend tests
    @test_throws ArgumentError Base.to_index(A, Block(1))

    A = PseudoBlockArray(collect(1:6), 1:3)
    @test view(A, Block(2)) == [2,3]
    view(A, Block(2))[2] = -1
    @test A[3] == -1

    # backend tests
    @test_throws ArgumentError Base.to_index(A, Block(1))

    A = BlockArray(reshape(collect(1:(6*12)),6,12), 1:3, 3:5)
    V = view(A, Block(2), Block(3))
    @test size(V) == (2, 5)
    V[1,1] = -1
    @test A[2,8] == -1

    V = view(A, Block(3, 2))
    @test size(V) == (3, 4)
    V[2,1] = -2
    @test A[5,4] == -2

    # test mixed blocks and other indices
    @test view(A, Block(2), 2) == [8,9]
    @test similar(A, (Base.OneTo(5), axes(A,2))) isa BlockArray{Int}
    @test view(A, Block(2), :) == A[2:3,:]

    @test view(A, 2, Block(1)) == [2,8,14]
    @test view(A, :, Block(1)) == A[:,1:3]

    @test view(V, Block(1, 1)) ≡ V

    @test_throws BlockBoundsError view(V, Block(1,2))
    @test_throws BlockBoundsError view(V, Block(2,1))

    A = BlockArray(reshape(collect(1:(6^3)),6,6,6), 1:3, 1:3, 1:3)
    V = view(A, Block(2), Block(3), Block(1))
    @test size(V) == (2, 3, 1)
    V[1,1,1] = -3
    @test A[2,4,1] == -3

    V = view(A,Block(1,1,1))
    @test size(V) == (1,1,1)
    V[1,1,1] = -4
    @test A[1,1,1] == -4

    # blocks mimic CartesianIndex in views
    V = view(A,Block(1,1),Block(2))
    @test size(V) == (1,1,2)
    V[1,1,1] = -5
    @test A[1,1,2] == -5

    V = view(A,Block(2),Block(1,1))
    @test size(V) == (2,1,1)
    V[1,1,1] = -6
    @test A[2,1,1] == -6

    # test mixed blocks and other indices
    @test view(A, Block(1), Block(2), 1) == A[1:1,2:3,1]
    @test view(A, Block(1,2), 1) == A[1:1,2:3,1]
    @test view(A, Block(1), 2, 1) == A[1:1,2,1]
    @test view(A, 1, Block(2), 1) == A[1,2:3,1]
    @test view(A, 1, 2, Block(2)) == A[1,2,2:3]
end

@testset "block view pointers" begin
    A = BlockArray(reshape(Vector{Float64}(1:(6^2)),6,6), 1:3, 1:3)

    V = view(A, Block(1,2))
    @test Base.unsafe_convert(Ptr{Float64}, V) == Base.unsafe_convert(Ptr{Float64}, A.blocks[1,2])
    @test unsafe_load(pointer(V)) == V[1,1]

    A = PseudoBlockArray(reshape(Vector{Float64}(1:(6^2)),6,6), 1:3, 1:3)

    V = view(A, Block(1,2))
    @test Base.unsafe_convert(Ptr{Float64}, V) == Base.unsafe_convert(Ptr{Float64}, view(A.blocks, 1:1, 2:3))
    @test unsafe_load(pointer(V)) == V[1,1]

    V = view(A, Block(2), 2:3)
    @test Base.unsafe_convert(Ptr{Float64}, V) == Base.unsafe_convert(Ptr{Float64}, view(A.blocks, 2:3, 2:3))
    @test unsafe_load(pointer(V)) == V[1,1]

    V = view(A, 2:3, Block(2))
    @test Base.unsafe_convert(Ptr{Float64}, V) == Base.unsafe_convert(Ptr{Float64}, view(A.blocks, 2:3, 2:3))
    @test unsafe_load(pointer(V)) == V[1,1]
end

@testset "block indx range of block range" begin
    A = PseudoBlockArray(collect(1:6), 1:3)
    V = view(A, Block.(1:2))
    @test V == 1:3
    @test axes(V,1) isa BlockArrays.BlockedUnitRange
    @test blockaxes(V,1) == Block.(1:2)
    @test view(V, Block(2)[1:2]) == [2,3]
    V = view(A, Block.(2:3))
    @test V == 2:6
    @test view(V, Block(2)[1:2]) == [4,5]
end

@testset "subarray implements block interface" begin
    A = PseudoBlockArray(reshape(Vector{Float64}(1:(6^2)),6,6), 1:3, 1:3)

    V = view(A, Block(2,3))
    @test PseudoBlockArray(V) isa PseudoBlockArray
    @test BlockArray(V) isa BlockArray
    @test PseudoBlockArray(V) == BlockArray(V) == V

    V = view(A, Block(2), Block.(2:3))
    @test PseudoBlockArray(V) isa PseudoBlockArray
    @test BlockArray(V) isa BlockArray
    @test PseudoBlockArray(V) == BlockArray(V) == V
    @test blocksize(V) == (1,2)

    V = view(A, Block.(2:3), Block(3))
    @test PseudoBlockArray(V) isa PseudoBlockArray
    @test BlockArray(V) isa BlockArray
    @test PseudoBlockArray(V) == BlockArray(V) == V
    @test blocksize(V) == (2,1)

    V = view(A, Block.(2:3), Block.(1:2))
    @test PseudoBlockArray(V) isa PseudoBlockArray
    @test BlockArray(V) isa BlockArray
    @test PseudoBlockArray(V) == BlockArray(V) == V
    @test blocksize(V) == (2,2)
end


text:
# https://selalib.github.io/vlasov-maxwell.html
# see reference here for conditions

@info "starting Weibel model"

e = species.e

function Weibel(vth, k, β, α, Tr)

    function P(x,v)
        if !(v isa Array)
            v = [v]
        end

        if !(x isa Array)
            x = [x]
        end

        v = sum(v .^2)
        x = sum(x .^2)

        1/(π*vth*sqrt(Tr)) * exp(-(v/Tr) / vth^2) * (1 + α*cos(k*x))
    end
end

D_e = Distribution(Weibel(0.02, 1.25, 10e-4, 10e-4, 48), e)
G = Geometry()

plasma = ElectrostaticPlasma([D_e], G)

sol = Plasma.solve(plasma, dim=3, GPU=false, time_ub = 50.0, ub=2*π)

@save "weibel.bson" sol

Plasma.plot(sol, 0.01)


text:
@generated function _object_and_preserve(x::T) where {T}
  # q = if sizeof(T) > 64
    # Expr(:block, :(rx = $(ManualMemory.Reference)(x)), :((rx, rx)))
  # elseif isbitstype(T)
  q = isbitstype(T) ? :((x, nothing)) : :((x, x))
  Expr(:block, Expr(:meta,:inline), q)
end
@inline object_and_preserve(x) = _object_and_preserve(x)

@inline object_and_preserve(x::String) = (x, x)
@inline object_and_preserve(A::AbstractArray{T}) where {T<:NativeTypes} = array_object_and_preserve(ArrayInterface.device(A), A)
@inline object_and_preserve(A::AbstractArray) = (A, A)
@inline array_object_and_preserve(::ArrayInterface.CPUPointer, A::AbstractArray) = (PtrArray(A), preserve_buffer(A))
@inline array_object_and_preserve(_, A::AbstractArray) = _object_and_preserve(A)

@inline object_and_preserve(r::Base.RefValue{T}) where {T} = (Base.unsafe_convert(Ptr{T}, r), r)

@inline function object_and_preserve(t::Tuple{T}) where {T}
  o, p = object_and_preserve(only(t))
  (o,), (p,)
end
@inline function object_and_preserve(t::Tuple{T1,T2,Vararg}) where {T1,T2}
  oh, ph = object_and_preserve(first(t))
  ot, pt = object_and_preserve(Base.tail(t))
  (oh, ot...), (ph, pt...)
end


text:
vector(x::Vector{T}) where {T} = x
vector(x::Array{T,N}) where {T,N} = vec(x)

function vector(x::NTuple{N,T}) where {N,T}
     result = Vector{T}(undef, N)
     for i=1:N
         result[i] = x[i]
     end
     return result
end

"""
     coordvecs(seq_of_coord, ::Type{T}=Float64)

```
vec_of_coord_tuples = [(1, 2), (2, 4)]
coordvecs(vec_of_coord_tuples) == [[1.0, 2.0], [2.0, 4.0]]

vec_of_coord_vectors = [[1, 2], [2, 4]]
coordvecs(vec_of_coord_vectors) == [[1.0, 2.0], [2.0, 4.0]]

tuple_of_coord_tuples = ((1, 2), (2, 4))
coordvecs(tuple_of_coord_tuples) == ([1.0, 2.0], [2.0, 4.0])
coordvecs(tuple_of_coord_tuples, Float32) == ([1.0f0, 2.0f0], [2.0f0, 4.0f0])
```
"""
coordvecs(coordseq::Vector{Vector{T}}) where {T<:Integer} = coordvecs(coordseq, Float64)
coordvecs(coordseq::Vector{NTuple{N,T}}) where {N, T<:Integer} = coordvecs(coordseq, Float64)
coordvecs(coordseq::NTuple{N1, NTuple{N2,T}}) where {N1, N2, T<:Integer} = coordvecs(coordseq, Float64)
coordvecs(coordseq::NTuple{N, Vector{T}}) where {N, T<:Integer} = coordvecs(coordseq, Float64)

coordvecs(coordseq::Vector{Vector{T}}) where {T<:AbstractFloat} = coordvecs(coordseq, T)
coordvecs(coordseq::Vector{NTuple{N,T}}) where {N, T<:AbstractFloat} = coordvecs(coordseq, T)
coordvecs(coordseq::NTuple{N1, NTuple{N2,T}}) where {N1, N2, T<:AbstractFloat} = coordvecs(coordseq, T)
coordvecs(coordseq::NTuple{N, Vector{T}}) where {N, T<:AbstractFloat} = coordvecs(coordseq, T)

coordvecs(coordseq::Vector{Vector{T1}}, ::Type{T2}) where {T1<:Real, T2<:AbstractFloat} = map(x->T2.(x), coordseq)
coordvecs(coordseq::Vector{NTuple{N,T1}}, ::Type{T2}) where {N, T1<:Real, T2<:AbstractFloat} = map(x->[T2.(x)...,], coordseq)
coordvecs(coordseq::NTuple{N1, NTuple{N2,T1}}, ::Type{T2}) where {N1, N2, T1<:Real, T2<:AbstractFloat} = map(x->[T2.(x)...,], coordseq)
coordvecs(coordseq::NTuple{N, Vector{T1}}, ::Type{T2}) where {N, T1<:Real, T2<:AbstractFloat} = map(x->[T2.(x)...,], coordseq)

vectorofcoordvecs(coordseq) = vector(coordvecs(coordseq))
matrixfromcoords(coordseq) = permutedims(reduce(hcat, vectorofcoordvecs(coordseq)))


text:
# Various algorithms for computing quantile

function quantile_bisect(d::ContinuousUnivariateDistribution, p::Float64,
                         lx::Float64, rx::Float64, tol::Float64)

    # find quantile using bisect algorithm
    cl = cdf(d, lx)
    cr = cdf(d, rx)
    @assert cl <= p <= cr
    while rx - lx > tol
        m = 0.5 * (lx + rx)
        c = cdf(d, m)
        if p > c
            cl = c
            lx = m
        else
            cr = c
            rx = m
        end
    end
    return 0.5 * (lx + rx)
end

quantile_bisect(d::ContinuousUnivariateDistribution, p::Float64) =
    quantile_bisect(d, p, minimum(d), maximum(d), 1.0e-12)

# if starting at mode, Newton is convergent for any unimodal continuous distribution, see:
#   Göknur Giner, Gordon K. Smyth (2014)
#   A Monotonically Convergent Newton Iteration for the Quantiles of any Unimodal
#   Distribution, with Application to the Inverse Gaussian Distribution
#   http://www.statsci.org/smyth/pubs/qinvgaussPreprint.pdf

function quantile_newton(d::ContinuousUnivariateDistribution, p::Float64, xs::Float64=mode(d), tol::Float64=1e-12)
    if 0.0 < p < 1.0
        while true
            x = xs + (p - cdf(d, xs)) / pdf(d, xs)
            abs(x-xs) >= max(abs(x),abs(xs))*tol || return x
            xs = x
        end
    elseif p == 0.0
        return minimum(d)
    elseif p == 1.0
        return maximum(d)
    else
        return NaN
    end
end

function cquantile_newton(d::ContinuousUnivariateDistribution, p::Float64, xs::Float64=mode(d), tol::Float64=1e-12)
    if 0.0 < p < 1.0
        while true
            x = xs + (ccdf(d, xs)-p) / pdf(d, xs)
            abs(x-xs) >= max(abs(x),abs(xs))*tol || return x
            xs = x
        end
    elseif p == 1.0
        return minimum(d)
    elseif p == 0.0
        return maximum(d)
    else
        return NaN
    end
end

function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Float64, xs::Float64=mode(d), tol::Float64=1e-12)
    if -Inf < lp < 0.0
        if lp < logcdf(d,xs)
            while true
                x = xs - exp(lp - logpdf(d,xs) + logexpm1(max(logcdf(d,xs)-lp,0.0)))
                abs(x-xs) >= max(abs(x),abs(xs))*tol || return x
                xs = x
            end
        else
            while true
                x = xs + exp(lp - logpdf(d,xs) + log1mexp(min(logcdf(d,xs)-lp,0.0)))
                abs(x-xs) >= max(abs(x),abs(xs))*tol || return x
                xs = x
            end
        end
    elseif lp == -Inf
        return minimum(d)
    elseif lp == 0.0
        return maximum(d)
    else
        return NaN
    end
end

function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Float64, xs::Float64=mode(d), tol::Float64=1e-12)
    if -Inf < lp < 0.0
        if lp < logccdf(d,xs)
            while true
                x = xs + exp(lp - logpdf(d,xs) + logexpm1(max(logccdf(d,xs)-lp,0.0)))
                abs(x-xs) >= max(abs(x),abs(xs))*tol || return x
                xs = x
            end
        else
            while true
                x = xs - exp(lp - logpdf(d,xs) + log1mexp(min(logccdf(d,xs)-lp,0.0)))
                abs(x-xs) >= max(abs(x),abs(xs))*tol || return x
                xs = x
            end
        end
    elseif lp == -Inf
        return maximum(d)
    elseif lp == 0.0
        return minimum(d)
    else
        return NaN
    end
end

# A macro: specify that the quantile (and friends) of distribution D
# is computed using the newton method
macro quantile_newton(D)
    esc(quote
        quantile(d::$D, p::Float64) = quantile_newton(d,p)
        cquantile(d::$D, p::Float64) = cquantile_newton(d,p)
        invlogcdf(d::$D, lp::Float64) = invlogcdf_newton(d,lp)
        invlogccdf(d::$D, lp::Float64) = invlogccdf_newton(d,lp)
    end)
end


text:
@testset "reverse_view" begin

    a = randn((3,4,3,5))
    f(a, dims) = @test reverse(a, dims=dims) == reverse_view(a, dims=dims)

    @test reverse(a) == reverse_view(a)
    for i in [1,2,4, (1,2), (3,4), (1,2,3,4), (1,4)]
        f(a, i)
    end
end


text:
# This file was generated, do not modify it. # hide
y1 = log.(y0)
y = transform(fit!(machine(UnivariateStandardizer(), y1),
                   verbosity=0), y1);

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
# Example of the piecewise linear regration with Mosek, MosekTools, and JuMP
# See Mosek manual 9.2.3 Convex piecewise linear regression

using JuMP, Mosek, MosekTools
using RDatasets, Plots

# prepare data
data = dataset("datasets", "pressure");
Y = data.Temperature ./1000;    Y = Y[8:end]
X = data.Pressure ./1000;       X = X[8:end]
n = length(X)
k = 3
(M_lb, M_ub) = (-5, 5)

# build and solve the model
m = Model(Mosek.Optimizer)
@variable(m, t)
@variable(m, M_lb <= f[1:n] <= M_ub)
@variable(m, a[1:k])
@variable(m, b[1:k])
@variable(m, z[1:n,1:k], Bin)
@objective(m, Min,  t)
@constraint(m, [t;Y - f] in SecondOrderCone())
@constraint(m, cut_one[i=1:n,j=1:k], (a[j]*X[i]+b[j]) - f[i] <= M_ub * (1 - z[i,j]))
@constraint(m, cut_two[i=1:n,j=1:k], (a[j]*X[i]+b[j]) - f[i] >= M_lb * (1 - z[i,j]))
@constraint(m, sum_bin[i=1:n], sum(z[i,:]) == 1)
optimize!(m)
status=termination_status(model)

@show a_ = JuMP.value.(a)
@show b_ = JuMP.value.(b)

# make a plot
Plots.gr(ylims=(0,1))
Plots.scatter(X,Y,legend=false)
for j in 1:k
  f_z(x) = a_[j]*x + b_[j]
  Plots.plot!(X, f_z)
end
current()