text:
## ================================================================================ ##

const UnaryOps = (
    :precision, :exponent_mask, :signficand_mask, :exponent, :significand, :signbit,
    :iszero, :isone, :isfinite, :isinf, :isnan, :issubnormal, :isinteger)

const BinaryOps = (:(<), :(<=), :(>=), :(>), :(!=), :(==), :isless, :isequal, :cmp)

const UnaryOps_oftype = (
    :zero, :one, :ceil, :floor, :trunc, :round,
    :(+), :(-), :sign, :abs, :inv, :sqrt, :cbrt,
    :log, :log1p, :log2, :log10, :exp, :expm1, :exp2, :exp10,
    :sin, :cos, :tan, :csc, :sec, :cot, :sinpi, :cospi,
    :sincos,
    :asin, :acos, :atan, :acsc, :asec, :acot,
    :sinh, :cosh, :tanh, :csch, :sech, :coth,
    :asinh, :acosh, :atanh, :acsch, :asech, :acoth)

const BinaryOps_oftype = (
    :(+), :(-), :(*), :(/), :(\), :hypot, :flipsign, :copysign,
    :mod, :mod1, :fld, :fld1, :div, :rem, :cld)

const TrinaryOps_oftype = (
    :clamp, :muladd, :fma)

const UnaryVectorOps_oftype = (
    :norm, :sum, :prod, :normalize
)

const UnaryMatrixOps_oftype = (
    :tr,
    :sqrt, :cbrt,
    :log, :log1p, :log2, :log10, :exp, :expm1, :exp2, :exp10,
    :sin, :cos, :tan, :csc, :sec, :cot, :sinpi, :cospi,
    :sincos,
    :asin, :acos, :atan, :acsc, :asec, :acot,
    :sinh, :cosh, :tanh, :csch, :sech, :coth,
    :asinh, :acosh, :atanh, :acsch, :asech, :acoth
)

const MatrixOperations = (
    :lu, :qr, :factorize,
)

#  ================================================================================  #

for (XT, FT) in ((:XFloat16, :Float32), (:XFloat32, :Float64))
  for Op in UnaryOps
    @eval $Op(x::$XT) = $Op(reinterpret($FT, x))
  end
  for Op in BinaryOps
    @eval $Op(x::$XT, y::$XT) = $Op(reinterpret($FT, x), reinterpret($FT, y))
  end
  for Op in UnaryOps_oftype
    @eval $Op(x::$XT) = reinterpret($XT, $Op(reinterpret($FT, x)))
  end
  for Op in BinaryOps_oftype
    @eval $Op(x::$XT, y::$XT) = reinterpret($XT, $Op(reinterpret($FT, x), reinterpret($FT, y)))
  end
  for Op in TrinaryOps_oftype
    @eval $Op(x::$XT, y::$XT, z::$XT) = reinterpret($XT, $Op(reinterpret($FT, x), reinterpret($FT, y), reinterpret($FT, z)))
  end

  for Op in UnaryVectorOps_oftype
    @eval $Op(x::Vector{$XT}) = reinterpret($XT, $Op(reinterpret($FT, x)))
  end

  for Op in UnaryMatrixOps_oftype
    @eval $Op(x::Matrix{$XT}) = reinterpret($XT, $Op(reinterpret($FT, x)))
  end

  (+)(x::Matrix{XFloat16}, y::Matrix{XFloat16}) =
      reinterpret(XFloat16, :(*)(reinterpret(Float32,x), reinterpret(Float32,y)))
  (+)(x::Matrix{XFloat32}, y::Matrix{XFloat32}) =
      reinterpret(XFloat32, :(*)(reinterpret(Float64,x), reinterpret(Float64,y)))
  (-)(x::Matrix{XFloat16}, y::Matrix{XFloat16}) =
      reinterpret(XFloat16, :(*)(reinterpret(Float32,x), reinterpret(Float32,y)))
  (-)(x::Matrix{XFloat32}, y::Matrix{XFloat32}) =
      reinterpret(XFloat32, :(*)(reinterpret(Float64,x), reinterpret(Float64,y)))
end

Base.:(^)(x::XFloat16, y::XFloat16) = XFloat16(Float32(x)^(Float32(y)))
Base.:(^)(x::XFloat32, y::XFloat32) = XFloat32(Float64(x)^(Float64(y)))
Base.:(^)(x::XFloat16, y::Rational) = XFloat16(Float32(x)^(y))
Base.:(^)(x::XFloat32, y::Rational) = XFloat32(Float64(x)^(y))
Base.:(^)(x::Rational, y::XFloat16) = XFloat16(x^Float32(y))
Base.:(^)(x::Rational, y::XFloat32) = XFloat32(x^Float64(y))

Base.:(<)(x::XFloat16, y::Rational)  = Float32(x) < y
Base.:(<=)(x::XFloat16, y::Rational) = Float32(x) <= y
Base.:(>)(x::XFloat16, y::Rational)  = Float32(x) > y
Base.:(>=)(x::XFloat16, y::Rational) = Float32(x) >= y
Base.:(==)(x::XFloat16, y::Rational) = Float32(x) == y
Base.:(!=)(x::XFloat16, y::Rational) = Float32(x) != y

Base.:(<)(x::XFloat32, y::Rational)  = Float64(x) < y
Base.:(<=)(x::XFloat32, y::Rational) = Float64(x) <= y
Base.:(>)(x::XFloat32, y::Rational)  = Float64(x) > y
Base.:(>=)(x::XFloat32, y::Rational) = Float64(x) >= y
Base.:(==)(x::XFloat32, y::Rational) = Float64(x) == y
Base.:(!=)(x::XFloat32, y::Rational) = Float64(x) != y

Base.:(<)(x::Rational, y::XFloat16)  = x < Float32(x)
Base.:(<=)(x::Rational, y::XFloat16)  = x <= Float32(x)
Base.:(>)(x::Rational, y::XFloat16)  = x >= Float32(x)
Base.:(>=)(x::Rational, y::XFloat16)  = x >= Float32(x)
Base.:(==)(x::Rational, y::XFloat16)  = x == Float32(x)
Base.:(!=)(x::Rational, y::XFloat16)  = x != Float32(x)

Base.:(<)(x::Rational, y::XFloat32)  = x < Float64(x)
Base.:(<=)(x::Rational, y::XFloat32)  = x <= Float64(x)
Base.:(>)(x::Rational, y::XFloat32)  = x > Float64(x)
Base.:(>=)(x::Rational, y::XFloat32)  = x >= Float64(x)
Base.:(==)(x::Rational, y::XFloat32)  = x == Float64(x)
Base.:(!=)(x::Rational, y::XFloat32)  = x != Float64(x)

Base.:(^)(x::XFloat16, y::AbstractFloat) = XFloat16(Float32(x)^(y))
Base.:(^)(x::XFloat32, y::AbstractFloat) = XFloat32(Float64(x)^(y))
Base.:(^)(x::AbstractFloat, y::XFloat16) = XFloat16(x^Float32(y))
Base.:(^)(x::AbstractFloat, y::XFloat32) = XFloat32(x^Float64(y))

Base.:(<)(x::XFloat16, y::AbstractFloat)  = Float32(x) < y
Base.:(<=)(x::XFloat16, y::AbstractFloat) = Float32(x) <= y
Base.:(>)(x::XFloat16, y::AbstractFloat)  = Float32(x) > y
Base.:(>=)(x::XFloat16, y::AbstractFloat) = Float32(x) >= y
Base.:(==)(x::XFloat16, y::AbstractFloat) = Float32(x) == y
Base.:(!=)(x::XFloat16, y::AbstractFloat) = Float32(x) != y

Base.:(<)(x::XFloat32, y::AbstractFloat)  = Float64(x) < y
Base.:(<=)(x::XFloat32, y::AbstractFloat) = Float64(x) <= y
Base.:(>)(x::XFloat32, y::AbstractFloat)  = Float64(x) > y
Base.:(>=)(x::XFloat32, y::AbstractFloat) = Float64(x) >= y
Base.:(==)(x::XFloat32, y::AbstractFloat) = Float64(x) == y
Base.:(!=)(x::XFloat32, y::AbstractFloat) = Float64(x) != y

Base.:(<)(x::AbstractFloat, y::XFloat16)  = x < Float32(x)
Base.:(<=)(x::AbstractFloat, y::XFloat16)  = x <= Float32(x)
Base.:(>)(x::AbstractFloat, y::XFloat16)  = x >= Float32(x)
Base.:(>=)(x::AbstractFloat, y::XFloat16)  = x >= Float32(x)
Base.:(==)(x::AbstractFloat, y::XFloat16)  = x == Float32(x)
Base.:(!=)(x::AbstractFloat, y::XFloat16)  = x != Float32(x)

Base.:(<)(x::AbstractFloat, y::XFloat32)  = x < Float64(x)
Base.:(<=)(x::AbstractFloat, y::XFloat32)  = x <= Float64(x)
Base.:(>)(x::AbstractFloat, y::XFloat32)  = x > Float64(x)
Base.:(>=)(x::AbstractFloat, y::XFloat32)  = x >= Float64(x)
Base.:(==)(x::AbstractFloat, y::XFloat32)  = x == Float64(x)
Base.:(!=)(x::AbstractFloat, y::XFloat32)  = x != Float64(x)

Base.:(^)(x::XFloat16, y::Real) = XFloat16(Float32(x)^(y))
Base.:(^)(x::XFloat32, y::Real) = XFloat32(Float64(x)^(y))
Base.:(^)(x::Real, y::XFloat16) = XFloat16(x^Float32(y))
Base.:(^)(x::Real, y::XFloat32) = XFloat32(x^Float64(y))

Base.:(<)(x::XFloat16, y::Real)  = Float32(x) < y
Base.:(<=)(x::XFloat16, y::Real) = Float32(x) <= y
Base.:(>)(x::XFloat16, y::Real)  = Float32(x) > y
Base.:(>=)(x::XFloat16, y::Real) = Float32(x) >= y
Base.:(==)(x::XFloat16, y::Real) = Float32(x) == y
Base.:(!=)(x::XFloat16, y::Real) = Float32(x) != y

Base.:(<)(x::XFloat32, y::Real)  = Float64(x) < y
Base.:(<=)(x::XFloat32, y::Real) = Float64(x) <= y
Base.:(>)(x::XFloat32, y::Real)  = Float64(x) > y
Base.:(>=)(x::XFloat32, y::Real) = Float64(x) >= y
Base.:(==)(x::XFloat32, y::Real) = Float64(x) == y
Base.:(!=)(x::XFloat32, y::Real) = Float64(x) != y

Base.:(<)(x::Real, y::XFloat16)  = x < Float32(x)
Base.:(<=)(x::Real, y::XFloat16)  = x <= Float32(x)
Base.:(>)(x::Real, y::XFloat16)  = x >= Float32(x)
Base.:(>=)(x::Real, y::XFloat16)  = x >= Float32(x)
Base.:(==)(x::Real, y::XFloat16)  = x == Float32(x)
Base.:(!=)(x::Real, y::XFloat16)  = x != Float32(x)

Base.:(<)(x::Real, y::XFloat32)  = x < Float64(x)
Base.:(<=)(x::Real, y::XFloat32)  = x <= Float64(x)
Base.:(>)(x::Real, y::XFloat32)  = x > Float64(x)
Base.:(>=)(x::Real, y::XFloat32)  = x >= Float64(x)
Base.:(==)(x::Real, y::XFloat32)  = x == Float64(x)
Base.:(!=)(x::Real, y::XFloat32)  = x != Float64(x)

#= conservative

Base.typemin(::Type{XFloat16}) = XFloat16(Base.typemin(Float16))
Base.typemin(::Type{XFloat32}) = XFloat32(Base.typemin(Float32))
Base.typemax(::Type{XFloat16}) = XFloat16(Base.typemax(Float16))
Base.typemax(::Type{XFloat32}) = XFloat32(Base.typemax(Float32))

Base.floatmin(::Type{XFloat16}) = XFloat16(Base.floatmin(Float16))
Base.floatmin(::Type{XFloat32}) = XFloat32(Base.floatmin(Float32))
Base.floatmax(::Type{XFloat16}) = XFloat16(Base.floatmax(Float16))
Base.floatmax(::Type{XFloat32}) = XFloat32(Base.floatmax(Float32))

Base.eps(::Type{XFloat16}) = XFloat16(Base.eps(Float16))
Base.eps(::Type{XFloat32}) = XFloat32(Base.eps(Float32))

Base.round(x::XFloat16, r::RoundingMode{:Up}) = XFloat16(Base.round(Float16(x), RoundUp))
Base.round(x::XFloat16, r::RoundingMode{:Down}) = XFloat16(Base.round(Float16(x), RoundDown))
Base.round(x::XFloat16, r::RoundingMode{:ToZero}) = XFloat16(Base.round(Float16(x), RoundToZero))
Base.round(x::XFloat16, r::RoundingMode{:Nearest}) = XFloat16(Base.round(Float16(x), RoundNearest))

Base.round(x::XFloat32, r::RoundingMode{:Up}) = XFloat32(Base.round(Float32(x), RoundUp))
Base.round(x::XFloat32, r::RoundingMode{:Down}) = XFloat32(Base.round(Float32(x), RoundDown))
Base.round(x::XFloat32, r::RoundingMode{:ToZero}) = XFloat32(Base.round(Float32(x), RoundToZero))
Base.round(x::XFloat32, r::RoundingMode{:Nearest}) = XFloat32(Base.round(Float32(x), RoundNearest))

=#

#= aggressive =#

Base.typemin(::Type{XFloat16}) = XFloat16(Base.typemin(Float32))
Base.typemin(::Type{XFloat32}) = XFloat32(Base.typemin(Float64))
Base.typemax(::Type{XFloat16}) = XFloat16(Base.typemax(Float32))
Base.typemax(::Type{XFloat32}) = XFloat32(Base.typemax(Float64))

Base.floatmin(::Type{XFloat16}) = XFloat16(Base.floatmin(Float32))
Base.floatmin(::Type{XFloat32}) = XFloat32(Base.floatmin(Float64))
Base.floatmax(::Type{XFloat16}) = XFloat16(Base.floatmax(Float32))
Base.floatmax(::Type{XFloat32}) = XFloat32(Base.floatmax(Float64))

