text:
using Colors

import Escher: @d

getdoc(fn) =
    Escher.escher_meta[fn]

const badge_themes = @d(
    :default => fillcolor("#e1f0e5"),
    :curried => fillcolor("#72b5a1"),
    :typedkwarg => fillcolor("#444A95"),
    :kwarg => fillcolor("#DC7F38"),
    :typedarg => fillcolor("#DCA938"),
)

badge(x, theme=:default) =
    fontsize(0.94em, x) |>
    pad([left, right], 0.4em) |>
    badge_themes[theme] |>
    pad(0.2em)

#=
@doc """
format an argument metadata for use in function signature.

an optional argument is formatted as:

    [arg::T=value, ]

a curried argument is formatted as:

    «arg::T»

The Type will be of two colors. A green color means that
convert(::Type{T}, x) is called on the argument (x). While
a red color means that the argument must be of the specified type.

""" ->
=#
function argsig(arg)
    name = string(arg[:name])
    typ = string(arg[:type])
    typ_color = arg[:coerced] ? "#49796B" : "#A45A52" # Hooker green and redwood

    sig = [name, fontcolor("#aaa", "::"), fontcolor(typ_color, typ)]
    if haskey(arg, :default)
        sig = [sig, "=", fontcolor("#888", print_val(arg[:default]))]
        if !arg[:kwarg]
            sig = ["[ ", sig, " ]"]
        end
    end

    if arg[:curried]
        sig = ["« ", abbr("Curried", sig), " »"]
    end

    hbox(sig)
end

