text:
# =============================================================================
# =
# = FlatZinc model.
# =
# =============================================================================

mutable struct VariableInfo
    index::MOI.VariableIndex
    name::String
    set::MOI.AbstractSet
end

mutable struct ConstraintInfo
    index::MOI.ConstraintIndex
    # FlatZinc does not allow constraint names.
    f::MOI.AbstractFunction
    s::MOI.AbstractSet
    output_as_part_of_variable::Bool
end

mutable struct Model <: MOI.AbstractOptimizer
    # A mapping from the MOI.VariableIndex to the variable object.
    # VariableInfo also stores some additional fields like the type of variable.
    variable_info::CleverDicts.CleverDict{MOI.VariableIndex, VariableInfo}

    # A mapping from the MOI.ConstraintIndex to the variable object.
    # ConstraintInfo also stores some additional fields like the type of
    # constraint. Deletion of constraints is not supported (having a vector
    # ensures that the order in which constraints are added is respected when
    # outputting the model, which makes testing easier).
    constraint_info::Vector{ConstraintInfo}

    # For FlatZinc-level sets/arrays, map the constraints that are defined on
    # sets/arrays to the assigned name of the set/array. This structure is
    # built when outputting the FZN file.
    sets_id::Dict{MOI.ConstraintIndex, String}
    arrs_id::Dict{MOI.ConstraintIndex, String}

    # Memorise the objective sense and the function separately.
    # The function can only be a single variable, as per FlatZinc limitations.
    objective_sense::MOI.OptimizationSense
    objective_function::Union{Nothing, MOI.VariableIndex}

    # Map names to MOI variable indices.
    name_to_var::Dict{String, MOI.VariableIndex}

    """
        Model()

    Create a new Model object.
    """
    function Model()
        model = new()
        model.variable_info =
            CleverDicts.CleverDict{MOI.VariableIndex, VariableInfo}()
        model.constraint_info = ConstraintInfo[]

        model.objective_sense = MOI.FEASIBILITY_SENSE
        model.objective_function = nothing

        model.sets_id = Dict{MOI.ConstraintIndex, String}()
        model.arrs_id = Dict{MOI.ConstraintIndex, String}()

        model.name_to_var = Dict{String, MOI.VariableIndex}()

        MOI.empty!(model)
        return model
    end
end

function Base.show(io::IO, ::Model)
    print(io, "A FlatZinc (fzn) model")
    return
end

function MOI.empty!(model::Model)
    empty!(model.variable_info)
    empty!(model.constraint_info)

    model.objective_sense = MOI.FEASIBILITY_SENSE
    model.objective_function = nothing

    model.sets_id = Dict{MOI.ConstraintIndex, String}()
    model.arrs_id = Dict{MOI.ConstraintIndex, String}()

    model.name_to_var = Dict{String, MOI.VariableIndex}()

    return
end

function MOI.is_empty(model::Model)
    !isempty(model.variable_info) && return false
    !isempty(model.constraint_info) && return false
    model.objective_sense != MOI.FEASIBILITY_SENSE && return false
    model.objective_function !== nothing && return false
    return true
end

# Set the objective.

function MOI.get(model::Model, ::MOI.ObjectiveFunction{MOI.VariableIndex})
    return model.objective_function
end

function MOI.set(
    model::Model,
    ::MOI.ObjectiveFunction{MOI.VariableIndex},
    f::MOI.VariableIndex,
)
    model.objective_function = f
    return nothing
end

function MOI.get(model::Model, ::MOI.ObjectiveSense)
    return model.objective_sense
end

function MOI.set(
    model::Model,
    ::MOI.ObjectiveSense,
    s::MOI.OptimizationSense,
)
    model.objective_sense = s
    return nothing
end

# Helpers.

function _create_variable(
    model::Model,
    set::Union{MOI.AbstractScalarSet, MOI.Reals},
)
    index = CleverDicts.add_item(
        model.variable_info,
        VariableInfo(MOI.VariableIndex(0), "", set),
    )
    model.variable_info[index].index = index
    return index
end

function _create_constraint(
    model::Model,
    f::F,
    set::S,
    as_part_of_variable::Bool,
) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    index = MOI.ConstraintIndex{F, S}(length(model.constraint_info) + 1)
    push!(
        model.constraint_info,
        ConstraintInfo(index, f, set, as_part_of_variable),
    )
    return index
end

# Fallback for copying operations.
function MOI.copy_to(dest::Model, src::MOI.ModelLike)
    return MOIU.default_copy_to(dest, src)
end

function MOI.supports_incremental_interface(::Model)
    return true
end

# Names.
# No support for constraint names in fzn, hence no ConstraintName.

function MOI.get(model::Model, ::MOI.VariableName, v::MOI.VariableIndex)
    return model.variable_info[v].name
end

function MOI.get(model::Model, ::Type{MOI.VariableIndex}, n::String)::MOI.VariableIndex
    # TODO: test this.
    # TODO: not terribly efficient, and used in FZN for parsing a solution (finding a variable by its name). Rather build a dict to retrieve the indices faster.
    for i in keys(model.variable_info)
        if model.variable_info[i].name == n
            return i
        end
    end
    return nothing
end

function MOI.set(
    model::Model,
    ::MOI.VariableName,
    v::MOI.VariableIndex,
    name::AbstractString,
)
    model.variable_info[v].name = string(name)
    return
end

# Variables.

function MOI.supports_add_constrained_variable(
    ::Model,
    ::Type{F},
) where {
    F <: Union{
        MOI.EqualTo{Float64},
        MOI.LessThan{Float64},
        MOI.Interval{Float64},
        MOI.EqualTo{Int},
        MOI.LessThan{Int},
        MOI.Interval{Int},
        MOI.EqualTo{Bool},
        MOI.ZeroOne,
        MOI.Integer,
    },
}
    return true
end

function MOI.supports_add_constrained_variables(
    ::Model,
    ::Type{F},
) where {
    F <: Union{
        MOI.EqualTo{Float64},
        MOI.LessThan{Float64},
        MOI.Interval{Float64},
        MOI.EqualTo{Int},
        MOI.LessThan{Int},
        MOI.Interval{Int},
        MOI.EqualTo{Bool},
        MOI.ZeroOne,
        MOI.Integer,
    },
}
    return true
end

function MOI.add_variable(model::Model)
    return _create_variable(model, MOI.Reals(1))
end

function MOI.add_constrained_variables(
    model::Model,
    sets::AbstractVector{<:MOI.AbstractScalarSet},
)
    # TODO: memorise that these variables are part of the same call, so that
    # the generated FlatZinc is shorter (array of variables)? This would
    # require that all sets are identical, though.
    vidx = [_create_variable(model, sets[i]) for i in 1:length(sets)]
    cidx = [
        _create_constraint(model, vidx[i], sets[i], true) for i in 1:length(sets)
    ]
    return vidx, cidx
end

function MOI.add_constrained_variable(
    model::Model,
    set::MOI.AbstractScalarSet,
)
    # Unlike FZN, MOI does not assume that a variable >= 0.0 is any different
    # from a >= 0 (int).
    vidx = _create_variable(model, set)
    cidx = _create_constraint(model, vidx, set, true)
    return vidx, cidx
end

function MOI.is_valid(model::Model, v::MOI.VariableIndex)
    return haskey(model.variable_info, v)
end

function MOI.get(model::Model, ::MOI.NumberOfVariables)
    return length(model.variable_info)
end

function MOI.get(model::Model, ::MOI.ListOfVariableIndices)
    return collect(keys(model.variable_info))
end

# Constraints.

function MOI.is_valid(
    model::Model,
    c::MOI.ConstraintIndex{F, S},
) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    info = get(model.constraint_info, c.value, nothing)
    return info !== nothing && typeof(info.s) == S
end

function MOI.add_constraint(
    model::Model,
    f::F,
    s::S,
) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    index = MOI.ConstraintIndex{F, S}(length(model.constraint_info) + 1)
    push!(model.constraint_info, ConstraintInfo(index, f, s, false))
    return index
end

function MOI.get(
    model::Model,
    ::MOI.ConstraintFunction,
    c::MOI.ConstraintIndex{F, S},
) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    return model.constraint_info[c.value].f
end

function MOI.get(
    model::Model,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{F, S},
) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    return model.constraint_info[c.value].s
end