Base.eps(::Type{XFloat16}) = XFloat16(Base.eps(Float32))
Base.eps(::Type{XFloat32}) = XFloat32(Base.eps(Float64))

Base.round(x::XFloat16, r::RoundingMode{:Up}) = XFloat16(Base.round(Float32(x), RoundUp))
Base.round(x::XFloat16, r::RoundingMode{:Down}) = XFloat16(Base.round(Float32(x), RoundDown))
Base.round(x::XFloat16, r::RoundingMode{:ToZero}) = XFloat16(Base.round(Float32(x), RoundToZero))
Base.round(x::XFloat16, r::RoundingMode{:Nearest}) = XFloat16(Base.round(Float32(x), RoundNearest))

Base.round(x::XFloat32, r::RoundingMode{:Up}) = XFloat32(Base.round(Float64(x), RoundUp))
Base.round(x::XFloat32, r::RoundingMode{:Down}) = XFloat32(Base.round(Float64(x), RoundDown))
Base.round(x::XFloat32, r::RoundingMode{:ToZero}) = XFloat32(Base.round(Float64(x), RoundToZero))
Base.round(x::XFloat32, r::RoundingMode{:Nearest}) = XFloat32(Base.round(Float64(x), RoundNearest))

#= =#


text:
using Flux, Flux.Tracker
using Distributions

struct ModelParam
  m::TrackedArray
  log_s::TrackedArray
  support::String # real, unit, simplex, positive
  size::Integer

  function ModelParam(K::Integer, support::String)
    @assert(support in ["real", "unit", "simplex", "positive"])
    return new(param(randn(K)), param(randn(K)), support, K)
  end
end

"""
Get variational parameters
"""
function vp(mp::ModelParam)
  if mp.support in ["unit", "simplex"]
    m = sigmoid.(mp.m) .* 20.0 .- 10.0
    s = sigmoid.(mp.log_s) * 10.0
  elseif mp.support == "positive"
    m = exp.(mp.m)
    s = exp.(mp.log_s)
  else
    m = mp.m
    s = exp.(mp.log_s)
  end

  return (m, s)
end

"""
Reparameterized sampling from variational distribution
"""
function rsample(mp::ModelParam)
  m, s = vp(mp)
  return randn(mp.size) .* s .+ m
end

#= TEST
mp = ModelParam(3, "unit")
vp(mp)
rsample(mp)
=#

ext: jl


text:
"""
    generate([repo/s], url/id, [shortname]; show_failures=false)

Generates a DataDeps code block.
The only reuired parameter is the url/id.

 - `url/id` The identifier for the dataset
     - a URL for a landing page is normally best
     - Other IDs like a DOIs also work.
 - `repo/s` either a single repository/API or a list of such
     - this takes on of the `DataRepo` types exported by this package.
         - E.g `CKAN()`, or `Figshare()`.
     - If not provided, this defaults to checking all of them.
     - If only one repo is provided, and it fails, the error will be thrown.
     - If multiple repos are provided, them the metadata from all of them is combined.
 - `shortname`, the name to use in the generated DataDep
     - if not provided will use the dataset's title, but these are often very long.
 - `show_failures`, weather or not to list all the `repos` that fail and why.
     - You generally do not want to turn this on, unless you are debugging.
     - It is fine and expected for most repos to fail (after all the data is probably only on one of them)
     - If *all* repos fail, then the failure list will be shown, regardless of if this is set or not.
"""
function generate(repo::DataRepo, dataname, shortname=nothing; kwargs...)
    metadata = find_metadata(repo, dataname, shortname)
    format_codeblock(metadata)
end

function generate(dataname, shortname=nothing; kwargs...)
    all_repos = [T() for T in leaf_subtypes(DataRepo)]
    generate(all_repos, dataname, shortname; kwargs...)
end

function generate(repos::Vector, dataname, shortname = nothing; show_failures=false)
    retrieved_metadatas_ch = Channel{Any}(128)
    failures_ch = Channel{Tuple{DataRepo, Exception}}(128)

    # Get all the metadata we can
    @sync for repo in repos
        @async try
            metadata = find_metadata(repo, dataname, shortname)
            push!(retrieved_metadatas_ch, metadata)
        catch err
            push!(failures_ch, (repo, err))
        end
    end
    close(retrieved_metadatas_ch)
    close(failures_ch)

    retrieved_metadatas = collect(retrieved_metadatas_ch)
    failures = collect(failures_ch)
    # ============ Handle errors ===========

    # Display errors
    if length(retrieved_metadatas) == 0 || show_failures
        for (repo, err) in failures
            @warn("$repo failed", exception=err)
        end
    end

    # ============= Produce final result ===========
    merged_metadata = aggregate(retrieved_metadatas)
    format_codeblock(merged_metadata)
end

function find_metadata(repo, dataname, shortname)
    mainpage, url = mainpage_url(repo, dataname)
    fullname = data_fullname(repo, mainpage)
    shortname = data_shortname(repo, shortname, fullname)
    Metadata(
        shortname,
        fullname,
        website(repo, url, mainpage),
        description(repo, mainpage),
        author(repo, mainpage),
        maintainer(repo, mainpage),
        license(repo, mainpage),
        published_date(repo, mainpage),
        create_date(repo, mainpage),
        modified_date(repo, mainpage),
        paper_cite(repo, mainpage),
        dataset_cite(repo, mainpage),
        get_urls(repo, mainpage),
        get_checksums(repo, mainpage)
    )
end

function data_shortname(repo, shortname::Nothing, fullname)
    # Remove any characters not allowed in a file path
    reduce((s,r)->replace(s, r => ""), ['\\', '/', ':', '*', '?', '<', '>', '|']; init=fullname)
end

data_shortname(repo, shortname, fullname) = shortname

website(::DataRepo, url, mainpage) = url
description(::DataRepo, mainpage) = missing
author(::DataRepo, mainpage) = missing
maintainer(::DataRepo, mainpage) = missing
license(::DataRepo, mainpage) = missing
published_date(::DataRepo, mainpage) = missing
create_date(::DataRepo, mainpage) = missing
modified_date(::DataRepo, mainpage) = missing
paper_cite(::DataRepo, mainpage) = missing
dataset_cite(::DataRepo, mainpage) = missing
get_checksums(::DataRepo, mainpage) = missing
get_urls(::DataRepo, mainpage) = missing

function mainpage_url(repo::DataRepo, dataname)
    if startswith(dataname, "http")
        url = dataname
        #TODO: This isn't going to take https/http differences.
        dataname = first(split(replace(url, base_url(repo) => ""), "/"))
    else # not a URL
        url = joinpath(base_url(repo), dataname)
    end
    getpage(url), url
end


text:
using Test
using JustSayIt
using JustSayIt.API
using PyCall
import JustSayIt: DEFAULT_MODEL_NAME, TYPE_MODEL_NAME, MODELDIR_PREFIX, DEFAULT_NOISES
import JustSayIt: init_jsi, finalize_jsi

commands = Dict("help"      => Help.help,
                "type"      => Keyboard.type)
modeldirs = Dict(DEFAULT_MODEL_NAME => joinpath(MODELDIR_PREFIX, "vosk-model-small-en-us-0.15"),
                 TYPE_MODEL_NAME    => joinpath(MODELDIR_PREFIX, "vosk-model-small-en-us-0.15"))

init_jsi(commands, modeldirs, DEFAULT_NOISES)

@testset "$(basename(@__FILE__))" begin
    @testset "1. finalizing" begin
        finalize_jsi()
        @test true # meaning reached that point.
    end;
end;


text:
module SyntheticDatasetsDashboard
using Dash
using DashCoreComponents
using DashHtmlComponents
using SyntheticDatasets
using StatsPlots

app = dash( suppress_callback_exceptions = true );

app.title = "Synthetic Datasets Dashboard";

app.index_string = """<!DOCTYPE html>

        {%metas%}
        <title>{%title%}</title>
        {%css%}

        {%app_entry%}

            {%config%}
            {%scripts%}
            {%renderer%}

"""

include("functions.jl")
include("index.jl")

include("../layouts/blobs.jl")
include("../layouts/moons.jl")
include("../layouts/gaussian.jl")
include("../layouts/two_spirals.jl")
include("../layouts/circles.jl")

end # module


text:
include("autodiff.jl")
include("utils.jl")
include("numbers.jl")
include("groups/group_utils.jl")
include("sized_abstract_array.jl")
include("orthonormal_bases.jl")

# starting with tests of simple manifolds
include("circle.jl")
include("cholesky_space.jl")
include("euclidean.jl")
include("fixed_rank.jl")
include("grassmann.jl")
include("hyperbolic.jl")
include("sphere.jl")
include("stiefel.jl")
include("symmetric.jl")
include("rotations.jl")
include("symmetric_positive_definite.jl")
include("torus.jl")

#meta manifolds
include("product_manifold.jl")
include("power_manifold.jl")
include("vector_bundle.jl")
include("graph.jl")

include("metric.jl")
include("statistics.jl")

# Lie groups and actions
include("groups/groups_general.jl")
include("groups/array_manifold.jl")
include("groups/circle_group.jl")
include("groups/translation_group.jl")
include("groups/special_orthogonal.jl")
include("groups/product_group.jl")
include("groups/semidirect_product_group.jl")
include("groups/special_euclidean.jl")
include("groups/group_operation_action.jl")
include("groups/rotation_action.jl")
include("groups/translation_action.jl")


text:
module PopupProxies

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, template

using Stipple

export popup_proxy, PopupProxy

function __init__()
  Genie.Renderer.Html.register_normal_element("q__popup__proxy", context = Genie.Renderer.Html)
end