#= @doc """
Show the signature of a function.
""" -> =#
function signature(meta)
    fn = string(meta[:name])
    args   = map(argsig, filter(x -> !x[:kwarg], meta[:args]))
    kwargs = map(argsig, filter(x -> x[:kwarg], meta[:args]))

    sep = [",", hskip(0.5em)]
    sig = [ fontweight(bold, fn), "(", intersperse(sep, args)... ]
    if !isempty(kwargs)
        sig = [ sig, "; ", hskip(0.5em), intersperse(sep, kwargs) ]
    end
    sig = [ sig, ")" ]
    hbox(sig..., hskip(0.5em), "→", hskip(0.5em), rettype(meta[:type]) |> fontcolor("#777"))
end

rettype(typ) =
    typ.head == :curly ?
        string(typ.args[2]) :
        typ.head == :comparison ?
            abbr(string(typ.args[1]), string(typ.args[3])) :
            string(typ)

rettype(typ::Symbol) = string(typ)

print_val(x::Symbol) = ":$x"
print_val(x::Expr) = sprint(io -> print(io, x))
print_val(x) = sprint(io -> show(io, x))

function argrow(arg)
    notes = Any[]
    if !arg[:coerced]
        push!(notes, badge(hbox("type ", hskip(0.5em), code(string(arg[:type]))), :typedarg))
    end
    if arg[:kwarg]
        push!(notes, badge(abbr("Keyword argument", "kwarg"), :kwarg))
    end
    if arg[:curried]
        push!(notes, badge(abbr("Curried argument: if you leave this argument out, you will get a function which takes this argument", "curried"), :curried))
    end
    if haskey(arg, :default)
        push!(notes, badge("default=" * print_val(arg[:default])))
    end
    Any[code(string(arg[:name])), vskip(2em), hskip(1em), hbox(notes) |> wrap, hskip(1em), arg[:doc]]'
end

function argstable(args)
    table(
        Union{}[],
        reduce(vcat, map(argrow, args))
    )
end

band(t, bg="#f1f1f1", fg="#000") =
      fontweight(500, t) |> fontcolor(fg) |> fontstyle(italic)

import Escher:@api, render

@api table => (Table <: Tile) begin
    arg(headers::AbstractVector)
    arg(columns::AbstractMatrix)
end

render_cell(x, state) = Elem(:td, render(x, state))
render_row(xs, state) = Elem(:tr, map(x -> render_cell(x, state), xs))

render(t::Table, state) =
    Elem(:table, [
        isempty(t.headers) ? [] : Elem(:thead, Elem(:th, map(h -> render_cell(h, state), t.headers))),
        Elem(:tbody, map(row -> render_row(row, state), Any[sub(t.columns, i, :) for i in 1:size(t.columns, 1)]))
    ])

function showdoc(fn)
    d = getdoc(fn)
    docmd = d[:doc]

    vbox(
        signature(d) |> wrap |> fonttype(monospace),
        vbox(
             vskip(1em),
             pad([left], 1em, docmd),
             vskip(1em),
             band("Arguments"),
             argstable(d[:args]) |> pad([left, right], 1em) |> pad([top, bottom], 0.5em)
        ) |> pad([left], 1em)
    )
end

showdocs(fns) =
   vbox(intersperse([vskip(2em),
        border([bottom], solid, 1Escher.px, colorant"#ddd", empty),
        vskip(2em)],
        map(showdoc, fns)))


text:
__precompile__()

module Cosmology

using QuadGK, Unitful
import Unitful: km, s
using UnitfulAstro: Mpc, Gpc, Gyr

export cosmology,
       age,
       angular_diameter_dist,
       comoving_radial_dist,
       comoving_transverse_dist,
       comoving_volume,
       comoving_volume_element,
       distmod,
       H,
       hubble_dist,
       hubble_time,
       luminosity_dist,
       lookback_time,
       scale_factor

abstract type AbstractCosmology end
abstract type AbstractClosedCosmology <: AbstractCosmology end
abstract type AbstractFlatCosmology <: AbstractCosmology end
abstract type AbstractOpenCosmology <: AbstractCosmology end

struct FlatLCDM{T<:Real} <: AbstractFlatCosmology
    h::T
    Ω_Λ::T
    Ω_m::T
    Ω_r::T
end
FlatLCDM(h::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real) =
    FlatLCDM(promote(float(h), float(Ω_Λ), float(Ω_m), float(Ω_r))...)

a2E(c::FlatLCDM, a) = sqrt(c.Ω_r + c.Ω_m*a + c.Ω_Λ*a^4)

struct ClosedLCDM{T<:Real} <: AbstractClosedCosmology
    h::T
    Ω_k::T
    Ω_Λ::T
    Ω_m::T
    Ω_r::T
end
ClosedLCDM(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real) =
    ClosedLCDM(promote(float(h), float(Ω_k), float(Ω_Λ), float(Ω_m),
                       float(Ω_r))...)

struct OpenLCDM{T<:Real} <: AbstractOpenCosmology
    h::T
    Ω_k::T
    Ω_Λ::T
    Ω_m::T
    Ω_r::T
end
OpenLCDM(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real) =
    OpenLCDM(promote(float(h), float(Ω_k), float(Ω_Λ), float(Ω_m),
                     float(Ω_r))...)

function a2E(c::Union{ClosedLCDM,OpenLCDM}, a)
    a2 = a*a
    sqrt(c.Ω_r + c.Ω_m*a + (c.Ω_k + c.Ω_Λ*a2)*a2)
end

for c in ("Flat", "Open", "Closed")
    name = Symbol("$(c)WCDM")
    @eval begin
        struct $(name){T<:Real} <: $(Symbol("Abstract$(c)Cosmology"))
            h::T
            Ω_k::T
            Ω_Λ::T
            Ω_m::T
            Ω_r::T
            w0::T
            wa::T
        end
        function $(name)(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real,
                         w0::Real, wa::Real)
            $(name)(promote(float(h), float(Ω_k), float(Ω_Λ), float(Ω_m),
                            float(Ω_r), float(w0), float(wa))...)
        end
    end
end

function WCDM(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real, w0::Real, wa::Real)
    if Ω_k < 0
        ClosedWCDM(h, Ω_k, Ω_Λ, Ω_m, Ω_r, w0, wa)
    elseif Ω_k > 0
        OpenWCDM(h, Ω_k, Ω_Λ, Ω_m, Ω_r, w0, wa)
    else
        FlatWCDM(h, Ω_k, Ω_Λ, Ω_m, Ω_r, w0, wa)
    end
end

function a2E(c::Union{FlatWCDM,ClosedWCDM,OpenWCDM}, a)
    ade = exp((1 - 3*(c.w0 + c.wa))*log(a) + 3*c.wa*(a - 1))
    sqrt(c.Ω_r + (c.Ω_m + c.Ω_k*a)*a + c.Ω_Λ*ade)
end

function cosmology(;h=0.69,
                   Neff=3.04,
                   OmegaK=0,
                   OmegaM=0.29,
                   OmegaR=nothing,
                   Tcmb=2.7255,
                   w0=-1,
                   wa=0)

    if OmegaR === nothing
        OmegaG = 4.48131e-7*Tcmb^4/h^2
        OmegaN = Neff*OmegaG*(7/8)*(4/11)^(4/3)
        OmegaR = OmegaG + OmegaN
    end

    OmegaL = 1 - OmegaK - OmegaM - OmegaR

    if !(w0 == -1 && wa == 0)
        return WCDM(h, OmegaK, OmegaL, OmegaM, OmegaR, w0, wa)
    end

    if OmegaK < 0
        return ClosedLCDM(h, OmegaK, OmegaL, OmegaM, OmegaR)
    elseif OmegaK > 0
        return OpenLCDM(h, OmegaK, OmegaL, OmegaM, OmegaR)
    else
        return FlatLCDM(h, OmegaL, OmegaM, OmegaR)
    end
end

# hubble rate

scale_factor(z) = 1/(1 + z)
E(c::AbstractCosmology, z) = (a = scale_factor(z); a2E(c,a)/a^2)
H(c::AbstractCosmology, z) = 100 * c.h * E(c, z) * km / s / Mpc

hubble_dist0(c::AbstractCosmology) = 2997.92458/c.h * Mpc
hubble_dist(c::AbstractCosmology, z) = hubble_dist0(c)/E(c,z)

hubble_time0(c::AbstractCosmology) = 9.77814/c.h * Gyr
hubble_time(c::AbstractCosmology, z) = hubble_time0(c)/E(c,z)

# distances

Z(c::AbstractCosmology, z::Real; kws...) =
    QuadGK.quadgk(a->1/a2E(c,a), scale_factor(z), 1; kws...)[1]
Z(c::AbstractCosmology, z₁::Real, z₂::Real; kws...) =
    QuadGK.quadgk(a->1/a2E(c,a), scale_factor(z₂), scale_factor(z₁); kws...)[1]

comoving_radial_dist(c::AbstractCosmology, z; kws...) = hubble_dist0(c)*Z(c, z; kws...)
comoving_radial_dist(c::AbstractCosmology, z₁, z₂; kws...) = hubble_dist0(c)*Z(c, z₁, z₂; kws...)

comoving_transverse_dist(c::AbstractFlatCosmology, z; kws...) =
    comoving_radial_dist(c, z; kws...)
comoving_transverse_dist(c::AbstractFlatCosmology, z₁, z₂; kws...) =
    comoving_radial_dist(c, z₁, z₂; kws...)
function comoving_transverse_dist(c::AbstractOpenCosmology, z; kws...)
    sqrtok = sqrt(c.Ω_k)
    hubble_dist0(c)*sinh(sqrtok*Z(c, z; kws...))/sqrtok
end
function comoving_transverse_dist(c::AbstractOpenCosmology, z₁, z₂; kws...)
    sqrtok = sqrt(c.Ω_k)
    hubble_dist0(c)*sinh(sqrtok*Z(c, z₁, z₂; kws...))/sqrtok
end
function comoving_transverse_dist(c::AbstractClosedCosmology, z; kws...)
    sqrtok = sqrt(abs(c.Ω_k))
    hubble_dist0(c)*sin(sqrtok*Z(c,z; kws...))/sqrtok
end
function comoving_transverse_dist(c::AbstractClosedCosmology, z₁, z₂; kws...)
    sqrtok = sqrt(abs(c.Ω_k))
    hubble_dist0(c)*sin(sqrtok*Z(c, z₁, z₂; kws...))/sqrtok
end

angular_diameter_dist(c::AbstractCosmology, z; kws...) =
    comoving_transverse_dist(c, z; kws...)/(1 + z)
angular_diameter_dist(c::AbstractCosmology, z₁, z₂; kws...) =
    comoving_transverse_dist(c, z₁, z₂; kws...)/(1 + z₂)

luminosity_dist(c::AbstractCosmology, z; kws...) =
    comoving_transverse_dist(c, z; kws...)*(1 + z)

distmod(c::AbstractCosmology, z; kws...) =
    5 * log10(luminosity_dist(c, z; kws...) / Mpc) + 25

# volumes

comoving_volume(c::AbstractFlatCosmology, z; kws...) =
    (4pi/3)*(comoving_radial_dist(Gpc, c, z; kws...))^3
function comoving_volume(c::AbstractOpenCosmology, z; kws...)
    DH = hubble_dist0(Gpc, c)
    x = comoving_transverse_dist(Gpc, c, z; kws...)/DH
    sqrtok = sqrt(c.Ω_k)
    2pi*(DH)^3*(x*sqrt(1 + c.Ω_k*x^2) - asinh(sqrtok*x)/sqrtok)/c.Ω_k
end
function comoving_volume(c::AbstractClosedCosmology, z; kws...)
    DH = hubble_dist0(Gpc, c)
    x = comoving_transverse_dist(Gpc, c,z; kws...)/DH
    sqrtok = sqrt(abs(c.Ω_k))
    2pi*(DH)^3*(x*sqrt(1 + c.Ω_k*x^2) - asin(sqrtok*x)/sqrtok)/c.Ω_k
end

comoving_volume_element(c::AbstractCosmology, z; kws...) =
    hubble_dist0(Gpc, c)*angular_diameter_dist(Gpc, c, z; kws...)^2/a2E(c, scale_factor(z))

# times

T(c::AbstractCosmology, a0, a1; kws...) = QuadGK.quadgk(x->x/a2E(c,x), a0, a1; kws...)[1]
age(c::AbstractCosmology, z; kws...) = hubble_time0(c)*T(c, 0, scale_factor(z); kws...)
lookback_time(c::AbstractCosmology, z; kws...) = hubble_time0(c)*T(c, scale_factor(z), 1; kws...)

# Easily select a different unit
for f in (:hubble_dist0, :hubble_dist, :hubble_time0, :hubble_time, :comoving_radial_dist,
          :comoving_transverse_dist, :angular_diameter_dist, :luminosity_dist,
          :comoving_volume, :comoving_volume_element, :age, :lookback_time)
    @eval $f(u::Unitful.Unitlike, args...; kws...) = uconvert(u, $f(args...; kws...))
end

###############
# Deprecations
#

@deprecate hubble_dist_mpc0(c::AbstractCosmology) ustrip(hubble_dist0(c::AbstractCosmology))
@deprecate hubble_dist_mpc(c::AbstractCosmology, z) ustrip(hubble_dist(c::AbstractCosmology, z))
@deprecate hubble_time_gyr0(c::AbstractCosmology) ustrip(hubble_time0(c::AbstractCosmology))
@deprecate hubble_time_gyr(c::AbstractCosmology, z) ustrip(hubble_time(c::AbstractCosmology, z))
@deprecate comoving_radial_dist_mpc(c::AbstractCosmology, z; kws...) ustrip(comoving_radial_dist(c::AbstractCosmology, z; kws...))
@deprecate comoving_transverse_dist_mpc(c::AbstractCosmology, z; kws...) ustrip(comoving_transverse_dist(c::AbstractCosmology, z; kws...))
@deprecate angular_diameter_dist_mpc(c::AbstractCosmology, z; kws...) ustrip(angular_diameter_dist(c::AbstractCosmology, z; kws...))
@deprecate luminosity_dist_mpc(c::AbstractCosmology, z; kws...) ustrip(luminosity_dist(c::AbstractCosmology, z; kws...))
@deprecate comoving_volume_gpc3(c::AbstractCosmology, z; kws...) ustrip(comoving_volume(c::AbstractCosmology, z; kws...))
@deprecate comoving_volume_element_gpc3(c::AbstractCosmology, z; kws...) ustrip(comoving_volume_element(c::AbstractCosmology, z; kws...))
@deprecate age_gyr(c::AbstractCosmology, z; kws...) ustrip(age(c::AbstractCosmology, z; kws...))
@deprecate lookback_time_gyr(c::AbstractCosmology, z; kws...) ustrip(lookback_time(c::AbstractCosmology, z; kws...))

end # module

ext: jl


text:
import UUIDs: UUID, uuid1

mutable struct Notebook
    path::String

    "Cells are ordered in a `Notebook`, and this order can be changed by the user. Cells will always have a constant UUID."
    cells::Array{Cell,1}

    notebook_id::UUID
    combined_funcdefs::Dict{Vector{Symbol}, SymbolsState}

    # buffer will contain all unfetched updates - must be big enough
    pendingupdates::Channel

    executetoken::Channel
end
# We can keep 128 updates pending. After this, any put! calls (i.e. calls that push an update to the notebook) will simply block, which is fine.
# This does mean that the Notebook can't be used if nothing is clearing the update channel.
Notebook(path::String, cells::Array{Cell,1}, notebooID) = let
    et = Channel{Nothing}(1)
    put!(et, nothing)
    Notebook(path, cells, notebooID, Dict{Vector{Symbol}, SymbolsState}(), Channel(1024), et)
end
Notebook(path::String, cells::Array{Cell,1}) = Notebook(path, cells, uuid1())
Notebook(cells::Array{Cell,1}) = Notebook(tempname() * ".jl", cells)

function cellindex_fromID(notebook::Notebook, cell_id::UUID)::Union{Int,Nothing}
    findfirst(c->c.cell_id == cell_id, notebook.cells)
end

# We use a creative delimiter to avoid accidental use in code
const _cell_id_delimiter = "# ╔═╡ "
const _order_delimiter = "# ╠═"
const _order_delimiter_folded = "# ╟─"
const _cell_suffix = "\n\n"

emptynotebook(path) = Notebook(path, [Cell("")])
emptynotebook() = emptynotebook(tempname() * ".jl")

function save_notebook(io, notebook::Notebook)
    println(io, "### A Pluto.jl notebook ###")
    println(io, "# ", PLUTO_VERSION_STR)
    # Anything between the version string and the first UUID delimiter will be ignored by the notebook loader.
    println(io, "")
    println(io, "using Markdown")
    # Super Advanced Code Analysis™ to add the @bind macro to the saved file if it's used somewhere.
    if any(occursin("@bind", c.code) for c in notebook.cells)
        println(io, PlutoRunner.fake_bind)
    end
    println(io)

    # TODO: this can be optimised by caching the topological order:
    # maintain cache with ordered UUIDs
    # whenever a run_reactive is done, move the found cells **up** until they are in one group, and order them topologcally within that group. Errable cells go to the bottom.

    # the next call took 2ms for a small-medium sized notebook: (so not too bad)
    celltopology = topological_order(notebook, notebook.cells)

    cells_ordered = union(celltopology.runnable, keys(celltopology.errable))

    for c in cells_ordered
        println(io, _cell_id_delimiter, string(c.cell_id))
        print(io, c.code)
        print(io, _cell_suffix)
    end

    println(io, _cell_id_delimiter, "Cell order:")
    for c in notebook.cells
        delim = c.code_folded ? _order_delimiter_folded : _order_delimiter
        println(io, delim, string(c.cell_id))
    end
end

function save_notebook(notebook::Notebook, path::String)
    open(path, "w") do io
        save_notebook(io, notebook)
    end
end

save_notebook(notebook::Notebook) = save_notebook(notebook, notebook.path)

"Load a notebook without saving it or creating a backup; returns a `Notebook`. REMEMBER TO CHANGE THE NOTEBOOK PATH after loading it to prevent it from autosaving and overwriting the original file."
function load_notebook_nobackup(io, path)::Notebook
    firstline = String(readline(io))

    if firstline != "### A Pluto.jl notebook ###"
        @error "File is not a Pluto.jl notebook"
    end

    file_VERSION_STR = readline(io)[3:end]
    if file_VERSION_STR != PLUTO_VERSION_STR
        # @info "Loading a notebook saved with Pluto $(file_VERSION_STR). This is Pluto $(PLUTO_VERSION_STR)."
    end

    collected_cells = Dict()

    # ignore first bits of file
    readuntil(io, _cell_id_delimiter)

    last_read = ""
    while !eof(io)
        cell_id_str = String(readline(io))
        if cell_id_str == "Cell order:"
            break
        else
            cell_id = UUID(cell_id_str)
            code_raw = String(readuntil(io, _cell_id_delimiter))
            # change Windows line endings to Linux
            code_normalised = replace(code_raw, "\r\n" => "\n")
            # remove the cell appendix
            code = code_normalised[1 : prevind(code_normalised, end, length(_cell_suffix))]

            read_cell = Cell(cell_id, code)
            collected_cells[cell_id] = read_cell
        end
    end

    ordered_cells = Cell[]
    while !eof(io)
        cell_id_str = String(readline(io))
        o, c = startswith(cell_id_str, _order_delimiter),
        if length(cell_id_str) >= 36
            cell_id = let
                UUID(cell_id_str[end-35:end])
            end
            next_cell = collected_cells[cell_id]
            next_cell.code_folded = startswith(cell_id_str, _order_delimiter_folded)
            push!(ordered_cells, next_cell)
        end
    end

    Notebook(path, ordered_cells)
end

function load_notebook_nobackup(path::String)::Notebook
    local loaded
    open(path, "r") do io
        loaded = load_notebook_nobackup(io, path)
    end
    loaded
end

"Create a backup of the given file, load the file as a .jl Pluto notebook, save the loaded notebook, compare the two files, and delete the backup of the newly saved file is equal to the backup."
function load_notebook(path::String)::Notebook
    local backupNum = 1
    backupPath = path
    while isfile(backupPath)
        backupPath = path * ".backup" * string(backupNum)
        backupNum += 1
    end
    cp(path, backupPath)

    loaded = load_notebook_nobackup(path)
    # Analyze cells so that the initial save is in topological order
    update_caches!(loaded, loaded.cells)
    save_notebook(loaded)
    # Clear symstates if autorun/autofun is disabled. Otherwise running a single cell for the first time will also run downstream cells.
    if ENV["PLUTO_RUN_NOTEBOOK_ON_LOAD"] != "true"
        for cell in loaded.cells
            cell.symstate = SymbolsState()
        end
    end

    if only_versions_or_lineorder_differ(path, backupPath)
        rm(backupPath)
    else
        @warn "Old Pluto notebook might not have loaded correctly. Backup saved to: " backupPath
    end

    loaded
end

function move_notebook(notebook::Notebook, newpath::String)
    # Will throw exception and return if anything goes wrong, so at least one file is guaranteed to exist.
    oldpath = notebook.path
    save_notebook(notebook, oldpath)
    save_notebook(notebook, newpath)
    @assert only_versions_differ(oldpath, newpath)
    notebook.path = newpath
    rm(oldpath)
end

"Check if two savefiles are identical, up to their version numbers and a possible line shuffle.

If a notebook has not yet had all of its cells run, we can't deduce the topological cell order."
function only_versions_or_lineorder_differ(pathA::AbstractString, pathB::AbstractString)::Bool
    Set(readlines(pathA)[3:end]) == Set(readlines(pathB)[3:end])
end

function only_versions_differ(pathA::AbstractString, pathB::AbstractString)::Bool
    readlines(pathA)[3:end] == readlines(pathB)[3:end]
end

function tryexpanduser(path)
	try
		expanduser(path)
	catch ex
		path
	end
end

tamepath = abspath ∘ tryexpanduser


text:
using SuperStiff.PeriodizeSC
using SuperStiff.Stiffness

macro assertion(ex, text)
    :($ex ? nothing : error("Assertion failed: ", $(text)))
end

# Reading parameters from JSON file
paramsfile = "params.json"
params = JSON.parsefile(paramsfile)
println(typeof(params))

# Defining own dictionnary for convenience
block_params = Dict{String,Any}(
	"path_to_files" => params["path_to_files"],
	"data_loop" => params["data_loop"],
	"pattern" => params["pattern"], # Pattern of the files to be loaded
	"beta" => params["beta"],
	"print_mu_dop" => params["Print_mu_dop"], # Parameter used to create stiffness.dat file (file used with SE_G_fig_producer.py)
	"w_discretization" => params["w_discretization"], # Matsubara frequency grid
	"AFM_SC_NOCOEX" => params["AFM_SC_NOCOEX"],
	"cumulant" => params["cumulant"],
	"Periodization" => params["Periodization"],
	"fout_name" => params["fout_name"]
	)

# Check if appropriate input was entered
if block_params["cumulant"] in [0,1]; nothing else; throw(ErrorException("Cumulant takes in 0 or 1 as arguments")) end
if block_params["Periodization"] in [0,1]; nothing else; throw(ErrorException("Periodization takes in 0 or 1 as arguments")) end
if block_params["AFM_SC_NOCOEX"] in [0,1]; nothing else; throw(ErrorException("AFM_SC_NOCOEX takes in 0 or 1 as arguments")) end

# Check out if appriopriate inputs entered (Debugging)
match_data_loop = match(r"^(.*?)(?=/)",block_params["data_loop"]).match
match_data_loop_arr = Array{String,1}([match_data_loop])
M_tol = -0.01

# Decide the name of the file to be created
fout_name = block_params["fout_name"]

super_data_M, super_datap, list_of_files = Stiffness.read_data_loop(block_params,block_params["print_mu_dop"],block_params["pattern"])

@assertion(length(list_of_files) != 0, "Length of list_of_files must differ from 0. Check if the path_to_files variable has been given the right path!")
@assertion(length(super_datap) != 0, "Length of super_datap must differ from 0. Check your file containing the chemical potential and the doping and make sure they are respectively marked \'mu\' and \'ave_mu\' in the header")

if Stiffness.verbose > 0
    println("Length list_of_files is : ", length(list_of_files),"\n")
    println("Size super_datap is: ", size(super_datap),"\n")
end

length(list_of_files) != length(super_datap[:,1]) && throw(Stiffness.File_Input_Error(:list_file))
list_of_files_and_datap = hcat(list_of_files,super_datap)

# Building list_modulevec defining the instance of modulevec
super_list_modulevec = Stiffness.gen_modulevec_args(block_params,list_of_files,super_datap[:,1],params)

if Stiffness.verbose > 0
    println("typeof list_modulevec: ",typeof(super_list_modulevec))
end

println(Stiffness.ts(super_list_modulevec[1]))
if Stiffness.verbose > 0
    println("Size of sEvec_c in StiffnessArray: ", size(super_list_modulevec[1].data_), " | typeof StiffnessArray: ", typeof(super_list_modulevec[1]))
end

########################################### Beginning of the superfluid stiffness calculations ###########################################
# Computing the superfluid stiffness
#@time list_stiff, list_kIntegral_DOS, list_kIntegral_DOS_nocoex = Stiffness.calc_stiff_funct(super_data_M, super_list_modulevec, match_data_loop_arr, M_tol)

beta = super_list_modulevec[1].vals_["beta"]
println("beta = $(@sprintf("%.2f", beta))")
w_discretization = super_list_modulevec[1].vals_["w_discretization"]
zvec = 1.0im*[(2*n+1)*pi/beta for n in 1:w_discretization]

list_stiff = Array{Complex{Float64},1}()
@time for (l,modulevec_el) in enumerate(super_list_modulevec)
    arr_w_kIntegrated = Matrix{Complex{Float64}}(0,0)
    println("Iteration number: ",l)
    if match_data_loop_arr[1] == "COEX"
        super_data_M_el = super_data_M[l]
    end
    println("Treating $(match_data_loop_arr[1]) data")
    zvec = 1.0im*[(2*n+1)*pi/beta for n in 1:w_discretization]
    modelvec = PeriodizeSC.ModelVector(Stiffness.ts(modulevec_el)[1],Stiffness.ts(modulevec_el)[2],Stiffness.ts(modulevec_el)[3],Stiffness.ts(modulevec_el)[4],zvec[1:w_discretization], modulevec_el.data_[1,1:w_discretization,:,:])
    arr_w_kIntegrated = Matrix{Complex{Float64}}(0,0)
    if match_data_loop_arr[1] == "COEX"
        @time arr_w_kIntegrated = Stiffness.calc_stiff_funct_COEX(super_data_M_el, modelvec, modulevec_el, M_tol)
    elseif match_data_loop_arr[1] == "NOCOEX"
        @time arr_w_kIntegrated = Stiffness.calc_stiff_funct_NOCOEX(modelvec, modulevec_el)
    end
    dop_stiff = 2.0/beta*sum(arr_w_kIntegrated[:,2])
    f = open(block_params["fout_name"], "a")
    if l == 1
        write(f, "#M_tol "*"$(M_tol)"*"\n")
    end
    write(f, "$(super_datap[:,1][l])"*"\t\t"*"$(super_datap[:,2][l])"*"\t\t"*"$(real(dop_stiff))"*"\n")
    close(f)
    push!(list_stiff,dop_stiff)
    println(dop_stiff)
end
println("Superfluid Stiffness: ", list_stiff)

println("Program exited")

max_stars_repo_name: oliviersimard/SuperStiff

max_issues_repo_name: oliviersimard/SuperStiff

max_forks_repo_name: oliviersimard/SuperStiff


text:
abstract type Material end
hard(m::Material)=mat数据[typeof(m).name.name][:hard]
function ltheme(m::Material)
	lt=mat数据[typeof(m).name.name][:lt]
	return RGB{N0f8}(lt[1],lt[2],lt[3])
end
function theme(m::Material)
	t=mat数据[typeof(m).name.name][:lt]
	return RGB{N0f8}(t[1],t[2],t[3])
end
function rtheme(m::Material)
	rt=mat数据[typeof(m).name.name][:lt]
	return RGB{N0f8}(rt[1],rt[2],rt[3])
end
macro material_types(s::Symbol...)
	for i in s
		eval(:(struct $i <: Material end))
	end
end
@material_types 苍穹木 石 粘土 陶瓷 铁 金 铜


text:
# Monetary policy switches
const dmp_switches = (:dmpintay, :dmptay, :dmptlr, :dmpalt, :dmpgen, :dmpex, :dmprr)

# Fiscal policy switches
const dfp_switches = (:dfpex, :dfpsrp, :dfpdbt)

"""
    set_mp!(data, switch[, range])

Set the monetary policy function over the given range. The given switch is set
to 1 and all other monetary policy switches are set to 0. If range is not given,
the switch is set over the full range of the data.

The set of valid switch values is in `dmp_switches`.
"""
function set_mp!(data::SimData, switch, range=mitrange(data))
    sw = Symbol(switch)
    if sw ∉ dmp_switches
        throw(ArgumentError("Not a valid Monetary Policy switch: $switch"))
    end
    data[range, dmp_switches] = 0.0
    data[range, sw] = 1.0
    return data
end

"""
    set_fp!(data, switch[, range])

Set the fiscal policy function over the given range. The given switch is set
to 1 and all other fiscal policy switches are set to 0. If range is not given,
the switch is set over the full range of the data.

The set of valid switch values is in `dfp_switches`.
"""
function set_fp!(data::SimData, switch, range=mitrange(data))
    sw = Symbol(switch)
    if sw ∉ dfp_switches
        throw(ArgumentError("Not a valid Fiscal Policy switch: $switch"))
    end
    data[range, dfp_switches] = 0.0
    data[range, sw] = 1.0
    return data
end

nothing


text:
using .MeasureTheory
using ..TupleVectors

function Base.rand(ω::Hypercube, d::ParameterizedMeasure)
    Dists.quantile(distproxy(d), rand(ω))
end

function Base.rand(ω::Hypercube, d::ProductMeasure)
    [rand(ω, dj) for dj in d.data]
end


text:
#=
This file contains code for calculating intervals for plotting purposes.
These functions should, at minimum, take in some form of data input and
return a tuple (min, max) of the limits corresponding to the function
=#

using Statistics: mean, median, std

"""
    zscale(input::AbstractArray,
        nsamples=1000;
        contrast=0.25,
        max_reject=0.5,
        min_npixels=5,
        k_rej=2.5,
        max_iterations=5)

Implementation of the `zscale` IRAF function for finding
colorbar limits of `input`, which showcase data near the median.
This is useful for data with extreme outliers, such as astronomical images.

## Keyword arguments
* `nsamples` - The number of samples to use from `input`. If fewer than `nsamples` are present, will use the full input
* `contrast` - The desired contrast
* `k_rej` - The number of standard deviations above which data is rejected
* `max_iteration` - The number of iterations used for fitting samples
* `max_reject` - The maximum number of pixels to reject during the iterative fitting
* `min_npixels` - The minimum number of pixels to calculate the limits after the iterative fitting

See the extended help (in REPL, `??zscale`) for technical details.

## Examples
```jldoctest
julia> img = 0:9999

julia> zscale(img)
(0, 9990)
```

## Extended help
## Description

The zscale algorithm is designed to display the image values
near the median image value, without the time consuming process
of computing a full image histogram.

This is particularly useful for astronomical images, which generally
have a very peaked histogram corresponding to the background sky
in direct imaging, or the continuum in a two dimensional spectrum.

A subset of the image is examined, and approximately `nsamples`
pixels are sampled evenly over the image. The number of lines is a
user parameter, `nsample_lines`.

The pixels are ranked in brightness to form the function `I(i)`,
where `i` is the rank of the pixel, and `I` is its value. Generally,
the midpoint of this function (the median) is very near the peak
of the image histogram.  There is a well defined slope about the
midpoint, which is related to the width of the histogram.

At the ends of the `I(i)` function, there are a few very bright and
dark pixels due to objects and defects in the field. To determine
the slope, a linear function is fit with iterative rejection:

```
I(i) = intercept + slope * (i - midpoint)
```

If more than half of the points are rejected, then there is no well
defined slope, and the full range of the sample defines `z1` and
`z2`. Otherwise, the endpoints of the linear function are used
(provided they are within the original range of the sample):

```
z1 = I(midpoint) + (slope / contrast) * (1 - midpoint)
z2 = I(midpoint) + (slope / contrast) * (npoints - midpoint)
```
"""
function zscale(
    input::AbstractArray,
    nsamples::Int = 1000;
    contrast = 0.25,
    max_reject = 0.5,
    min_npixels = 5,
    k_rej = 2.5,
    max_iterations = 5,
)

    # get samples from finite values of input
    values = float(filter(isfinite, input))
    stride = max(1, round(Int, length(values) / nsamples))
    samples = @view values[1:stride:end][1:min(nsamples, end)]
    sort!(samples)

    N = length(samples)
    vmin = first(samples)
    vmax = last(samples)

    # fit a line to the sorted samples
    min_pix = max(min_npixels, round(Int, N * max_reject))
    x = 0:(N - 1)

    ngood = N
    last_good = N + 1

    # good pixel mask (array bool faster than bitarray)
    mask = ones(Bool, N)

    # get number of neighbors to mask if a pixel is bad
    ngrow = max(1, round(Int, N / 100) ÷ 2)

    local β # line slope
    # iteratively fit samples and reject sigma-clipped outliers
    for _ in 1:max_iterations
        (ngood ≥ last_good || ngood < min_pix) && break

        # linear fit using mask
        x_ = @view x[mask]
        y = @view samples[mask]
        α, β = fit_line(x_, y)
        flat = @. samples - (α + β * x)

        # k-sigma rejection threshold
        threshold = k_rej * std(flat[mask])

        # detect and reject outliers based on threshold
        @. mask[!(-threshold ≤ flat ≤ threshold)] = false

        # dilate mask
        mask = dilate_mask(mask, ngrow)

        last_good = ngood
        ngood = count(mask)
    end

    if ngood ≥ min_pix
        slope = contrast > 0 ? β / contrast : β
        center = (N - 1) ÷ 2
        m = median(samples)
        vmin = max(vmin, m - (center - 1) * slope)
        vmax = min(vmax, m + (N - center) * slope)
    end

    return vmin, vmax
end

"""
    dilate_mask(mask, ngrow)

Takes a mask and dilates each `false` with `ngrow` `false`s on
either side.  This is equivalent to boolean "convolution".
"""
function dilate_mask(mask, ngrow)
    out = similar(mask)
    idxs = CartesianIndices(mask)
    mindx = idxs[1].I[1]
    maxdx = idxs[end].I[1]
    @inbounds for idx in idxs
        lower = max(mindx, idx.I[1] - ngrow)
        upper = min(maxdx, idx.I[1] + ngrow)
        # output will only be true if there are no falses in the input section
        out[idx] = all(mask[lower:upper])
    end
    return out
end

"""
    fit_line(x, y)

Simple linear regression, returns `(intercept, slope)`.
https://en.wikipedia.org/wiki/Simple_linear_regression
"""
function fit_line(x, y)
    mx = mean(x)
    my = mean(y)
    SSxy = sum(zip(x, y)) do (xi, yi)
        (xi - mx) * (yi - my)
    end
    SSx = sum(x) do xi
        (xi - mx)^2
    end
    β = SSxy / SSx
    α = my - β * mx
    return α, β
end


text:
using Bifurcations
using Bifurcations.Examples: Calcium

solver = init(Calcium.prob)
@time solve!(solver)
sol = solver.sol
display(sol)

using Plots
plot = Bifurcations.plot  # workaround
plt1 = plot(sol)
display(plt1)

using Bifurcations: BifurcationProblem, special_intervals
using Setfield: @lens

point_list = sort!(special_intervals(solver), by=p->p.u0[end])
point = point_list[1]
@show point

codim2_prob = BifurcationProblem(
    point,
    solver,
    (@lens _.gca),
    (0.0, 8.0),
)
codim2_solver = init(codim2_prob)
solve!(codim2_solver)

plt2 = plot(codim2_solver.sol)
display(plt2)


text:
module SchwarzOpt

using Printf
using DataStructures
using LinearAlgebra
import Base.@kwdef
# import LinearAlgebra: BLAS

using MathOptInterface, Plasmo, JuMP

include("utils.jl")

include("optimizer.jl")

include("moi.jl")

end


text:
const dot_signs = [0b10 0b01]

const dot_decode = Array{Char}(undef, 5)
dot_decode[0b00 + 1] = ' '
dot_decode[0b01 + 1] = '.'
dot_decode[0b10 + 1] = '\''
dot_decode[0b11 + 1] = ':'

"""
Similar to the `AsciiCanvas`, the `DotCanvas` only uses
ASCII characters to draw it's content. Naturally,
it doesn't look quite as nice as the Unicode-based
ones. However, in some situations it might yield
better results. Printing plots to a file is one
of those situations.

The DotCanvas is best utilized in combination
with `scatterplot`.
For `lineplot` we suggest to use the `AsciiCanvas`
instead.
"""
struct DotCanvas <: LookupCanvas
    grid::Array{UInt8,2}
    colors::Array{UInt8,2}
    pixel_width::Int
    pixel_height::Int
    origin_x::Float64
    origin_y::Float64
    width::Float64
    height::Float64
end

@inline x_pixel_per_char(::Type{DotCanvas}) = 1
@inline y_pixel_per_char(::Type{DotCanvas}) = 2

@inline lookup_encode(c::DotCanvas) = dot_signs
@inline lookup_decode(c::DotCanvas) = dot_decode

function DotCanvas(args...; nargs...)
    CreateLookupCanvas(DotCanvas, args...; nargs...)
end


text:
module TestUnordered

import SplittablesBase

raw_examples = dict_examples = """
Dict{Symbol,Int}()
Dict(:a => 1)
Dict(:a => 1, :b => 2)
Dict(:a => 1, :b => 2, :c => 3)
Dict(:a => 1, :b => 2, :c => 3, :d => 4)
Dict(:a => 1, :b => 2, :c => 3, :d => 4, :e => 5)
Dict(zip('a':'z', 1:26))
"""

function wrapdict(f)
    lines = map(split(dict_examples, "\n", keepempty = false)) do code
        "$f($code)"
    end
    return join(lines, "\n") * "\n"
end

raw_examples *= wrapdict("keys") * wrapdict("values")

raw_examples *= """
Set(1:0)
Set(1:1)
Set(1:2)
Set(1:3)
Set(1:4)
Set(1:5)
Set(1:100)
"""

# An array of `(label = ..., data = ...)`
examples = map(split(raw_examples, "\n", keepempty = false)) do code
    (label = code, data = Base.include_string(@__MODULE__, code))
end

SplittablesBase.Testing.test_unordered(examples)

end  # module


text:
#############################################################################
#Copyright (c) 2016 Charles Le Losq
#
#The MIT License (MIT)
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the #Software without restriction, including without limitation the rights to use, copy, #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, #and to permit persons to whom the Software is furnished to do so, subject to the #following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#
# diffusion.jl contains several functions focused on the treatment of infrared spectra along diffusion profiles in minerals
#
#
#############################################################################
"""
The 'IRdataprep' function allows to obtain the frequency, distance and y vectors in a region of interest, ready to input in an optimisation algorithm.

   IRdataprep(data,distance_step,start_sp,stop_sp,low_x,high_x,norm_low_x,norm_high_x)
with
data (array, Float64): an array containing in the first column the x values and in the subsequent columns the correcponding intensities;
distance_step (Float64): the steps between each spectrum along the diffusion profile, in metres;
start_sp and stop_sp (Int): the starting and stopping spectra of interest (warning: number of data columns - 1, as data also contains the x axis such that size(data)[2] = number of spectras + 1);
low_x,high_x (Float64): frequencies between which the signal is of interest, defined as low_x < high_x;
norm_low_x and norm_high_x (Float64): frequencies between which the signal must be integrated for mormalising the final signals.

This function first checks that the x values are increasing, and correct that if this is not the case.

Then, it selects the data in the region of interest, delimited by frequencies defined as low_x < high_x.
Normalisation to the area between the norm_low_x and norm_high_x frequency is also performed (norm_low_x < norm_high_x)

The code returns the x and the y arrays, plus arrays for optimisation: x_input_fit contains the frequency-distance couples and y_input_fit their corresponding y values.
"""
function IRdataprep(data::Array{Float64},distance_step::Float64,spectra_numbers::Array{Int64}, portion_interest::Array{Float64}, normalization_region::Array{Float64},integration_regions::Array{Float64},polarizer_orientation::AbstractString;fig_select::Int64=1,bas_switch::AbstractString = "No",roi::Array{Float64} = [3000. 3100.; 3425. 3530.; 3630. 3650.],baseline_type = "gcvspline", smoothing_coef = 10.0, Spline_Order = 3,noise_estimation = "No", noise_calculation_portion::Array{Float64} = [3000. 3100.])

    start_sp::Int = spectra_numbers[1]
    stop_sp::Int = spectra_numbers[2]
    low_x::Float64 = portion_interest[1]
    high_x::Float64 = portion_interest[2]
    norm_low_x::Float64 = normalization_region[1]
    norm_high_x::Float64 = normalization_region[2]
    ese_low_x::Float64 = noise_calculation_portion[1]
    ese_high_x::Float64 = noise_calculation_portion[2]

    if data[end,1] <= data[1,1]
        data=data[end:-1:1,:]
        data[data[:,:] .== 0]= 1e-20# to avoid pure 0 values
    end

    x::Array{Float64} = data[find(low_x .< data[:,1] .< high_x),1]
    x_sili::Array{Float64} = data[find(norm_low_x .< data[:,1] .< norm_high_x),1]

    if stop_sp .>= size(data,2)
        y::Array{Float64} = data[find(low_x .< data[:,1] .< high_x),start_sp:end]
        y_sili::Array{Float64} = data[find(norm_low_x .< data[:,1] .< norm_high_x),start_sp:end]
        if noise_estimation == "Yes"
            ese::Array{Float64} =  std(data[find(ese_low_x .< data[:,1] .< ese_high_x),start_sp:end],1) #relative error
            y_ese_r::Array{Float64} = ese[1,:]./y[:,:]
        end
    else
        y = data[find(low_x .< data[:,1] .< high_x),start_sp:stop_sp]
        y_sili = data[find(norm_low_x .< data[:,1] .< norm_high_x),start_sp:stop_sp]
		if noise_estimation == "Yes"
            ese =  std(data[find(ese_low_x .< data[:,1] .< ese_high_x),start_sp:stop_sp],1) #relative error
        	y_ese_r = ese[1,:]./y[:,:]
        end

    end

	nb_exp = size(y,2)
	nb_points = size(y,1)

   # integration of silicate bands for normalisation depending on the polariser orientation, we use different thickness determination
    if polarizer_orientation == "A"
        K1 = 0.8345
    elseif polarizer_orientation == "C"
        K1 = 0.5620
    end

    for i = 1:nb_exp
		y_sili_corr,~ = baseline(x_sili[:,1],y_sili[:,i],[norm_low_x norm_low_x+5;norm_high_x-5 norm_high_x],"poly",p=1.0)
        y[:,i] = y[:,i]./(trapz(x_sili[:,1],y_sili_corr[:,1])./K1.*1e-4) # integration of silicate bands for normalisation, thickness is in cm
    end

    if (bas_switch == "Yes") && (roi[1,1] != 0)
        data_bas = ones(size(y))
        base = ones(size(y))
        for i = 1:nb_exp
			data_bas[:,i], base[:,i] = baseline(x, y[:,i],roi,baseline_type,p=smoothing_coef,SplOrder=Spline_Order)
			#for k = 1:size(roi,1)-1
			#	data_bas[roi[k,1] .<= x .<= roi[k+1,2],i], base[roi[k,1] .<= x .<= roi[k+1,2],i] = baseline(x[roi[k,1] .<= x .<= roi[k+1,2]], y[roi[k,1] .<= x .<= roi[k+1,2],i],roi[k:k+1,:],baseline_type,p=smoothing_coef)
			#end
        end

        figure()
        plot(x,y[:,fig_select],"k-",label="Data")
        plot(x,base[:,fig_select+1],"b-",label="Baseline")
        y = data_bas
        plot(x,y[:,fig_select],"g-",label="Corrected")
		legend(loc="best")
    end

   # constructing the good distance vector + X and Y associated values
    steps::Float64 = 0.0
	integrale_matrix::Array{Float64} = zeros(nb_exp,size(integration_regions,1))
	distance_vector::Array{Float64} = zeros(nb_exp)
    y_interest::Array{Float64} = zeros(nb_points,nb_exp)
    ese_y_interest::Array{Float64} = zeros(length(x),1)

    for i = 1:nb_exp

		distance_vector[i] = steps # individual distance vector
		for k = 1:size(integration_regions,1)
			integrale_matrix[i,k] = trapz(x[integration_regions[k,1] .<= x .<= integration_regions[k,2]],y[integration_regions[k,1] .<= x .<= integration_regions[k,2],i])
		end

        steps = steps + distance_step
    end

	if noise_estimation == "Yes"
        #return x, y, y_ese_r, distance_vector, integrale_matrix
    else
        return x, y, distance_vector, integrale_matrix
	end
end

"""
The `peak_diffusion` function allows to fit/plot gaussian peaks for which the amplitude is defined by a 1D diffusion law
Call as:
    peak_diffusion(g_c0::Array{Float64,1},g_c1::Array{Float64,1},g_D::Array{Float64,1}, g_freq::Array{Float64,1}, g_hwhm::Array{Float64,1},x::Array{Float64,2},time::Float64)

where:
g_c0 is the vector of concentration imposed, at the diffusive bondary (a.u., float64 number);
g_c1 is the vector of initial concentration in the host (a.u., float64 number);
g_D is the vector of diffusivity coefficient, in log units
g_freq is the vector containing the frequencies fo the peaks
g_hwhm is the vector containing the hwhm of the peaks
x is an array containing the distances and associated frequencies
time is the duration of the run in seconds

The function returns a vector of the values of the peak intensity at the input frequencies.

The current version assumes the peak being a Gaussian peak. Updates will integrate Lorentzian and Pseudo-Voigt peaks.
"""
function peak_diffusion(g_c0::Array{Float64,1},g_c1::Array{Float64,1},g_D::Array{Float64,1}, g_freq::Array{Float64,1}, g_hwhm::Array{Float64,1},x::Array{Float64,2},time::Float64)
    segments = zeros(size(x)[1],size(g_c0)[1])
    for i = 1:size(g_c0)[1]
        segments[:,i] = ((g_c1[i] - g_c0[i]) .* erfc(x[:,1]./(2. .* sqrt( 10.0 .^g_D[i] .* time))) + g_c0[i]) .*exp(-log(2) .* ((x[:,2]-g_freq[i]) ./g_hwhm[i]) .^2)
    end
    return sum(segments,2), segments
end

"""
The global model for using with the Julia.optim package. Not necessarily up to date because I prefer using the JuMP way.
   model(p,x_fit,temps;compo_out = 0)
"""
function model(p::Array{Float64,1},x_fit::Array{Float64,2},temps::Float64;compo_out::Int = 0)
    number_gauss = length(p)/5
    components = zeros(length(x_fit[:,2]),round(Int8,number_gauss))
    compteur = 1
    for i = 1:5:length(p)
        components[:,compteur] = peak_diffusion(x_fit[:,1],temps,x_fit[:,2],p[i],p[i+1],p[i+2],p[i+3],p[i+4])
        compteur = compteur + 1
    end
    if compo_out == 0
        return sum(components,2)
    else
        return components
    end
end


text:
@with_kw struct AbstractKktRhs{T <: AbstractFloat}
    dx::Vector{T}
    dy::Vector{T}
    dz::Vector{T}
    ds::Vector{T}
    dτ::T
    dκ::T
    @assert length(ds) == length(dz)
end

function get_aff_dir!(iter::Iterate{T}, resid::Residual{T}, λ::Vector{T}, prob::Problem{T},
                      kkt::KktCache{T}, scaling::CompositeConeScaling{T},
                      settings::Settings{T}, stats::Stats{T}) where T <: AbstractFloat
    # Form the The RHS of (44a-b), i.e., (dx, dy, dz, dτ, ds, dκ).
    ds = jordan(λ, λ, prob.C)
    dκ = iter.κ * iter.τ

    # Solve KKT.
    rhs = AbstractKktRhs{T}(resid.rx, resid.ry, resid.rz, ds, resid.rτ, dκ)
    aff_dir, ir_iters = solve_abstract_kkt(iter, λ, prob, rhs, kkt, scaling)
    register_aff_dir_stats!(stats, ir_iters)
    return aff_dir
end

function get_comb_dir!(iter::Iterate{T}, resid::Residual{T}, λ::Vector{T}, prob::Problem{T},
                       kkt::KktCache{T}, scaling::CompositeConeScaling{T}, aff_dir::StepDirection{T},
                       settings::Settings{T}, σ::T, stats::Stats{T}) where T <: AbstractFloat
    # These parameters map (38a-b) to (44a-b), by defining (dx, dy, dz, dτ, ds, dκ).
    μ_hat = (iter.s' * iter.z + iter.κ * iter.τ) / (degree(prob.C) + 1)
    dx = (1 - σ) .* resid.rx
    dy = (1 - σ) .* resid.ry
    dz = (1 - σ) .* resid.rz
    dτ = (1 - σ) .* resid.rτ

    ds = (jordan(λ, λ, prob.C)
            + jordan(scaling \ aff_dir.Δs, scaling * aff_dir.Δz, prob.C)
            - σ * μ_hat * get_id(prob.C))
    dκ = iter.κ * iter.τ + aff_dir.Δκ * aff_dir.Δτ - σ*μ_hat

    # Solve KKT.
    rhs = AbstractKktRhs(dx, dy, dz, ds, dτ, dκ)
    comb_dir, ir_iters = solve_abstract_kkt(iter, λ, prob, rhs, kkt, scaling)
    register_comb_dir_stats!(stats, ir_iters, μ_hat)
    return comb_dir
end

function solve_abstract_kkt(iter::Iterate{T}, λ::Vector{T}, prob::Problem{T}, rhs::AbstractKktRhs{T},
                            kkt::KktCache{T}, scaling::CompositeConeScaling{T}) where T <: AbstractFloat
    W = scaling
    temp_z = W * jordan_inv(λ, rhs.ds, prob.C)
    new_rhs = [rhs.dx, -rhs.dy, temp_z - rhs.dz, zeros(2*kkt.q)]
    (x2, y2, z2), best_ref_ratio, last_ref_ratio, num_iters, corr_ratio = solve_and_refine(
                     kkt, new_rhs; max_iters=3)

    Wz1 = W * kkt.z1
    Δτ = ((rhs.dτ - rhs.dκ / iter.τ + prob.c'*x2 + prob.b'*y2 + prob.h'*z2)
                / (iter.κ ./ iter.τ - prob.c'*kkt.x1 - prob.b'*kkt.y1 - prob.h'*kkt.z1))
    Δx = x2 + Δτ * kkt.x1
    Δy = y2 + Δτ * kkt.y1
    Δz = z2 + Δτ * kkt.z1
    Δs = -temp_z - W*(W*Δz)
    Δκ = -(rhs.dκ + iter.κ*Δτ)/iter.τ
    return StepDirection(Δx, Δy, Δz, Δτ, Δs, Δκ), num_iters
end


text:
@testset "Testing matching_tools.jl" begin

    @testset "random_prefs for one-to-one" begin
        nums = (8, 6)
        prefs_arrays = random_prefs(nums..., allow_unmatched=false)
        prefs_arrays_allowed = random_prefs(nums..., allow_unmatched=true)
        prefs_arrays_all = tuple(prefs_arrays..., prefs_arrays_allowed...)

        @testset "size" begin
            for (i, prefs_array) in enumerate(prefs_arrays_all)
                @test size(prefs_array) == (nums[i%2+1]+1, nums[(i+1)%2+1])
            end
        end

        @testset "permutation" begin
            sorted_arrays = [repeat(0:nums[i%2+1], 1, nums[i]) for i in 1:2]
            for (i, prefs_array) in enumerate(prefs_arrays_all)
                @test sort(prefs_array; dims=1) == sorted_arrays[(i+1)%2+1]
            end
        end

        @testset "unmatched not allowed" begin
            for prefs_array in prefs_arrays
                @test all(prefs_array[end, :] .== 0)
            end
        end

        @testset "unmatched not most preferred" begin
            for prefs_array in prefs_arrays_allowed
                @test all(prefs_array[1, :] .!= 0)
            end
        end
    end

    @testset "random_prefs for many-to-one" begin
        nums = (8, 6)
        s_prefs, c_prefs, caps =
            random_prefs(nums..., ReturnCaps, allow_unmatched=false)
        s_prefs_allowed, c_prefs_allowed, caps_allowed =
            random_prefs(nums..., ReturnCaps, allow_unmatched=true)
        prefs_arrays_all = (s_prefs, c_prefs, s_prefs_allowed, c_prefs_allowed)

        @testset "permutation" begin
            sorted_arrays = [repeat(0:nums[i%2+1], 1, nums[i]) for i in 1:2]
            for (i, prefs_array) in enumerate(prefs_arrays_all)
                @test sort(prefs_array; dims=1) == sorted_arrays[(i+1)%2+1]
            end
        end

        @testset "permutation" begin
            sorted_arrays = [repeat(0:nums[i%2+1], 1, nums[i]) for i in 1:2]
            for (i, prefs_array) in enumerate(prefs_arrays_all)
                @test sort(prefs_array; dims=1) == sorted_arrays[(i+1)%2+1]
            end
        end

        @testset "unmatched not allowed" begin
            for prefs_array in (s_prefs, c_prefs)
                @test all(prefs_array[end, :] .== 0)
            end
        end

        @testset "unmatched not most preferred" begin
            for prefs_array in (s_prefs_allowed, c_prefs_allowed)
                @test all(prefs_array[1, :] .!= 0)
            end
        end

        @testset "caps positive" begin
            for x in (caps, caps_allowed)
                @test all(x .> 0)
            end
        end

        @testset "caps not exceed rankings of unmatched" begin
            for (prefs_array, x) in ((c_prefs, caps),
                                     (c_prefs_allowed, caps_allowed))
                I = findall(prefs_array .== 0)
                rankings_unmatched = getindex.(I, 1)
                @test all(x .<= rankings_unmatched)
            end
        end
    end

end


text:
using Random

using Dates, TimeZones
using StatsBase, Impute
using JuliaDB, DataFrames, DataFramesMeta
using Mise
using MicroLogging

function run_model()

    @info "Start OU Process Prediction"
    flush(stdout); flush(stderr)

    seoul_codes = [
        111121,111123,111131,111141,111142,
        111151,111152,111161,111171,111181,
        111191,111201,111212,111221,111231,
        111241,111251,111261,111262,111273,
        111274,111281,111291,111291,111301,
        111301,111311]
    seoul_names = [
        "중구","종로구","용산구","광진구","성동구",
        "중랑구","동대문구","성북구","도봉구","은평구",
        "서대문구","마포구","강서구","구로구","영등포구",
        "동작구","관악구","강남구","서초구","송파구",
        "강동구","금천구","강북구","강북구","양천구",
        "양천구","노원구"]
    # construct named tuple
    seoul_stations = (; zip(Symbol.(seoul_names), seoul_codes)...)

    # NamedTuple (::Symbol -> ::DataFrame)
    rawdf = load_data_DNN("/input/jongro_seoul.csv", seoul_stations)
    #=
    first(df, 5)
    │ Row │ stationCode │ date                      │ lat      │ long      │ SO2     │ CO      │ O3      │ NO2     │ PM10    │ PM25    │ temp    │ u       │ v       │ pres    │ prep     │ snow     │
    │     │ Int64       │ TimeZones.ZonedDateTime   │ Float64  │ Float64   │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64  │ Float64  │
    ├─────┼─────────────┼───────────────────────────┼──────────┼───────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┤─────────┤──────────┤──────────┤
    │ 1   │ 111123      │ 2015-01-01T01:00:00+09:00 │ 37.572   │ 127.005   │ 0.004   │ 0.2     │ 0.02    │ 0.009   │ 57.0    │ 0.0     │ -7.4    │ 4.41656 │ 1.60749 │ 1011.8  │ missing  │ missing  │
    │ 2   │ 111123      │ 2015-01-01T02:00:00+09:00 │ 37.572   │ 127.005   │ 0.005   │ 0.2     │ 0.019   │ 0.008   │ 70.0    │ 3.0     │ -8.0    │ 4.22862 │ 1.53909 │ 1011.7  │ missing  │ missing  │
    │ 3   │ 111123      │ 2015-01-01T03:00:00+09:00 │ 37.572   │ 127.005   │ 0.005   │ 0.2     │ 0.02    │ 0.006   │ 92.0    │ 5.0     │ -8.4    │ 3.57083 │ 1.29968 │ 1012.1  │ missing  │ missing  │
    │ 4   │ 111123      │ 2015-01-01T04:00:00+09:00 │ 37.572   │ 127.005   │ 0.004   │ 0.2     │ 0.019   │ 0.005   │ 111.0   │ 2.0     │ -8.8    │ 4.60449 │ 1.6759  │ 1012.3  │ missing  │ missing  │
    │ 5   │ 111123      │ 2015-01-01T05:00:00+09:00 │ 37.572   │ 127.005   │ 0.005   │ 0.2     │ 0.019   │ 0.006   │ 127.0   │ 5.0     │ -9.1    │ 5.35625 │ 1.94951 │ 1011.8  │ missing  │ missing  │
    =#

    #===== start of parameter zone =====#
    total_fdate, total_tdate = get_date_range(rawdf)
    train_fdate = ZonedDateTime(2012, 1, 1, 0, tz"Asia/Seoul")
    train_tdate = ZonedDateTime(2017, 12, 31, 23, tz"Asia/Seoul")
    test_fdate = ZonedDateTime(2018, 1, 1, 1, tz"Asia/Seoul")
    test_tdate = ZonedDateTime(2018, 12, 31, 23, tz"Asia/Seoul")

    # due to seasonality, train date range must be over 1 year
    @assert train_tdate - train_fdate >= Dates.Day(365)

    # stations
    #=
    train_stn_names = [
        :중구,:종로구,:용산구,:광진구,:성동구,
        :중랑구,:동대문구,:성북구,:도봉구,:은평구,
        :서대문구,:마포구,:강서구,:구로구,:영등포구,
        :동작구,:관악구,:강남구,:서초구,:송파구,
        :강동구,:금천구,:강북구,:강북구,:양천구,
        :양천구,:노원구]
    =#
    train_stn_names = [:종로구, :강서구, :송파구, :강남구]
    train_stn_names = [:종로구]
    #=
    train_stn_names = [
        :중구,:종로구,:용산구,:광진구,:성동구,
        :중랑구,:동대문구,:성북구,:도봉구,:은평구,
        :서대문구,:마포구,:강서구,:구로구,:영등포구]
    =#

    # test set is Jongro only
    test_stn_names = [:종로구]

    df = filter_raw_data(rawdf, train_fdate, train_tdate, test_fdate, test_tdate)
    @show first(df, 10)

    features = [:SO2, :CO, :O3, :NO2, :PM10, :PM25, :temp, :u, :v, :pres, :humid, :prep, :snow]
    # If you want exclude some features, modify train_features
    # exclude :PM10, :PM25 temporarily for log transform
    train_features = [:SO2, :CO, :O3, :NO2, :PM10, :PM25, :temp, :u, :v, :pres, :humid, :prep, :snow]
    target_features = [:PM10, :PM25]

    # For GPU, change precision of Floating numbers
    _eltype::DataType = Float32
    scaling_method = :zscore

    sample_size = 24
    output_size = 24
    epoch_size = 300
    batch_size = 32
    input_size = sample_size * length(train_features)
    #===== end of parameter zone =====#

    log_prefix = "log_"
    inv_prefix = "inv_"
    scaled_prefix = "scaled_"
    scaled_features = [Symbol(scaled_prefix, f) for f in features]
    scaled_train_features = [Symbol(scaled_prefix, f) for f in train_features]
    scaled_target_features = [Symbol(scaled_prefix, f) for f in target_features]

    # Preprocessing
    # 0. imputation
    DataFrames.allowmissing!(df, train_features)
    for ycol in train_features
        df[!, ycol] = Impute.srs(df[!, ycol])
    end
    DataFrames.disallowmissing!(df, train_features)

    # 2. Scaling
    if scaling_method == :zscore
        # scaling except :PM10, :PM25
        zscore!(df, train_features, scaled_train_features)

        statvals = extract_col_statvals(df, train_features)
    elseif scaling_method == :minmax
        minmax_scaling!(df, train_features, scaled_train_features, 0.0, 10.0)

        statvals = extract_col_statvals(df, train_features)
    elseif scaling_method == :logzscore
        # zscore except targets
        _train_features = setdiff(train_features, target_features)
        _scaled_train_features = setdiff(scaled_train_features, scaled_target_features)
        zscore!(df, _train_features, _scaled_train_features)

        # Log transform :PM10, :PM25
        # add 10.0 not to log zero values
        # Y = zscore(log(X + 10))
        log_target_features = [Symbol(log_prefix, f) for f in target_features]
        for (target, log_target) in zip(target_features, log_target_features)
            df[!, log_target] = log.(df[!, target] .+ 10.0)
        end

        zscore!(df, log_target_features, scaled_target_features)

        statvals = extract_col_statvals(df, vcat(train_features, log_target_features))
    elseif scaling_method == :invzscore
        # zscore except targets
        _train_features = setdiff(train_features, target_features)
        _scaled_train_features = setdiff(scaled_train_features, scaled_target_features)
        zscore!(df, _train_features, _scaled_train_features)

        # Inverse transform :PM10, :PM25
        # add 10.0 not to divide zero value
        # Y = zscore(1.0 / (X + 10))
        inv_target_features = [Symbol(inv_prefix, f) for f in target_features]
        for (target, inv_target) in zip(target_features, inv_target_features)
            df[!, inv_target] = 1.0 ./ (df[!, target] .+ 10.0)
        end

        zscore!(df, inv_target_features, scaled_target_features)

        statvals = extract_col_statvals(df, vcat(train_features, inv_target_features))
    elseif scaling_method == :logminmax
        # minmax except targets
        _train_features = setdiff(train_features, target_features)
        _scaled_train_features = setdiff(scaled_train_features, scaled_target_features)
        minmax_scaling!(df, _train_features, _scaled_train_features, 0.0, 10.0)

        # Log transform :PM10, :PM25
        # add 10.0 not to log zero values
        log_target_features = [Symbol(log_prefix, f) for f in target_features]
        for (target, log_target) in zip(target_features, log_target_features)
            df[!, log_target] = log.(df[!, target] .+ 10.0)
        end

        # then minmax
        minmax_scaling!(df, log_target_features, scaled_target_features, 0.0, 10.0)

        statvals = extract_col_statvals(df, vcat(train_features, log_target_features))
    end

    @show first(df, 5)
    # convert Float types
    for feature in features
        df[!, feature] = _eltype.(df[!, feature])
    end

    # plot histogram regardless to station
    for target in target_features
        plot_histogram(df, target, "/mnt/")
        plot_histogram(df, Symbol(:scaled_, target), "/mnt/")

        if scaling_method == :logzscore || scaling_method == :logminmax
            plot_histogram(df, Symbol(:log_, target), "/mnt/")
        end

        if scaling_method == :invzscore
            plot_histogram(df, Symbol(:inv_, target), "/mnt/")
        end
    end

    # line plot of first train station
    plot_lineplot_total(filter_station(df, seoul_stations[train_stn_names[1]]), :PM10, "/mnt/")
    plot_lineplot_total(filter_station(df, seoul_stations[train_stn_names[1]]), :PM25, "/mnt/")
    #plot_pcorr(df, scaled_features, features, "/mnt/")

    # windowed dataframe for train/valid (input + output)
    test_wd = []

    # random permutation train_valid_wd itself for splitting train/valid set
    rng = MersenneTwister()

    # windows dataframe for test (only for 종로구)
    for name in test_stn_names
        code = seoul_stations[name]
        stn_df = filter_station_DNN(df, code)

        # to add zero(0) offset for X₀
        push!(test_wd,
            window_df(stn_df, sample_size, output_size + 1, test_fdate, test_tdate))
    end
    # Flatten
    test_wd = collect(Base.Iterators.flatten(test_wd))
    test_size = length(test_wd)

    # simply collect dates, determine exact date for prediction (for 1h, 24h, and so on) later
    test_dates = collect(test_fdate + Hour(sample_size - 1):Hour(1):test_tdate - Hour(output_size))

    for target in target_features
        @info "$(string(target)) Training..."
        flush(stdout); flush(stderr)

        @info "training feature : " train_features
        @info "sizes (sample, output, epoch, batch) : ", sample_size, output_size, epoch_size, batch_size

        for name in test_stn_names
            code = seoul_stations[name]
            stn_df = filter_station_DNN(df, code)

            Base.Filesystem.mkpath("/mnt/$(name)/")

            max_time = input_size * 15
            acf = StatsBase.autocor(stn_df[!, Symbol(scaled_prefix, target)], 0:max_time)

            # free minibatch after training because of memory usage
            ou_df = evolve_OU(test_wd, target, Symbol(scaled_prefix, target),
                acf, statvals, _eltype, sample_size, output_size, test_dates, "/mnt/$(name)/")

            # create directory per each time
            for i = 1:output_size
                i_pad = lpad(i, 2, '0')
                Base.Filesystem.mkpath("/mnt/$(name)/$(i_pad)/")
            end

            dfs_out = export_CSV(DateTime.(test_dates, Local), ou_df, target, output_size, "/mnt/$(name)/", String(target))
            plot_DNN_scatter(dfs_out, target, output_size, "/mnt/$(name)/", String(target))
            plot_DNN_histogram(dfs_out, target, output_size, "/mnt/$(name)/", String(target))

            plot_DNN_lineplot(dfs_out, target, output_size, "/mnt/$(name)/", String(target))

            df_corr = compute_corr(dfs_out, output_size, "/mnt/$(name)/", String(target))
            plot_corr(df_corr, output_size, "/mnt/$(name)/", String(target))
        end
    end
end

run_model()


text:
export ProductAngularQuadrature
export angular_quadrature

# Angular quadrature defined on the unit sphere octant in the upper right, closest to the viewer.
# The angles and weights are transformed to the other octants by symmetry.
#     +----+----+
#    /    /    /|
#   +----+----+ |
#  /    /    /| +
# +----+----+ |/|
# |    |this| + |
# |    | one|/| +
# +----+----+ |/
# |    |    | +
# |    |    |/
# +----+----+
#
# The spherical coordinates are defined in the following manner
# Ω̂ = (Ω_i, Ω_j, Ω_k) = (cos(θ),  sin(θ)cos(γ),  sin(θ)sin(γ))
#                     = (     μ, √(1-μ²)cos(γ), √(1-μ²)sin(γ))
#
#        j
#        ^
#        |   θ is the polar angle about the i-axis (x-direction)
#        |   γ is the azimuthal angle in the j-k plane, from the j-axis
#        |
#        |
#        |
#       /|
#      (γ|
#       \|--------------------> i
#       / \θ)
#      /   \
#     /     \
#    /       \ Ω̂
#   /         v
#  𝘷
#  k

# TODO: Function to convert product into general quadrature
#struct GeneralAngularQuadrature{T <: AbstractFloat} <: AngularQuadrature
#    Ω̂::Tuple{Point{T}} # Points on the unit sphere satisfying θ ∈ (0, π/2), γ ∈ (0, π)
#    w::Tuple{T} # Weights for each point
#end

struct ProductAngularQuadrature{nγ, nθ, T}
    wγ::SVector{nγ,T}   # Weights for the azimuthal angles
    γ::SVector{nγ,T}    # Azimuthal angles, γ ∈ (0, π)
    wθ::SVector{nθ,T}   # Weights for the polar angles
    θ::SVector{nθ,T}    # Polar angles, θ ∈ (0, π/2)
end

function angular_quadrature(azimuthal_form::Symbol,
                            azimuthal_degree::Integer,
                            polar_form::Symbol,
                            polar_degree::Integer,
                            type::Type{T}) where {T}
    if azimuthal_form === :chebyshev
        azi_weights_half, azi_angles_half = chebyshev_angular_quadrature(azimuthal_degree, T)
        azi_weights = vcat(azi_weights_half, azi_weights_half)
        azi_angles = vcat(azi_angles_half, reverse(π .- azi_angles_half))
    else
        error("Cannot identify azimuthal quadrature.")
    end
    if polar_form === :chebyshev
        pol_weights, pol_angles = chebyshev_angular_quadrature(polar_degree, T)
    else
        error("Cannot identify azimuthal quadrature.")
    end
    return ProductAngularQuadrature{2azimuthal_degree, polar_degree, T}(
                                       azi_weights, azi_angles,
                                       pol_weights, pol_angles
                                      )
end

function chebyshev_angular_quadrature(M::Integer, type::Type{T}) where {T}
    # A Chebyshev-type quadrature for a given weight function is a quadrature formula
    # with equal weights. This function produces evenly spaced angles with equal weights.
    weights = SVector(ntuple(m->T(1)/M,           M))
    angles  = SVector(ntuple(m->π*(2T(m) - 1)/4M, M))
    return weights, angles
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
include("../src/PowerFlowNetworks.jl")

using .PowerFlowNetworks
using .PowerFlowNetworks: get_features_graph, get_features_opf, _get_edges_distance, get_nv, get_ne
using Test
using Graphs
using SQLite, DataFrames
using Dates

const db_path = "test/data/TEST_PowerFlowNetworks_SQLite.sqlite"
isfile(db_path) && rm(db_path)
const testdir = dirname(@__FILE__)

tests = [
   "core",
   "read_features",
   "io/read_matpower",
   "io/read_go",
   "graphs/graphs",
   "graphs/operations",
   "db/setup_db",
   "db/inserts",
   "db/operations",
   "db/infos",
   "utils/ischordal",
   "utils/clique_cliquetree"
]

try
    @testset "PowerFlowNetwork" begin
        for t in tests
            tp = joinpath(testdir, "$(t).jl")
            include(tp)
        end
    end
finally
    if any(x -> occursin("db", x), tests)
        rm(db_path)
    end
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
#!/usr/bin/julia

println("""
WARNING: this is just a proof of concept, not a fully featured terminal application!
         TTY->serial works only per line (instead of per character)

""")
println("Compiling, please wait…")
using PySerial

function serial2tty(s)
  while true
   w = inWaiting(ser)
   if w > 0
     s = read(ser,w)
     print(s)
   end
   sleep(0.001)
  end
end

function tty2serial(s)
  while !eof(STDIN)
    c = read(STDIN, Uint8) # it seems to be buffered until newline, that's not so good for a general terminal application, but ok...
    write(s,[c])
  end
end

port = length(ARGS) > 0 ? ARGS[1] : "/dev/ttyUSB0"
baudrate = length(ARGS) > 1 ? int(ARGS[2]) : 9600

ser = Serial(port, baudrate=baudrate)
println("== Welcome to juliaterm.jl! Port: $port, baudrate: $baudrate ==")

@async serial2tty(ser)
tty2serial(ser)


text:
import Base: ∘

abstract type AbstractMobiusTransformation{T <: Number} <: Function end

"""
    MobiusTransformation(a,b,c,d)

Möbius transformation, a conformal automorphism in the Riemann sphere
\$z \\mapsto \\frac{az+b}{cz+d}\$
with \$ad-bc\\neq 0\$.
"""
struct MobiusTransformation{T} <: AbstractMobiusTransformation{T}
  a::T
  b::T
  c::T
  d::T

  function MobiusTransformation{T}(a0::T, b0::T, c0::T, d0::T) where T <: Number
    #@assert a0 * d0 - b0 * c0 != 0 "Invalid parameters for Möbius transformation."
    new(a0,b0,c0,d0)
  end

end

MobiusTransformation(a0::T, b0::T, c0::T, d0::T) where {T <: Number} =
MobiusTransformation{T}(a0,b0,c0,d0)

MobiusTransformation(a0::Number, b0::Number, c0::Number, d0::Number) =
MobiusTransformation(promote(a0,b0,c0,d0)...)

"""
    LinearTransformation(a)

Linear transformation (homotecy and rotation) in the Riemann sphere
\$z \\mapsto a z\$.
"""
struct LinearTransformation{T} <: AbstractMobiusTransformation{T}
  a::T

  function LinearTransformation{T}(a0::T) where T <: Number
    #@assert a0 == 0 "Invalid parameters for Möbius transformation."
    new(a0)
  end
end

LinearTransformation(a0::T) where {T <: Number} = LinearTransformation{T}(a0)

MobiusTransformation(a0::Number) = LinearTransformation(a0)

"""
    Translation(b)

Translation in the Riemann sphere
\$z \\mapsto z + b\$.
"""
struct Translation{T} <: AbstractMobiusTransformation{T}
  b::T

  function Translation{T}(b0::T) where T <:Number
    new(b0)
  end
end

Translation(b0::T) where {T <: Number} = Translation{T}(b0)

"""
    AffineTransformation(a,b)

Affine transformation in the Riemann sphere
\$z \\mapsto a z + b\$.
"""
struct AffineTransformation{T} <: AbstractMobiusTransformation{T}
  a::T
  b::T

  function AffineTransformation{T}(a0::T, b0::T) where T <:Number
    new(a0,b0)
  end
end

AffineTransformation(a0::T, b0::T) where {T <: Number} = AffineTransformation{T}(a0,b0)

AffineTransformation(a0::Number, b0::Number) = AffineTransformation(promote(a0,b0)...)

MobiusTransformation(a0::Number, b0::Number) = AffineTransformation(promote(a0,b0)...)

"""
    InversionReflection(b)

Inversion reflection in the Riemann sphere
\$z \\mapsto \\frac{b}{z}\$.
"""
struct InversionReflection{T} <: AbstractMobiusTransformation{T}
  b::T

  function InversionReflection{T}(b0::T) where T <:Number
    new(b0)
  end
end

InversionReflection(b0::T) where {T <: Number} = InversionReflection{T}(b0)

Base.show(io::IO, f::MobiusTransformation{T}) where T =
  print(io, "MöbiusTransformation{$T}: z -> ((", f.a, ")z + (", f.b, ")) / ((", f.c, ")z + (", f.d, "))")

Base.show(io::IO, f::LinearTransformation{T}) where T  =
  print(io, "LinearTransformation{$T}: z -> (", f.a, ")z")

Base.show(io::IO, f::Translation{T}) where T  =
  print(io, "Translation{$T}: z -> z + (", f.b, ")")

Base.show(io::IO, f::AffineTransformation{T}) where T  =
  print(io, "AffineTransformation{$T}: z -> (", f.a, ")z + (", f.b, ")")

Base.show(io::IO, f::InversionReflection{T}) where T  =
  print(io, "InversionReflection{$T}: z -> (", f.b, ") / z")

function (f::MobiusTransformation)(z::Number)
  # ToDo! Create a method to eval infinity case
  if isinf(z)
    return f.a / f.c
  end

  (f.a * z + f.b) / (f.c * z + f.d)
end

function (f::LinearTransformation)(z::Number)
  f.a * z
end

function (f::Translation)(z::Number)
  z + f.b
end

function (f::AffineTransformation)(z::Number)
  f.a * z + f.b
end

function (f::InversionReflection)(z::Number)
  # ToDo! Create a method to eval infinity case
  if isinf(z)
    return 0
  end

  f.b / z
end

"""
Möbius transformation parameters accesors.
"""
a(f::MobiusTransformation) = f.a
b(f::MobiusTransformation) = f.b
c(f::MobiusTransformation) = f.c
d(f::MobiusTransformation) = f.d

a(f::LinearTransformation) = f.a
b(f::LinearTransformation) = 0
c(f::LinearTransformation) = 0
d(f::LinearTransformation) = 1

a(f::Translation) = 1
b(f::Translation) = f.b
c(f::Translation) = 0
d(f::Translation) = 1

a(f::AffineTransformation) = f.a
b(f::AffineTransformation) = f.b
c(f::AffineTransformation) = 0
d(f::AffineTransformation) = 1

a(f::InversionReflection) = 0
b(f::InversionReflection) = f.b
c(f::InversionReflection) = 1
d(f::InversionReflection) = 0

"""
Trace

\$tr(T)=a+d\$
"""
tr(f::AbstractMobiusTransformation) = a(f) + d(f)

tr(f::MobiusTransformation) = f.a + f.d

tr(f::LinearTransformation) = f.a + 1

tr(f::Translation) = 1

tr(f::AffineTransformation) = f.a + 1

tr(f::InversionReflection) = 0

"""
Determinant.

\$det(T)=ad-bc\$
"""
det(f::AbstractMobiusTransformation) = a(f)*d(f) - b(f)*c(f)

det(f::MobiusTransformation) = f.a * f.d - f.b * f.c

det(f::LinearTransformation) = f.a

det(f::Translation) = 1

det(f::AffineTransformation) = f.a

det(f::InversionReflection) = -f.b

#function normalize(f::MobiusTransformation)
#  dt = det(f)
#  MobiusTransformation(f.a/dt, f.b/det, f.c/det, f.d/det)
#end

#function normalize(f::LinearTransformation)
#  f
#end

"""
Inverse.

\$T^{-1}(z)=\\frac{d z - b}{-c z + a}\$
"""
inverse(f::MobiusTransformation) = MobiusTransformation(f.d, -f.b, -f.c, f.a)

inverse(f::LinearTransformation) = LinearTransformation(1.0 / f.a)

inverse(f::Translation) = Translation(-f.b)

inverse(f::AffineTransformation) = AffineTransformation(1.0 / f.a, -f.b / f.a)

inverse(f::InversionReflection) = InversionReflection(f.b)

"""
Derivative.

\$T'(z)=\\frac{ad-bc}{(c z + d)^2}\$
"""
function derivative(f::MobiusTransformation)
  function der(z::Number)
    (f.a * f.d - f.b * f.c) / ((f.c * z + f.d)^2)
  end
end

function derivative(f::LinearTransformation)
  function der(z::Number)
    f.a
  end
end

function derivative(f::Translation)
  function der(z::Number)
    1
  end
end

function derivative(f::AffineTransformation)
  function der(z::Number)
    f.a
  end
end

function derivative(f::InversionReflection)
  function der(z::Number)
    -f.b / (z^2)
  end
end

"""
    fixedpoints(f)

Returns the fixed points of a Möbius transformation. Always returns two points,
even in the parabolic case returns the one fixed point duplicated.
"""
function fixedpoints(f::MobiusTransformation{T}) where T <: Real
  if f.c == 0
    return f.b / (f.d - f.a), Inf
  end
  discr = (f.a + f.d)^2 - 4 * (f.a * f.d - f.b * f.c)
  rd = discr < 0 ? sqrt(complex(discr)) : sqrt(discr)
  (f.a - f.d + rd)/(2*f.c), (f.a - f.d - rd)/(2*f.c)
end

function fixedpoints(f::MobiusTransformation)
  if f.c == 0
    return f.b / (f.d - f.a), Inf
  end
  rd = sqrt((f.a + f.d)^2 - 4 * (f.a * f.d - f.b * f.c))
  (f.a - f.d + rd)/(2*f.c), (f.a - f.d - rd)/(2*f.c)
end

function fixedpoints(f::LinearTransformation)
  0, Inf
end

function fixedpoints(f::Translation)
  Inf, Inf
end

function fixedpoints(f::AffineTransformation)
  if f.a == 1
    return Inf, Inf
  end
  f.b / (1 - f.a), Inf
end

function fixedpoints(f::InversionReflection{T}) where T <: Real
  rd = f.b < 0 ? sqrt(complex(f.b)) : sqrt(f.b)
  rd, -rd
end

function fixedpoints(f::InversionReflection)
  rd = sqrt(f.b)
  rd, -rd
end

"""
Kind of Möbius transformation:
- `:parabolic` if \$t^2=4\$.
- `:elliptic` if \$t^2\\in [0,4)\$.
- `:hyperbolic` if \$t^2\\in (-\\infty,0)\\cup (4,\\infty)\$.
- `:loxodromic` otherwise.
where \$t=tr(T)/det(T)\$.
"""
function kind(f::AbstractMobiusTransformation)
  t = tr(f) / det(f)
  t2 = t^2
  if imag(t2) ≈ 0
    t2 = real(t2)
    if t2 == 4
      return :parabolic
    end

    if 0 <= t2 < 4
      return :elliptic
    end

    return :hyperbolic
  end

  :loxodromic
end

kind(f::Translation) = :parabolic

"""
Compose two Möbius transformations.
"""
compose(f::AbstractMobiusTransformation, g::AbstractMobiusTransformation) =
  MobiusTransformation(a(f)*a(g)+b(f)*c(g), a(f)*b(g)+b(f)*d(g),
    c(f)*a(g)+d(f)*c(g), c(f)*b(g)+c(f)*d(g))

compose(f::MobiusTransformation, g::MobiusTransformation) =
  MobiusTransformation(f.a*g.a+f.b*g.c, f.a*g.b+f.b*g.d,
    f.c*g.a+f.d*g.c, f.c*g.b+f.d*g.d)

compose(f::LinearTransformation, g::LinearTransformation) =
  LinearTransformation(f.a * g.a)

compose(f::Translation, g::Translation) =
  Translation(f.b + g.b)

compose(f::AffineTransformation, g::AffineTransformation) =
  AffineTransformation(f.a * g.a, f.a * g.b + f.b)

compose(f::InversionReflection, g::InversionReflection) =
  InversionReflection(f.b / g.b)

# LinearTransformation vs ...
compose(f::LinearTransformation, g::MobiusTransformation) =
  MobiusTransformation(f.a*g.a, f.a*g.b, g.c, g.d)

compose(f::MobiusTransformation, g::LinearTransformation) =
  MobiusTransformation(f.a*g.a, f.b, f.c*g.a, f.d)

compose(f::LinearTransformation, g::Translation) =
  AffineTransformation(f.a, f.a * g.b)

compose(f::Translation, g::LinearTransformation) =
  AffineTransformation(g.a, f.b)

compose(f::LinearTransformation, g::AffineTransformation) =
  AffineTransformation(f.a * g.a, f.a * g.b)

compose(f::AffineTransformation, g::LinearTransformation) =
  AffineTransformation(f.a * g.a, f.b)

# Translation vs ...
compose(f::Translation, g::MobiusTransformation) =
  MobiusTransformation(g.a+f.b*g.c, g.b+f.b*g.d, g.c, g.d)

compose(f::MobiusTransformation, g::Translation) =
  MobiusTransformation(f.a, f.a*g.b+f.b, f.c, f.c*g.b+f.d)

compose(f::Translation, g::AffineTransformation) =
  AffineTransformation(g.a, g.b+f.b)

compose(f::AffineTransformation, g::Translation) =
  AffineTransformation(f.a, f.a*g.b+f.b)

# AffineTransformation vs ...
compose(f::AffineTransformation, g::MobiusTransformation) =
  MobiusTransformation(f.a*g.a+f.b*g.c, f.a*g.b+f.b*g.d, g.c, g.d)

compose(f::MobiusTransformation, g::AffineTransformation) =
  MobiusTransformation(f.a*g.a, f.a*g.b+f.b, f.c*g.a, f.c*g.b+f.d)

# InversionReflection vs ...
# result in a MobiusTransformation, fall back to AbstractMobiusTransformation implementation

"""
Binary operator to compose two Möbius transformations.
"""
∘(f::AbstractMobiusTransformation, g::AbstractMobiusTransformation) = compose(f,g)

"""
    maptozerooneinf(z1,z2,z3)

Create the Möbius transformation such that \$z_1 \\mapsto 0\$,
\$z_2 \\mapsto 1\$ and \$z_3 \\mapsto \\infty\$.
"""
maptozerooneinf(z1::Number, z2::Number, z3::Number) =
  MobiusTransformation(z2-z3, z1*(z3-z2), z2-z1, z3*(z1-z2))

"""
    mapfromzerooneinf(z1,z2,z3)

Create the Möbius transformation such that \$0 \\mapsto z_1\$,
  \$1 \\mapsto z_2\$ and \$\\infty \\mapsto z_3\$.
"""
mapfromzerooneinf(z1::Number, z2::Number, z3::Number) =
  MobiusTransformation(z3*(z1-z2), z1*(z2-z3), z1-z2, z2-z3)

"""
Create the Möbius transformation such that \$z_1 \\mapsto w_1\$,
  \$z_2 \\mapsto w_2\$ and \$z_3 \\mapsto w_3\$.
"""
function MobiusTransformation(z1::Number, z2::Number, z3::Number,
  w1::Number, w2::Number, w3::Number)
  f = maptozerooneinf(z1,z2,z3)
  g = mapfromzerooneinf(w1,w2,w3)
  g∘f
end

#=
function (f::MobiusTransformation)(c::CLine)
  CLine(
   abs2(f.d) * c.A - 2 * real( conj(f.c) * f.d * c.B ) + abs2(f.c) * c.C,
   conj(f.a) * (f.d * c.B - f.c * c.C) - conj(f.b) * (f.c * conj(c.B) - f.d * c.A),
   abs2(f.b) * c.A - 2 * real( conj(f.a) * f.b * c.B ) + abs2(f.a) * c.C
  )
end

# ToDO: Methods!!!
=#

"""
Given a Möbius transformation \$f(z)=\\frac{az+b}{cz+d}\$ and a circle
\$C: |center-z|=radius\$, the circle  \$f(C)\$ has center

\$ w = f(center - \\frac{radius^2}{\\overline{\\frac{d}{c}+center}})

and radius

\$|w - f(center+radius)|\$.
"""
function (f::MobiusTransformation)(c::Circle)
  if f.c ≈ 0.0
    return Circle( f(c.center), abs(f.a/f.d)*c.radius )
  elseif abs(c.center + f.d/f.c) ≈ c.radius
    z = c.center + f.d/f.c
    w1 = f(c.center + z)
    w2 = f(c.center + im*z)
    return Line(w1, angle(w1-w2))
  end

  # "Grandma's recipe" from "Indra's Pearls" to calculate image of a circle
  z = c.center - c.radius*c.radius/conj(f.d/f.c + c.center)
  Circle(f(z), abs(z - f(c.center + c.radius)))
end

function (f::LinearTransformation)(c::Circle)
  Circle( f(c.center), abs(f.a)*c.radius )
end

function (f::Translation)(c::Circle)
  Circle( f(c.center), c.radius )
end

function (f::AffineTransformation)(c::Circle)
  Circle( f(c.center), abs(f.a)*c.radius )
end

function (f::InversionReflection)(c::Circle)
  if abs(c.center) ≈ c.radius
    w1 = f(2*c.center)
    w2 = f(c.center + im*c.center)
    return Line(w1, angle(w1-w2))
  end

  # "Grandma's recipe" from "Indra's Pearls" to calculate image of a circle
  z = c.center - c.radius*c.radius/conj(c.center)
  Circle(f(z), abs(z - f(c.center + c.radius)))
end

function (f::MobiusTransformation)(c::Arc)
  if f.c ≈ 0.0
    return Arc( f(c.center), f(c.p1), f(c.p2) )
  elseif abs(c.center + f.d/f.c) ≈ radius(c)
    #ToDo!: Case 2 rays!
    return LineSegment(f(c.p1), f(c.p2))
  end

  # "Grandma's recipe" from "Indra's Pearls" to calculate image of a circle
  z = c.center - radius2(c)/conj(f.d/f.c + c.center)
  Arc(f(z), f(c.p1), f(c.p2))
end

function (f::LinearTransformation)(c::Arc)
  Arc(f(c.center), f(c.p1), f(c.p2))
end

function (f::Translation)(c::Arc)
  Arc(f(c.center), f(c.p1), f(c.p2))
end

function (f::AffineTransformation)(c::Arc)
  Arc(f(c.center), f(c.p1), f(c.p2))
end

function (f::InversionReflection)(c::Arc)
  if abs(c.center) ≈ c.radius
    #ToDo!: Case 2 rays!
    return LineSegment(f(c.p1), f(c.p2))
  end

  # "Grandma's recipe" from "Indra's Pearls" to calculate image of a circle
  z = c.center - radius2(c)/conj(c.center)
  Arc(f(z), f(c.p1), f(c.p2))
end

function (f::MobiusTransformation)(l::Line)
  if f.c ≈ 0.0
    w1 = f(l.base)
    return Line( w1, angle(w1-f(l.base+exp(l.angle*im))) )
  elseif l.base ≈ -f.d/f.c
    w = f(Inf)
    return Line(w, angle(w-f(Inf)))
  else
    a = angle(l.base + f.d/f.c)
    if a ≈ l.angle || a ≈ (l.angle + pi)
      w = f(l.base)
      return Line(w, angle(w-f(Inf)))
    end
  end

  v = exp(im*l.angle)
  w1 = f(l.base)
  w2 = f(l.base + v)
  w3 = f(l.base - v)
  Circle(circumcenter(w1,w2,w3), circumradius(w1,w2,w3))
end

function (f::LinearTransformation)(l::Line)
  Line(f(l.base), l.angle + angle(f.a))
end

function (f::Translation)(l::Line)
  Line(f(l.base), l.angle)
end

function (f::AffineTransformation)(l::Line)
  Line(f(l.base), l.angle + angle(f.a))
end

function (f::InversionReflection)(l::Line)
  if l.base ≈ 0
    return Line(l.base, l.angle + angle(f.b))
  else
    a = angle(l.base)
    if a ≈ l.angle || a ≈ (l.angle + pi)
      return Line(l.base, l.angle + angle(f.b))
    end
  end

  v = exp(im*l.angle)
  w1 = f(l.base)
  w2 = f(l.base + v)
  w3 = f(l.base - v)
  Circle(circumcenter(w1,w2,w3), circumradius(w1,w2,w3))
end

function (f::MobiusTransformation)(l::LineSegment)
  if f.c ≈ 0.0
    return LineSegment( f(l.p1), f(l.p2) )
  elseif l.p1 ≈ -f.d/f.c
    # ToDo! Ray
  elseif l.p2 ≈ -f.d/f.c
    # ToDo! Ray
  else
    # ToDo! Two Rays!
    # ToDo! Line Segment!
  end

  w1 = f(l.z1)
  w2 = f(l.z2)
  w3 = f((l.z1+l.z2)/2)
  Arc(circumcenter(w1,w2,w3), w1, w2)
end

function (f::LinearTransformation)(l::LineSegment)
  LineSegment(f(l.p1), f(l.p2))
end

function (f::Translation)(l::LineSegment)
  LineSegment(f(l.p1), f(l.p2))
end

function (f::AffineTransformation)(l::LineSegment)
  LineSegment(f(l.p1), f(l.p2))
end

function (f::InversionReflection)(l::LineSegment)
  if l.p1 ≈ 0
    # ToDo! Ray
  elseif l.p2 ≈ 0
    # ToDo! Ray
  else
    # ToDo! Two Rays!
    # ToDo! Line Segment!
  end

  w1 = f(l.z1)
  w2 = f(l.z2)
  w3 = f((l.z1+l.z2)/2)
  Arc(circumcenter(w1,w2,w3), w1, w2)
end


text:
using MooreJuliaMultigrid
using LinearAlgebra
using Test

@testset "MooreJuliaMultigrid.jl" begin
    # Matrix system to test
    n = 128
    h = 1/(n-1)
    A = Tridiagonal(-1*ones(n-2),2*ones(n-1),-1*ones(n-2))/h^2
    b = zeros(n-1)
    x0 = [.2*sin(3*pi*x/n)+.7*sin(4*pi*x/n)+.4*sin(10*pi*x/n) for x = 1:(n-1)]
    prolongOp = BuildLinearInterpOp(n-1)
    restrictOp = 0.5*prolongOp'

    @test norm(Jacobi(A,b,x0,numiters=10,weight=2/3)) < 10^(-2)
    @test norm(BuildLinearInterpOp(7) - 0.5*[1 0 0;
                                         2 0 0;
                                         1 1 0;
                                         0 2 0;
                                         0 1 1;
                                         0 0 2;
                                         0 0 1]) == 0
    @test norm(Vcycle(A,b,x0,Jacobi,(numiters=2,weight=2/3),(numiters=2,weight=2/3),prolongOp,restrictOp)) < 10^(-2)

end


text:
# This file is for the attractor of Sierpinski

using FractalTools
using Plots
using BenchmarkTools
using Plots
# Construct ifs
ifs = Sierpinski()
initset = [rand(2)]
atr = attractor(ifs, initset, alg=RandAlg(), numiter=Int(1e4), allocated=true)

# Plot the attractor
scatter(getindex.(atr.set,1), getindex.(atr.set,2), markersize=2)

# ch = Channel(0)
# task = @async FractalTools._randalg_sequential(ch, only(initset), ifs.ws, 1000, ifs.probs)


text:
## %PREDICTMD_GENERATED_BY%

using PredictMDExtra
PredictMDExtra.import_all()

using PredictMD
PredictMD.import_all()

using CSVFiles
using CategoricalArrays
using DataFrames
using DecisionTree
using Distributions
using FileIO
using GLM
using IterTools
using Knet
using LIBSVM
using LinearAlgebra
using PredictMD
using PredictMDAPI
using PredictMDExtra
using RDatasets
using Random
using StatsModels
using Test
using Unitful

const Schema = StatsModels.Schema

# PREDICTMD IF INCLUDE TEST STATEMENTS
logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
if isnothing(logger)
    logger_stream = devnull
else
    logger_stream = logger.stream
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### Begin project-specific settings

DIRECTORY_CONTAINING_THIS_FILE = @__DIR__
PROJECT_DIRECTORY = dirname(
    joinpath(splitpath(DIRECTORY_CONTAINING_THIS_FILE)...)
    )
PROJECT_OUTPUT_DIRECTORY = joinpath(
    PROJECT_DIRECTORY,
    "output",
    )
mkpath(PROJECT_OUTPUT_DIRECTORY)
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "data"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "models"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "plots"))

# PREDICTMD IF INCLUDE TEST STATEMENTS
if PredictMD.is_travis_ci()
    PredictMD.cache_to_path!(
        ;
        from = ["cpu_examples", "breast_cancer_biopsy", "output",],
        to = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

# PREDICTMD IF INCLUDE TEST STATEMENTS
@debug("PROJECT_OUTPUT_DIRECTORY: ", PROJECT_OUTPUT_DIRECTORY,)
if PredictMD.is_travis_ci()
    PredictMD.cache_to_path!(
        ;
        from = ["cpu_examples", "breast_cancer_biopsy", "output",],
        to = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### End project-specific settings

### Begin model output code

Kernel = LIBSVM.Kernel

Random.seed!(999)

trainingandtuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_features_df.csv",
    )
trainingandtuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_labels_df.csv",
    )
testing_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_features_df.csv",
    )
testing_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_labels_df.csv",
    )
training_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_features_df.csv",
    )
training_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_labels_df.csv",
    )
tuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_features_df.csv",
    )
tuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_labels_df.csv",
    )
trainingandtuning_features_df = DataFrames.DataFrame(
    FileIO.load(
        trainingandtuning_features_df_filename;
        type_detect_rows = 100,
        )
    )
trainingandtuning_labels_df = DataFrames.DataFrame(
    FileIO.load(
        trainingandtuning_labels_df_filename;
        type_detect_rows = 100,
        )
    )
testing_features_df = DataFrames.DataFrame(
    FileIO.load(
        testing_features_df_filename;
        type_detect_rows = 100,
        )
    )
testing_labels_df = DataFrames.DataFrame(
    FileIO.load(
        testing_labels_df_filename;
        type_detect_rows = 100,
        )
    )
training_features_df = DataFrames.DataFrame(
    FileIO.load(
        training_features_df_filename;
        type_detect_rows = 100,
        )
    )
training_labels_df = DataFrames.DataFrame(
    FileIO.load(
        training_labels_df_filename;
        type_detect_rows = 100,
        )
    )
tuning_features_df = DataFrames.DataFrame(
    FileIO.load(
        tuning_features_df_filename;
        type_detect_rows = 100,
        )
    )
tuning_labels_df = DataFrames.DataFrame(
    FileIO.load(
        tuning_labels_df_filename;
        type_detect_rows = 100,
        )
    )

smoted_training_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "smoted_training_features_df.csv",
    )
smoted_training_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "smoted_training_labels_df.csv",
    )
smoted_training_features_df = DataFrames.DataFrame(
    FileIO.load(
        smoted_training_features_df_filename;
        type_detect_rows = 100,
        )
    )
smoted_training_labels_df = DataFrames.DataFrame(
    FileIO.load(
        smoted_training_labels_df_filename;
        type_detect_rows = 100,
        )
    )

logistic_classifier_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "logistic_classifier.jld2",
    )
random_forest_classifier_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "random_forest_classifier.jld2",
    )
c_svc_svm_classifier_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "c_svc_svm_classifier.jld2",
    )
nu_svc_svm_classifier_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "nu_svc_svm_classifier.jld2",
    )
knet_mlp_classifier_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "knet_mlp_classifier.jld2",
    )

# PREDICTMD IF INCLUDE TEST STATEMENTS
logistic_classifier = nothing
Test.@test isnothing(logistic_classifier)

random_forest_classifier = nothing
Test.@test isnothing(random_forest_classifier)

c_svc_svm_classifier = nothing
Test.@test isnothing(c_svc_svm_classifier)

nu_svc_svm_classifier = nothing
Test.@test isnothing(nu_svc_svm_classifier)

knet_mlp_classifier = nothing
Test.@test isnothing(knet_mlp_classifier)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

logistic_classifier =
    PredictMD.load_model(logistic_classifier_filename)
random_forest_classifier =
    PredictMD.load_model(random_forest_classifier_filename)
c_svc_svm_classifier =
    PredictMD.load_model(c_svc_svm_classifier_filename)
nu_svc_svm_classifier =
    PredictMD.load_model(nu_svc_svm_classifier_filename)
knet_mlp_classifier =
    PredictMD.load_model(knet_mlp_classifier_filename)

PredictMD.parse_functions!(logistic_classifier)
PredictMD.parse_functions!(random_forest_classifier)
PredictMD.parse_functions!(c_svc_svm_classifier)
PredictMD.parse_functions!(nu_svc_svm_classifier)
PredictMD.parse_functions!(knet_mlp_classifier)

# PREDICTMD IF INCLUDE TEST STATEMENTS
PredictMD.parse_functions!(logistic_classifier)
PredictMD.parse_functions!(random_forest_classifier)
PredictMD.parse_functions!(c_svc_svm_classifier)
PredictMD.parse_functions!(nu_svc_svm_classifier)
PredictMD.parse_functions!(knet_mlp_classifier)

PredictMD.parse_functions!(logistic_classifier)
PredictMD.parse_functions!(random_forest_classifier)
PredictMD.parse_functions!(c_svc_svm_classifier)
PredictMD.parse_functions!(nu_svc_svm_classifier)
PredictMD.parse_functions!(knet_mlp_classifier)

PredictMD.parse_functions!(logistic_classifier)
PredictMD.parse_functions!(random_forest_classifier)
PredictMD.parse_functions!(c_svc_svm_classifier)
PredictMD.parse_functions!(nu_svc_svm_classifier)
PredictMD.parse_functions!(knet_mlp_classifier)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

PredictMD.predict_proba(logistic_classifier, smoted_training_features_df)
PredictMD.predict_proba(random_forest_classifier, smoted_training_features_df)
PredictMD.predict_proba(c_svc_svm_classifier, smoted_training_features_df)
PredictMD.predict_proba(nu_svc_svm_classifier, smoted_training_features_df)
PredictMD.predict_proba(knet_mlp_classifier, smoted_training_features_df)

PredictMD.predict_proba(logistic_classifier,testing_features_df)
PredictMD.predict_proba(random_forest_classifier,testing_features_df)
PredictMD.predict_proba(c_svc_svm_classifier,testing_features_df)
PredictMD.predict_proba(nu_svc_svm_classifier,testing_features_df)
PredictMD.predict_proba(knet_mlp_classifier,testing_features_df)

PredictMD.predict(logistic_classifier,smoted_training_features_df)
PredictMD.predict(random_forest_classifier,smoted_training_features_df)
PredictMD.predict(c_svc_svm_classifier,smoted_training_features_df)
PredictMD.predict(nu_svc_svm_classifier,smoted_training_features_df)
PredictMD.predict(knet_mlp_classifier,smoted_training_features_df)

PredictMD.predict(logistic_classifier,testing_features_df)
PredictMD.predict(random_forest_classifier,testing_features_df)
PredictMD.predict(c_svc_svm_classifier,testing_features_df)
PredictMD.predict(nu_svc_svm_classifier,testing_features_df)
PredictMD.predict(knet_mlp_classifier,testing_features_df)

single_label_name = :Class
negative_class = "benign"
positive_class = "malignant"

PredictMD.predict(logistic_classifier,smoted_training_features_df, positive_class, 0.3)
PredictMD.predict(random_forest_classifier,smoted_training_features_df, positive_class, 0.3)
PredictMD.predict(c_svc_svm_classifier,smoted_training_features_df, positive_class, 0.3)
PredictMD.predict(nu_svc_svm_classifier,smoted_training_features_df, positive_class, 0.3)
PredictMD.predict(knet_mlp_classifier,smoted_training_features_df, positive_class, 0.3)

PredictMD.predict(logistic_classifier,testing_features_df, positive_class, 0.3)
PredictMD.predict(random_forest_classifier,testing_features_df, positive_class, 0.3)
PredictMD.predict(c_svc_svm_classifier,testing_features_df, positive_class, 0.3)
PredictMD.predict(nu_svc_svm_classifier,testing_features_df, positive_class, 0.3)
PredictMD.predict(knet_mlp_classifier,testing_features_df, positive_class, 0.3)

### End model output code

# PREDICTMD IF INCLUDE TEST STATEMENTS
PredictMD.get_underlying(logistic_classifier)
PredictMD.get_underlying(random_forest_classifier)
PredictMD.get_underlying(c_svc_svm_classifier)
PredictMD.get_underlying(nu_svc_svm_classifier)
PredictMD.get_underlying(knet_mlp_classifier)

PredictMD.get_history(logistic_classifier)
PredictMD.get_history(random_forest_classifier)
PredictMD.get_history(c_svc_svm_classifier)
PredictMD.get_history(nu_svc_svm_classifier)
PredictMD.get_history(knet_mlp_classifier)

PredictMD.parse_functions!(logistic_classifier)
PredictMD.parse_functions!(random_forest_classifier)
PredictMD.parse_functions!(c_svc_svm_classifier)
PredictMD.parse_functions!(nu_svc_svm_classifier)
PredictMD.parse_functions!(knet_mlp_classifier)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

# PREDICTMD IF INCLUDE TEST STATEMENTS
if PredictMD.is_travis_ci()
    PredictMD.path_to_cache!(
        ;
        to = ["cpu_examples", "breast_cancer_biopsy", "output",],
        from = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

## %PREDICTMD_GENERATED_BY%