function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.VariableIndex},
    ::Type{MOI.LessThan{Int}},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.VariableIndex},
    ::Type{MOI.LessThan{Float64}},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.VariableIndex},
    ::Type{CP.Strictly{MOI.LessThan{Float64}}},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.VariableIndex},
    ::Type{CP.Domain{Int}},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.VariableIndex},
    ::Type{MOI.Interval{Float64}},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.VectorOfVariables},
    ::Type{CP.Element{Int}},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.VectorOfVariables},
    ::Type{CP.Element{Bool}},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.VectorOfVariables},
    ::Type{CP.Element{Float64}},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.VectorOfVariables},
    ::Type{CP.MaximumAmong},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.VectorOfVariables},
    ::Type{CP.MinimumAmong},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.ScalarAffineFunction{Int}},
    ::Type{MOI.EqualTo{Int}},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.ScalarAffineFunction{Int}},
    ::Type{MOI.LessThan{Int}},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.ScalarAffineFunction{Int}},
    ::Type{CP.DifferentFrom{Int}},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{MOI.EqualTo{Float64}},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{MOI.LessThan{Float64}},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{CP.Strictly{MOI.LessThan{Float64}}},
)
    return true
end
function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{CP.DifferentFrom{Float64}},
)
    return true
end

function MOI.get(
    model::Model,
    ::MOI.ListOfConstraintIndices{F, S},
) where {F, S}
    return [
        c.index for
        c in model.constraint_info if typeof(c.f) == F && typeof(c.s) == S
    ]
end

function MOI.get(model::Model, ::MOI.NumberOfConstraints{F, S}) where {F, S}
    count = 0
    for c in model.constraint_info
        if typeof(c.f) == F && typeof(c.s) == S
            count += 1
        end
    end
    return count
end

function MOI.get(model::Model, ::MOI.ListOfConstraintTypesPresent)
    types = Set{Tuple{Any, Any}}()
    for info in model.constraint_info
        push!(types, (typeof(info.f), typeof(info.s)))
    end
    return collect(types)
end


text:
"""
`File` has following public fields:

* `path` - path to the file
* `metadata` - dictionary parsed from JSON-sidecar
* `entities` - dictionary parsed from key-value filename

The `File` can be initialized by specifying only `path`. Other behavior can also be
tweaked accordingly with the optional parameters:

* `load_metadata::Bool` - defaults to `true`. If `true`, this expects a JSON-sidecar to
  be present alongside with every BIDS files.
* `require_modality::Bool` - defaults to `true`. If `true` this expects a modality
  presents in every file name, e.g. sub-subtest_ses-1_run-001_**T1w**.nii.gz
* `strict::Bool` - defaults to `true`. If `true`, `BIDSTools` will throw an error on
  invalid BIDS filename, this can be turned to `false` to not parse those file names
  and display a warning instead. This will result in empty dictionary in `entities`.
* `extract_from_full_path::Bool` - defaults to `true`. If `true`, will try to extract
  `sub` and `ses` from full path and append to `entities` if they can't be found in the
  parsed filename.

# Example

```julia-repl
julia> file = File("/path/to/bids/root/sub-Subtest/ses-1/anat/sub-Subtest-1_run-001_T1w.nii.gz")
File:
    path = /path/to/bids/root/sub-Subtest/ses-1/anat/sub-Amelll_ses-1_run-001_T1w.nii.gz
    metadata_exist = true
```
"""
struct File
    path::String
    # This is the metadata from JSON-sidecar (doesn't exists in BIDS-like directory)
    metadata::OrderedDict{String, Any}
    # Entities are the key-value from filename
    entities::OrderedDict{String, String}
end

function File(
    path::AbstractString;
    load_metadata::Bool=true,
    require_modality::Bool=true,
    strict::Bool=true,
    extract_from_full_path::Bool=true
)
    entities = parse_path(path, require_modality=require_modality, strict=strict)
    # Try extracting sub and ses from full path if not exists in parsed filename
    if extract_from_full_path
        if !haskey(entities, "sub")
            entities["sub"] = get_sub(
                                path,
                                from_fname=false,
                                require_modality=require_modality,
                                strict=strict
                              )
        end
        if !haskey(entities, "ses")
            entities["ses"] = get_ses(
                                path,
                                from_fname=false,
                                require_modality=require_modality,
                                strict=strict
                              )
        end
    end
    metadata = !load_metadata ? OrderedDict{String,Any}() :
               JSON.parsefile(
                   get_metadata_path(path), dicttype=OrderedDict{String,Any}
               )
    File(path, metadata, entities)
end

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

function Base.show(io::IO, file::File)
    if length(file.metadata) == 0
        print(io, "File($(repr(file.path)), no_metadata=true)")
    else
        print(io, "File($(repr(file.path)))")
    end
end

"""
    function get_metadata_path(path)

Get file path of metadata file (json sidecar) for a BIDS `path` which can be a string or
`File`.
"""
function get_metadata_path(file::File)
    get_metadata_path(file.path)
end

function get_metadata_path(path::AbstractString)
    fname_without_ext = split(basename(path), ".")[1]
    this_dirname = dirname(path)

    metadata_path = isfile(joinpath(this_dirname, fname_without_ext*".json")) ?
                    joinpath(this_dirname, fname_without_ext*".json") : nothing
    return metadata_path
end

"""
    function parse_path(path; require_modality::Bool=true, strict::Bool=true)

Function to parse path of a File to Dictionary. The argument `path` can be
a `File` or path (any `String` object implementing `AbstractString`). This function
implements `parse_fname`.

The following keyword arguments can be passed:

* `require_modality::Bool` - defaults to `true`. If `true` this expects a modality
  presents in every file name, e.g. sub-subtest_ses-1_run-001_**T1w**.nii.gz
* `strict::Bool` - defaults to `true`. If `true`, `BIDSTools` will throw an error on
  invalid BIDS filename, this can be turned to `false` to not parse those file names
  and display a warning instead. This will result in empty dictionary in `entities`.
"""
parse_path(file::File; kws...) = parse_path(file.path; kws...)

function parse_path(
    path::AbstractString; require_modality::Bool=true, strict::Bool=true
)
    fname = basename(path)
    parse_fname(fname, require_modality=require_modality, strict=strict)
end

"""
    function parse_fname(
        fname::AbstractString; require_modality::Bool=true, strict::Bool=true
    )

Function to parse filename keys and values with this following structure
`<k1>-<v1>_<k2>-<v2>_..._<kn>-<vn>`

The following keyword arguments can be passed:

* `require_modality::Bool` - defaults to `true`. If `true` this expects a modality
  presents in every file name, e.g. sub-subtest_ses-1_run-001_**T1w**.nii.gz
* `strict::Bool` - defaults to `true`. If `true`, `BIDSTools` will throw an error on
  invalid BIDS filename, this can be turned to `false` to not parse those file names
  and display a warning instead. This will result in empty dictionary in `entities`.

This function returns a dictionary containing those keys and values.
"""
function parse_fname(
    fname::AbstractString; require_modality::Bool=true, strict::Bool=true
)
    fname_without_ext = split(fname, ".")[1]
    parts = [split(part, '-') for part in split(fname_without_ext, '_')]
    if require_modality
        # Add modality entitiy
        modality = pop!(parts)
        @assert length(modality) == 1 "Got unexpected modality $(modality)"
        push!(parts, ["modality", modality[1]])
    end
    d = OrderedDict{String,String}()
    for (idx, part) in enumerate(parts)
        if length(part) != 2
            msg = """Invalid BIDS file name $fname
                     (part $(join(part, '-')) should have exactly one '-')"""
            if strict
                error(msg)
            else
                @warn msg
                return OrderedDict{String, String}()
            end
        end
        k,v = part
        if haskey(d, k)
            msg = "Invalid BIDS file name (key $k occurs twice)"
            if strict
                error(msg)
            else
                @warn msg
                return OrderedDict{String, String}()
            end
        elseif isempty(k)
            msg = "Empty key in pair $k-$v"
            if strict
                error(msg)
            else
                @warn msg
                return OrderedDict{String, String}()
            end
        end
        d[k] = v
    end
    d
end

"""
    function check_entities_meta(file::File; kws...)

Private function to check whether keywords argument in metadata or in entities. Return
false if the key could not be found anywhere or the value of given key is not the same.
"""
function check_entities_meta(file::File; kws...)
    for (k,v) in kws
        keystr = string(k)
        meta_val = get(file.entities, keystr) do
            get(file.metadata, keystr, nothing)
        end
        if isnothing(meta_val) || meta_val != v
            return false
        end
    end
    return true
end

"""
    function get_files(
        files::Vector{File}; path::Union{String, Regex, Nothing}=nothing, kws...
    )

Function to query files based on their `path`, `entities`, and `metadata`. `path` is
optional and can be passed as either `String` or `Regex`, while entities and metadata
can be passed as keyword args, i.e. `key="value"`.

In addition to querying vector of File, this can also be used to obtain desired files
in a `Layout`, `Subject`, or `Session` by simply replacing Vector{File} into desired
object.

# Example

```julia
filtered_files = get_files(files, path="anat", run="002", modality="T1w")

# Filter from layout only from entities and metadata
filtered_files = get_files(layout, run="002", modality="T1w")
```
"""
function get_files(
    files::Vector{File}; path::Union{String, Regex, Nothing}=nothing, kws...
)
    if !isnothing(path)
        filter!(x->occursin(path, x.path), files)
    end
    result = filter(files) do f
        check_entities_meta(f; kws...)
    end
    return result
end

"""
    function get_sub(
        path;
        from_fname::Bool=true,
        require_modality::Bool=true,
        strict::Bool=true
    )

Function to get subject_id from path or File object.

The following keyword arguments can be passed:

* `from_fname`::Bool - detauls to true. If `true`, only looks the `subject_id` from
  filename, otherwise, looks into full path if `subject_id` can't be found in filename
* `require_modality::Bool` - defaults to `true`. If `true` this expects a modality
  presents in every file name, e.g. sub-subtest_ses-1_run-001_**T1w**.nii.gz
* `strict::Bool` - defaults to `true`. If `true`, `BIDSTools` will throw an error on
  invalid BIDS filename, this can be turned to `false` to not parse those file names
  and display a warning instead. This will result in empty dictionary in `entities`.

Returns nothing if no subject ID found.
"""
get_sub(file::File; kws...) = get_sub(file.path; kws...)

function get_sub(
    path::AbstractString;
    from_fname::Bool=true,
    require_modality::Bool=true,
    strict::Bool=true
)
    sub_rgx =  r"[\\/]sub-(.+?)[\\/]"
    sub_match = get(
        parse_fname(basename(path), require_modality=require_modality, strict=strict),
        "sub",
        nothing
    )
    if isnothing(sub_match) && !from_fname
        sub_match = isnothing(match(sub_rgx, path)) ? nothing : match(sub_rgx, path)[1]
    end
    !isnothing(sub_match) && return sub_match
    nothing
end

"""
    function get_ses(
        path;
        from_fname::Bool=true,
        require_modality::Bool=true,
        strict::Bool=true
    )

Function to get session_id from path or File object.

The following keyword arguments can be passed:

* `from_fname`::Bool - defaults to true. If `true`, only looks the `session_id` from
  filename, otherwise, looks into full path if `session_id` can't be found in filename
* `require_modality::Bool` - defaults to `true`. If `true` this expects a modality
  presents in every file name, e.g. sub-subtest_ses-1_run-001_**T1w**.nii.gz
* `strict::Bool` - defaults to `true`. If `true`, `BIDSTools` will throw an error on
  invalid BIDS filename, this can be turned to `false` to not parse those file names
  and display a warning instead. This will result in empty dictionary in `entities`.

Returns nothing if no session ID found.
"""
get_ses(file::File; kws...) = get_ses(file.path; kws...)

function get_ses(
    path::AbstractString;
    from_fname::Bool=true,
    require_modality::Bool=true,
    strict::Bool=true
)
    ses_rgx =  r"[\\/]ses-(.+?)[\\/]"
    ses_match = get(
        parse_fname(basename(path), require_modality=require_modality, strict=strict),
        "ses",
        nothing
    )
    if isnothing(ses_match) && !from_fname
        ses_match = isnothing(match(ses_rgx, path)) ? nothing : match(ses_rgx, path)[1]
    end
    !isnothing(ses_match) && return ses_match
    nothing
end

"""
    function construct_fname(entities::AbstractDict; ext::Union{String,Nothing}=nothing)

Function to construct BIDS filename from `entities`. It is recommended to use
`OrderedDict` for this purpose to retain the order of the elements. To supply modality,
e.g. `_T1w`, use `modality` key, i.e. "modality"=>"T1w".
"""
function construct_fname(entities::AbstractDict; ext::Union{String,Nothing}=nothing)
    result_fname = ""
    for (k,v) in entities
        k == "modality" && continue
        isnothing(v) && continue
        !occursin(r"[-_]", k) ||
            throw(ArgumentError("Cannot have - or _ in bids key $k"))
        !occursin(r"[-_]", v) ||
            throw(ArgumentError("Cannot have - or _ in bids key $v"))
        result_fname = result_fname == "" ? "$k-$v" :
                       join([result_fname, "$k-$v"], "_")
    end
    result_fname = !haskey(entities, "modality") ? result_fname :
                   join([result_fname, entities["modality"]], "_")
    result_fname = isnothing(ext) ? result_fname : "$(result_fname).$ext"
    result_fname
end

ext: jl


text:
"""
    PShere

Set of points x ∈ ℝ^n in such that ||x||_p = r.
"""
struct PSphere <: Manifold{ℝ}
    p::Float64
    r::Float64
    n::Int64
end

==(M::PSphere, N::PSphere) = M.p == N.p && M.r == N.r && M.n == N.n

<(::PSphere, ::PSphere) = false
<(M1::PSphere, M2::Euclidean) = representation_size(M1) == representation_size(M2)
<(::Euclidean, ::PSphere) = false

copy(M::PSphere) = PSphere(M.p, M.r, M.n)

function representation_size(M::PSphere)
    return (M.n,)
end

function check_manifold_point(M::PSphere, x; kwargs...)
    if size(x) != (M.n,)
        return DomainError(size(x), "x should be vector of size $(M.n).")
    end
    norm_x = norm(x, M.p)
    if !isapprox(norm_x, M.r)
        return DomainError(
            norm_x - M.r,
            "p-norm of x supposed to be r, residual is $norm_x.",
        )
    end
    return nothing
end

function check_tangent_vector(M::PSphere, x, ξ; check_base_point = true, kwargs...)
    if check_base_point
        mpe = check_manifold_point(M, x; kwargs...)
        mpe === nothing || return mpe
    end

    normalvector = x .^ (M.p - 1) ./ norm(x, M.p)^(1 - M.p)
    tangent_component = dot(x, ξ)
    if !isapprox(tangent_component, 0.0)
        return DomainError(
            tangent_component,
            "p-norm of x supposed to be r, residual is $tangent_component.",
        )
    end
    return nothing
end

embedding_dimension(M::PSphere) = M.n

function name(M::PSphere; short = true)
    return short ? "𝕊($(M.p), $(M.r))" : "$(M.p)-sphere of radius $(M.r)"
end

function Base.show(io::IO, M::PSphere)
    return print(io, name(M))
end

# distance(::PSphere, p, q) = norm(p - q)

# # Euclidean to riemannian gradients, hessians at vectors.
# egrad_to_rgrad!(M::PSphere, gradf_x, x, ∇f_x) = project!(M, gradf_x, x, ∇f_x)  # ! this should factor out
# egrad_to_rgrad(M::PSphere, x, ∇f_x) = project(M, x, ∇f_x)                      # ! this should factor out

# function ehess_to_rhess!(M::PSphere, Hessf_xξ, x, ∇f_x, ∇²f_ξ, ξ)
#     return project!(M, Hessf_xξ, x, ∇²f_ξ)
# end
# function ehess_to_rhess(M::PSphere, x, ∇f_x, ∇²f_ξ, ξ)              # ! this should factor out
#     Hessf_xξ = zeros(size(x))
#     return ehess_to_rhess!(M, Hessf_xξ, x, ∇f_x, ∇²f_ξ, ξ)
# end

# # embed!(::PSphere, X, x) = copyto!(X, x)
# # embed!(::PSphere, Xi, x, ξ) = copyto!(Xi, ξ)

# embed(::PSphere, x) = x
# embed(::PSphere, x, ξ) = ξ

# exp!(::PSphere, y, x, ξ) = (@. y = x + ξ)

# @inline inner(::PSphere, x, ξ, η) = dot(ξ, η)

# @inline log!(M::PSphere, ξ, x, y) = (ξ .= y .- x)

# manifold_dimension(::PSphere{n}) where {n} = n

# function name(M::PSphere; short = true)
#     return short ? "l1-$(sum(M.nnz_coords))/$(length(M.nnz_coords))" :
#            "PSphere with $(sum(M.nnz_coords)) nnz"
# end

# norm(::PSphere, x, ξ) = norm(ξ)

# function project!(M::PSphere, ξ, x, X)
#     (@. ξ = X * M.nnz_coords)
#     return ξ
# end
# project!(M::PSphere, x, X) = (@. x = X * M.nnz_coords)

# function randomMPoint(M::PSphere)
#     n = length(M.nnz_coords)
#     return (2 * rand(n) .- 1) .* M.nnz_coords
# end

# function randomTVector(M::PSphere, x)
#     pt = randomMPoint(M)
#     pt /= norm(pt)
#     return pt
# end

# representation_size(::PSphere{n}) where {n} = (n)

# retract!(::PSphere, y, x, ξ) = (@. y = x + ξ)

# function Base.show(io::IO, M::PSphere{n}) where {n}
#     return print(io, name(M))
# end

# # zero_tangent_vector(::PSphere, ::Any)
# zero_tangent_vector!(M::PSphere, v, p) = fill!(v, 0)

# # const l1MPoint = Vector
# # const l1TVector = Vector

# # copy(m::PSphere) = PSphere(m.nnz_coords)

# # ==(m1::PSphere, m2::PSphere) = m1.nnz_coords == m2.nnz_coords

# # # Functions
# # # ---
# # manifold_dimension(M::PSphere) = sum(M.nnz_coords)

# # # metric
# # inner(M::PSphere, x::l1MPoint, ξ::l1TVector, ν::l1TVector) = dot(ξ, ν)
# # norm(M::PSphere, x::l1MPoint, ξ::l1TVector) = inner(M, x, ξ, ξ)

# # # Tangent space
# # function project!(M::PSphere, ξ::l1TVector, x::l1MPoint, q::AbstractArray)
# #     @. ξ = q * M.nnz_coords
# #     return ξ
# # end

# # function project(M::PSphere, x::l1MPoint, q::AbstractArray)                  # ! This should factor out, with `zro_tangent_vector`
# #     ξ = zeros(size(x))
# #     return project!(M, ξ, x, q)
# # end

# # function retract!(M::PSphere, y::l1MPoint, x::l1MPoint, ξ::l1TVector)
# #     return @. y = x + ξ
# # end
# # function retract(M::PSphere, x::l1MPoint, ξ::l1TVector)
# #     y = zeros(size(x))
# #     return retract!(M, y, x, ξ)
# # end

# # # Display
# # # ---
# # show(io::IO, M::PSphere) = print(io, M.abbreviation)


text:
"""
Approximates the area under the curve using the trapezoidal rule
Arguments:
 - f: function for the
 - x_start: starting value for x
 - x_end: ending value for x
 - steps: steps taken while integrating.
"""

function trapazoidal_area(f, x_start, x_end, steps)
	x1 = x_start
	fx1 = f(x_start)
	area = 0.0

	for _ in 1:steps
		x2 = ((x_end - x_start) / steps) + x1
		fx2 = f(x2)
		area += abs(fx2 + fx1) * (x2 - x1) / 2
		x1 = x2
		fx1 = fx2
	end

	return area
end

max_stars_repo_name: madisonestabrook/Julia

max_issues_repo_name: madisonestabrook/Julia

max_forks_repo_name: madisonestabrook/Julia


text:
module AutoMLPipeline

greet() = print("Hello World!")
export fit!, transform!, fit_transform!

include("abstracttypes.jl")
using .AbsTypes

include("utils.jl")
using .Utils

include("basefilters.jl")
using .BaseFilters

include("pipelines.jl")
using .Pipelines
export @pipeline
export @pipelinex

include("featureselector.jl")
using .FeatureSelectors

include("skpreprocessor.jl")
using .SKPreprocessors

include("decisiontree.jl")
using .DecisionTreeLearners

include("ensemble.jl")
using .EnsembleMethods

include("crossvalidator.jl")
using .CrossValidators

include("skcrossvalidator.jl")
using .SKCrossValidators

end # module


text:
using Documenter, NumericalAnalysis

makedocs(
    modules = [NumericalAnalysis],
    clean = false,
    sitename = "NumericalAnalysis.jl",
    pages = Any[
        "Home" => "index.md",
        "Functions" => Any[
            "Basic function" => "m/basic.md",
            "SEq in one variable" => "m/SEq1.md",
            "Interpolation and the Lagrange Polynomial" => "m/Polynomial.md",
            "Numerical Differentiation and Integration" => "m/NCalculus.md",
        ],
    ],
)


text:
using Test
using SubModule2

@test SubModule2.f() == 1

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
function device_timeseries_ub(ps_m::CanonicalModel,
                              ts_data::Tuple{Vector{String}, Vector{Vector{Float64}}},
                              time_range::UnitRange{Int64},
                              cons_name::Symbol,
                              var_name::Symbol)

    ps_m.constraints[cons_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, ts_data[1], time_range)

    for t in time_range, (ix, name) in enumerate(ts_data[1])

        ps_m.constraints[cons_name][name, t] = JuMP.@constraint(ps_m.JuMPmodel, 0.0 <= ps_m.variables[var_name][name, t] <= ts_data[2][ix][t])

    end

    return

end

function device_timeseries_lb(ps_m::CanonicalModel,
                              ts_data::Tuple{Vector{String}, Vector{Vector{Float64}}},
                              time_range::UnitRange{Int64},
                              cons_name::Symbol,
                              var_name::Symbol)

    ps_m.constraints[cons_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, ts_data[1], time_range)

    for t in time_range, (ix, name) in enumerate(ts_data[1])

        ps_m.constraints[cons_name][name, t] = JuMP.@constraint(ps_m.JuMPmodel, ts_data[2][ix][t] <= ps_m.variables[var_name][name, t])

    end

    return

end

function device_timeseries_param_ub(ps_m::CanonicalModel,
                                    ts_data::Tuple{Vector{String}, Vector{Vector{Float64}}},
                                    time_range::UnitRange{Int64},
                                    cons_name::Symbol,
                                    param_name::Symbol,
                                    var_name::Symbol)

    ps_m.parameters[param_name] = JuMP.Containers.DenseAxisArray{PJ.ParameterRef}(undef, ts_data[1], time_range)
    ps_m.constraints[cons_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, ts_data[1], time_range)

    for t in time_range, (ix, name) in enumerate(ts_data[1])

        ps_m.parameters[param_name][name, t] = PJ.add_parameter(ps_m.JuMPmodel, ts_data[2][ix][t]);
                                               JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][name, t] >= 0.0)
        ps_m.constraints[cons_name][name, t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][name, t] <= ps_m.parameters[param_name][name, t])

    end

    return

end

function device_timeseries_param_lb(ps_m::CanonicalModel,
                                    ts_data::Tuple{Vector{String}, Vector{Vector{Float64}}},
                                    time_range::UnitRange{Int64},
                                    cons_name::Symbol,
                                    param_name::Symbol,
                                    var_name::Symbol)

ps_m.parameters[param_name] = JuMP.Containers.DenseAxisArray{PJ.ParameterRef}(undef, ts_data[1], time_range)
ps_m.constraints[cons_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, ts_data[1], time_range)

    for t in time_range, (ix, name) in enumerate(ts_data[1])

        ps_m.parameters[param_name][name, t] = PJ.add_parameter(ps_m.JuMPmodel, ts_data[2][ix][t])
        ps_m.constraints[cons_name][name, t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.parameters[param_name][name, t] <= ps_m.variables[var_name][name, t])

    end

    return

end


text:
# @inline faces(grid::TreeGrid; filter::Function = face -> true, min_level::Int = 0, max_level::Int = typemax(Int)) = FullyThreadedTree.all_faces(grid.tree, filter = filter, min_level = min_level, max_level = max_level)
# @inline faces(grid::TreeGrid, level::Int; filter::Function = face -> true) = FullyThreadedTree.all_faces(grid.tree, level, filter = filter)

# @inline cells(grid::TreeGrid; filter::Function = cell -> true, min_level = 0, max_level = typemax(Int)) = FullyThreadedTree.cells(grid.tree; filter = filter, min_level = min_level, max_level = max_level)
# @inline cells(grid::TreeGrid, level::Int; filter::Function = cell -> true) = FullyThreadedTree.cells(grid.tree, level, filter = filter)

# TODO iterate cells/faces per level (they are sorted according to their level)

function show(grid::TreeGrid)

end


text:
using LibSodium
using Test

@testset "Public-key signatures - Ed25519 - Signing" begin

    sodium_init()

    string_to_sign = "test"

    seed_hex = "3ba564b9eb1a7e4549da027dc6dea3fb8ba0c2970f65cdcca13bd3a8dc8e7de7"
    seed = hex2bytes(seed_hex)

    # Generate public & secret key
    pk = Vector{Cuchar}(undef, crypto_sign_ed25519_PUBLICKEYBYTES)
    sk = Vector{Cuchar}(undef, crypto_sign_ed25519_SECRETKEYBYTES)
    crypto_sign_ed25519_seed_keypair(pk, sk, seed)

    @test bytes2hex(sk) == "3ba564b9eb1a7e4549da027dc6dea3fb8ba0c2970f65cdcca13bd3a8dc8e7de7e539bc6e0faa5df63176178f8809bfcbd4a6ac3178803a89e697ebe309384879"
    @test sizeof(sk) == 64

    signed_message = Vector{Cuchar}(undef, crypto_sign_ed25519_BYTES + sizeof(string_to_sign))
    signed_message_len = Ref{Culonglong}(0)

    signature = Vector{Cuchar}(undef, crypto_sign_ed25519_BYTES)
    signature_len = Ref{Culonglong}(0)

    crypto_sign_ed25519(signed_message, signed_message_len, string_to_sign, sizeof(string_to_sign), sk)
    crypto_sign_ed25519_detached(signature, signature_len, string_to_sign, sizeof(string_to_sign), sk)

    @test bytes2hex(signed_message) == "13247659998dd51b975a273d39f646bbbf1e48b61168605b6c727803cef30c6ef5370e9f4cb8faf2cb773fc82886bffe5d5603165293b62c83edc34a86af0d0774657374"
    @test signed_message_len[] == 68
    @test bytes2hex(signature) == "13247659998dd51b975a273d39f646bbbf1e48b61168605b6c727803cef30c6ef5370e9f4cb8faf2cb773fc82886bffe5d5603165293b62c83edc34a86af0d07"
    @test signature_len[] == 64

end


text:
#
#  Initialise.jl
#  CollagenPolymer
#
#  Created by Christopher Revell on 01/05/2020.
#
#

module Initialise

# Import Julia packages
using Random
using Distributions
using LinearAlgebra
using StaticArrays

# Import local program modules

function initialise(nMonomers,nDomains,nParticles,domainLength,boxSize)

    # Allocate system arrays
    pos            = Array{SVector{3,Float64}}(undef,nParticles)            # xyz positions of all particles
    F              = Array{SVector{3,Float64}}(undef,nParticles)            # xyz dimensions of all forces applied to particles
    fill!(F,@SVector zeros(3))
    magsF          = zeros(Float64,nParticles)                              # Vector of force magnitudes for all particules
    W              = Array{SVector{3,Float64}}(undef,nParticles)            # xyz values of stochastic Wiener process for all particles
    magsW          = zeros(Float64,nParticles)                              # Vector of force magnitudes for all particules
    pairsList      = Tuple{Int64, Int64}[]                                  # Array storing tuple of particle interaction pairs eg pairsList[2]=(1,5) => 2nd element of array shows that particles 1 and 5 are in interaction range
    boundaryList   = Tuple{Int64,Int64,Int64}[]                             # Array storing list of particles in boundary cells, with 2nd and 3rd components of tuples storing which part of boundary cell is at
    dxMatrix       = SMatrix{3,3}(Matrix(1I, 3, 3))                         # Identity matrix for later calculations
    neighbourCells = MVector{13}(Vector{Tuple{Int64,Int64,Int64}}(undef, 13))     # Vector storing 13 neighbouring cells for a given cell

    # Arrays to prevent reallocation
    initialx = zeros(3)
    dx       = zeros(3)

    # Create each monomer
    for ii=0:nMonomers-1
        notFound = true
        # Loop to find random monomer position and orientation such that it will fit within system box
        while notFound
            initialx .= rand(Uniform(-boxSize/2.0,boxSize/2.0),3)
            dx       .= rand(Float64,3).-0.5
            normalize!(dx)
            if false in (-boxSize/2.0 .< (initialx .+ dx.*(nDomains-1)*domainLength) .< boxSize/2.0)
                # Repeat
            else
                notFound = false
            end
        end
        # Once position and orientation is found, initialise all particles within monomer
        for jj=1:nDomains
            pos[ii*nDomains+jj] = SVector{3}(initialx .+ (jj-1)*domainLength.*dx)
        end
    end

    @info "Initialised $nMonomers monomers with $(nDomains*nMonomers) particles"

    # Create random number generators for each thread
    # threadRNG = Vector{Random.MersenneTwister}(undef, nthreads())
    # for i in 1:nthreads()
    #     threadRNG[i] = MersenneTwister()
    # end

    threadRNG = MersenneTwister()

    # Allocate variables to reuse in calculations and prevent memory reallocations
    AA     = @SVector zeros(3)
    BB     = @SVector zeros(3)
    CC     = @SVector zeros(3)
    DD     = @SVector zeros(3)
    EE     = @SVector zeros(3)

    return pos,F,W,magsF,magsW,pairsList,boundaryList,dxMatrix,neighbourCells,threadRNG,AA,BB,CC,DD,EE

end

export initialise

end

max_stars_repo_name: chris-revell/collagen-model-jl

max_issues_repo_name: chris-revell/collagen-model-jl

max_forks_repo_name: chris-revell/collagen-model-jl


text:
module TestLazyInverses
using LinearAlgebra
using LazyInverses
using Test

@testset "inverse" begin
    n, k = 3, 2
    element_types = (Float32, Float64, ComplexF32, ComplexF64)
    for elty in element_types
        @testset "eltype $elty" begin
            R = randn(elty, n, n)
            matrix_types = ["generic", "pos. def.", "diagonal"]
            matrices = [R'R + I, R + I, Diagonal(randn(elty, n))]  # loop through different matrix types
            for (A, Aty) in zip(matrices, matrix_types)
                @testset "$Aty matrix" begin
                    Inv = inverse(A)
                    @test Inv isa Inverse{<:elty}
                    @test A*Inv ≈ (one(elty)*I)(n)
                    @test Inv*A ≈ (one(elty)*I)(n)

                    M = inv(A)
                    @test Matrix(adjoint(Inv)) ≈ adjoint(M)
                    @test Matrix(transpose(Inv)) ≈ transpose(M)

                    @test Matrix(Inv') ≈ inv(A)'
                    @test AbstractMatrix(Inv') ≈ inv(A)'
                    @test diag(Inv) ≈ diag(inv(A))
                    @test diag(Inv') ≈ diag(inv(A)')

                    isposdef_A = isposdef(A)
                    @test isposdef(Inv) == isposdef_A
                    @test ishermitian(Inv) == ishermitian(A)
                    @test issymmetric(Inv) == issymmetric(A)

                    # determinant
                    @test det(Inv) ≈ det(M)
                    if isposdef_A # log det errors on non-positive definite systems
                        @test logdet(Inv) ≈ logdet(M)
                    end
                    @test all(logabsdet(Inv) .≈ logabsdet(M))

                    # factorize
                    @test factorize(Inv) ≡ Inv # no-op

                    # algebra - multiplication and division with scalars, vectors, and matrices
                    # first, test scaling by numbers
                    @testset "scalar algebra" begin
                        x = randn(elty)
                        # multiplication of Inverse with scalar
                        @test Inv * x isa Inverse
                        @test Matrix(Inv * x) ≈ M * x
                        @test x * Inv isa Inverse
                        @test Matrix(x * Inv) ≈ x * M
                        # multiplication of adjoint of Inverse with scalar
                        @test Matrix(Inv' * x) ≈ M' * x
                        @test Matrix(x * Inv') ≈ x * M'

                        # division of Inverse by scalar
                        @test Inv / x isa Inverse
                        @test Matrix(Inv / x) ≈ M / x
                        @test x \ Inv isa Inverse
                        @test Matrix(x \ Inv) ≈ x \ M
                        # division of adjoint of Inverse by scalar
                        @test Matrix(Inv' / x) ≈ M' / x
                        @test Matrix(x \ Inv') ≈ x \ M'
                    end

                    @testset "vector algebra" begin
                        x = randn(elty, n)
                        # Inverse, multiplication with vector
                        @test Inv * x ≈ A \ x
                        @test x' * Inv ≈ x' / A
                        # Inverse, multiplication with vector
                        @test Inv \ x ≈ A * x
                        @test x' / Inv ≈ x' * A

                        # tests with adjoint of Inverse
                        @test Inv' * x ≈ A' \ x
                        @test Inv' \ x ≈ A' * x
                        @test x' * Inv' ≈ x' / A'
                        @test x' / Inv' ≈ x' * A'
                        y = randn(elty, n)
                        @test dot(x, Inv, y) ≈ dot(x, A \ y)
                    end

                    @testset "matrix algebra" begin
                        X = randn(elty, n, n)
                        @test Inv * X ≈ A \ X
                        @test Inv \ X ≈ A * X
                        @test X * Inv ≈ X / A
                        @test X / Inv ≈ X * A
                        @test Inv * X' ≈ A \ X'
                        @test Inv \ X' ≈ A * X'
                        @test X' * Inv ≈ X' / A
                        @test X' / Inv ≈ X' * A
                    end

                    @testset "in-place algebra" begin
                        X = randn(elty, n, n)
                        Y = zero(X)
                        InvFact = Inverse(factorize(A))
                        @test mul!(Y, InvFact, X) ≈ InvFact * X
                        # 4 and 5 arg mul!
                        α, β = randn(2)
                        @test mul!(Y, InvFact, X, α) ≈ α * (InvFact * X)
                        Y = randn(elty, n, n)
                        Z = copy(Y)
                        @test mul!(Y, InvFact, X, α, β) ≈ α * (InvFact * X) + β * Z
                        Z = copy(X)
                        @test lmul!(InvFact, Z) ≈ Inv * X
                        Z = copy(X')
                        @test rmul!(Z, InvFact) ≈ X' * Inv

                        Z = copy(Y)
                        @test mul!(Z, X, InvFact) ≈ X * InvFact
                        α, β = randn(2)
                        @test mul!(Y, X, InvFact, α) ≈ α * (X * InvFact)
                        Z = copy(Y)
                        @test mul!(Y, X, InvFact, α, β) ≈ α * (X * InvFact) + β * Z

                        @test ldiv!(Y, Inv, X) ≈ Inv \ X
                        Z = copy(X)
                        @test ldiv!(Inv, Z) ≈ Inv \ X
                        Z = copy(X')
                        @test rdiv!(Z, Inv) ≈ X' / Inv
                    end

                    # inv
                    @test inv(Inv) isa AbstractMatrix
                    @test inv(Inv) ≈ A
                    @test inv(Inverse(Inv)) isa AbstractMatrix
                    @test inv(Inverse(Inv)) ≈ inv(A)
                    @test AbstractMatrix(Inv) ≈ inv(A)
                end # matrices testset
            end # loop over matrices

            # diagonal case
            D = Diagonal(rand(elty, n))
            @test inv(Inverse(D)) isa AbstractMatrix
            @test inv(Inverse(D)) ≈ D
            @test AbstractMatrix(Inverse(D)) isa Diagonal
            @test Matrix(Inverse(D)) isa Matrix

            # 1 x 1 matrix
            x = rand(elty, (1, 1))
            @test inverse(x) isa elty
            @test inverse(x) ≈ inv(x)[1] # smart pseudo-constructor returns scalar on 1x1 matrix

            # scalar
            x = rand(elty)
            @test Inverse(x) * x ≈ 1
            @test x * Inverse(x) ≈ 1
            @test inverse(x) ≈ 1/x
            @test AbstractMatrix(Inverse(x)) ≈ fill(1/x, (1, 1))
        end # elty testset
    end # loop over elty
end

end # module


text:
struct Mock{TAG} end

tagof(x::Sym{Mock}) = DEFAULT_TAG
tagof(x::Sym{Mock{TAG}}) where {TAG} = TAG

ismocking(x) = false
ismocking(x::Sym{<:Mock}) = true

##################################################
# Function mocking
##################################################

function Cassette.execute(ctx::CassetteCtx, f, args...)
    if Cassette.canoverdub(ctx, f, args...)
        context, dict = ctx.metadata
        return apply(f, (get(dict, arg, arg) for arg in args)...; context=context)
    else
        return f(args...)
    end
end

function mock(f, xs::Sym{<:Mock}...; context::Context = @__context__)
    tags = map(tagof, xs)
    ps = map(AbstractInstances.oftype, tags)
    metadata = (
        context = context,
        dict = IdDict{Any,Sym}(key => val for (key, val) in zip(ps, xs))
    )
    ctx = CassetteCtx(metadata = metadata)
    result = Cassette.overdub(ctx, invoke, f, Tuple{tags...}, ps...)
    return get(ctx.metadata.dict, result, result)
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ d2e14b19-98d7-4fba-8bd9-bbbb47e002f6
begin
	bigbreak = HTML("<br>"^4);
	using PlutoUI
	using Test
end

# ╔═╡ ce3c8fde-802a-11ec-1b1a-31be58ae3474
md"""
# Conditionals
with the if keyword
In Julia, the syntax

```julia
if *condition 1*
    *option 1*
elseif *condition 2*
    *option 2*
else
    *option 3*
end
```
allows us to conditionally evaluate one of our options.

For example, we might want to implement the FizzBuzz test:
given a number, N, print
-  "Fizz" if N is divisible by 3,
-  "Buzz" if N is divisible by 5, and
-  "FizzBuzz" if N is divisible by 3 and 5.

Otherwise just print the number itself! Enter your choice for N here:
"""

# ╔═╡ fd2ee11d-1734-457b-b587-4128b474b772
@bind N PlutoUI.Slider(1:20, show_value = true)

# ╔═╡ 7b619e85-abb8-4115-a9df-b5fedd04a8e4
md"N = $N"

# ╔═╡ da9796ce-4572-4d2d-bb96-5949f52d46d0
with_terminal() do

	if (N % 3 == 0) && (N % 5 == 0)
		# `&&` means "AND"; % computes the remainder after division
	    println("FizzBuzz")
	elseif N % 3 == 0
	    println("Fizz")
	elseif N % 5 == 0
	    println("Buzz")
	else
	    println(N)
	end

end

# ╔═╡ aca5d158-2902-4b7a-bf50-6f7a722f5ea7
md"""
with ternary operators
For this last block, we could instead use the ternary operator with the syntax
```julia
a ? b : c
```
which equates to
```julia
if a
    b
else
    c
end
```
"""

# ╔═╡ 8a3c8101-0c2a-438f-ba58-b80e9b14bc10
md"
Now let's say we want to return the larger of two numbers. Give x and y values here: "

# ╔═╡ 5732d2dc-c411-4c56-9720-ce6d979fc60f
x = 10

# ╔═╡ d993ffa0-c613-4b4b-baeb-a223ce408217
y = 30

# ╔═╡ 782a5034-9184-4a55-b877-bd0bb4e4aebb
md"Using the if and else keywords, we might write:
"

# ╔═╡ 14d69119-9a68-428b-a857-56298b42aee5
if x > y
    x
else
    y
end

# ╔═╡ 46533172-0fe1-406b-a52d-23305282d43b
md"and as a ternary operator, the conditional looks like this:"

# ╔═╡ b13c9cc4-3c4e-4416-896f-d5cbdc2d51c0
(x > y) ? x : y

# ╔═╡ 8b81523f-79dc-46ba-8e4e-5dbc19d5d434
true ? 5 : 2 * 6

# ╔═╡ 7925dfa8-63b1-4fc2-beae-922fa93aa081
( false  ?  1  : 2 )

# ╔═╡ fdfb1f2e-6105-43dd-92cf-4b78e1db8e73
bigbreak

# ╔═╡ 61aa9056-ebfe-4c06-a752-d081ece68167
md"""
## Exercises

5.1

Write a conditional statement that prints "negative", "zero" or "odd" according to the input number `N`.

"""

# ╔═╡ 75efd803-be55-4d0f-b8c8-6b63dd90d68c
md"""M = $(@bind M PlutoUI.Slider(-100:100, show_value = true)) """

# ╔═╡ 8f1298ab-309f-4d4d-8165-4ee522f3b1c8
with_terminal() do
	# modify this cell to return "negative", "zero" and "positive".
	println(M)
end

# ╔═╡ 602b3c77-e46e-4592-a787-4aeed77eb059
bigbreak

# ╔═╡ c7c1ec65-49ed-4ae9-b9a9-0d82a346af81
md"""
5.2

Use the ternary operator to construct a reasonable answer from Bob in the dialog below.
"""

# ╔═╡ 283c45af-47a3-4446-b873-205d8aa568e7
md"""
**Alice:** _Is the number $M even?_

**Bob:**   _I don't know._
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[compat]
PlutoUI = "~0.7.32"
"""

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

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

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

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

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

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

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

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

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

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "0b5cfbb704034b5b4c1869e36634438a047df065"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.1"

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

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "ae6145ca68947569058866e443df69587acc1806"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.32"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╟─d2e14b19-98d7-4fba-8bd9-bbbb47e002f6
# ╟─ce3c8fde-802a-11ec-1b1a-31be58ae3474
# ╟─fd2ee11d-1734-457b-b587-4128b474b772
# ╟─7b619e85-abb8-4115-a9df-b5fedd04a8e4
# ╠═da9796ce-4572-4d2d-bb96-5949f52d46d0
# ╟─aca5d158-2902-4b7a-bf50-6f7a722f5ea7
# ╟─8a3c8101-0c2a-438f-ba58-b80e9b14bc10
# ╠═5732d2dc-c411-4c56-9720-ce6d979fc60f
# ╠═d993ffa0-c613-4b4b-baeb-a223ce408217
# ╟─782a5034-9184-4a55-b877-bd0bb4e4aebb
# ╠═14d69119-9a68-428b-a857-56298b42aee5
# ╟─46533172-0fe1-406b-a52d-23305282d43b
# ╠═b13c9cc4-3c4e-4416-896f-d5cbdc2d51c0
# ╠═8b81523f-79dc-46ba-8e4e-5dbc19d5d434
# ╠═7925dfa8-63b1-4fc2-beae-922fa93aa081
# ╟─fdfb1f2e-6105-43dd-92cf-4b78e1db8e73
# ╟─61aa9056-ebfe-4c06-a752-d081ece68167
# ╟─75efd803-be55-4d0f-b8c8-6b63dd90d68c
# ╠═8f1298ab-309f-4d4d-8165-4ee522f3b1c8
# ╟─602b3c77-e46e-4592-a787-4aeed77eb059
# ╟─c7c1ec65-49ed-4ae9-b9a9-0d82a346af81
# ╠═283c45af-47a3-4446-b873-205d8aa568e7
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: slides/4 - Conditionals.jl

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_path: slides/4 - Conditionals.jl

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_path: slides/4 - Conditionals.jl

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
# ndgrid-oneto.jl

using LazyGrids: ndgrid, ndgrid_array, GridOT
using Test: @test, @testset, @test_throws, @inferred

@testset "oneto" begin
    (L, M, N) = (5, 6, 7)
    (xa, ya, za) = @inferred ndgrid_array(1:L, 1:M, 1:N)
    (xl, yl, zl) = @inferred ndgrid(L, M, N)

    @test eltype(xl) === Int
    @test size(xl) === (L, M, N)
    @test xl isa GridOT
    @test yl isa GridOT
    @test zl isa GridOT
    @test xl == xa
    @test yl == ya
    @test zl == za
    @test xl[2,3,4] == xa[2,3,4] == 2
    @test yl[2,3,4] == ya[2,3,4] == 3
    @test zl[2,3,4] == za[2,3,4] == 4
end


text:
module Day06
export both_parts

mutable struct Node
    name     :: String
    parent   :: String
    children :: Vector{String}

    orbits_count :: Int
end

Node(name::String) = Node(name, "COM", [], -1)

mutable struct Tree
    nodes        :: Dict{String, Node}
    orbits_count :: Int
end

Base.getindex(t::Tree, name::String) = get(t.nodes, name, Node(name))
Base.setindex!(t::Tree, val, name::String) = t.nodes[name]=val
function Base.show(io::IO, t::Tree)
    for (name, node) in t.nodes
        print(io, "$(node.orbits_count) $name:\t")
        for child in node.children
            print(io, "  $child")
        end
        println(io)
    end
end

Tree() = Tree(Dict(), 0)
function Tree(edges)
    t = Tree()
    for (name1, name2) in edges
        node1 = t[name1]
        node2 = t[name2]

        node2.parent = name1
        push!(node1.children, name2)

        t[name1] = node1
        t[name2] = node2
    end
    t
end

function count!(tree)
    todo = ["COM"]
    while length(todo) > 0
        name = pop!(todo)
        node = tree[name]
        node.orbits_count  = tree[node.parent].orbits_count + 1
        tree.orbits_count += node.orbits_count
        for c in node.children
            push!(todo, c)
        end
    end
end

function read_input(fname)
    open(fname) do f
        map(readlines(f)) do line
            (name1, name2) = split(line, ")")
            (String(name1), String(name2))
        end
    end
end

function ancestors(tree, name)
    res = String[]
    node = tree[name]
    while node.name != "COM"
        parent = node.parent
        push!(res, parent)
        node = tree[parent]
    end
    res
end

function orbital_transfers(tree, n1, n2)
    a1 = ancestors(tree, n1)
    a2 = ancestors(tree, n2)
    @inbounds for i in 1:length(a1)
        j = findnext(x->x==a1[i], a2, 1)
        if j != nothing
            return i+j-2
        end
    end
end

function both_parts()
    tree = Tree(read_input("input06.dat"))
    count!(tree)
    (tree.orbits_count,
     orbital_transfers(tree, "SAN", "YOU"))
end

end


text:
"""
optimalorder(hc::Hclust, dm::Array{Float64,2})

Given a hierarchical cluster, use fast algorithm to determine optimal leaf order
minimizing the distance between adjacent leaves. This is done using a heuristic
where, when combining multi-leaf sub branches, only the outermost leaves are
compared (a maximum of 4 comparisons per intersection).

Sub branches are flipped if necessary to minimize the distance between adjacent
nodes, and then the combined branches are treated as a block for future
comparisons.

Based on:
[Bar-Joseph et. al. "Fast optimal leaf ordering for hierarchical clustering." _Bioinformatics_. (2001)](https://doi.org/10.1093/bioinformatics/17.suppl_1.S22)
"""
function optimalorder(hc::Hclust, dm::Array{Float64,2})
    ord = deepcopy(hc)
    optimalorder!(ord, dm)
    return ord
end

function optimalorder!(hc::Hclust, dm::Array{Float64,2})
    ord = hc.order
    orderleaves!(ord, hc, dm)
end

function orderleaves!(order::Vector{Int}, hcl::Hclust, dm::Array{Float64,2})
    extents = Tuple{Int,Int}[]
    for (vl, vr) in zip(hcl.merge[:,1], hcl.merge[:,2])
        (u, m, uidx, midx) = leaflocs(vl, order, extents)
        (k, w, kidx, widx) = leaflocs(vr, order, extents)
        if vl < 0 && vr < 0
            # Nothing needs to be done
        elseif vl < 0
            flp = flip1(m, k, w, dm)
            flp == 2 && reverse!(order, kidx, widx)
        elseif vr < 0
            flp = flip1(k, m, u, dm)
            flp == 2 && reverse!(order, uidx, midx)
        elseif vl > 0 && vr > 0
            flp = flip2(u, m, k, w, dm)
            (flp == 2 || flp == 4) && reverse!(order, uidx, midx)
            (flp == 3 || flp == 4) && reverse!(order, kidx, widx)
        else
            error("invalid 'merge' order in Hclust: ($vl, $vr) ")
        end
        push!(extents, (uidx, widx))
    end
end

function leaflocs(v::Int, order::Vector{Int}, extents::Vector{Tuple{Int,Int}})
    ##################
    if v < 0
        leftextent = findfirst(abs(v) .== order)
        rightextent = leftextent
    elseif v > 0
        leftextent = extents[v][1]
        rightextent = extents[v][2]
    else
        error("leaf position cannot be zero")
    end
    left = order[leftextent]
    right = order[rightextent]
    return left, right, leftextent, rightextent
end

"""
For 1 multi-leaf branch and a leaf, determine if flipping branch is required
1 = do not flip
2 = flip right
"""
function flip1(m::Int, k::Int, w::Int, dm::Array{Float64,2})
    dm[m,k] <= dm[m,w] ? 1 : 2
end

"""
For 2 multi-leaf branches, determine if one or two flips is required
1 = do not flip
2 = flip left
3 = flip right
4 = flip both
"""
function flip2(u::Int, m::Int, k::Int, w::Int, dm::Array{Float64,2})
    argmin([dm[m,k], dm[u,k], dm[m,w], dm[u,w]])
end


text:
const SYSTEM_UUIDS = Dict{Int,UUID}()

function system_uuid_fallback()
    get!(SYSTEM_UUIDS, myid()) do
        uuid_file = joinpath(tempdir(), "dagger-system-uuid")
        if !isfile(uuid_file)
            uuid = uuid4()
            open(uuid_file, "w") do io
                write(io, string(uuid))
            end
            return uuid
        else
            return parse(UUID, read(uuid_file, String))
        end
    end
    SYSTEM_UUIDS[myid()]
end

if Sys.islinux()
    system_uuid() = system_uuid_fallback()
    # TODO: dmidecode
else
    system_uuid() = system_uuid_fallback()
end

"Gets the unique UUID for the server hosting worker `wid`."
function system_uuid(wid::Integer)
    get!(SYSTEM_UUIDS, wid) do
        remotecall_fetch(system_uuid, wid)
    end
end


text:
function GasOutflow(Depth,Thp,GasPVT, Qmax; di=2.99, Ts=80, ∇T=1, ϵ=0.0006, Qn=10, Dn=50, tol=0.05, sg=0.6, Inc=0)

    # Create the range of Depth and Rates
    DepthRange=range(0, stop=Depth, length=Dn)
    if size(Qmax,1)==1
        Qrange=range(100, stop=Qmax, length=Qn)
    else
        Qrange=Qmax
        Qn=size(Qmax,1)
    end

    #Size of sensibilities
    Thpn=size(Thp,1)
    din=size(di,1)

        ## Gas
    Z=LinearInterpolation(GasPVT.P.*1, GasPVT.Z.*1)
    ρg=LinearInterpolation(GasPVT.P.*1, GasPVT.Rho.*1)
    βg=LinearInterpolation(GasPVT.P.*1, GasPVT.Bg.*1)
    μg=LinearInterpolation(GasPVT.P.*1, GasPVT.Mug.*1)
    Cg=LinearInterpolation(GasPVT.P.*1, GasPVT.Cg.*1)

    Ma=GasPVT.Rho[1]*10.73*(GasPVT.T[1]+460)/GasPVT.P[1]
    sg=Ma/28.96

        #Temperature
    T=LinearInterpolation([0, Depth], [Ts+460, 460+Ts+Depth*(∇T/100)])

    #create variables
    Pwf=zeros(Thpn*din,Qn)
    Pprofile=zeros(Dn,Qn,Thpn*din)
    angle=90 .- fill(Inc,Dn)
    leg=Array{String,1}(undef,Thpn*din)
       count=0
    for k=1:Thpn  #Loop for Thp
        for n=1:din #Loop for Diameter

            count=count+1
             P=zeros(Dn,Qn)      #Matrix of Pressures
            ∇P=zeros(Dn,Qn)      #Matriz of perssure gradient
            ∇Pnew=0.0
            P[1,:].=Thp[k]       #Initial pressure is Thp
            ∇P[1,:].=0.012       #Initial PRessure gradient

                #General Loop for each Rate
                for j=1:Qn

                    #General Loop for each depth
                    for i=2:Dn
                        e=tol+0.01
                        ∇Pguess=∇P[i-1,j]

                        while e>=tol

                            Pguess=∇Pguess*(DepthRange[i]-DepthRange[i-1])+P[i-1,j]
                            Nre=(4*28.97*sg*Qrange[j]*14.7)/(π*di[n]*μg(Pguess)*10.73*520)
                            f=(1/(-4*log((ϵ/3.7065)-(5.0452/Nre)*log((ϵ^1.1098/2.8257)+(7.149/Nre)^0.8981))))^2   #Friction factor
                            s=(-0.0375*sg*sin(angle[i]*π/180)*(DepthRange[i]-DepthRange[i-1]))/(Z(Pguess)*T[DepthRange[i]])
                            Pnew=sqrt(exp(-s)*P[i-1,j]^2-2.685e-3*((f*(Z(Pguess)*T[DepthRange[i]]*Qrange[j])^2 * (1-exp(-s)))/(sin(angle[i]*π/180)*di[n]^5)))
                            ∇Pnew=(Pnew-P[i-1,j])/(DepthRange[i]-DepthRange[i-1])
                            e=  abs(∇Pguess-∇Pnew)/∇Pnew
                            ∇Pguess=∇Pnew

                        end #end While
                    ∇P[i,j]=∇Pnew
                    P[i,j]=∇P[i,j]*(DepthRange[i]-DepthRange[i-1])+P[i-1,j]
                    end #End Loop Depth

                end #End Loop Rate
             Pwf[count,:].=P[end,:]
             Pprofile[:,:,count]=P
             leg[count]="Thp=$(Thp[k]); di=$(di[n])"
         end #End Loop Diameter

    end #End Loop Thp
    return DepthRange, Qrange, Pwf, leg, Pprofile
      end #End Function


text:
using Genie, Genie.Renderer.Html, Stipple, StipplePlotly

Genie.config.log_requests = false

xx = ["start", "1d", "1w", "2w", "6w"]
spaces = "       "

# Data:

y1 = [98.5, 96.0, 95.1, 94.4, 94.0]
pb1 = PlotData(
    x = xx, y = y1, xaxis="x", yaxis="y",
    error_y = ErrorBar(3 .* ones(Float64,size(y1)); color="rgba(0,0,0,0.6)"),
    plot = StipplePlotly.Charts.PLOT_TYPE_BAR, name = "Product 1",
    text = spaces .* string.(y1), textposition="auto"
)

y2 = [99.5, 94.0, 93.1, 92.4, 92.0]
pb2 = PlotData(
    x = xx, y = y2, xaxis="x", yaxis="y",
    error_y = ErrorBar(3 .* ones(Float64,size(y1)); color="rgba(0,0,0,0.6)"),
    plot = StipplePlotly.Charts.PLOT_TYPE_BAR, name = "Product 2",
    text = spaces .* string.(y2), textposition="auto"
)

y3 = [97.5, 92.0, 92.1, 91.4, 91.0]
pb3 = PlotData(
    x = xx, y = y3, xaxis="x", yaxis="y",
    error_y = ErrorBar(3 .* ones(Float64,size(y1)); color="rgba(0,0,0,0.6)"),
    plot = StipplePlotly.Charts.PLOT_TYPE_BAR, name = "Product 3",
    text = spaces .* string.(y3), textposition="auto"
)

y4 = [92.5, 88.1, 87.1, 85.9, 84.0]
pb4 = PlotData(
    x = xx, y = y4, xaxis="x", yaxis="y",
    error_y = ErrorBar(3 .* ones(Float64,size(y1)); color="rgba(0,0,0,0.6)"),
    plot = StipplePlotly.Charts.PLOT_TYPE_BAR, name = "Product 4",
    text = spaces .* string.(y4), textposition="auto"
)

plotdata = [pb1, pb2, pb3, pb4]

layout = PlotLayout(barmode="group", font=Font(16),
    title = PlotLayoutTitle(text="Product comparison", font=Font(24)),
    xaxis = [PlotLayoutAxis(xy = "x", index = 1, title="duration", font=Font(size=24), showline = true, zeroline = false)],
    yaxis = [PlotLayoutAxis(xy = "y", index = 1, ticks = "outside", title="ratio (%)", font=Font(size=24), autorange=false, range=[70.0,100.0], showline = true, zeroline = false)]
)

Base.@kwdef mutable struct Model <: ReactiveModel
  data::R{Vector{PlotData}} = plotdata, READONLY
  layout::R{PlotLayout} = layout, READONLY
  config::R{PlotConfig} = PlotConfig(), READONLY
end

model = Stipple.init(Model())

function ui()
  page(
    vm(model), class="container", [
        plot(:data, layout = :layout, config = :config)
    ]
  ) |> html
end

route("/", ui)

up()