"""
    popup_proxy()

Renders a popup.
"""
function popup_proxy(
                      fieldname::Union{Symbol,Nothing} = nothing,
                      args...;
                      content::Union{String,Vector,Function} = "",
                      kwargs...)

  Genie.Renderer.Html.q__popup__proxy([isa(content, Function) ? content() : join(content)],
                    args...;
                    attributes(
                              [
                                :fieldname => fieldname,
                                kwargs...
                              ], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
end

popup_proxy(content::Union{Vector,Function},
            fieldname::Union{Symbol,Nothing} = nothing,
            args...;
            kwargs...) = popup_proxy(fieldname, args...; content, kwargs...)

mutable struct PopupProxy
  fieldname
  args
  content
  kwargs

  PopupProxy( fieldname::Union{Symbol,Nothing} = nothing,
              args...;
              content::Union{String,Vector,Function} = "",
              kwargs...) = new(fieldname, args, content, kwargs)
end

Base.string(pp::PopupProxy) = popup_proxy(pp.fieldname, pp.args...; pp.content, pp.kwargs...)

end


text:
"""
PetroPhysics(args...)

Calculate basic petrophysics properties according the WellLog
It is requiered that the WellLog be a DataFrame type with the columns of the new properties already created
They can be created with the function ``AddColPetro``

The next table show the list of variables allowed:

|PropertyName|Args|Default|Input|Description
|---|---|---|---|---|
|Log|Mandatory|--|DataFrame|Well Logs dataframe file|
|DepthFrom|Mandatory| -- | Number |Set the Depth upper limit to estimate properties|
|DepthTo|Mandatory| -- | Number |  Set the Depth lower limit to estimate properties|
|Vsh|Optional|false|Vsh=[GrSand, GrShale]|Calculates the Vsh of the given interval|
|Phie|Optional|false|Phie=true|Calculates the Phie of the given interval.It is requiered Vsh, DenPhi, and Neutron Logs|
|Sw|Optional|false|Sw=[a,m,n,Rw]|Calculates the Sw of the given interval with archie parameters. It is requiered Phie and DeepRes logs|
|Perm|Optional|false|Perm=[Fluid, Author]|Calculates the Permeability of the given interval with using Author and Fluids coeficients|
|PayFlag|Optional|false|PayFlag=[VshCutoff,PhieCutOff,SwCutoff,KCutOff]|Calculates the Pay of the given interval using the CutOff specified|
|Kh|Optional|false|Kh=true|Calculates the Flow Capacity and Normalized Cumulative Flow Capacity for picking perfs of the given interval|
"""
function PetroPhysics(Log,DepthFrom,DepthTo;    #range of depth to calculate Petrophysics
                     Vsh=false,                #If caculate Vshale   [GrSand, GrShale]
                     DenPhi=false,              #If caculate Porosity from Density Log  [RhoMatrix, RhoFluid]
                     Phie=false,               #If caculate Efective Porosity. It is requiered Vsh, DenPhi, and Neutron
                     Sw=false,                 #If caculate Water Saturation. It is requiered Phie and DeepRes [a,m,n,Rw]
                     Perm=false,               #If caculate Permeability).    Phie, Sw, Fluid, Author [Fluid, Author]
                     PayFlag=false,          #If calculate PayFlag. It is requiered Vsh,Phie,Sw,Perm [VshCutoff,PhieCutOff,SwCutoff,KCutOff]
                     Kh=false,               #If requiered Flow Capacity percentage
                     DepthType= :Md)

#Get the index of the range depth
frr=findall(x->x==DepthFrom,Log[!,DepthType])
too=findall(x->x==DepthTo,Log[!,DepthType])

# convert the index from Array{Int64} to Int64
fr=frr[1]
to=too[1]

ColumnList=names(Log)

####################################################
# VSHALE
####################################################

if Vsh != false
  #Vsh[1]= GrSand
  #Vsh[2]= GrShale
    if (:Vsh in ColumnList) == false
      Log.Vsh=zeros(size(Log,1)).+1
    end

        for i=fr:to
         Log.Vsh[i]=Vshale(Log[!,Vsh[1]][i],Vsh[2],Vsh[3])
         if Log.Vsh[i]>1
            Log.Vsh[i]=1
        end
        end

end

####################################################
# DenPhi
####################################################

    if DenPhi != false
    #DenPhi[2]= RhoMatrix
    #DenPhi[3]= RhoFluid
    if (:DenPhi in ColumnList) == false

       Log.DenPhi=zeros(size(Log,1))
    end

    for i=fr:to
           Log.DenPhi[i]= RhoPorosity(Log[!,DenPhi[1]][i],DenPhi[2],DenPhi[3])
        if Log.DenPhi[i]<0
            Log.DenPhi[i]=0
        end
        end

    end
####################################################
# Phie
####################################################
    if Phie != false
        if (:Phie in ColumnList) == false
          Log.Phie=zeros(size(Log,1))
        end

        if size(Phie,1)==3
            for i=fr:to
                Log.Phie[i]=Phiefec(Log[!,Phie[1]][i],Log[!,Phie[2]][i],Log[!,Phie[3]][i])
            end
        elseif size(Phie,1)==2
            for i=fr:to
                Log.Phie[i]=Log[!,Phie[1]][i].*(1 .-Log[!,Phie[2]][i] )
            end

        end
    end
####################################################
# Sw
####################################################
    if Sw != false
     #Sw[1]=a
    #Sw[2]=m
    #Sw[3]=n
    #Sw[4]=Rw
    #Sw[5]=Phie
    #Sw[6]=DeepRes
    #Sw[7]=Shale Resistivity
    #Sw[8]=Vshale
    #Sw[9]=Method [:Archie, :Smdx, :Indo]
    if (:Sw in ColumnList) == false
      Log.Sw=zeros(size(Log,1)).+2
    end

      for i=fr:to
            try
                Log.Sw[i]=SwFunction(Sw[1],Sw[2],Sw[3],Sw[4],Log[!,Sw[5]][i],Log[!,Sw[6]][i],
                             Rsh=Sw[7],Vsh=Log[!,Sw[8]][i],α=Sw[10],Method=Sw[9])
            catch e
                Log.Sw[i]=1
            end
          if Log.Sw[i]>1
            Log.Sw[i]=1
            end
      end

    end

####################################################
# Perm
####################################################
  if Perm != false
      #Perm[1]=phi
      #Perm[2]=Sw
      #Perm[3]=Fluid "Oil" or "Gas"
      #Perm[4]=Author "Timur" or "Morris"
  if (:Perm in ColumnList) == false
    Log.Perm=zeros(size(Log,1))
  end

    for i=fr:to
     Log.Perm[i]=PermWR(Log[!,Perm[1]][i],Log[!,Sw[2]][i],Fluid=Perm[3],Author=Perm[4])
    end

  end

 ####################################################
# PayFlag
####################################################
  if PayFlag != false

 #PayFlag[1]=VshCutoff
 #PayFlag[2]=PhieCutOff
 #PayFlag[3]=SwCutoff
 #PayFlag[4]=KCutOff
 if (:PayFlag in ColumnList) == false
   Log.PayFlag=zeros(size(Log,1))
 end

    for i=fr:to
     Log.PayFlag[i]=(Log[!,PayFlag[1]][i].<=PayFlag[2]).*
            (Log[!,PayFlag[3]][i].>=PayFlag[4]).*
            (Log[!,PayFlag[5]][i].<=PayFlag[6]).*
            (Log[!,PayFlag[7]][i].>=PayFlag[8]).*1
    end

  end

 ####################################################
# Flow Capacity Kh
####################################################
  if Kh != false
      #md
      #Perm
      #Payflag
      if (:Kh in ColumnList) == false
        Log.Kh=zeros(size(Log,1))
      end

      if (:KhCum in ColumnList) == false
        Log.KhCum=zeros(size(Log,1))
      end

     KH=zeros(to-fr+2)
     KHCum=zeros(to-fr+2)
    for i=fr:to
        j=i-fr+2
        KH[j]=Log[!,Kh[2]][i].*(Log[!,Kh[1]][i].-Log[!,Kh[1]][i-1]).*Log[!,Kh[3]][i]
        Log.Kh[i]=KH[j]
        KHCum[j]=KHCum[j-1]+KH[j]
    end
     KhTotal=sum(KH)
    if KhTotal > 0
        for i=fr:to
            j=i-fr+2
            Log.KhCum[i]=1-(KHCum[j]/KhTotal)
        end
    end

  end

return Log

end


text:
#
# CartesianBoxes.jl -
#
# Extends CartesianIndices.
#
#-------------------------------------------------------------------------------
#
# This file is part of the `CartesianBoxes.jl` package which is licensed under
# the MIT "Expat" License.
#
# Copyright (c) 2017-2022 Éric Thiébaut.
#

__precompile__(true)

"""
`CartesianBoxes` extends `CartesianIndices`.
"""
module CartesianBoxes

export
    CartesianBox,
    boundingbox,
    isnonemptypartof,
    isnonzero

using Base: tail, @propagate_inbounds

import Base:
    intersect, issubset,
    simd_outer_range, simd_inner_length, simd_index

# Starting with Julia-1.6.0-beta1, Cartesian indices can have non-unit step.
const IndexRange{I} = @static if VERSION < v"1.6.0-beta1"
    AbstractUnitRange{I}
else
    OrdinalRange{I,I}
end

# Method to convert an argument (of CartesianBox constructor) to an index
# range.
to_index_range(x::IndexRange{Int}) = x
to_index_range(x::IndexRange{<:Integer}) = convert(IndexRange{Int}, x)
to_index_range(x::Integer) = Base.OneTo{Int}(x)

"""

`CartesianBox{N}` defines a rectangular region of `N`-dimensional indices and
can be constructed by:

    CartesianBox(A)
    CartesianBox(axes(A))
    CartesianBox(CartesianIndex(istart,jstart,...), CartesianIndex(istop,jstop,...))
    CartesianBox((istart,jstart,...), (istop,jstop,...))
    CartesianBox((istart:[istep:]istop, jstart:[jstep:]jstop, ...))
    CartesianBox(istart:[istep:]istop, jstart:[jstep:]jstop, ...)

where `A` is an array (to define a region consisting in all the indices in the
array), `istart`, `istop`, etc. are integers (to define a region from
`(istart,jstart,...)` to `(istop,jstop,...)`.  Starting with Julia 1.6, a
non-unit step `istep`, `jstep`, etc. may also be specified.

`CartesianBox(size(A))` yields a Cartesian region of the same size as `A` but
whose indices all start at 1.  In most cases, `CartesianBox(axes(A))` is more
likely to do the right thing.  In fact, `CartesianBox(A)` is equivalent to
`CartesianBox(axes(A))`.

An instance of `CartesianBox` can also be constructed from an instance of
`CartesianIndices` and conversely.  The conversion is lossless in the sense of
the following example:

    B = CartesianBox(...)
    R = CartesianIndices(B)
    CartesianBox(R) === B # always true

An instance of `CartesianBox` can be efficiently used in a loop as follows:

    for i in CartesianBox(...)
       A[i] = ...
    end

where `i` will be set to a `CartesianIndex` with all the multi-dimensional
indices of the rectangular region defined by `B`.

When at least one of `A` or `B` is a Cartesian box, the expression `A ∩ B`, or
equivalently `intersect(A,B)`, yields the Cartesian box contining all indices
in `A` and `B`.  This may be used to write safe loops like:

    A = ...               # some array
    B = CartesianBox(...) # some region of interest
    @inbounds for i in B ∩ A
        A[i] = ...
    end

to operate on the indices of the Cartesian box `B` that are valid for `A`.

When at least one of `A` or `B` is a Cartesian box, the expression `A ⊆ B`, or
equivalently `intersect(A,B)`, yields whether all Cartesian indices defined by
`A` are also indices of `B`.  This may be used as:

    A = ...               # some array
    B = CartesianBox(...) # some region of interest
    if B ⊆ A
        @inbounds for i in B
            A[i] = ...
        end
    end

to only access the indices of the Cartesian box `B` if they are all valid for
`A`.

See also: [`boundingbox`](@ref), `CartesianIndices`, `CartesianIndex`,
[`intersect`](@ref).

"""
struct CartesianBox{N,I<:NTuple{N,IndexRange{Int}}} <: AbstractArray{CartesianIndex{N},N}
    indices::I
end

"""
    indices(A) -> inds

yields an `N`-tuple of index ranges if `A` is an instance of
[`CartesianBox{N}`](@ref) or of `CartesianIndices{N}`; yields an `N`-tuple of
indices if `A` is an instance of `CartesianIndex{N}` or an `N`-tuple of
integers.

"""
indices(R::CartesianBox) = getfield(R, :indices)
indices(R::CartesianIndices) = getfield(R, :indices)
indices(I::CartesianIndex) = Tuple(I)
indices(I::NTuple{N,Integer}) where {N} = I

Base.Tuple(R::CartesianBox) = indices(R)

# Fast conversion between CartesianIndices and CartesianBox.
Base.CartesianIndices(R::CartesianBox{N,I}) where {N,I} =
    CartesianIndices{N,I}(indices(R))
CartesianBox(R::CartesianIndices{N,I}) where {N,I} =
    CartesianBox{N,I}(indices(R))

# Other constructors.
CartesianBox(B::CartesianBox) = B
CartesianBox(A::AbstractArray) = CartesianBox(axes(A))
CartesianBox(inds::Vararg{Union{Integer,IndexRange{<:Integer}}}) =
        CartesianBox(inds)
CartesianBox(first::CartesianIndex{N}, last::CartesianIndex{N}) where {N} =
    CartesianBox(indices(first), indices(last))
CartesianBox(first::NTuple{N,Integer}, last::NTuple{N,Integer}) where {N} =
    CartesianBox(map((i,j) -> i:j, first, last))
CartesianBox(inds::Tuple{Vararg{Union{Integer,IndexRange{<:Integer}}}}) =
        CartesianBox(map(to_index_range, inds))

@deprecate ranges(B::CartesianBox) indices(B) false

Base.first(B::CartesianBox) = CartesianIndex(map(first, indices(B)))
Base.last(B::CartesianBox) = CartesianIndex(map(last, indices(B)))
Base.size(B::CartesianBox) = map(length, indices(B))
Base.axes(B::CartesianBox) = map(Base.axes1, indices(B))

Base.show(io::IO, ::MIME"text/plain", B::CartesianBox) = show(io, B)
Base.show(io::IO, B::CartesianBox) = begin
    print(io, "CartesianBox(")
    print(io, indices(B))
    print(io, ")")
end

Base.view(A::AbstractArray{<:Any,N}, B::CartesianBox{N}) where {N} =
    view(A, indices(B)...)

function Base.getindex(A::AbstractArray{T,N}, B::CartesianBox{N}) where {T,N}
    empty = isempty(B)
    empty || isnonemptypartof(B, A) ||
        throw(BoundsError("sub-region is not part of array"))
    C = Array{T}(undef, size(B))
    if ! empty
        Ifirst = first(B)
        if any(i -> i != 1, indices(Ifirst))
            off = CartesianIndex(map(i -> i - 1, indices(Ifirst)))
            @inbounds @simd for i in B
                C[i - off] = A[i]
            end
        else
            @inbounds @simd for i in B
                C[i] = A[i]
            end
        end
    end
    return C
end

Base.setindex!(A::AbstractArray{<:Any,N}, x, B::CartesianBox{N}) where {N} =
    A[indices(B)...] = x

Base.fill!(A::AbstractArray{T,N}, B::CartesianBox{N}, x) where {T,N} =
    fill!(A, B, convert(T, x)::T)
function Base.fill!(A::AbstractArray{T,N},
                    B::CartesianBox{N}, x::T) where {T,N}
    if ! isempty(B)
        isnonemptypartof(B, A) ||
            throw(BoundsError("sub-region is not part of array"))
        @inbounds @simd for i in B
            A[i] = x
        end
    end
    return A
end

Base.IndexStyle(::Type{<:CartesianBox{N,R}}) where {N,R} =
    IndexStyle(CartesianIndices{N,R})
@inline @propagate_inbounds Base.getindex(B::CartesianBox, I...) =
    getindex(CartesianIndices(B), I...)

"""
    CartesianBoxable{N}

is a union of types (other than `CartesianBox{N}`) which can be automatically
converted into a `CartesianBox{N}`.

!!! note
    Although the constructor `CartesianBox` can be also applied to any instance
    of `AbstractArray`, this abstract type does not belong to the union
    `CartesianBoxable` as it is considered that such a conversion cannot be
    automatic.

"""
const CartesianBoxable{N} = Union{CartesianIndices{N},
                                  NTuple{N,IndexRange{<:Integer}},
                                  NTuple{N,Integer}}

# Extend eachindex() method.
Base.eachindex(::IndexCartesian, B::CartesianBox) = B

# Make CartesianBox iterable.
Base.iterate(iter::CartesianBox) = iterate(CartesianIndices(iter))
Base.iterate(iter::CartesianBox, state) = iterate(CartesianIndices(iter), state)

# Extend methods for fast SIMD iterations.
simd_outer_range(iter::CartesianBox{0}) = iter
simd_outer_range(iter::CartesianBox) = CartesianBox(tail(indices(iter)))
simd_inner_length(iter::CartesianBox{0}, ::CartesianIndex) = 1
simd_inner_length(iter::CartesianBox, I::CartesianIndex) = length(indices(iter)[1])
simd_index(iter::CartesianBox{0}, ::CartesianIndex, I1::Int) = first(iter)
@inline @propagate_inbounds function simd_index(iter::CartesianBox,
                                                Ilast::CartesianIndex, I1::Int)
    # For maximum portability, delegate work to do to the embedded
    # CartesianIndices.
    simd_index(CartesianIndices(iter), Ilast, I1)
end

# Increment a range.
incr(rng::AbstractUnitRange, adj::Number) =
    (first(rng) + adj):(last(rng) + adj)
incr(rng::AbstractRange, adj::Number) =
    (first(rng) + adj):step(rng):(last(rng) + adj)

# Decrement a range.
decr(rng::AbstractUnitRange, adj::Number) =
    (first(rng) - adj):(last(rng) - adj)
decr(rng::AbstractRange, adj::Number) =
    (first(rng) - adj):step(rng):(last(rng) - adj)

# Shifting of a CartesianBox.
function Base.:(+)(R::CartesianBox{N},
                   I::Union{NTuple{N,Integer},CartesianIndex{N}}) where {N}
    CartesianBox(map((r,i) -> incr(r, i), indices(R), indices(I)))
end
function Base.:(-)(R::CartesianBox{N},
                   I::Union{NTuple{N,Integer},CartesianIndex{N}}) where {N}
    CartesianBox(map((r,i) -> decr(r, i), indices(R), indices(I)))
end

Base.isempty(B::CartesianBox) = prod(map(isempty, indices(B)))

# Extend ∩ operator (\cap-tab) for Cartesian boxes.
Base.intersect(A::CartesianBox{N}, B::CartesianBox{N}) where {N} =
    CartesianBox(map(intersect, indices(A), indices(B)))

# Extend ⊆ operator (\subseteq-tab) for Cartesian boxes.
Base.issubset(A::CartesianBox{N}, B::CartesianBox{N}) where {N} =
    prod(map(issubset, indices(A), indices(B)))

# Override ∩(A,B) and ⊆(A,B) when at least one of A or B is a Cartesian box,
# the other being boxable or an abstract array.
for f in (:intersect, :issubset)
    @eval begin
        $f(A::CartesianBox{N}, B::CartesianBoxable{N}) where {N} =
            $f(A, CartesianBox(B))
        $f(A::CartesianBox{N}, B::AbstractArray{<:Any,N}) where {N} =
            $f(A, CartesianBox(B))
        $f(A::CartesianBoxable{N}, B::CartesianBox{N}) where {N} =
            $f(CartesianBox(A), B)
        $f(A::AbstractArray{<:Any,N}, B::CartesianBox{N}) where {N} =
            $f(CartesianBox(A), B)
    end
end

"""
    intersect(CartesianBox, A, B)

yields the Cartesian box given by the intersection of the Cartesian regions
defined by `A` and `B`.  In this context, a Cartesian region can specified by a
Cartesian box, a list of integer valued ranges, a list of dimensions, or an
instance of `CartesianIndices`.

This method is equivalent to `intersect(A,B)`, or `A ∩ B` for short, when at
least one of `A` or `B` is a Cartesian box.

See also: [`CartesianBox`](@ref), [`isnonemptypartof`](@ref).

""" intersect

@deprecate intersection(A, B) intersect(CartesianBox, A, B) false

intersect(::Type{CartesianBox}, A::CartesianBox{N}, B::CartesianBox{N}) where {N} = A ∩ B
intersect(::Type{CartesianBox}, A::CartesianBox{N}, B::CartesianBoxable{N}) where {N} =
    A ∩ CartesianBox(B)
intersect(::Type{CartesianBox}, A::CartesianBoxable{N}, B::CartesianBox{N}) where {N} =
    CartesianBox(A) ∩ B
intersect(::Type{CartesianBox}, A::CartesianBoxable{N}, B::CartesianBoxable{N}) where {N} =
    CartesianBox(A) ∩ CartesianBox(B)

"""
    isnonemptypartof(A, B)

yields whether the region defined by `A` is nonempty and a valid part of the
region defined by `B` or of the contents of `B` if `B` is an array.  This is
equivalent to:

    !isempty(CartesianBox(A)) && (CartesianBox(A) ⊆ CartesianBox(B))

except that `A` may not be an array.

See also: [`CartesianBox`](@ref), [`issubset`](@ref).

"""
function isnonemptypartof(A::Union{CartesianBoxable{N},CartesianBox{N}},
                          B::Union{CartesianBoxable{N},
                                   AbstractArray{<:Any,N}}) where {N}
    A_box = CartesianBox(A)
    if isempty(A_box)
        return false
    else
        return A_box ⊆ CartesianBox(B)
    end
end

isnonemptypartof(A::CartesianBox{N}, B::CartesianBox{N}) where {N} =
    !isemppty(A) && issubset(A, B)

"""
    boundingbox([pred,] A [, B])

yields the bounding-box of the elements in array `A` for which the predicate
function `pred` is true.  If the predicate function `pred` is omitted, the
result is the bounding-box of non-zero values in array `A` or of the `true`
values in `A` if its elements are of type `Bool`.  Optional argument `B` is to
only consider a sub-region `B` of `A` (`B` can be a `CartesianBox`, a
`CartesianIndices`, or a tuple of integer valued unit ranges).

!!! warning
    The algorithm is pretty silly for now and could be made faster than
    `O(length(A))`.

See also: [`CartesianBox`](@ref), [`intersect`](@ref), [`isnonzero`](@ref).

"""
boundingbox(A::AbstractArray{Bool}) = boundingbox(identity, A)
boundingbox(A::AbstractArray) = boundingbox(isnonzero, A)
boundingbox(A::AbstractArray{Bool,N}, B::CartesianBox{N}) where {N} =
    boundingbox(identity, A, B)
boundingbox(A::AbstractArray{<:Any,N}, B::CartesianBox{N}) where {N} =
    boundingbox(isnonzero, A, B)
boundingbox(A::AbstractArray{<:Any,N}, B::CartesianBoxable{N}) where {N} =
    boundingbox(A, CartesianBox(B))
boundingbox(pred, A::AbstractArray{<:Any,N}, B::CartesianBoxable{N}) where {N} =
    boundingbox(pred, A, CartesianBox(B))

function boundingbox(pred,
                     A::AbstractArray{T,N}) where {T,N}
    Imin = typemax_(CartesianIndex{N})
    Imax = typemin_(CartesianIndex{N})
    @inbounds for I in CartesianBox(A)
        if pred(A[I])
            Imin = min(Imin, I)
            Imax = max(Imax, I)
        end
    end
    if isempty_(Imin, Imax)
        Imin = one_(CartesianIndex{N})
        Imax = zero(CartesianIndex{N})
    end
    return CartesianBox(Imin, Imax)
end

function boundingbox(pred,
                     A::AbstractArray{T,N},
                     B::CartesianBox{N}) where {T,N}
    Imin = typemax_(CartesianIndex{N})
    Imax = typemin_(CartesianIndex{N})
    @inbounds for I in CartesianBox(A) ∩ B
        if pred(A[I])
            Imin = min(Imin, I)
            Imax = max(Imax, I)
        end
    end
    if isempty_(Imin, Imax)
        Imin = one_(CartesianIndex{N})
        Imax = zero(CartesianIndex{N})
    end
    return CartesianBox(Imin, Imax)
end

@inline isempty_(first::CartesianIndex{N}, last::CartesianIndex{N}) where {N} =
    any(map(isless, indices(last), indices(first)))

# Provide equivalent methods for `typemin`, `typemax`, and `one` applied
# to a Cartesian index with different names to avoid type piracy.
@inline typemin_(::Type{<:CartesianIndex{N}}) where {N} =
    CartesianIndex(ntuple(i -> typemin(Int), Val(N)))
@inline typemax_(::Type{<:CartesianIndex{N}}) where {N} =
    CartesianIndex(ntuple(i -> typemax(Int), Val(N)))

one_(I::CartesianIndex) = one_(typeof(I))
@static if VERSION < v"1.1.0-rc1"
    one_(T::Type{<:CartesianIndex}) = one(T)
else
    one_(T::Type{<:CartesianIndex}) = oneunit(T)
end

"""
    isnonzero(x)

yields whether `x` is non-zero, the negation of `iszero(x)`.

"""
isnonzero(x) = !iszero(x)

end # module


text:
include("../src/JsonCsvConverter.jl")
using Base.Test

c = JsonCsvConverter.Converter()

function _output_file(filename::String)
    Base.source_dir()*"/output/$filename"
end

function _input_file(filename::String)
    Base.source_dir()*"/input/$filename"
end

function _expected_file(filename::String)
    Base.source_dir()*"/expected/$filename"
end

function set_up()
    files = [
        "out.csv",
        "out.json",
    ]

    for filename in files
        filepath = _output_file(filename)
        if isfile(filepath)
            rm(filepath)
        end
    end
end

function _test(
    input_filename::String,
    output_filename::String,
    expected_output_filename::String,
    parse_nested_json::Bool=false,
)
    set_up()
    input_filepath = _input_file(input_filename)
    output_filepath = _output_file(output_filename)
    expected_filepath = _expected_file(expected_output_filename)
    c.convert(input_filepath, output_filepath, parse_nested_json=parse_nested_json)
    @test readstring(output_filepath) == readstring(expected_filepath)
end

function _test_double(
    input_filename::String,
    first_output_filename::String,
    second_output_filename::String,
    parse_nested_json::Bool=false,
)
    set_up()
    input_filepath = _input_file(input_filename)
    first_output_filepath = _output_file(first_output_filename)
    second_output_filepath = _output_file(second_output_filename)

    c.convert(input_filepath, first_output_filepath, parse_nested_json=parse_nested_json)
    c.convert(first_output_filepath, second_output_filepath, parse_nested_json=parse_nested_json)

    @test readstring(second_output_filepath) == readstring(input_filepath)
end

function test_flat_json_to_csv()
    _test("flat.json", "out.csv", "flat_json_to_csv.csv")
end

function test_flat_quoted_csv_to_json()
    _test("flat_quoted.csv", "out.json", "flat_csv_to_json.json")
end

function test_flat_unquoted_csv_to_json()
    _test("flat_unquoted.csv", "out.json", "flat_csv_to_json.json")
end

function test_flat_json_to_csv_and_back()
    _test_double("flat.json", "out.csv", "out.json")
end

function test_flat_csv_to_json_and_back()
    _test_double("flat_quoted.csv", "out.json", "out.csv")
end

function test_nested_json_to_csv()
    _test("nested.json", "out.csv", "nested_json_to_csv.csv")
end

function test_nested_csv_to_nested_json()
    _test("nested.csv", "out.json", "nested_csv_to_json.json", true)
end

function test_nested_csv_to_flat_json()
end

function test_nested_json_to_csv_and_back()
    _test_double("nested.json", "out.csv", "out.json", true)
end

function test_nested_csv_to_json_and_back()
    _test_double("nested.csv", "out.json", "out.csv", true)
end

function test_na_handling()
    _test("flat_unquoted_nas.csv", "out.json", "na_handling.json")
end

test_flat_json_to_csv()
test_flat_quoted_csv_to_json()
test_flat_unquoted_csv_to_json()
test_flat_json_to_csv_and_back()
test_flat_csv_to_json_and_back()

test_nested_json_to_csv()
test_nested_csv_to_nested_json()
test_nested_csv_to_flat_json()
test_nested_json_to_csv_and_back()
test_nested_csv_to_json_and_back()

test_na_handling()


text:
"""
    Logistic(μ,θ)

The *Logistic distribution* with location `μ` and scale `θ` has probability density function

```math
f(x; \\mu, \\theta) = \\frac{1}{4 \\theta} \\mathrm{sech}^2
\\left( \\frac{x - \\mu}{2 \\theta} \\right)
```

```julia
Logistic()       # Logistic distribution with zero location and unit scale, i.e. Logistic(0, 1)
Logistic(u)      # Logistic distribution with location u and unit scale, i.e. Logistic(u, 1)
Logistic(u, b)   # Logistic distribution with location u ans scale b

params(d)       # Get the parameters, i.e. (u, b)
location(d)     # Get the location parameter, i.e. u
scale(d)        # Get the scale parameter, i.e. b
```

External links

* [Logistic distribution on Wikipedia](http://en.wikipedia.org/wiki/Logistic_distribution)

"""
struct Logistic{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    θ::T
    Logistic{T}(µ::T, θ::T) where {T} = new{T}(µ, θ)
end

function Logistic(μ::T, θ::T; check_args=true) where {T <: Real}
    check_args && @check_args(Logistic, θ > zero(θ))
    return Logistic{T}(μ, θ)
end

Logistic(μ::Real, θ::Real) = Logistic(promote(μ, θ)...)
Logistic(μ::Integer, θ::Integer) = Logistic(float(μ), float(θ))
Logistic(μ::T) where {T <: Real} = Logistic(μ, one(T))
Logistic() = Logistic(0.0, 1.0, check_args=false)

@distr_support Logistic -Inf Inf

#### Conversions
function convert(::Type{Logistic{T}}, μ::S, θ::S) where {T <: Real, S <: Real}
    Logistic(T(μ), T(θ))
end
function convert(::Type{Logistic{T}}, d::Logistic{S}) where {T <: Real, S <: Real}
    Logistic(T(d.μ), T(d.θ), check_args=false)
end

#### Parameters

location(d::Logistic) = d.μ
scale(d::Logistic) = d.θ

params(d::Logistic) = (d.μ, d.θ)
@inline partype(d::Logistic{T}) where {T<:Real} = T

#### Statistics

mean(d::Logistic) = d.μ
median(d::Logistic) = d.μ
mode(d::Logistic) = d.μ

std(d::Logistic) = π * d.θ / sqrt3
var(d::Logistic) = (π * d.θ)^2 / 3
skewness(d::Logistic{T}) where {T<:Real} = zero(T)
kurtosis(d::Logistic{T}) where {T<:Real} = T(6)/5

entropy(d::Logistic) = log(d.θ) + 2

#### Evaluation

zval(d::Logistic, x::Real) = (x - d.μ) / d.θ
xval(d::Logistic, z::Real) = d.μ + z * d.θ

pdf(d::Logistic, x::Real) = (lz = logistic(-abs(zval(d, x))); lz*(1-lz)/d.θ)
logpdf(d::Logistic, x::Real) = (u = -abs(zval(d, x)); u - 2*log1pexp(u) - log(d.θ))

cdf(d::Logistic, x::Real) = logistic(zval(d, x))
ccdf(d::Logistic, x::Real) = logistic(-zval(d, x))
logcdf(d::Logistic, x::Real) = -log1pexp(-zval(d, x))
logccdf(d::Logistic, x::Real) = -log1pexp(zval(d, x))

quantile(d::Logistic, p::Real) = xval(d, logit(p))
cquantile(d::Logistic, p::Real) = xval(d, -logit(p))
invlogcdf(d::Logistic, lp::Real) = xval(d, -logexpm1(-lp))
invlogccdf(d::Logistic, lp::Real) = xval(d, logexpm1(-lp))

function gradlogpdf(d::Logistic, x::Real)
    e = exp(-zval(d, x))
    ((2e) / (1 + e) - 1) / d.θ
end

mgf(d::Logistic, t::Real) = exp(t * d.μ) / sinc(d.θ * t)

function cf(d::Logistic, t::Real)
    a = (π * t) * d.θ
    a == zero(a) ? complex(one(a)) : cis(t * d.μ) * (a / sinh(a))
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

module UnitCommitment

include("instance/structs.jl")
include("model/formulations/base/structs.jl")
include("solution/structs.jl")

include("model/formulations/ArrCon2000/structs.jl")
include("model/formulations/CarArr2006/structs.jl")
include("model/formulations/DamKucRajAta2016/structs.jl")
include("model/formulations/Gar1962/structs.jl")
include("model/formulations/KnuOstWat2018/structs.jl")
include("model/formulations/MorLatRam2013/structs.jl")
include("model/formulations/PanGua2016/structs.jl")
include("solution/methods/XavQiuWanThi2019/structs.jl")

include("import/egret.jl")
include("instance/read.jl")
include("model/build.jl")
include("model/formulations/ArrCon2000/ramp.jl")
include("model/formulations/base/bus.jl")
include("model/formulations/base/line.jl")
include("model/formulations/base/psload.jl")
include("model/formulations/base/sensitivity.jl")
include("model/formulations/base/system.jl")
include("model/formulations/base/unit.jl")
include("model/formulations/CarArr2006/pwlcosts.jl")
include("model/formulations/DamKucRajAta2016/ramp.jl")
include("model/formulations/Gar1962/pwlcosts.jl")
include("model/formulations/Gar1962/status.jl")
include("model/formulations/Gar1962/prod.jl")
include("model/formulations/KnuOstWat2018/pwlcosts.jl")
include("model/formulations/MorLatRam2013/ramp.jl")
include("model/formulations/MorLatRam2013/scosts.jl")
include("model/formulations/PanGua2016/ramp.jl")
include("model/jumpext.jl")
include("solution/fix.jl")
include("solution/methods/XavQiuWanThi2019/enforce.jl")
include("solution/methods/XavQiuWanThi2019/filter.jl")
include("solution/methods/XavQiuWanThi2019/find.jl")
include("solution/methods/XavQiuWanThi2019/optimize.jl")
include("solution/optimize.jl")
include("solution/solution.jl")
include("solution/warmstart.jl")
include("solution/write.jl")
include("transform/initcond.jl")
include("transform/slice.jl")
include("transform/randomize/XavQiuAhm2021.jl")
include("utils/log.jl")
include("utils/benchmark.jl")
include("validation/repair.jl")
include("validation/validate.jl")

end


text:
using OrdinaryDiffEq, Test, DiffEqDevTools
using LinearAlgebra, Random

# Linear exponential solvers
A = DiffEqArrayOperator([2.0 -1.0; -1.0 2.0])
u0 = ones(2)
prob = ODEProblem(A,u0,(0.0,1.0))
sol1 = solve(prob, LinearExponential(krylov=:off))(1.0)
sol2 = solve(prob, LinearExponential(krylov=:simple))(1.0)
sol3 = solve(prob, LinearExponential(krylov=:adaptive))(1.0)
sol_analytic = exp(1.0 * Matrix(A)) * u0

@test isapprox(sol1, sol_analytic, rtol=1e-10)
@test isapprox(sol2, sol_analytic, rtol=1e-10)
@test isapprox(sol3, sol_analytic, rtol=1e-10)

# u' = A(t)u solvers

function update_func(A,u,p,t)
    A[1,1] = cos(t)
    A[2,1] = sin(t)
    A[1,2] = cos(t)*sin(t)
    A[2,2] = sin(t)^2
end
A = DiffEqArrayOperator(ones(2,2),update_func=update_func)
prob = ODEProblem(A, ones(2), (0.0, 5.))
dts = 1 ./2 .^(10:-1:1)
sol  = solve(prob,OrdinaryDiffEq.MagnusMidpoint(),dt=1/4)

dts = 1 ./2 .^(10:-1:1)
test_setup = Dict(:alg=>Vern9(),:reltol=>1e-14,:abstol=>1e-14)
sim = analyticless_test_convergence(dts,prob,MagnusMidpoint(),test_setup)
@test sim.𝒪est[:l2] ≈ 2 atol=0.2
sim = analyticless_test_convergence(dts,prob,MagnusMidpoint(krylov=true),test_setup)
@test sim.𝒪est[:l2] ≈ 2 atol=0.2

A = DiffEqArrayOperator(ones(2,2),update_func=update_func)
prob = ODEProblem(A, ones(2), (0.0, 5.))
dts = 1 ./2 .^(10:-1:1)
sol  = solve(prob,OrdinaryDiffEq.MagnusLeapfrog(),dt=1/4)

dts = 1 ./2 .^(10:-1:1)
test_setup = Dict(:alg=>Vern9(),:reltol=>1e-14,:abstol=>1e-14)
sim = analyticless_test_convergence(dts,prob,MagnusLeapfrog(),test_setup)
@test sim.𝒪est[:l2] ≈ 2 atol=0.2
sim = analyticless_test_convergence(dts,prob,MagnusLeapfrog(krylov=true),test_setup)
@test sim.𝒪est[:l2] ≈ 2 atol=0.2

function update_func(A,u,p,t)
    A[1,1] = cos(t)*u[1]
    A[2,1] = sin(t)*u[2]
    A[1,2] = -cos(t)*sin(t)
    A[2,2] = -(sin(t)^2) *u[1] * u[2]
end
A = DiffEqArrayOperator(ones(2,2),update_func=update_func)
prob = ODEProblem(A, ones(2), (0.5, 5.))
dts = 1 ./2 .^(10:-1:1)
sol  = solve(prob,OrdinaryDiffEq.LieEuler(),dt=1/4)

dts = 1 ./2 .^(10:-1:1)
test_setup = Dict(:alg=>Vern9(),:reltol=>1e-14,:abstol=>1e-14)
sim = analyticless_test_convergence(dts,prob,LieEuler(),test_setup)
@test sim.𝒪est[:l2] ≈ 1 atol=0.2
sim = analyticless_test_convergence(dts,prob,LieEuler(krylov=true),test_setup)
@test sim.𝒪est[:l2] ≈ 1 atol=0.2

function B(y::AbstractMatrix)
    b = similar(y)
    N = size(b, 1)
    for l=1:N
        for k=1:l
            if k < l
                @inbounds b[k,l] = y[k, l-1] - y[k+1, l]
            else
                @inbounds b[k,l] = 0
            end
        end
    end
    for l=1:N
        for k=l+1:N
            @inbounds b[k,l] = -b[l,k]
        end
    end

    return b
end

function update_func(A, u, p, t)
    A .= B(u)
    return nothing
end

η = diagm([1.,2,3,4,5])
A = DiffEqArrayOperator(Matrix{eltype(η)}(I(size(η,1))), update_func=update_func)
dts = 1 ./2 .^(10:-1:2)
tspan = (0., 20.)

# IIP
f = SplitFunction(A, (du,u,p,t)->du.=-u*B(u), _func_cache=similar(η))
prob = SplitODEProblem(f, η, tspan)
sol = solve(prob, CayleyEuler(), dt=1/10)

@test sol.retcode == :Success
eig_err = [norm(eigvals(sol[i])-eigvals(η)) for i in eachindex(sol)]
@test all(≈(e,0, atol=1e-13) for e in eig_err)

test_setup = Dict(:alg=>Vern9(),:reltol=>1e-14,:abstol=>1e-14)

sim = analyticless_test_convergence(dts,prob,CayleyEuler(),test_setup)
@test sim.𝒪est[:l2] ≈ 1 atol=0.2

# OOP
f = SplitFunction(A, (u,p,t)->-u*B(u), _func_cache=similar(η))
prob = SplitODEProblem(f, η, tspan)
sol = solve(prob, CayleyEuler(), dt=1/10)

@test sol.retcode == :Success
eig_err = [norm(eigvals(sol[i])-eigvals(η)) for i in eachindex(sol)]
@test all(≈(e,0, atol=1e-13) for e in eig_err)

test_setup = Dict(:alg=>Vern9(),:reltol=>1e-14,:abstol=>1e-14)

sim = analyticless_test_convergence(dts,prob,CayleyEuler(),test_setup)
@test sim.𝒪est[:l2] ≈ 1 atol=0.2


text:
using Dates
using NCDatasets
using OrderedCollections
using Test
using ClimateMachine.Writers

@testset "Writers" begin
    odims = OrderedDict(
        "x" => (collect(1:5), Dict()),
        "y" => (collect(1:5), Dict()),
        "z" => (collect(1010:10:1050), Dict()),
    )
    ovartypes = OrderedDict(
        "v1" => (("x", "y", "z"), Float64, Dict()),
        "v2" => (("x", "y", "z"), Float64, Dict()),
    )
    vals1 = rand(5, 5, 5)
    vals2 = rand(5, 5, 5)

    nc = NetCDFWriter()
    nfn, _ = mktemp()
    nfull = full_name(nc, nfn)

    init_data(nc, nfn, odims, ovartypes)
    append_data(nc, OrderedDict("v1" => vals1, "v2" => vals2), 2.0)

    NCDataset(nfull, "r") do nds
        xdim = nds["x"][:]
        ydim = nds["y"][:]
        zdim = nds["z"][:]
        @test xdim == odims["x"][1]
        @test ydim == odims["y"][1]
        @test zdim == odims["z"][1]
        @test try
            adim = nds["a"][:]
            adim == ones(5)
            false
        catch e
            true
        end
        t = nds["time"][:]
        v1 = nds["v1"][:]
        v2 = nds["v2"][:]
        @test length(t) == 1
        @test t[1] == DateTime(1900, 1, 1, 0, 0, 2)
        @test v1[:, :, :, 1] == vals1
        @test v2[:, :, :, 1] == vals2
    end
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using Pkg
Pkg.activate("./GraphBLASJuliaCvdlab/")

include("../../src/GraphBLAS_bfs_cvd.jl")
include("../../src/GraphBLAS_sparse_square_matrix_cvd.jl")

@testset "GraphBLASJuliaCvdlab" begin

### algorithm correctness test
    mat = [0 1 0 1 0 0 0; 0 0 0 0 1 0 1; 0 0 0 0 0 1 0; 1 0 1 0 0 0 0 ; 0 0 0 0 0 1 0; 0 0 1 0 0 0 0; 0 0 1 1 1 0 0]
    A = from_matrix(mat)

    sol = from_lists([1, 2, 3, 4, 5, 6, 7], [1, 2, 3, 2, 3, 4, 3])

    @test GraphBLAS_bfs_cvd(A, 1) == sol

### exception throw tests

    # Input matrix is not square
    B = from_matrix([1 0 1 0 ; 0 0 1 0 ; 0 0 0 1 ; 1 0 0 0 ; 0 0 1 0])
    @test_throws DimensionMismatch GraphBLAS_bfs_cvd(B, 1)

    # source is a not a positive Integer
    @test_throws ArgumentError GraphBLAS_bfs_cvd(A, -1)

    # source is not an Integer
    @test_throws MethodError GraphBLAS_bfs_cvd(A, 1.2)

    # Input matrix is not a GBMatrix
    @test_throws MethodError GraphBLAS_bfs_cvd(mat, 1)

end


text:
function allornothing(network::AbstractNetwork,r,trips::Vector,costs::Vector)
    #nlinks = numlinks(network)
    zoneids = id.(zones(network))
    flows = zero(costs)

    _, parentvec = dijkstra(network, r, costs)
    for j in zoneids
        t = trips[j]
        if isapprox(t, 0.) | (r == j)
            continue
        end

        # Finding the shortest path
        p = pathto(parentvec, j)

        # Assigning flows
        for l in p
            flows[idx(network, l)] += t
        end

    end
    return flows
end

function allornothing(network::AbstractNetwork,
                      trips::Matrix{U},
                      costs::Vector{T};
                      basedon::Symbol = :link) where {T<:Real, U<:Real}
    # Initialization
    nnodes, nzones = numnodes(network), numzones(network)
    nlinks = length(costs)
    zonelist = zones(network)
    zoneids = id.(zonelist)

    flows = (basedon == :link) ? zeros(nlinks) : (basedon == :origin) ? zeros(nlinks, nzones) : error("Invalid input!")

    for (i,z) in enumerate(zoneids)
        if basedon == :link
            flows += allornothing(network, z, trips[z,:], costs)
        else
            flows[:,i] += allornothing(network, z, trips[z,:], costs)
        end
    end

    return flows
end

allornothing(network::AbstractNetwork, trips::AbstractMatrix{T}, costfn::CostFunction; basedon::Symbol = :link) where {T<:Real} = allornothing(network, trips, costfn(zeros(numlinks(network))).costs; basedon=basedon)


text:
SQ._collect{R<:Relation}(_src::Tuple{R}, q::SQ.Node{SQ.Select}) =
    Relations._collect(first(_src), q)

function pre_apply{S<:Attributed,T}(r::Relation{S,T}, q::SQ.Node{SQ.Select})
    _cols = Vector{Any}()
    _attrs = [ d.res_field for d in SQ.dos(q) ]
    fs = [ d.f for d in SQ.dos(q) ]
    nrow = length(r)
    for f in fs
        U = Core.Inference.return_type(f, Tuple{T})
        if isleaftype(U)
            U <: Nullable ? push!(_cols, NullableArray{eltype(U)}(nrow)) :
                            push!(_cols, Array{U}(nrow))
        # else TODO: make this a real code path for when inferred type is non-concrete
        end
    end
    return _cols, _attrs, fs
end

function _collect{S<:Attributed,T}(src::Relation{S,T}, q::SQ.Node{SQ.Select})
    # @time _cols, _attrs, fs = pre_apply(q, T)
    # @time res = Relation(_cols, _attrs)
    # @time apply!(res, tuple(fs...), src)
    _cols, _attrs, fs = pre_apply(src, q)
    res = Relation(tuple(_attrs...), _cols...)
    apply!(res, tuple(fs...), src)
    return res
end

@generated function apply!(res, fs, src)
    n = length(fs.parameters)
    v_exs, setindex!_exs = Expr(:block), Expr(:block)
    for j in 1:n
        push!(v_exs.args,
              Expr(:(=), Symbol("v_$j"),
              Expr(:call, Expr(:call, :getfield, :fs, j), :row)))
        push!(setindex!_exs.args,
              Expr(:call, :setindex!,
                   Expr(:call, :getfield, :(res.src), j), Symbol("v_$j"), :i))
    end
    res = quote
        for (i, row) in enumerate(src)
            $v_exs; $setindex!_exs
        end
        return
    end
end


text:
# Julia 1.0
using Combinatorics
function make(str, test )
  filter(test, collect( permutations(split(str))) )
end

men = make("danish english german norwegian swedish",
           x -> "norwegian" == x[1])

drinks = make("beer coffee milk tea water", x -> "milk" == x[3])

#Julia 1.0 compatible
colors = make("blue green red white yellow",
              x -> 1 == findfirst(c -> c == "white", x) - findfirst(c -> c == "green",x))

pets = make("birds cats dog horse zebra")

smokes = make("blend blue-master dunhill pall-mall prince")

function eq(x, xs, y, ys)
  findfirst(xs, x) == findfirst(ys, y)
end

function adj(x, xs, y, ys)
  1 == abs(findfirst(xs, x) - findfirst(ys, y))
end

function print_houses(n, pet, nationality, colors, drink, smokes)
  println("$n, $pet,    $nationality       $colors       $drink    $smokes")
end

for m = men, c = colors
  if eq("red",c, "english",m) && adj("norwegian",m, "blue",c)
    for d = drinks
      if eq("danish",m, "tea",d) && eq("coffee",d,"green",c)
        for s = smokes
          if eq("yellow",c,"dunhill",s) &&
             eq("blue-master",s,"beer",d) &&
             eq("german",m,"prince",s)
            for p = pets
              if eq("birds",p,"pall-mall",s) &&
                 eq("swedish",m,"dog",p) &&
                 adj("blend",s,"cats",p) &&
                 adj("horse",p,"dunhill",s)
                println("Zebra is owned by ", m[findfirst(c -> c == "zebra", p)])
                println("Houses:")
                for house_num in 1:5
                  print_houses(house_num,p[house_num],m[house_num],c[house_num],d[house_num],s[house_num])
                end
              end
            end
          end
        end
      end
    end
  end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using PyPlot, MAT
;

using AD4SM
;

mean(x) = sum(x)/length(x)

#=
makeHexaModels is a helper function to create a N1xN2
hexagonal lattice with side length L0
=#
function makeHexaModels(;N1 = 5, N2=5, L0=1, dTol=1e-6)

  # unit cell and periodic directions
  a1 = sqrt(3)L0.*(cos(π/6), sin(π/6))
  a2 = sqrt(3)L0.*(cos(π/6), -sin(π/6))
  p0 = [[cos(θ)L0, sin(θ)L0] for θ = range(0, stop=5π/3, length=6)]

  # replicate the nodes
  nodes = []
  for jj in  range(1, stop=N2)
    for ii in range(0, stop=N1-1)
      Δ = a1.*floor(ii/2) .+ a2.*ceil(ii/2) .+ (jj-1).*(a1.-a2)
      for item in p0
        push!(nodes, [item[1]+Δ[1], item[2]+Δ[2]])
      end
    end
  end

  #remove duplicated nodes
  global jj = 1
  while jj < length(nodes)
    rmidx = []
    for (ii, p0) in enumerate(nodes[jj+1:end])
      dd = sqrt((p0[1]-nodes[jj][1])^2 + (p0[2]-nodes[jj][2])^2)
      (dd <= dTol) && push!(rmidx, ii)
    end
    deleteat!(nodes, jj.+rmidx)
    jj+=1
  end

  # add beams
  jj    = 1
  beams = []
  while jj < length(nodes)
    for (ii, p0) in enumerate(nodes[jj+1:end])
      dd = sqrt((p0[1]-nodes[jj][1])^2 + (p0[2]-nodes[jj][2])^2)
      (abs(dd-L0) <= dTol) && push!(beams, [jj, jj+ii])
    end
    jj+=1
  end

  return (nodes, beams)
end
;

@time (nodes, beams) = makeHexaModels(N1 = 21, N2 = 21, L0 = 1.2)

nNodes   = length(nodes)
t, w, Es = .1, .1, 5
A        = t*w
elems    = [Elements.Beam(beam, nodes[beam], t, w,
    mat=Materials.Hooke(Es, 0.3)) for beam in beams]
;

minx = minimum([node[1] for node in nodes])
maxx = maximum([node[1] for node in nodes])
miny = minimum([node[2] for node in nodes])
maxy = maximum([node[2] for node in nodes])

idbtm = findall(abs.([node[2] for node in nodes].-miny) .<1e-6)
idtop = findall(abs.([node[2] for node in nodes].-maxy) .<1e-6)

LY    = maxy-miny
;

ax = getproperty(figure(),:add_subplot)(1,1,1)
for item in beams
    n1 = nodes[item[1]]
    n2 = nodes[item[2]]
    PyPlot.plot(
        [n1[1], n2[1]],
        [n1[2], n2[2]],
        markersize=1, color=:b)
end
getproperty(ax, :set_aspect)("equal")
title("undeformed model")
;

u     = fill!(zeros(3, nNodes), NaN)
ΔY    = 0.75LY

u[:,idbtm] .= 0.0
u[:,idtop] .= 0.0
u[2,idtop] .= ΔY
;

bfreeu = isnan.(u)
uold   = zeros(size(u))
unew   = zeros(size(u))
unew[.!bfreeu] .= u[.!bfreeu]*1e-4

fe    = zeros(length(unew))
;

allus = Elements.solve(elems, u, N=20, bprogress=false, bechoi=false, becho=true, ballus=true)
;

unew = allus[end][1]
ax   = getproperty(figure(),:add_subplot)(1,1,1)

for beam in beams
    n1, n2 = nodes[beam]
    u1 = unew[:,beam[1]]
    u2 = unew[:,beam[2]]

    PyPlot.plot(
        [n1[1], n2[1]],
        [n1[2], n2[2]],
        markersize=1, color=:r, alpha=0.25)
    PyPlot.plot(
        [n1[1]+u1[1], n2[1]+u2[1]],
        [n1[2]+u1[2], n2[2]+u2[2]],
        markersize=1, color=:b)
end
getproperty(ax, :set_aspect)("equal")
title("deformed  model")
;

idd = LinearIndices(u)[2,idtop]

RY  = [sum(u[2][idd]) for u in allus];
Δu  = [mean(u[1][idd]) for u in allus]

PyPlot.plot(Δu/LY, RY/Es/A)
;

matwrite("HexaLattice.mat", Dict(
  "Ry"     => RY,
  "nodes"  => hcat(nodes[:]...)|> transpose |> collect,
  "beams"  => hcat(beams[:]...)|> transpose |> collect,
  "LY"     => LY,
  "DeltaY" => Δu,
  "u"      => unew,
  "Es"     => Float64(Es),
  "A"      => A))


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()

## ============================================================================
@info "load packages"
using CovidSurvey
using Random, Dates, Turing
using Base.Threads
using StatsPlots
using PrettyTables
using DataFrames
using BSON
using ColorSchemes
colors = ColorSchemes.tableau_10.colors

using Memoization
using ReverseDiff
setadbackend(:reversediff)
Turing.setrdcache(true)

plotlyjs()
@info "number of threads available: $(nthreads())"
## ============================================================================

ps = (
    warmup = 100,
    model = "deaths",
    steps  = 100,
    seed   = nothing,
    observ = "2021-03-01",#"2021-02-06",#"2021-03-25"
    chains = 1,
    preds = "CS,MG,MA",#"CF,CC,CR,CS"
    semipara=true,
    rwstep = 7,
    include= false,
    std = false,
)
## ============================================================================
# load data
@info "load data"

data_params = (
      observationsend  = ps.observ
    , predictors        = ps.preds|> CovidSurvey.parse_predictors
    # , cases_start       = ps.cases
    # , hospitmodel       = Regional.HospitInit1(obs_stop="2020-07-01")
    , casemodel         = Regional.CaseInit(obs_start="2020-06-01")
    , seromodel         = Regional.SeroInit(delay=0, std=1.)
    , rwstep            = ps.rwstep
    , epidemicstart    = 20
    , numimpute        = 6
    , include_early_dynamic = ps.include
    , link              = KLogistic(4.5)
    , invlink           = KLogit(4.5)
    , predictive        = false
    , covariates_kwargs = Dict(
      :semiparametric => ps.semipara,
      :fname          => [projectdir("data/contacts/dk/", "averaged_contact_rates_region=$(r).csv" ) for r in 1:Regional.nregions],#[projectdir("data", "mean_contact_rates_region=$r.csv" ) for r in 1:Regional.nregions],
      :shift          => -1,
      :startdate      => "2020-11-10", # >= 1000 cases / day
      :datecol        => "dates",
      :aggregation    => ps.rwstep,
      :standartize    => ps.std,
      :mobility       => [projectdir("data/mobility/mobility_region=$(Regional.regions[i]).csv") for i in 1:Regional.nregions]
      )
    )
data = Regional.load_data(; data_params... )
d = turing_data = data.turing_data;

gethospits(d,i) = DataFrame( :hospits => d.hospits[i] )
getdeaths(d,i) = DataFrame( :deaths => d.deaths[i] )
getcases(d,i) = DataFrame( :cases => d.cases[i] )
getdates(d,i) = DataFrame( :date => d.dates[i] )

function getrw(d, i)
    no = d.num_observations[i]
    cs = d.covariates_starts[i]
    ri = d.rt_step_indices[i]
    li = d.lockdown_indices[i]

    rw = ones(Int64, no)
    rw[li+1:cs-1] = ri .+ 1
    rw[cs:no] .=  rw[cs-1]
    return DataFrame(:rw => rw)
end

function getcovariates(d, i)
    no = d.num_observations[i]
    cs = d.covariates_starts[i]
    names = data.predictors
    co = d.covariates[i]
    np = length(names)

    dat = zeros(no, np)
    dat[cs:no, :] = co
    return DataFrame(dat, names)
end

function epidata(d, i, r)
    df = hcat(
          getdates(d,i)
        , getcases(d,i)
        , gethospits(d,i)
        , getdeaths(d,i)
        , getcovariates(d, i)
        , getrw(d, i)
    )
    df[!,:region] .= r
    return df
end

function epidata(data)
    return vcat( [epidata(data.turing_data, i, r) for (i, r) in enumerate(data.regions)]... )
end

epid = epidata(data)
CSV.write(projectdir("data/epidemia", "denmark_predictors=CS-MG-MA.csv"), epid)

plotlyjs()
let
    p = plot()
    for (i, r) in enumerate(data.regions)
        df = epidata(d, i, Regional.regions[i])
        plot!(df.date, cumsum(df.hospit), hover=["$i, cases: $n" for (i,n) in enumerate(cumsum(df.hospits))], c=colors[i])
        vline!([df.date[20]], c=colors[i])
    end
    hline!([30], c=colors[end])
    p
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
# using Base.Tests
using Test
using SemanticModels
using SemanticModels.Unitful
import SemanticModels.Unitful: DimensionError, uconvert, NoUnits, s
using DifferentialEquations
using Distributions: Uniform
using GLM
using DataFrames
using Plots

tests = ["parse.jl",
         "cassette.jl",
         "transform/ode.jl"]

for test in tests
  @testset "Running $test" begin
    include(test)
  end
end

examples = ["agentbased.jl",
            "agentgraft.jl",
            "modelmacro.jl",
            "monomial_regression.jl",
            "multivariate_regression.jl",
            "polynomial_regression.jl",
            "workflow.jl",
            "pseudo_polynomial_regression.jl",
            "odegraft.jl",
            ]
@testset "Test all examples" begin
  for ex in examples
      @info "Running example: " file=ex
      try
        include(ex)
        @test true == true
      catch err
        println(err)
        @info "Error running " file=ex
        @test true == false
      end
  end
end


text:
using Distributions, PyCall, PyPlot
@pyimport matplotlib.gridspec as gridspec

include(joinpath(dirname(dirname(@__FILE__)), "test", "common.jl"))
complex = !isempty(ARGS) && ARGS[1] == "complex"

function bench(fns, input)
    [t=>begin
        times = Array(Vector{Float64}, length(fns))
        for ifn = 1:length(fns)
            fn = fns[ifn]
            inp = input[t][ifn]
            fn(inp...)
            gc()
            nrep = max(iceil(2/(@elapsed (gc_disable(); fn(inp...); gc_enable(); gc()))), 3)
            println("Running $nrep reps of $fn($t)")
            @time times[ifn] = [begin
                gc()
                gc_disable()
                time = @elapsed fn(inp...)
                gc_enable()
                time
            end for i = 1:nrep]
            # println((mean(times[ifn]), std(times[ifn])))
        end
        times
    end for t in types]
end

function ratioci(y, x, alpha=0.05)
    tq² = abs2(quantile(TDist(length(x) + length(y) - 2), alpha))
    μx = mean(x)
    σx² = varm(x, μx)
    μy = mean(y)
    σy² = varm(y, μy)
    a = sqrt((μx*μy)^2 - (μx^2 - tq²*σx²)*(μy^2 - tq²*σy²))
    b = μx^2 - tq²*σx²
    (((μx*μy) - a)/b, ((μx*μy) + a)/b)
end

# First generate some random data and test functions in Base on it
const NVALS = 1_000_000
base_unary = complex ? base_unary_complex : base_unary_real
base_binary = complex ? base_binary_complex : base_binary_real
types = complex ? (Complex64, Complex128) : (Float32, Float64)
input = [t=>[[(randindomain(t, NVALS, domain),) for (fn, domain) in base_unary];
             [(randindomain(t, NVALS, domain1), randindomain(t, NVALS, domain2))
              for (fn, domain1, domain2) in base_binary];
             (randindomain(t, NVALS, (0, 100)), randindomain(t, 1, (-1, 20))[1])]
            for t in types]
fns = [[x[1] for x in base_unary]; [x[1] for x in base_binary]; (complex ? [] : .^)]

builtin = bench(fns, input)

# Now with VML
using VML

vml = bench(fns, input)

# Print ratio
clf()
colors = ["r", "y"]
for itype = 1:length(types)
    builtint = builtin[types[itype]]
    vmlt = vml[types[itype]]
    μ = vec(map(mean, builtint)./map(mean, vmlt))
    ci = zeros(Float64, 2, length(fns))
    for ifn = 1:length(builtint)
        lower, upper = ratioci(builtint[ifn], vmlt[ifn])
        ci[1, ifn] = μ[ifn] - lower
        ci[2, ifn] = upper - μ[ifn]
    end
    bar(0.2+(0.4*itype):length(fns), μ, 0.4, yerr=ci, color=colors[itype], ecolor="k")
end
ax = gca()
ax[:set_xlim](0, length(fns)+1)
fname = [string(fn.env.name) for fn in fns]
if !complex
    fname[end-1] = "A.^B"
    fname[end] = "A.^b"
end
xticks(1:length(fns)+1, fname, rotation=70, fontsize=10)
title("VML Performance")
ylabel("Relative Speed (Base/VML)")
legend([string(x) for x in types])
ax[:axhline](1; color="black", linestyle="--")
savefig("performance$(complex ? "_complex" : "").png")


text:
#
# Description of the input and the output structures.
#

# Canonical singleton and empty types.
typealias Unit Void
typealias Zero Union{}

# The type has one or no elements.
isunit(T::Type) = T <: Void || T <: Tuple{} || T <: Nullable{Union{}} || T <: Vector{Union{}}
iszero(T::Type) = T <: Union{}

#
# Structure of input.
#

# Maps a tuple of parameter names to a tuple of types.
typealias Params Pair{Tuple{Vararg{Symbol}}, Tuple{Vararg{Type}}}
const NoParams = Params((), ())

immutable InputMode
    # Depends on the preceeding and subsequent input values.
    relative::Bool
    # List of parameters with types.
    params::Params
end

InputMode(;relative::Bool=false, params::Params=NoParams) = InputMode(relative, params)

InputMode(mode::InputMode; relative=nothing, params=nothing) =
    InputMode(
        relative != nothing ? relative : mode.relative,
        params != nothing ? params : mode.params)

# Predicates and properties.
isfree(mode::InputMode) = !mode.relative && isempty(mode.params.first)
isrelative(mode::InputMode) = mode.relative
params(mode::InputMode) = mode.params

# Lattice on input structures.
typemin(::Type{InputMode}) = InputMode(false, NoParams)
max(mode::InputMode) = mode
max(mode1::InputMode, mode2::InputMode) =
    begin
        relative = max(mode1.relative, mode2.relative)
        if mode2.params == NoParams || mode1.params == mode2.params
            params = mode1.params
        elseif mode1.params == NoParams
            params = mode2.params
        else
            # Merging parameters.
            p1 = Dict{Symbol,Type}(zip(mode1.params.first, mode1.params.second))
            p2 = Dict{Symbol,Type}(zip(mode2.params.first, mode2.params.second))
            names = Symbol[]
            types = Type[]
            for name in sort(unique([collect(keys(p1)); collect(keys(p2))]))
                T = !haskey(p1, name) ? p2[name] :
                    !haskey(p2, name) ? p1[name] : typeintersect(p1[name], p2[name])
                push!(names, name)
                push!(types, T)
            end
            params = Params((names...), (types...))
        end
        return InputMode(relative, params)
    end

# Type and structure of input.
immutable Input
    domain::Type
    mode::InputMode
end

Input(domain::Type; relative::Bool=false, params::Params=NoParams) =
    Input(domain, InputMode(relative, params))

Input(input::Input; domain=nothing, relative=nothing, params=nothing) =
    Input(
        domain != nothing ? domain : input.domain,
        InputMode(input.mode; relative=relative, params=params))

show(io::IO, input::Input) =
    begin
        mode = input.mode
        if !isempty(mode.params.first)
            print(io, "{")
        end
        print(io, input.domain)
        if mode.relative
            print(io, "...")
        end
        for (name, T) in zip(mode.params.first, mode.params.second)
            print(io, ", $name => $T")
        end
        if !isempty(mode.params.first)
            print(io, "}")
        end
    end

domain(input::Input) = input.domain
mode(input::Input) = input.mode

typemin(::Type{Input}) = Input(Any, typemin(InputMode))
max(input::Input) = input
max(input1::Input, input2::Input) =
    Input(typeintersect(input1.domain, input2.domain), max(input1.mode, input2.mode))

# Predicates and properties.
isfree(input::Input) = isfree(input.mode)
isrelative(input::Input) = isrelative(input.mode)
params(input::Input) = params(input.mode)

# How the value is represented in the pipeline.
kind(input::Input) =
    let T = input.domain,
        mode = input.mode
        if !isempty(mode.params.first)
            Ns = Val{mode.params.first}
            Vs = Tuple{mode.params.second...}
            input.mode.relative ?
                EnvRel{Ns,Vs,T} : Env{Ns,Vs,T}
        elseif input.mode.relative
            Rel{T}
        else
            Iso{T}
        end
    end

data(input::Input) = data(kind(input))

convert(input::Input, args...; params...) =
    convert(kind(input), args...; params...)

#
# Structure of output.
#

immutable OutputMode
    # At most one output value for each input value.
    lunique::Bool
    # At least one output for each input value.
    ltotal::Bool
    # At most one input for each output value.
    runique::Bool
    # At least one input for each output value.
    rtotal::Bool
end

# Default is a plain function.
OutputMode(; lunique::Bool=true, ltotal::Bool=true, runique::Bool=false, rtotal::Bool=false) =
    OutputMode(lunique, ltotal, runique, rtotal)

OutputMode(mode::OutputMode; lunique=nothing, ltotal=nothing, runique=nothing, rtotal=nothing) =
    OutputMode(
        lunique != nothing ? lunique : mode.lunique,
        ltotal != nothing ? ltotal : mode.ltotal,
        runique != nothing ? runique : mode.runique,
        rtotal != nothing ? rtotal : mode.rtotal)

# Predicates.
isplain(mode::OutputMode) = mode.lunique && mode.ltotal
ispartial(mode::OutputMode) = mode.lunique && !mode.ltotal
isplural(mode::OutputMode) = !mode.lunique
issingular(mode::OutputMode) = mode.lunique
isnonempty(mode::OutputMode) = mode.ltotal
ismonic(mode::OutputMode) = mode.runique
iscovering(mode::OutputMode) = mode.rtotal

# Lattice.
typemin(::Type{OutputMode}) = OutputMode(true, true, true, true)
max(mode::OutputMode) = mode
max(mode1::OutputMode, mode2::OutputMode) =
    OutputMode(
        min(mode1.lunique, mode2.lunique),
        min(mode1.ltotal, mode2.ltotal),
        min(mode1.runique, mode2.runique),
        min(mode1.rtotal, mode2.rtotal))

# Type and structure of output.
immutable Output
    domain::Type
    mode::OutputMode
end

Output(domain::Type; lunique::Bool=true, ltotal::Bool=true, runique::Bool=false, rtotal::Bool=false) =
    Output(domain, OutputMode(lunique, ltotal, runique, rtotal))

Output(output::Output; domain=nothing, lunique=nothing, ltotal=nothing, runique=nothing, rtotal=nothing) =
    Output(
        domain != nothing ? domain : output.domain,
        OutputMode(output.mode, lunique=lunique, ltotal=ltotal, runique=runique, rtotal=rtotal))

show(io::IO, output::Output) =
    begin
        T = output.domain
        props = []
        if isplural(output)
            T = Vector{T}
            if output.mode.ltotal
                push!(props, :nonempty)
            end
        elseif ispartial(output)
            T = Nullable{T}
        end
        if output.mode.runique
            push!(props, :unique)
        end
        if output.mode.rtotal
            push!(props, :covering)
        end
        print(io, T)
        if !isempty(props)
            print(io, " # $(join(props, ", "))")
        end
    end

domain(output::Output) = output.domain
mode(output::Output) = output.mode

typemin(::Type{Output}) = Output(Null, typemin(OutputMode))
max(output::Output) = output
max(output1::Output, output2::Output) =
    Output(typejoin(output1.domain, output2.domain), max(output1.mode, output2.mode))

# Predicates.
isplain(output::Output) = isplain(output.mode)
ispartial(output::Output) = ispartial(output.mode)
isplural(output::Output) = isplural(output.mode)
issingular(output::Output) = issingular(output.mode)
isnonempty(output::Output) = isnonempty(output.mode)
ismonic(output::Output) = ismonic(output.mode)
iscovering(output::Output) = iscovering(output.mode)

# How the value is represented in the pipeline.
kind(output::Output) =
    let T = output.domain,
        mode = output.mode
        mode.lunique && mode.ltotal ? Iso{T} : mode.lunique ? Opt{T} : Seq{T}
    end

data(output::Output) = data(kind(output))

convert(output::Output, args...) =
    convert(kind(output), args...)

#
# Properties of any mapping that defines `input` and `output` methods.
#

input(input::Input) = input
output(output::Output) = output
input(mapping::Pair{Input,Output}) = mapping.first
output(mapping::Pair{Input,Output}) = mapping.second
mapping(input::Input, output::Output) = input => output
mapping(mapping) = input(mapping) => output(mapping)
idomain(mapping) = domain(input(mapping))
odomain(mapping) = domain(output(mapping))
imode(mapping) = mode(input(mapping))
omode(mapping) = mode(output(mapping))
ikind(mapping) = kind(input(mapping))
okind(mapping) = kind(output(mapping))
idata(mapping) = data(input(mapping))
odata(mapping) = data(output(mapping))
isfree(mapping) = isfree(input(mapping))
isrelative(mapping) = isrelative(input(mapping))
params(mapping) = params(input(mapping))
isplain(mapping) = isplain(output(mapping))
ispartial(mapping) = ispartial(output(mapping))
isplural(mapping) = isplural(output(mapping))
issingular(mapping) = issingular(output(mapping))
isnonempty(mapping) = isnonempty(output(mapping))
ismonic(mapping) = ismonic(output(mapping))
iscovering(mapping) = iscovering(output(mapping))

show(io::IO, mapping::Pair{Input,Output}) =
    begin
        if mapping.first != Input(Unit)
            print(io, mapping.first, " -> ")
        end
        print(io, mapping.second)
    end

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 Cstr
using LightXML, Printf
using DifferentialEquations, Sundials
include("Utils.jl")
include("Constants.jl")
include("SurfaceReactions.jl")

struct ConstParams{T1}
    ρ_in # inlet density
    q_in    #inlet flow rate
    avg_molwt_in   #average molecular weight at the inlet
    As  #surface area
    V   #reactor volume
    mass_fracs_in::Array{T1,1}  #inlet mass fractions
    T   # temperature
    p   # pressure
    ρq_yk::Array{T1,1} #inlet mass flow rate ρ x q x y_k
end

function cstr(input_file::AbstractString, lib_dir::AbstractString, sens= false)

    local gasphase = Array{String,1}

    xmldoc = parse_file(input_file)

    xmlroot = root(xmldoc)

    #get the gasphase from xml input
    gasphase = get_collection_from_xml(xmlroot,"gasphase")

    #create the thermo object
    thermo_file = lib_dir*"therm.dat"
    thermo_obj = SurfaceReactions.IdealGas.create_thermo(gasphase,thermo_file)

    #get the molefractions from xml input
    mole_fracs = get_molefraction_from_xml(xmlroot,thermo_obj.molwt,gasphase)
    mass_fracs_in = zeros(length(gasphase))
    molefrac_to_massfrac!(mass_fracs_in,mole_fracs,thermo_obj.molwt)
    avg_molwt = average_molwt(mole_fracs,thermo_obj.molwt)
    #get the operating temperature
    local T = get_value_from_xml(xmlroot,"T")

    #get the initial pressure
    local p = get_value_from_xml(xmlroot,"p")

    #get the reactor volume
    local V = get_value_from_xml(xmlroot,"volume")

    #get the flow rate. The flow rate is expected in m3/s
    local q = get_value_from_xml(xmlroot,"flow-rate")

    #get the surface area
    local As = get_value_from_xml(xmlroot,"As")

    #get the mechanism input file
    local mech_file = get_text_from_xml(xmlroot,"surface_mech")
    mech_file = lib_dir*"/"*mech_file
    md = SurfaceReactions.compile_mech(mech_file,thermo_obj,gasphase)

    #get the integration time
    it = get_value_from_xml(xmlroot,"time")

    #density at the inlet
    ρ = density(mole_fracs,thermo_obj.molwt,T,p)
    #toal number of species
    n_species = length(gasphase) + length(md.sm.species)
    #definition of array for the storage of reaction rates
    rates = zeros(n_species)
    #solution vector
    sol = copy(mass_fracs_in)
    #initial coverage as specified in the mechanism
    covg = md.sm.si.ini_covg
    #append the coverage to the solution vector
    append!(sol,md.sm.si.ini_covg)

    #file output stream for saving the data
    g_stream = open("gas_profile.dat","w")
    s_stream = open("surf_profile.dat","w")
    o_streams = (g_stream, s_stream)
    create_header(g_stream,["t","T","p","rho"],gasphase)
    create_header(s_stream,"t","T",md.sm.species)

    #create state object
    state = SurfaceReactions.State(mole_fracs,covg,T,p,rates)
    #inlet conditions
    t_span = (0,it)
    #the following is the first term in the right hand side of gasphase species residual
    ρq_mass_fracs = mass_fracs_in * (ρ*q)
    cp = ConstParams(ρ,q,avg_molwt,As,V,mass_fracs_in, T, p,ρq_mass_fracs)
    params = (state, thermo_obj, md, cp, o_streams)
    prob = ODEProblem(residual!,sol,t_span,params)
    if sens == true
        return (params, prob, t_span)
    end
    cb = FunctionCallingCallback(save_data)
    soln = solve(prob,alg_hints=[:stiff] , reltol=1e-6, abstol = 1e-8, save_everystep=false, callback=cb)
    # println("Solver Integration: ", soln.retcode, "\t", typeof(soln.retcode))
    #close the files
    close(g_stream)
    close(s_stream)
    return soln.retcode
end

function residual!(du,u,p,t)
    state, thermo_obj, md, cp, o_streams = p
    ng = length(state.mole_frac)
    ns = length(md.sm.species)
    mass_fracs = u[1:ng]
    #update the state with latest mole fractions
    massfrac_to_molefrac!(state.mole_frac,mass_fracs,thermo_obj.molwt)
    #update the state with latest coverages
    state.covg = u[ng+1:ng+ns]
    #calculate the density
    ρ = density(state.mole_frac,thermo_obj.molwt,cp.T,cp.p)
    #calculate the molar production rates
    SurfaceReactions.calculate_molar_production_rates!(state,thermo_obj,md)
    #outlet volumetric flow rate
    q_out = cp.q_in*cp.avg_molwt_in/average_molwt(state.mole_frac,thermo_obj.molwt)
    #gasphase species residual
    spIn = cp.ρq_yk/cp.V
    spOut = q_out * ρ * mass_fracs /cp.V
    rgVec = (state.s_rate[1:ng] .* thermo_obj.molwt)*cp.As/cp.V
    du[1:ng] = (spIn - spOut + rgVec)/ρ
    #surface species residual
    du[ng+1:ng+ns] = (state.s_rate[ng+1:ng+ns] .* md.sm.si.site_coordination)/(md.sm.si.density*1e4)
end

function save_data(u,t,integrator)
    state = integrator.p[1]
    thermo_obj = integrator.p[2]
    g_stream, s_stream = integrator.p[5]
    d = density(state.mole_frac,thermo_obj.molwt,state.T,state.p)
    write_to_file(g_stream,t,state.T,state.p,d,state.mole_frac)
    write_to_file(s_stream,t,state.T,state.covg)
    @printf("%.4e\n",t)
end

end


text:
module Tree

#using StatsBase
include("../utils.jl")

export DecisionTree, MV, RandomForest

include("utils.jl")
include("DT.jl")
include("RF.jl")
include("MV.jl")

end


text:
# already contained in Core
# struct OverflowError <: Exception
#     msg::String
# end

struct UnderflowError <: Exception
    msg::String
end

function overflow_error(msg = "")
    throw( OverflowError(msg) )
end

function underflow_error(msg = "")
    throw( UnderflowError(msg) )
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# https://discourse.julialang.org/t/generation-of-documentation-fails-qt-qpa-xcb-could-not-connect-to-display/60988
ENV["GKSwstype"] = "100"

using Literate

files = Set(readdir())
exclude_files = Set(["output", "build.jl", "notebook.jl", "index.jl", "index.md", "README.md"])
target_files = setdiff(files, exclude_files)

function replace_includes(str)
    included = target_files

    # Here the path loads the files from their proper directory,
    # which may not be the directory of the `examples.jl` file!
    path = "./"

    for ex in included
        content = read(path*ex, String)
        str = replace(str, "include(\"$(ex)\")" => content)
    end
    return str
end

function update_file_name(content, file_name)
    content = replace(content, "FILE_NAME" => file_name)
    return content
end

cp("../Project.toml", "output/Project.toml", force=true)
cp("../Manifest.toml", "output/Manifest.toml", force=true)
Literate.notebook("notebook.jl", "output", config=Dict("execute"=>false), preprocess = replace_includes, name = "PlotsGallery")
rm("output/Project.toml")
rm("output/Manifest.toml")

# Insert Table of Contents
using Pkg
Pkg.add("JSON")
import JSON
json = JSON.parsefile("output/PlotsGallery.ipynb")
cells = json["cells"]
sources = []
for cell in cells
    if cell["cell_type"] == "markdown"
        append!(sources, cell["source"])
    end
end

toc = []
for item in sources
    m = match(r"#+ [^\n]*", item)
    if !isnothing(m)
        push!(toc, m.match)
    end
end

function to_toc(item)
    h, c = match(r"(#+) ([^\n]+)", item).captures
    link = replace(c, " " => "-")
    n = length(h) - 1
    return repeat(" ", 2 * n) * "1. [$(c)](#$(link))"
end

str = join(to_toc.(toc), "\n")
str = "# Table of Contents\n" * str
toc_cell = Dict("cell_type" => "markdown", "source" => [str], "metadata" => Dict())
pushfirst!(json["cells"], toc_cell)

open("output/PlotsGallery.ipynb", "w") do f
    print(f, JSON.json(json))
end

# Generate markdown for site
for file in target_files
    Literate.markdown(file, "output/src", config=Dict("documenter"=>false),
                      preprocess = x -> update_file_name(x, splitext(file)[1]))
end

# Generate index.md
include("index.jl")

# Generate figures
ENV["GKS_ENCODING"]="utf-8" # for scientific formatter
using Plots;
gr()

for file in target_files
    include(file)
    name = splitext(file)[1] * ".png"
    OUTPUT_DIR = "output/src/images/"
    savefig(OUTPUT_DIR * name)
    run(`convert $(OUTPUT_DIR * name) -resize 30% $(OUTPUT_DIR)small_$(name)`)
end


text:
# TODO we should have an instruction for defining new problems
# module NewTopOptProblems
# export NewPointLoadCantilever

using Ferrite
using TopOpt
using TopOpt.TopOptProblems: RectilinearGrid, Metadata
using TopOpt.TopOptProblems:
    left,
    right,
    bottom,
    middley,
    middlez,
    nnodespercell,
    nfacespercell,
    find_black_and_white,
    find_varind
using TopOpt.Utilities: @params

@params struct NewPointLoadCantilever{dim,T,N,M} <: StiffnessTopOptProblem{dim,T}
    rect_grid::RectilinearGrid{dim,T,N,M}
    E::T
    ν::T
    ch::ConstraintHandler{<:DofHandler{dim,<:Cell{dim,N,M},T},T}
    load_dict::Dict{Int,Vector{T}}
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    metadata::Metadata
end
# force::T
# force_dof::Integer

function NewPointLoadCantilever(
    ::Type{Val{CellType}},
    nels::NTuple{dim,Int},
    sizes::NTuple{dim},
    E=1.0,
    ν=0.3,
    force=1.0,
) where {dim,CellType}
    iseven(nels[2]) && (length(nels) < 3 || iseven(nels[3])) ||
        throw("Grid does not have an even number of elements along the y and/or z axes.")

    _T = promote_type(eltype(sizes), typeof(E), typeof(ν), typeof(force))
    if _T <: Integer
        T = Float64
    else
        T = _T
    end
    if CellType === :Linear || dim === 3
        rect_grid = RectilinearGrid(Val{:Linear}, nels, T.(sizes))
    else
        rect_grid = RectilinearGrid(Val{:Quadratic}, nels, T.(sizes))
    end

    if haskey(rect_grid.grid.facesets, "fixed_all")
        pop!(rect_grid.grid.facesets, "fixed_all")
    end
    #addfaceset!(rect_grid.grid, "fixed_all", x -> left(rect_grid, x));
    addnodeset!(rect_grid.grid, "fixed_all", x -> left(rect_grid, x))

    if haskey(rect_grid.grid.nodesets, "down_force")
        pop!(rect_grid.grid.nodesets, "down_force")
    end
    if dim == 3
        addnodeset!(
            rect_grid.grid, "down_force", x -> right(rect_grid, x) && bottom(rect_grid, x)
        )
        #  && middlez(rect_grid, x));
    else
        addnodeset!(
            rect_grid.grid,
            "down_force",
            x -> right(rect_grid, x) && right(rect_grid, x) && middley(rect_grid, x),
        )
    end

    # Create displacement field u
    dh = DofHandler(rect_grid.grid)
    if CellType === :Linear || dim === 3
        push!(dh, :u, dim) # Add a displacement field
    else
        ip = Lagrange{2,RefCube,2}()
        push!(dh, :u, dim, ip) # Add a displacement field
    end
    close!(dh)

    ch = ConstraintHandler(dh)

    #dbc = Dirichlet(:u, getfaceset(rect_grid.grid, "fixed_all"), (x,t) -> zeros(T, dim), collect(1:dim))
    dbc = Dirichlet(
        :u, getnodeset(rect_grid.grid, "fixed_all"), (x, t) -> zeros(T, dim), collect(1:dim)
    )
    add!(ch, dbc)
    close!(ch)
    t = T(0)
    Ferrite.update!(ch, t)

    metadata = Metadata(dh)
    load_dict = Dict{Int,Vector{T}}()
    for fnode in getnodeset(rect_grid.grid, "down_force")
        load_dict[fnode] = [0, -force, 0]
    end

    N = nnodespercell(rect_grid)
    M = nfacespercell(rect_grid)

    black, white = find_black_and_white(dh)
    varind = find_varind(black, white)

    return NewPointLoadCantilever(
        rect_grid, E, ν, ch, load_dict, black, white, varind, metadata
    )
end

# used in FEA to determine default quad order
# we don't assume the problem struct has `rect_grid` to define its grid
TopOptProblems.nnodespercell(p::NewPointLoadCantilever) = nnodespercell(p.rect_grid)

# ! important, used for specification!
function TopOptProblems.getcloaddict(p::NewPointLoadCantilever{dim,T}) where {dim,T}
    # f = T[0, -p.force, 0]
    # fnode = Tuple(getnodeset(p.rect_grid.grid, "down_force"))[1]
    # return Dict{Int, Vector{T}}(fnode => f)
    return p.load_dict
end

# end # end NewTopOptProblems module


text:
# Mohammad Asif Zaman
# Date: June, 2018

# A very simple function to display a progress bar for a loop.
# Can be easily integrated in jupyter notebook.

# Progress_bar function.
# To use as a separate file: Should be added as a progress_bar.jl file in the same folder as your code. In the main code, use: include("progress_bar.jl");
# It is necessary to have PyCall pakcage installed for the include command to work.
# To use in Jupyter notebook, just include the funciton in a cell and run it.

# Input arguments:
# current_step = current position in the loop. Use the loop index.
# total_step = total number of steps in your loop. Note that it is assumed that the loop starts from 1 and ends at total_step.
# Nbar = how many blocks you want to include in the progress bar.

# The code has been tested using Atom IDE and Jupyter notebook.

# A simple code to test the progress_bar function.

using PyCall

include("progress_bar.jl");

max_count = 40;    # define maximum count of the loop
Nblock = 20;       # number of blocks in the bar plot

println("Starting Program...");
tic()
for i = 1:max_count
    progress_bar(i,max_count,Nblock);

    for k = 1:5e6    end   # delay loop. sleep() command causes some unexpected issues.

end

println("\nFinished Running.");
toc();


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    veccross(vec1::Matrix{F},vec2::Matrix{F}
            )::Matrix{F} where {F<:AbstractFloat}

or

    veccross(vec1::Vector{F},vec2::Vector{F}
            )::Vector{F} where {F<:AbstractFloat}

Determines cross product of two vectors.

# Arguments:
- `vec1`, `vec2`: An array of column vectors represented by a matrix of size
  `(dim,nv)`, where `dim` is the dimension of the vector and `nv` the number
  of vectors.

# Output:
- `cross`: Array of column vectors, each element is corresponding to the cross
  product of the respective components in `vec1` and `vec2`.

# Examples:
```julia
julia> cross = veccross(vec2,vec2)
```
Determine the cross product of

    [2.3;3.4;5.6] and [1.2;4.5;1.2]
    [5.1;0.0;2.3] and [2.5;3.2;4.0]

```julia
julia> cross=veccross([2.3 5.1; 3.4 0.0; 5.6 2.3],[1.2 2.5; 4.5 3.2; 1.2 4.0])
```
"""
function veccross(vec1::Matrix{F},vec2::Matrix{F}
                 )::Matrix{F} where {F<:AbstractFloat}
if size(vec1,1) == 2 && size(vec2,1) == 2
  # 2D vector
  cross = zeros(F,(3,size(vec1,2)));
  cross[3,:] .= vec1[1,:] .* vec2[2,:] .- vec1[2,:] .* vec2[1,:];
  return cross;
elseif size(vec1,1) == 3 && size(vec2,1) == 3
  # 3D vector
  return [vec1[2:2,:] .* vec2[3:3,:] .- vec1[3:3,:] .* vec2[2:2,:];
          vec1[3:3,:] .* vec2[1:1,:] .- vec1[1:1,:] .* vec2[3:3,:];
          vec1[1:1,:] .* vec2[2:2,:] .- vec1[2:2,:] .* vec2[1:1,:]];
else
  throw(ArgumentError("Input has to be either 2D or 3D with equal dims!"));
end # if
end # veccross

using LinearAlgebra;

# TODO: look if this is even necessary for this toolbox (using LinearAlgebra)
# overloaded function for vector input to get a float instead of a vector
function veccross(vec1::Vector{F},vec2::Vector{F}
                 )::Vector{F} where {F<:AbstractFloat}
return cross(vec1,vec2);
end # veccross


text:
#function _create_mesh_from_elements(is3D::Bool,
#                                    name::String,
#                                    points::Vector{Point{3,T}},
#                                    element_vecs::Vector{Vector{UInt64}},
#                                    element_sets::Dict{String, BitSet}) where {T}
#    # Determine element lengths
#    element_lengths = Int64[]
#    for element in element_vecs
#        l = length(element)
#        if l ∉ element_lengths
#            push!(element_lengths, l)
#        end
#    end
#    sort!(element_lengths)
#    U = _select_UInt_type(length(points))
#    if !is3D # is2D
#        K = 2
#        # Verify all points are approximately the same z-coordinate
#        if any(i->ϵ_Point < abs(points[1][3] - points[i][3]), 2:length(points))
#            error("Points of 2D mesh do not lie in the same plane")
#        end
#        points2D = convert.(Point{2,T}, points) # convert to 2D
#        if all(x->x < 6, element_lengths) # Linear mesh
#            P = 1
#        else
#            P = 2
#        end
#        polytopes = [ Polytope{K,P,length(elem),U}(elem) for elem in element_vecs ]
#        # Verify convexity
#        return PolytopeVertexMesh(name, points2D, polytopes, element_sets)
#    else
#        K = 3
#        if all(x->x < 10, element_lengths) # Linear mesh
#            P = 1
#        else # Quadratic Mesh
#            P = 2
#        end
#        polytopes = [ Polytope{K,P,length(elem),U}(elem) for elem in element_vecs ]
#        # Verify convexity
#        return PolytopeVertexMesh(name, points, polytopes, element_sets)
#    end
#    error("Invalid mesh type")
#end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport