text:
using Markdown
"""
    struct Element
        atomic_number::Int
        symbol_name::String
        chinese_name::String
        english_name::String
        latin_name::String
        pinyin_name::String
    end # struct Element
元素结构体，包含元素的原子序数，元素符号，中文名，英文名，拉丁文名，拼音名。

构造函数：
```julia
Element(Z, sym, chinese, english, latin, pinyin)
Element(ele::Element)                         # identity
Element(Z::Integer)                           # 从原子序数构造
Element(name::Union{AbstractString, Symbol})  # 从任何名称构造，比如 "Fe", "铝", "Oxygen"
```
"""
struct Element
    atomic_number::Int
    symbol_name::String
    chinese_name::String
    english_name::String
    latin_name::String
    pinyin_name::String
    Element(Z, sym, chinese, english, latin, pinyin) = new(Z, sym, chinese, english, latin, pinyin)
end # struct Element

"不存在的元素"
const NoneElement = Element(-1, "", "", "", "", "")

"判断是否为不存在的元素"
is_none(ele::Element)::Bool = !(1 <= ele.atomic_number <= 118)

"元素周期表"
const element_table = [
    Element(1, "H", "氢", "Hydrogen", "Hydrogenium", "qīng"),
    Element(2, "He", "氦", "Helium", "Helium", "hài"),
    Element(3, "Li", "锂", "Lithium", "Lithium", "lǐ"),
    Element(4, "Be", "铍", "Beryllium", "Beryllium", "pí"),
    Element(5, "B", "硼", "Boron", "Borum", "péng"),
    Element(6, "C", "碳", "Carbon", "Сarbonium (Carboneum)", "tàn"),
    Element(7, "N", "氮", "Nitrogen", "Nitrogenium", "dàn"),
    Element(8, "O", "氧", "Oxygen", "Oxygenium", "yǎng"),
    Element(9, "F", "氟", "Fluorine", "Fluorum", "fú"),
    Element(10, "Ne", "氖", "Neon", "Neon", "nǎi"),
    Element(11, "Na", "钠", "Sodium", "Natrium", "nà"),
    Element(12, "Mg", "镁", "Magnesium", "Magnesium", "měi"),
    Element(13, "Al", "铝", "Aluminum", "Aluminium", "lǚ"),
    Element(14, "Si", "硅", "Silicon", "Silicium", "guī"),
    Element(15, "P", "磷", "Phosphorus", "Phosphorus", "lín"),
    Element(16, "S", "硫", "Sulfur", "Sulphuris", "liú"),
    Element(17, "Cl", "氯", "Chlorine", "Сhlorum", "lǜ"),
    Element(18, "Ar", "氩", "Argon", "Argon", "yà"),
    Element(19, "K", "钾", "Potassium", "Kalium", "jiǎ"),
    Element(20, "Ca", "钙", "Calcium", "Сalcium", "gài"),
    Element(21, "Sc", "钪", "Scandium", "Scandium", "kàng"),
    Element(22, "Ti", "钛", "Titanium", "Titanium", "tài"),
    Element(23, "V", "钒", "Vanadium", "Vanadium", "fán"),
    Element(24, "Cr", "铬", "Chromium", "Chromium", "gè"),
    Element(25, "Mn", "锰", "Manganese", "Manganum", "měng"),
    Element(26, "Fe", "铁", "Iron", "Ferrum", "tiě"),
    Element(27, "Co", "钴", "Cobalt", "Cobaltum", "gǔ"),
    Element(28, "Ni", "镍", "Nickel", "Niccolum", "niè"),
    Element(29, "Cu", "铜", "Copper", "Cuprum", "tóng"),
    Element(30, "Zn", "锌", "Zinc", "Zincum", "xīn"),
    Element(31, "Ga", "镓", "Gallium", "Gallium", "jiā"),
    Element(32, "Ge", "锗", "Germanium", "Germanium", "zhě"),
    Element(33, "As", "砷", "Arsenic", "Arsenicum", "shēn"),
    Element(34, "Se", "硒", "Selenium", "Selenium", "xī"),
    Element(35, "Br", "溴", "Bromine", "Bromum", "xiù"),
    Element(36, "Kr", "氪", "Krypton", "Krypton", "kè"),
    Element(37, "Rb", "铷", "Rubidium", "Rubidium", "rú"),
    Element(38, "Sr", "锶", "Strontium", "Strontium", "sī"),
    Element(39, "Y", "钇", "Yttrium", "Yttrium", "yǐ"),
    Element(40, "Zr", "锆", "Zirconium", "Zirconium", "gào"),
    Element(41, "Nb", "铌", "Niobium", "Niobium", "ní"),
    Element(42, "Mo", "钼", "Molybdenum", "Molybdaenum", "mù"),
    Element(43, "Tc", "锝", "Technetium", "Technetium", "dé"),
    Element(44, "Ru", "钌", "Ruthenium", "Ruthenium", "liǎo"),
    Element(45, "Rh", "铑", "Rhodium", "Rhodium", "lǎo"),
    Element(46, "Pd", "钯", "Palladium", "Palladium", "bǎ"),
    Element(47, "Ag", "银", "Silver", "Argentum", "yín"),
    Element(48, "Cd", "镉", "Cadmium", "Cadmium", "gé"),
    Element(49, "In", "铟", "Indium", "Indium", "yīn"),
    Element(50, "Sn", "锡", "Tin", "Stannum", "xī"),
    Element(51, "Sb", "锑", "Antimony", "Stibium", "tī"),
    Element(52, "Te", "碲", "Tellurium", "Tellurium", "dì"),
    Element(53, "I", "碘", "Iodine", "Iodium", "diǎn"),
    Element(54, "Xe", "氙", "Xenon", "Xenon", "xiān"),
    Element(55, "Cs", "铯", "Caesium", "Caesium", "sè"),
    Element(56, "Ba", "钡", "Barium", "Barium", "bèi"),
    Element(57, "La", "镧", "Lanthanum", "Lanthanum", "lán"),
    Element(58, "Ce", "铈", "Cerium", "Cerium", "shì"),
    Element(59, "Pr", "镨", "Praseodymium", "Praseodymium", "pǔ"),
    Element(60, "Nd", "钕", "Neodymium", "Neodymium", "nǚ"),
    Element(61, "Pm", "钷", "Promethium", "Promethium", "pǒ"),
    Element(62, "Sm", "钐", "Samarium", "Samarium", "shān"),
    Element(63, "Eu", "铕", "Europium", "Europium", "yǒu"),
    Element(64, "Gd", "钆", "Gadolinium", "Gadolinium", "gá"),
    Element(65, "Tb", "铽", "Terbium", "Terbium", "tè"),
    Element(66, "Dy", "镝", "Dysprosium", "Dysprosium", "dī"),
    Element(67, "Ho", "钬", "Holmium", "Holmium", "huǒ"),
    Element(68, "Er", "铒", "Erbium", "Erbium", "ěr"),
    Element(69, "Tm", "铥", "Thulium", "Thulium", "diū"),
    Element(70, "Yb", "镱", "Ytterbium", "Ytterbium", "yì"),
    Element(71, "Lu", "镥", "Lutetium", "Lutetium", "lǔ"),
    Element(72, "Hf", "铪", "Hafnium", "Hafnium", "hā"),
    Element(73, "Ta", "钽", "Tantalum", "Tantalum", "tǎn"),
    Element(74, "W", "钨", "Tungsten", "Wolframium", "wū"),
    Element(75, "Re", "铼", "Rhenium", "Rhenium", "lái"),
    Element(76, "Os", "锇", "Osmium", "Osmium", "é"),
    Element(77, "Ir", "铱", "Iridium", "Iridium", "yī"),
    Element(78, "Pt", "铂", "Platinum", "Platinum", "bó"),
    Element(79, "Au", "金", "Gold", "Aurum", "jīn"),
    Element(80, "Hg", "汞", "Mercury", "Hydrargyrum", "gǒng"),
    Element(81, "Tl", "铊", "Thallium", "Thallium", "tā"),
    Element(82, "Pb", "铅", "Lead", "Plumbum", "qiān"),
    Element(83, "Bi", "铋", "Bismuth", "Bisemutum (Bismuthum, Bismutum)", "bì"),
    Element(84, "Po", "钋", "Polonium", "Polonium", "pō"),
    Element(85, "At", "砹", "Astatine", "Astatum", "ài"),
    Element(86, "Rn", "氡", "Radon", "Radon", "dōng"),
    Element(87, "Fr", "钫", "Francium", "Francium", "fāng"),
    Element(88, "Ra", "镭", "Radium", "Radium", "léi"),
    Element(89, "Ac", "锕", "Actinium", "Actinium", "ā"),
    Element(90, "Th", "钍", "Thorium", "Thorium", "tǔ"),
    Element(91, "Pa", "镤", "Protactinium", "Protactinium", "pú"),
    Element(92, "U", "铀", "Uranium", "Uranium", "yóu"),
    Element(93, "Np", "镎", "Neptunium", "Neptunium", "ná"),
    Element(94, "Pu", "钚", "Plutonium", "Plutonium", "bù"),
    Element(95, "Am", "镅", "Americium", "Americium", "méi"),
    Element(96, "Cm", "锔", "Curium", "Curium", "jú"),
    Element(97, "Bk", "锫", "Berkelium", "Berkelium", "péi"),
    Element(98, "Cf", "锎", "Californium", "Californium", "kāi"),
    Element(99, "Es", "锿", "Einsteinium", "Einsteinium", "āi"),
    Element(100, "Fm", "镄", "Fermium", "Fermium", "fèi"),
    Element(101, "Md", "钔", "Mendelevium", "Mendelevium", "mén"),
    Element(102, "No", "锘", "Nobelium", "Nobelium", "nuò"),
    Element(103, "Lr", "铹", "Lawrencium", "Lawrencium", "láo"),
    Element(104, "Rf", "𬬻", "Rutherfordium", "Rutherfordium", "lú"),
    Element(105, "Db", "𬭊", "Dubnium", "Dubnium", "dù"),
    Element(106, "Sg", "𬭳", "Seaborgium", "Seaborgium", "xǐ"),
    Element(107, "Bh", "𬭛", "Bohrium", "Bohrium", "bō"),
    Element(108, "Hs", "𬭶", "Hassium", "Hassium", "hēi"),
    Element(109, "Mt", "鿏", "Meitnerium", "Meitnerium", "mài"),
    Element(110, "Ds", "𫟼", "Darmstadtium", "Darmstadtium", "dá"),
    Element(111, "Rg", "𬬭", "Roentgenium", "Roentgenium", "lún"),
    Element(112, "Cn", "鎶", "Copernicium", "Copernecium", "gē"),
    Element(113, "Nh", "鉨", "Nihonium", "Nihonium", "nǐ"),
    Element(114, "Fl", "𫓧", "Flerovium", "Flerovium", "fū"),
    Element(115, "Mc", "镆", "Moscovium", "Moscovium", "mò"),
    Element(116, "Lv", "鉝", "Livermorium", "Livermorium", "lì"),
    Element(117, "Ts", "石田", "Tennessine", "Tennessine", "tián"),
    Element(118, "Og", "气奥", "Oganesson", "Oganneson", "ào")
]

"通过原子序数查找元素"
find_element_with_Z(Z::Integer) = (1 <= Z <= 118) ? element_table[Z] : NoneElement

"通过元素符号查找元素"
function find_element_with_symbol(sym::AbstractString)
    pos = findfirst(ele -> ele.symbol_name == sym, element_table)
    pos === nothing ? NoneElement : element_table[pos]
end

"通过元素中文名查找元素"
function find_element_with_chinese(chinese::AbstractString)
    pos = findfirst(ele -> ele.chinese_name == chinese, element_table)
    pos === nothing ? NoneElement : element_table[pos]
end

"通过元素英文名查找元素"
function find_element_with_english(english::AbstractString)
    pos = findfirst(ele -> ele.english_name == english, element_table)
    pos === nothing ? NoneElement : element_table[pos]
end

"通过元素拉丁文名查找元素"
function find_element_with_latin(latin::AbstractString)
    pos = findfirst(ele -> match(Regex("\\b" * latin * "\\b"), ele.latin_name) !== nothing, element_table)
    pos === nothing ? NoneElement : element_table[pos]
end

"通过元素拼音查找元素"
function find_element_with_pinyin(pinyin::AbstractString)
    pos = findfirst(ele -> ele.pinyin_name == pinyin, element_table)
    pos === nothing ? NoneElement : element_table[pos]
end

"""
    find_element(name::Union{Integer, AbstractString})
通过任何名称查找元素，包括使用原子序数
"""
function find_element(name::Union{Integer, AbstractString})
    name isa Integer && return find_element_with_Z(name)
    ele = find_element_with_symbol(name)
    is_none(ele) && (ele = find_element_with_chinese(name))
    is_none(ele) && (ele = find_element_with_english(name))
    is_none(ele) && (ele = find_element_with_latin(name))
    is_none(ele) && (ele = find_element_with_pinyin(name))
    return ele
end

Element(ele::Element) = ele
Element(Z::Integer) = find_element_with_Z(Z)
Element(name::Union{AbstractString, Symbol}) = find_element(string(name))

const ElementConstructType = Union{Element, Integer, AbstractString, Symbol}
"获取原子序数"
getZ(ele::ElementConstructType) = Element(ele).atomic_number
"获取原子序数"
atomic_number(ele::ElementConstructType) = Element(ele).atomic_number
"获取元素符号"
symbol(ele::ElementConstructType) = Element(ele).symbol_name
"获取元素中文名"
chinese(ele::ElementConstructType) = Element(ele).chinese_name
"获取元素英文名"
english(ele::ElementConstructType) = Element(ele).english_name
"获取元素拉丁文名"
latin(ele::ElementConstructType) = Element(ele).latin_name
"获取元素中文名拼音"
pinyin(ele::ElementConstructType) = Element(ele).pinyin_name

"""
    show_element_table(type::AbstractString="symbol")
打印元素周期表，`type`可以是`"symbol", "chinese"`两种
"""
function show_element_table(type::AbstractString="symbol")
    choose_field = nothing
    use_space = ""
    if type == "symbol"
        choose_field = symbol
        use_space = "  "
    elseif type == "chinese"
        choose_field = chinese
        use_space = "\u3000"
    else
        throw(ArgumentError("invalid element table type: $type"))
    end
    # 第一行
    str = @sprintf("%-2s%s %-2s\n", choose_field(Element(1)), repeat(" " * use_space, 16), choose_field(Element(2)))
    # 第二、三行
    for line = 2:3
        start = (line - 2) * 8 + 3
        for z = start:start+1
            str = str * @sprintf("%-2s ", choose_field(Element(z)))
        end
        str = str * repeat(use_space * " ", 10)
        for z = start+2:start+6
            str = str * @sprintf("%-2s ", choose_field(Element(z)))
        end
        str = str * @sprintf("%-2s\n", choose_field(Element(start+7)))
    end
    # 第四、五行
    for line = 4:5
        start = (line - 4) * 18 + 19
        for z = start:start+16
            str = str * @sprintf("%-2s ", choose_field(Element(z)))
        end
        str = str * @sprintf("%-2s\n", choose_field(Element(start + 17)))
    end
    # 第六、七行
    for line = 6:7
        start = (line - 6) * 32 + 55
        for z = start:start+2
            str = str * @sprintf("%-2s ", choose_field(Element(z)))
        end
        for z = start+17:start+30
            str = str * @sprintf("%-2s ", choose_field(Element(z)))
        end
        str = str * @sprintf("%-2s\n", choose_field(Element(start + 31)))
    end
    str = str * repeat(use_space * " ", 17) * use_space * "\n"
    # 附加行
    for line = 6:7
        str = str * repeat(use_space * " ", 2)
        start = (line - 6) * 32 + 57
        for z = start:start+14
            str = str * @sprintf("%-2s ", choose_field(Element(z)))
        end
        str = str * repeat(use_space * " ", 7)
        str = str * use_space * "\n"
    end
    print(str)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule TempestModel_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("TempestModel")
JLLWrappers.@generate_main_file("TempestModel", UUID("2962c967-6e2c-5a30-9da7-7ff16ffb0b2a"))
end  # module TempestModel_jll

ext: jl


text:
module TestHCEstimator
using Test
using JuMP, Ipopt
using DataFrames, CSV

include("../src/HCEstimator.jl")
using .HCEstimator
include("../src/Estimator.jl")
using .Estimator

include("cases_no_v.jl")
using .Cases

function runtests()

    @testset "HCEstimator" begin
        for case in [case3_dist, case33, case3_dist_no_dgs]
            sys, name = case()
            @testset "$name" begin
                other = nl_pf(Model(), sys)
                other_w_obj = add_objective(other, sys)

                sot = build_model(Model(), sys)
                @test sprint(print, sot) == sprint(print, other_w_obj)

                other_sot = build_model(Model(), sys, false)
                @test sprint(print, sot) == sprint(print, other)
            end

        end
    end
end

end
TestHCEstimator.runtests()


text:
using AnalyticalEOR

swr = 0.2
sor = 0.2

krw0 = 0.5
kro0 = 1.0
nw = 3.0
no = 1.0

kr = RelPerms(swr, sor, krw0, kro0, nw, no)

μw = 1.0
μo = 3.0

si = 0.20
sj = 0.80

wf = solve_waterflooding(si, sj, kr, μw, μo)

wf.S̃

plot_fw(wf)

plot_sw_profile(wf, 0.1, 0.2, 0.3, 0.5)

animate_sw_profile(wf)

max_stars_repo_path: examples/Ex1 - Waterflooding.jl

max_issues_repo_path: examples/Ex1 - Waterflooding.jl

max_forks_repo_path: examples/Ex1 - Waterflooding.jl


text:
@export ϕtoθ(ϕ) = π - 2atan(√3 * cos(ϕ))
@export ϕtoδz(ϕ, l=DEFAULT_LENGTH) = √3/2 * l * sin(ϕ)

abstract type AbstractMoleculeConstructor end

"""
    Cumulene(N; θ=π/2*(N % 2), l=DEFAULT_LENGTH, orbitals=(orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z"))

Molecule constructor for Cumulene with relative twist angle, `θ`, of the
hydrogen positions at either end.  The orbitals of the Carbon atoms
can be restricted by providing an iterable of `orbitals` as a keyword
argument.

"""
@export struct Cumulene <: AbstractMoleculeConstructor
    N
    θ
    l_ch
    l_cc
    orbitals
    valence
    left_index::Int
    right_index::Int

    function Cumulene(
        ;
        N = 5,
        θ = π/2*(N % 2),
        l_ch = CH_BOND_LENGTH,
        l_cc = CC_TRIPLE_BOND_LENGTH,
        orbitals = (orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z"),
        valence = missing,
        left_index = 3,
        right_index = 3
    )
        return new(N, θ, l_ch, l_cc, orbitals, valence, left_index, right_index)
    end
end

function (m::Cumulene)(; kwargs...)
    N = get(kwargs, :N, m.N)
    θ = get(kwargs, :θ, m.θ)
    l_ch = get(kwargs, :l_ch, m.l_ch)
    l_cc = get(kwargs, :l_cc, m.l_cc)
    orbitals = get(kwargs, :orbitals, m.orbitals)
    valence = get(kwargs, :valence, m.valence)

    mol = Molecule()

    r0 = [0.0, 0.0, 0.0]
    r1 = [-√3/2, 0, -0.5] * l_ch
    r2 = [√3/2, 0, -0.5] * l_ch
    axes = makeaxes(0.0, 0.0)

    A = Hydrogen(r1, axes=axes)
    B = Hydrogen(r2, axes=axes)

    C = Carbon(r0, orbitals..., axes=axes, valence=valence)

    push!(mol, A, B, C)
    push!(mol, Bond(A, C), Bond(B, C))

    C′ = C
    for i in 1:(N-1)
        axes = makeaxes(N > 1 ? θ*i/(N-1) : 0.0, 0.0)
        r = [0, 0, i*l_cc]
        C = Carbon(r, orbitals..., axes=axes, valence=valence)
        push!(mol, C)
        i > 0 && push!(mol, Bond(C, C′))
        C′ = C
    end

    r1 = [-√3/2 * cos(θ), -√3/2 * sin(θ), (N - 0.5)] * l_ch
    r2 = [√3/2 * cos(θ), √3/2 * sin(θ), (N - 0.5)] * l_ch
    axes = makeaxes(θ, 0.0)

    A = Hydrogen(r1, axes=axes)
    B = Hydrogen(r2, axes=axes)

    push!(mol, A, B)
    push!(mol, Bond(A, C′), Bond(B, C′))
    return mol
end

"""
    Polyacetylene(N; ϕ=π/4, l=DEFAULT_LENGTH, handedness=1, δz=nothing, θ=nothing, r0=nothing)

Molecule constructor for (twisted) Polyacetylene with twist angle, `ϕ`.

"""
@export struct Polyacetylene <: AbstractMoleculeConstructor
    N
    θ
    l_ch
    l_cc
    handedness
    δz
    ϕ
    r0
    left_index::Int
    right_index::Int

    function Polyacetylene(
        ;
        N = 12,
        θ = π/2,
        l_ch = CH_BOND_LENGTH,
        l_cc = (CC_SINGLE_BOND_LENGTH + CC_DOUBLE_BOND_LENGTH)/2,
        handedness = 1,
        δz = nothing,
        ϕ = nothing,
        r0 = nothing,
        left_index = 2,
        right_index = 3
    )
        return new(N, θ, l_ch, l_cc, handedness, δz, ϕ, r0, left_index, right_index)
    end
end

function (m::Polyacetylene)(; kwargs...)
    N = get(kwargs, :N, m.N)
    θ = get(kwargs, :θ, m.θ)
    l_ch = get(kwargs, :l_ch, m.l_ch)
    l_cc = get(kwargs, :l_cc, m.l_cc)
    handedness = get(kwargs, :handedness, m.handedness)
    δz = get(kwargs, :δz, m.δz)
    ϕ = get(kwargs, :ϕ, m.ϕ)
    r0 = get(kwargs, :r0, m.r0)

    ρ = atan(2tan(θ/2)) + π * (θ > π)

    if δz !== nothing && θ !== nothing
        r0 = [l_cc, 0, 0]
    end

    if δz === nothing
        δz = √3*l_cc*sin(ρ)/2
    end
    if ϕ === nothing
        ϕ = π - 2atan(√3*cos(ρ))*handedness
    end

    if r0 === nothing
        r0 = [√(l_cc^2 - δz^2) / 2sin(ϕ/2), 0, 0]
    end
    r0::Vector

    mol = Molecule()

    rotation_matrix = LinearAlgebra.Givens(1, 2, cos(ϕ), -sin(ϕ))
    shift_vector = [0, 0, δz]
    r1 = shift_vector + rotation_matrix*r0

    axes = makeaxes(π/2, ρ*handedness)
    A′ = Hydrogen(r1 + (r0 - r1)*l_ch/l_cc, axes=axes)
    push!(mol, A′)

    for i in 1:N
        rotation_matrix = LinearAlgebra.Givens(1, 2, cos(i*ϕ), -sin(i*ϕ))
        shift_vector = [0, 0, i*δz]
        axes = makeaxes(π/2 + i*ϕ, ρ*handedness)
        r1 = shift_vector + rotation_matrix*r0
        r2 = shift_vector + rotation_matrix*(r0 + [l_ch, 0, 0])
        A = Carbon(r1, axes=axes)
        B = Hydrogen(r2, axes=axes)
        push!(mol, A, B)
        push!(mol, Bond(A, A′), Bond(A, B))
        A′ = A
    end
    rotation_matrix = LinearAlgebra.Givens(1, 2, cos((N+1)*ϕ), -sin((N+1)*ϕ))
    shift_vector = [0, 0, (N+1)*δz]
    axes = makeaxes(π/2 + (N+1)*ϕ, ρ*handedness)
    r3 = shift_vector + rotation_matrix*r0
    A = Hydrogen(r1 + (r3 - r1)*l_ch/l_cc, axes=axes)
    push!(mol, A)
    push!(mol, Bond(A, A′))
    return mol
end

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

"""
    Helicene(N; ϕ=π/4, l=DEFAULT_LENGTH, handedness=1, δz=nothing, θ=nothing, r0=nothing)

Molecule constructor for Helicene with twist angle, `ϕ`.

"""
@export struct Helicene <: AbstractMoleculeConstructor
    N
    l_ch
    l_cc
    δz
    left_index::Int
    right_index::Int

    function Helicene(
        ;
        N = 4,
        l_ch = CH_BOND_LENGTH,
        l_cc = (CC_SINGLE_BOND_LENGTH + CC_DOUBLE_BOND_LENGTH)/2,
        δz = 0.5,
        left_index = 4,
        right_index = 3
    )
        return new(N, l_ch, l_cc, δz, left_index, right_index)
    end
end

function (m::Helicene)(; kwargs...)
    N = get(kwargs, :N, m.N)
    δz = get(kwargs, :δz, m.δz)
    l_ch = get(kwargs, :l_ch, m.l_ch)
    l_cc = get(kwargs, :l_cc, m.l_cc)

    shift_vector = [0, 0, δz]
    θ = 2π/6
    ϕ = 0

    r0 = [l_cc, 0, 0]
    r_12 = [l_cc, 0, 0]
    r_23 = [l_cc/2, √3/2*l_cc, 0]
    r_34 = [-l_cc/2, √3/2*l_cc, 0]
    r_35 = [l_ch, 0, 0]
    r_46 = [l_ch/2, √3/2*l_ch, 0]

    r1 = r0
    r2 = r1 + r_12
    r3 = r2 + r_23
    r4 = r3 + r_34
    r5 = r3 + r_35
    r6 = r4 + r_46

    ϕ1 = atan(r1[2], r1[1])
    ϕ2 = atan(r2[2], r2[1])
    ϕ3 = atan(r3[2], r3[1])
    ϕ4 = atan(r4[2], r4[1])
    ϕ5 = atan(r5[2], r5[1])
    ϕ6 = atan(r6[2], r6[1])

    mol = Molecule()

    i = -1
    axes = makeaxes(π/2 + i*θ, 0)
    rotation_matrix = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ), -sin(i*θ + ϕ))

    r1_ = (i + 3ϕ1/π)*shift_vector + rotation_matrix*r1
    r4_ = (i + 3ϕ4/π)*shift_vector + rotation_matrix*r4
    A′ = Hydrogen(r1 + (r1_ - r1)*l_ch/l_cc, axes=axes)
    D′ = Hydrogen(r2 + (r4_ - r2)*l_ch/l_cc, axes=axes)
    push!(mol, A′)
    push!(mol, D′)
    for i in 0:(N-1)
        rotation_matrix = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ), -sin(i*θ + ϕ))
        axes = makeaxes(π/2 + i*θ, 0)

        r1_ = (i + 3ϕ1/π)*shift_vector + rotation_matrix*r1
        r2_ = (i + 3ϕ2/π)*shift_vector + rotation_matrix*r2
        r3_ = (i + 3ϕ3/π)*shift_vector + rotation_matrix*r3
        r4_ = (i + 3ϕ4/π)*shift_vector + rotation_matrix*r4
        r5_ = (i + 3ϕ5/π)*shift_vector + rotation_matrix*r5
        r6_ = (i + 3ϕ6/π)*shift_vector + rotation_matrix*r6

        A = Carbon(r1_, axes=axes)
        B = Carbon(r2_, axes=axes)
        C = Carbon(r3_, axes=axes)
        D = Carbon(r4_, axes=axes)
        E = Hydrogen(r5_, axes=axes)
        F = Hydrogen(r6_, axes=axes)
        push!(mol, A, B, C, D, E, F)
        push!(mol, Bond(A, A′), Bond(A, B), Bond(B, C), Bond(C, D), Bond(B, D′), Bond(C, E), Bond(D, F))
        A′ = A
        D′ = D
    end
    i = N
    rotation_matrix = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ), -sin(i*θ + ϕ))
    axes = makeaxes(π/2 + i*θ, 0)

    r1_ = (i + 3ϕ1/π)*shift_vector + rotation_matrix*r1
    r2_ = (i + 3ϕ2/π)*shift_vector + rotation_matrix*r2
    r3_ = (i + 3ϕ3/π)*shift_vector + rotation_matrix*r3

    A = Carbon(r1_, axes=axes)
    B = Carbon(r2_, axes=axes)
    C = Hydrogen(r2_ + (r3_ - r2_)*l_ch/l_cc, axes=axes)

    i = N + 1
    axes = makeaxes(π/2 + i*θ, 0)

    rotation_matrix = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ), -sin(i*θ + ϕ))
    r7_ = (i + 3ϕ1/π)*shift_vector + rotation_matrix*r1

    D = Hydrogen(r1_ + (r7_ - r1_)*l_ch/l_cc, axes=axes)

    push!(mol, A, B, C, D)
    push!(mol, Bond(A, A′), Bond(A, B), Bond(B, C), Bond(B, D′), Bond(A, D))

    return mol
end

function makehelicene2(N, a=0.5, l=1.4; handedness=1)
    shift_vector = [0, 0, a]
    θ = 2π/6#π - 2atan(√3 * cos(ϕ))*handedness
    ϕ1 = 0
    ϕ2 = 0
    ϕ3 = atan(√3, 5)
    ϕ4 = atan(√3, 2)
    ϕ5 = atan(√3, 7)
    ϕ6 = atan(3*√3, 5)
    r0 = [l, 0, 0]

    mol = Molecule()

    axes = makeaxes(π/2, 0.0)
    i = 0
    rotation_matrix1 = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ1), -sin(i*θ + ϕ1))
    rotation_matrix4 = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ4), -sin(i*θ + ϕ4))
    rotation_matrix5 = LinearAlgebra.Givens(1, 2, cos((i-1)*θ + ϕ1), -sin((i-1)*θ + ϕ1))

    r1 = (i + 3ϕ1/π)*shift_vector + rotation_matrix1*r0
    r4 = (i + 3ϕ4/π)*shift_vector + rotation_matrix4*(√(5 + √3)*r0)
    r5 = ((i-1) + 3ϕ1/π)*shift_vector + rotation_matrix5*r0
    A′ = Oxygen(vec(r1...), orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z", axes=axes)
    D′ = Hydrogen(vec(r4...), orbital"1s", axes=axes)
    F′ = Hydrogen(vec(r5...), orbital"1s", axes=axes)
    push!(mol, A′, D′, F′)
    push!(mol, Bond(A′, F′))

    for i in 1:N
        rotation_matrix1 = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ1), -sin(i*θ + ϕ1))
        rotation_matrix2 = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ2), -sin(i*θ + ϕ2))
        rotation_matrix3 = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ3), -sin(i*θ + ϕ3))
        rotation_matrix4 = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ4), -sin(i*θ + ϕ4))
        rotation_matrix5 = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ5), -sin(i*θ + ϕ5))
        rotation_matrix6 = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ6), -sin(i*θ + ϕ6))
        axes = makeaxes(π/2 + i*θ, 0.0)

        r1 = (i + 3ϕ1/π)*shift_vector + rotation_matrix1*r0
        r2 = (i + 3ϕ2/π)*shift_vector + rotation_matrix2*(2r0)
        r3 = (i + 3ϕ3/π)*shift_vector + rotation_matrix3*(√7*r0)
        r4 = (i + 3ϕ4/π)*shift_vector + rotation_matrix4*(√7*r0)
        r5 = (i + 3ϕ5/π)*shift_vector + rotation_matrix5*(√13*r0)
        r6 = (i + 3ϕ6/π)*shift_vector + rotation_matrix6*(√13*r0)

        A = Carbon(vec(r1...), orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z", axes=axes)
        B = Carbon(vec(r2...), orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z", axes=axes)
        if i == 2
            C = Nitrogen(vec(r3...), orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z", axes=axes)
            F = Carbon(vec(r6...), orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z", axes=axes)
        else
            C = Carbon(vec(r3...), orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z", axes=axes)
            F = Hydrogen(vec(r6...), orbital"1s", axes=axes)
        end

        if i == 3
            D = Nitrogen(vec(r4...), orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z", axes=axes)
            E = Carbon(vec(r5...), orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z", axes=axes)
        else
            D = Carbon(vec(r4...), orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z", axes=axes)
            E = Hydrogen(vec(r5...), orbital"1s", axes=axes)
        end

        push!(mol, A, B, C, D, E, F)
        push!(mol, Bond(A, A′), Bond(A, B), Bond(B, C), Bond(C, D), Bond(B, D′), Bond(C, E), Bond(D, F))

        if i == 3
            ϕ9 = atan(√3, 4)
            ϕ10 = -ϕ9
            rotation_matrix9 = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ9), -sin(i*θ + ϕ9))
            rotation_matrix10 = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ10), -sin(i*θ + ϕ10))
            r7 = (i + 3ϕ1/π)*shift_vector + rotation_matrix1*(4r0)
            r8 = (i + 3ϕ1/π)*shift_vector + rotation_matrix1*(5r0)
            r9 = (i + 3ϕ9/π)*shift_vector + rotation_matrix9*(√19*r0)
            r10 = (i + 3ϕ10/π)*shift_vector + rotation_matrix10*(√19*r0)

            G = Carbon(vec(r7...), orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z", axes=axes)
            H = Hydrogen(vec(r8...), orbital"1s", axes=axes)

            I = Hydrogen(vec(r9...), orbital"1s", axes=axes)
            J = Hydrogen(vec(r10...), orbital"1s", axes=axes)
            push!(mol, G, H, I, J)
            push!(mol, Bond(E, G), Bond(F′, G), Bond(G, H), Bond(E, I), Bond(F′, J))
        end

        A′ = A
        D′ = D
        F′ = F
    end

    i = N + 1
    rotation_matrix1 = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ1), -sin(i*θ + ϕ1))
    rotation_matrix2 = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ2), -sin(i*θ + ϕ2))
    rotation_matrix3 = LinearAlgebra.Givens(1, 2, cos(i*θ + ϕ3), -sin(i*θ + ϕ3))
    rotation_matrix4 = LinearAlgebra.Givens(1, 2, cos((i+1)*θ + ϕ1), -sin((i+1)*θ + ϕ1))
    rotation_matrix5 = LinearAlgebra.Givens(1, 2, cos((i+2)*θ + ϕ1), -sin((i+2)*θ + ϕ1))
    axes = makeaxes(π/2 + i*θ, 0.0)

    r1 = (i + 3ϕ1/π)*shift_vector + rotation_matrix1*r0
    r2 = (i + 3ϕ2/π)*shift_vector + rotation_matrix2*(r0 + [abs(l), 0, 0])
    r3 = (i + 3ϕ3/π)*shift_vector + rotation_matrix3*([abs(l)*√7, 0, 0])
    r4 = ((i+1) + 3ϕ1/π)*shift_vector + rotation_matrix4*r0
    r5 = ((i+2) + 3ϕ1/π)*shift_vector + rotation_matrix5*r0

    A = Carbon(vec(r1...), orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z", axes=axes)
    B = Carbon(vec(r2...), orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z", axes=axes)
    C = Hydrogen(vec(r3...), orbital"1s", axes=axes)
    D = Oxygen(vec(r4...), orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z", axes=axes)
    E = Hydrogen(vec(r5...), orbital"1s", axes=axes)
    push!(mol, A, B, C, D, E)
    push!(mol, Bond(A, A′), Bond(A, B), Bond(B, C), Bond(B, D′), Bond(A, D), Bond(D, E))
    return mol
end


text:
using Revise
using AdFem
using PyPlot
using LinearAlgebra
using Statistics
using MAT
SMALL_SIZE = 20
MEDIUM_SIZE = 20
BIGGER_SIZE = 20

matplotlib.rc("font", size=SMALL_SIZE)          # controls default text sizes
matplotlib.rc("axes", titlesize=SMALL_SIZE)     # fontsize of the axes title
matplotlib.rc("axes", labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
matplotlib.rc("xtick", labelsize=SMALL_SIZE)    # fontsize of the tick labels
matplotlib.rc("ytick", labelsize=SMALL_SIZE)    # fontsize of the tick labels
matplotlib.rc("legend", fontsize=SMALL_SIZE)    # legend fontsize
matplotlib.rc("figure", titlesize=BIGGER_SIZE)  # fontsize of the figure title

function f(x, y)
    1/(1+x^2) + x * y + y^2
end

mmesh = Mesh(50, 50, 1/50, degree=2)

E = eval_f_on_gauss_pts(f, mmesh)

E2 = matread("fenics/data2-inverse10001.mat")["E"]
E3 = matread("pixel/data2-inverse981.mat")["E"]

close("all")
figure(figsize=(30, 8))
subplot(131)
title("Exact")
visualize_scalar_on_gauss_points(E, mmesh)
gca().set_rasterized(true)
subplot(132)
title("Difference (DNN)")
visualize_scalar_on_gauss_points(abs.(E-E2), mmesh)
xlabel("x")
subplot(133)
title("Difference (Discretization on Gauss Points)")
visualize_scalar_on_gauss_points(abs.(E-E3), mmesh)
xlabel("x")
tight_layout()
savefig("LinearElasticity.png")
savefig("LinearElasticity.pdf")


text:
module Helmholtz #hide
# # Helmholtz equation with fancy names
# In this example, we use the SLEPc to find the eigenvalues of the following Helmholtz equation:
# ``u'' + \omega^2 u = 0`` associated to Dirichlet boundary conditions on the domain ``[0,1]``. Hence
# the theoritical eigenvalues are ``\omega = k \pi`` with ``k \in \mathbb{Z}^*``; and the associated
# eigenvectors are ``u(x) = \sin(k\pix)``.
# A centered finite difference scheme is used for the spatial discretization.
#
# The equation is written in matrix form ``Au = \alpha Bu`` where ``\alpha = \omega^2``.
#
# To run this example, simplfy excute `mpirun -n your_favourite_integer julia helmholtz_FD.jl`
#
# In this example, "fancy" names are use to interface with PETSc/SLEPc. For legacy method names check the
# previous example.
#
# Note that the way we achieve things in the document can be highly improved and the purpose of this example
# is only demonstrate some method calls to give an overview.
#
# Start by importing both `PetscWrap`, for the distributed matrices, and `SlepcWrap` for the eigenvalues.
using PetscWrap
using SlepcWrap

# Number of mesh points and mesh step
n = 21
Δx = 1. / (n - 1)

# Initialize SLEPc. Either without arguments, calling `SlepcInitialize()` or using "command-line" arguments.
# To do so, either provide the arguments as one string, for instance
# `SlepcInitialize("-eps_max_it 100 -eps_tol 1e-5")` or provide each argument in
# separate strings : `PetscInitialize(["-eps_max_it", "100", "-eps_tol", "1e-5")`.
# Here we ask for the five closest eigenvalues to ``0``, using a non-zero pivot for the LU factorization and a
# "shift-inverse" process.
SlepcInitialize("-eps_target 0 -eps_nev 5 -st_pc_factor_shift_type NONZERO -st_type sinvert")

# Create the problem matrices, set sizes and apply "command-line" options. Note that we should
# set the number of preallocated non-zeros to increase performance.
A = create_matrix(n, n)
B = create_matrix(n, n)
set_from_options!(A)
set_from_options!(B)
set_up!(A)
set_up!(B)

# Get rows handled by the local processor
A_rstart, A_rend = get_range(A)
B_rstart, B_rend = get_range(B)

# Fill matrix A  with second order derivative central scheme
for i in A_rstart:A_rend
    if (i == 1)
        A[1, 1:2] = [-2., 1] / Δx^2
    elseif (i == n)
        A[n, n-1:n] = [1., -2.] / Δx^2
    else
        A[i, i-1:i+1] = [1., -2., 1.] / Δx^2
    end
end

# Fill matrix B with identity matrix
for i in B_rstart:B_rend
    B[i,i] = -1.
end

# Set boundary conditions : u(0) = 0 and u(1) = 0. Only the processor handling the corresponding rows are playing a role here.
(A_rstart == 1) && (A[1, 1:2] = [1. 0.] )
(B_rstart == 1) && (B[1,   1] = 0.      )

(A_rend == n) && (A[n, n-1:n] = [0. 1.] )
(B_rend == n) && (B[n,     n] = 0.      )

# Assemble the matrices
assemble!(A)
assemble!(B)

# Now we set up the eigenvalue solver
eps = create_eps(A, B)
set_from_options!(eps)
set_up!(eps)

# Then we solve
solve!(eps)

# And finally we can inspect the solution. Let's first get the number of converged eigenvalues:
nconv = neigs(eps)

# Then we can get/display these eigenvalues (more precisely their square root, i.e ``\simeq \omega``)
for ieig in 1:nconv
    eig = get_eig(eps, ieig)
    @show √(real(eig))
end

# You can also get all the converged eigenvalues in one call
eigs = get_eigenvalues(eps)

# We can also play with eigen vectors.
for ieig in 1:nconv
    vpr, vpi, vecpr, vecpi = get_eigenpair(eps, ieig)

    ## At this point, you can call VecGetArray to obtain a Julia array (see PetscWrap examples).
    ## If you are on one processor, you can even plot the solution to check that you have a sinus
    ## solution. On multiple processors, this would require to "gather" the solution on one processor only.
end

# Finally, let's free the memory
destroy!(A)
destroy!(B)
destroy!(eps)

# And call finalize when you're done
SlepcFinalize()

end #hide


text:
#include("printhook.jl")
#
include("benderstest.jl")
#
#include("farmer.jl")
#
#include("examples_smoketest.jl")


text:
using DBInterface, Test

@test_throws MethodError DBInterface.connect(Int64)


text:
using PyPlot
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCellDG
include("../../../test/useful_routines.jl")
include("analytical-solver.jl")
PS = PlaneStrainSolver
include("transformation-elasticity-solver.jl")

function lame_lambda(k, m)
    return k - 2m / 3
end

function angular_position(points)
    cpoints = points[1, :] + im * points[2, :]
    return rad2deg.(angle.(cpoints))
end

polyorder = 3
nelmts = 17
penaltyfactor = 1e2

folderpath = "examples/potential/circular-bc/"
filename =
    "polyorder-" * string(polyorder) * "-nelmts-" * string(nelmts) * ".png"

width = 1.0              # mm
K1, K2 = 247.0, 192.0    # GPa
mu1, mu2 = 126.0, 87.0   # GPa
# K1, K2 = 247.0, 247.0
# mu1,mu2 = 126.0, 87.0
lambda1 = lame_lambda(K1, mu1)
lambda2 = lame_lambda(K2, mu2)
stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)

rho1 = 3.93e3           # Kg/m^3
rho2 = 3.68e3           # Kg/m^3
V01 = 1.0 / rho1
V02 = 1.0 / rho2

ΔG0Jmol = -14351.0  # J/mol
molarmass = 0.147
ΔG0 = ΔG0Jmol / molarmass

theta0 = -0.067
transfstress =
    CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

meshwidth = [width, width]
numqp = required_quadrature_order(polyorder) + 2

dx = width / nelmts
penalty = penaltyfactor / dx * 0.5 * (lambda1 + mu1 + lambda2 + mu2)

basis = TensorProductBasis(2, polyorder)
interfacecenter = [0.0, 0.0]
interfaceradius = 0.5
outerradius = 1.5
analyticalsolution = PS.CylindricalSolver(
    interfaceradius,
    outerradius,
    interfacecenter,
    lambda1,
    mu1,
    lambda2,
    mu2,
    theta0,
)

mesh, cellquads, facequads, interfacequads, levelset =
    construct_mesh_and_quadratures(
        meshwidth,
        nelmts,
        basis,
        interfacecenter,
        interfaceradius,
        numqp,
    )

nodaldisplacement = nodal_displacement(
    mesh,
    basis,
    cellquads,
    facequads,
    interfacequads,
    stiffness,
    theta0,
    analyticalsolution,
    penalty,
)

cutmesh = CutCellDG.background_mesh(mesh)
refseedpoints, refseedcellids =
    CutCellDG.seed_zero_levelset(2, levelset, cutmesh)
spatialpoints = CutCellDG.map_to_spatial(
    refseedpoints[1, :, :],
    refseedcellids[1, :],
    CutCellDG.background_mesh(mesh),
)

normals = CutCellDG.collect_interface_normals(interfacequads, mesh)

productreferencepoints = refseedpoints[1, :, :]
parentreferencepoints = refseedpoints[2, :, :]
productreferencecellids = refseedcellids[1, :]
parentreferencecellids = refseedcellids[2, :]

relspatialpoints = spatialpoints .- interfacecenter
angularposition = angular_position(relspatialpoints)
sortidx = sortperm(angularposition)

angularposition = angularposition[sortidx]
parentreferencepoints = parentreferencepoints[:, sortidx]
productreferencepoints = productreferencepoints[:, sortidx]
parentreferencecellids = parentreferencecellids[sortidx]
productreferencecellids = productreferencecellids[sortidx]
spatialpoints = spatialpoints[:, sortidx]
normals = normals[:, sortidx]

parentstrain = parent_strain_at_reference_points(
    nodaldisplacement,
    basis,
    parentreferencepoints,
    parentreferencecellids,
    mesh,
)
parentstress = parent_stress(parentstrain, stiffness)

productstrain = product_strain_at_reference_points(
    nodaldisplacement,
    basis,
    theta0,
    productreferencepoints,
    productreferencecellids,
    mesh,
)
productstress = product_stress(productstrain, stiffness, theta0)

parentstrainenergy = V02 * strain_energy(parentstress, parentstrain)
productstrainenergy = V01 * strain_energy(productstress, productstrain)

exactparentstrainenergy = PS.core_strain_energy(analyticalsolution, V02)
exactproductstrainenergy =
    PS.shell_strain_energy(analyticalsolution, interfaceradius, V01)

parentradialtraction = CutCellDG.traction_force_at_points(parentstress, normals)
parentsrr = CutCellDG.traction_component(parentradialtraction, normals)

productradialtraction =
    CutCellDG.traction_force_at_points(productstress, normals)
productsrr = CutCellDG.traction_component(productradialtraction, normals)

parentdilatation = CutCellDG.dilatation(parentstrain)
productdilatation = CutCellDG.dilatation(productstrain)

parentcompwork = V02 * (1 .+ parentdilatation) .* parentsrr
productcompwork = V01 * (1 .+ productdilatation) .* productsrr

exactparentcompwork = PS.core_compression_work(analyticalsolution, V02)
exactproductcompwork =
    PS.shell_compression_work(analyticalsolution, interfaceradius, V01)

parentpotential = parentstrainenergy - parentcompwork
productpotential = productstrainenergy - productcompwork

exactparentpotential = exactparentstrainenergy - exactparentcompwork
exactproductpotential = exactproductstrainenergy - exactproductcompwork

potentialdifference = (productpotential - parentpotential) * 1e9 / abs(ΔG0)

exactpotentialdifference =
    (exactproductpotential - exactparentpotential) * 1e9 / abs(ΔG0)

pderr =
    maximum(abs.(potentialdifference .- exactpotentialdifference)) /
    abs(exactpotentialdifference)

# Δylim = 0.2 * abs(exactpotentialdifference)
# ylim = (exactpotentialdifference - Δylim, exactpotentialdifference + Δylim)
#
# fig, ax = PyPlot.subplots()
# ax.plot(angularposition, potentialdifference, color = "black")
# ax.plot(
#     angularposition,
#     exactpotentialdifference * ones(length(angularposition)),
#     linestyle = "dotted",
#     color = "black",
# )
# ax.grid()
# ax.set_ylim(ylim...)
# ax.set_xlabel("Angular Position (deg)")
# ax.set_ylabel(L"([\Phi] - [G0])/[G0]")
# fig.tight_layout()
# fig
# fig.savefig(folderpath * filename)

# fig,ax = PyPlot.subplots()
# numplotpts = length(angularposition)
# ax.plot(angularposition,parentpotential)
# ax.plot(angularposition,exactparentpotential*ones(numplotpts))
# ax.grid()
# fig
#
# fig,ax = PyPlot.subplots()
# numplotpts = length(angularposition)
# ax.plot(angularposition,productpotential)
# ax.plot(angularposition,exactproductpotential*ones(numplotpts))
# ax.grid()
# fig

# fig,ax = PyPlot.subplots()
# numplotpts = length(angularposition)
# ax.plot(angularposition,productstrainenergy)
# ax.plot(angularposition,exactproductstrainenergy*ones(numplotpts))
# ax.grid()
# fig

# fig,ax = PyPlot.subplots()
# numplotpts = length(angularposition)
# ax.plot(angularposition,productcompwork)
# ax.plot(angularposition,exactproductcompwork*ones(numplotpts))
# ax.grid()
# fig

# productpressure = CutCellDG.pressure_at_points(productstress)
# exactproductpressure = PS.shell_pressure(analyticalsolution,interfaceradius)

########################################################################
# productpressure = pressure(productstress)
# parentpressure = pressure(parentstress)
# exactparentpressure = -1/3*sum(core_stress(analyticalsolution))
# exactproductpressure = -1/3*sum(shell_stress(analyticalsolution,interfaceradius))
# fig,ax = PyPlot.subplots(2,1)
# numplotpts = length(angularposition)
# ax[1].plot(angularposition,productpressure)
# ax[1].plot(angularposition,exactproductpressure*ones(numplotpts),"--")
# ax[1].grid()
# ax[1].set_ylim(-3.0,-2.5)
# ax[2].plot(angularposition,parentpressure)
# ax[2].plot(angularposition,exactparentpressure*ones(numplotpts),"--")
# ax[2].grid()
# ax[2].set_ylim(3.25,3.75)
# fig
########################################################################


text:
################################################################################
# Copyright (C) Marcus Greiff 2020
#
# @file example_FOF_attitude.jl
# @author Marcus Greiff
# @date July 2020
#
# This example demonstrates the filtered output feedback for the attitude
# dynamics with errors related to distances on SO(3), and a closed loop
# simulation, visualizing the state trajectories along with the observer states
# and the state trajectories. SImilar to the FSF attitude examples, it also
# plots the associated Lyapunov function and distances as computed in the
# controller. Contrary to the FSF controllers, this FOF controller has memory,
# by virtue of the integration of the observer dynamcis - consequently, we can
# only call the function once pwe time-step, which means the ay multistep-method
# has to be avoided in the simulation of the closed loop system in Julia.
# Consequently, the system is simulated with a constant step Euler method
# instead of Tsit5().
#
# Furthermore, the method does not have access to the true system states, so
# The lyapunov function and errors have to be evaluated externally in Julia,
# instead of in C. So contrary to the FSF feedback controllers, the estimated
# state (updated by the FOF controller) is stored in the State structure,
# instead of the actual system state.
#
################################################################################
using DifferentialEquations
using LinearAlgebra
using LaTeXStrings
using AerialVehicleControl

showPlot = [1,2];     # Set to 0 for no plot, or any number in (1,2)
savePlot = false; # Set to true if the plot is to be saved
namePlot = "continuous_SO3"

# Recompile the C-code
recompile()

# Wrapper for the calling of the FOF controller from C
function update_control(R::ref_state_qw_t, S::dyn_state_qw_t, C::con_state_qw_fof_t, dt::Cdouble)
    status = ccall((:update_attitude_FOF_SO3_continuous, AerialVehicleControl.CONT_LIB_PATH),
        Cint,
        (Ref{ref_state_qw_t}, Ref{dyn_state_qw_t}, Ref{con_state_qw_fof_t}, Cdouble,),
        R,
        S,
        C,
        dt)
    return status;
end

# Wrapper for the updating of the measurements in the controller structure
function update_measurements(
    C::con_state_qw_fof_t,
    y0m::matrix_double_t,
    y1m::matrix_double_t,
    y2m::matrix_double_t,
    y3m::matrix_double_t
)
    status = ccall((:update_attitude_FOF_SO3_continuous_measurements, AerialVehicleControl.CONT_LIB_PATH),
        Cint,
        (Ref{con_state_qw_fof_t}, Ref{matrix_double_t}, Ref{matrix_double_t},Ref{matrix_double_t}, Ref{matrix_double_t},),
        C,
        y0m,
        y1m,
        y2m,
        y3m)
    return status;
end

# Right-hand side in the ODE used for numerical integration. Here, we have
#
# x[1:4]   = q     (quaternion, integrated in Julia)
# x[5:7]   = w     (att. rates, integrated in Julia)
# x[8:11]  = qr    (quaternion reference, integrated in Julia)
# x[12:14] = wr    (att. rates reference, integrated in Julia)
# x[15:18] = qhat  (quaternion estimate, integrated in C using CG)
# x[19:21] = what  (att. rates estimate, integrated in C using CG)
# x[22:24] = taur  (torques reference, computed in Julia)
# x[25:27] = tau   (torques actuation, computed in C)
# x[28]    = prevt (time at the previous call of the OF-controller)
# x[29]    = lyap  (Lyapunov function, computed in Julia)
# x[30]    = Psi1  (Estimate attitude error, computed in Julia)
# x[31]    = Psi2  (Tracking attitude error, computed in Julia)
#
function odefun!(dx,x,params,t)
    # Extract the states
    q     = x[1:4];
    w     = x[5:7];
    qr    = x[8:11];
    wr    = x[12:14];
    tprev = x[28];

    S = params[1];
    C = params[2];

    # Extract the parameters
    J = ntuple_2_mat(C.inertia, 3, 3)

    # Reference torques driving the reference dynamics
    taur = [sin(t), sin(2*t), sin(3*t)];

    # Reference dynamics
    dqr, dwr = attitude_dynamics(qr, wr, J, taur)

    # Create state structs for ccall
    R = ref_state_qw_t(hcat(qr), hcat(wr), hcat(dwr), 0.1)

    # Compute control signal through a ccall to the FSF configured on SO3
    if tprev != t
        timestep = Cdouble(t - tprev);
        #println("Computed new control action at a time t=$t, dt = $timestep")

        # Update the measurements in the controller object
        v1 = [C.globalDirections[1][i] for i = 1:3];
        v2 = [C.globalDirections[2][i] for i = 1:3];
        v3 = [C.globalDirections[3][i] for i = 1:3];
        y0 = hcat(w);
        y1 = hcat(quat_2_SO3(hcat(q))'*v1);
        y2 = hcat(quat_2_SO3(hcat(q))'*v2);
        y3 = hcat(quat_2_SO3(hcat(q))'*v3);
        y0m = matrix_double_t(y0);
        y1m = matrix_double_t(y1);
        y2m = matrix_double_t(y2);
        y3m = matrix_double_t(y3);
        @assert(1 == update_measurements(C, y0m, y1m, y2m, y3m));

        # Update the controller
        @assert(1 == update_control(R, S, C, timestep));
    end

    # Actuate the dynamics with the computed control signal
    tau = [C.torques[i] for i in 1:3]
    # Physical dynamics
    dq , dw  = attitude_dynamics( q,  w, J, tau )

    # Externalize state derivatives (the observer dynamics are simulated in the FOF controller)
    dx[1:4]   = dq;
    dx[5:7]   = dw;
    dx[8:11]  = dqr;
    dx[12:14] = dwr;

    # Normalize the quaternions which are simulated by DifferentialEquations.jl
    x[1:4]    = q  / norm(q );
    x[8:11]   = qr / norm(qr);

    # Output the estimates
    qhat = [S.quaternion[i] for i = 1:4];
    what = [S.attrates[i] for i = 1:3];
    x[15:18] = qhat;
    x[19:21] = what;

    # Externalize other signals for visualization
    x[22:24]  = taur;
    x[25:27]  = tau;
    x[28]     = t;   # Set the current time so that we can compute the dt

    # Evaluate distances and the Lyapunov fucntion
    Rmat   = quat_2_SO3(hcat(q))
    Rhat   = quat_2_SO3(hcat(qhat))
    Rr     = quat_2_SO3(hcat(qr))
    Re     = Rr*Rmat';
    Rtilde = Rhat*Rmat';
    we     = wr   - w;
    wtilde = what - w;

    v1 = [C.globalDirections[1][i] for i = 1:3];
    v2 = [C.globalDirections[2][i] for i = 1:3];
    v3 = [C.globalDirections[3][i] for i = 1:3];

    erre1 = Re * Rtilde' * v1 - v1;
    erre2 = Re * Rtilde' * v2 - v2;
    erre3 = Re * Rtilde' * v3 - v3;

    errt1 = Rtilde' * v1 - v1;
    errt2 = Rtilde' * v2 - v2;
    errt3 = Rtilde' * v3 - v3;

    x[29] = (+(
               +erre1'*erre1*C.gain_ki[1]
               +erre2'*erre2*C.gain_ki[2]
               +erre3'*erre3*C.gain_ki[3]
             ) / 2.0
             +we'*J*we
             +(
               +errt1'*errt1*C.gain_ki[1]
               +errt2'*errt2*C.gain_ki[2]
               +errt3'*errt3*C.gain_ki[3]
             ) / 2.0
             +(1.0/(2.0 * C.gain_cw)) * wtilde'*J*wtilde)

    x[30] = tr(Matrix{Float64}(I, 3, 3) - Rhat'*Rmat) / 2.0
    x[31] = tr(Matrix{Float64}(I, 3, 3) - Rr'*Rmat) / 2.0
end

function initialize_example()
    # Randomize the initial system states
    q0 = randn(4,1); q0 = q0./norm(q0); w0 = 2*randn(3,1);

    # Randomize an initial reference of states
    qr0 = randn(4,1); qr0 = qr0./norm(qr0); wr0 = 0*rand(3,1); ar0 = 0.0*rand(3,1); fr0 = 0.5;
    R = ref_state_qw_t(hcat(qr0), hcat(wr0), hcat(ar0), fr0)

    # Randomize a set of observer states
    qhat0 = randn(4,1); qhat0 = qhat0./norm(qhat0); what0 = 2*randn(3,1);
    S = dyn_state_qw_t(hcat(qhat0), hcat(what0))

    # Randomize the controller struct
    J  = rand_PSD(3); J  = J ./norm(J );
    Kw = rand_PSD(3); Kw = Kw./norm(Kw); Kw = 2*Kw;
    ki = (rand(1,3) .+ 0.5); ki = ki./norm(ki); ki = 2*ki;
    Cw = rand_PSD(3); Cw = Cw./norm(Cw); Cw = 2*Cw;
    cw, cR = 2.0, 2.0;
    a, b, c, d = 1.0, 2.0, 3.0, 4.0;
    Y0 = rand(3,1) .- 0.5;
    Yi = rand(3,3) .- 0.5;
    Vi = rand(3,3) .- 0.5;
    Vi = Matrix{Float64}(I, 3, 3);
    # Normalize the directions in Vi
    for i = 1:3
        Vi[:,i] = Vi[:,i] / norm(Vi[:,i])
    end

    # Form controller object
    C  = con_state_qw_fof_t(J, Kw, ki, Cw, cw, cR, a, b, c, d, Y0, Yi, Vi);

    # Initial states for the debugging and logging of signals
    output = zeros(10,1);

    x0    = vcat(q0, w0, qr0, wr0, qhat0, what0, output)
    tspan = (0.0,20.0)
    return x0, tspan, R, S, C
end

# Compute the solution
x0, tspan, _, S, C = initialize_example()
prob               = ODEProblem(odefun!, x0, tspan, [S, C]);
sol                = solve(prob, Euler(), adaptive=false, dt = 1e-4)

# Visualization of the results
if 1 in showPlot
    using Plots
    pqqr = plot(sol, vars = (0,1:4),   color=:black, style=:solid, linewidth=2, xaxis="Time (t)",yaxis="Attitude [.]", label=[L"q(t)" nothing nothing nothing])
    plot!(      sol, vars = (0,8:11),  color=:blue,  style=:dash,  linewidth=2, xaxis="Time (t)",yaxis="Attitude [.]", label=[L"q_r(t)" nothing nothing nothing])
    plot!(      sol, vars = (0,15:18), color=:red,   style=:dot,   linewidth=2, xaxis="Time (t)",yaxis="Attitude [.]", label=[L"\hat{q}(t)" nothing nothing nothing])
    pwwr = plot(sol, vars = (0,5:7),   color=:black, style=:solid, linewidth=2, xaxis="Time (t)",yaxis="Attitude rate [rad/s]", label=[L"\omega(t)" nothing nothing])
    plot!(      sol, vars = (0,12:14), color=:blue,  style=:dash,  linewidth=2, xaxis="Time (t)",yaxis="Attitude rate [rad/s]", label=[L"\omega_{r}(t)" nothing nothing])
    plot!(      sol, vars = (0,19:21), color=:red,   style=:dot,   linewidth=2, xaxis="Time (t)",yaxis="Attitude rate [rad/s]", label=[L"\hat{\omega}(t)" nothing nothing])
    pttr = plot(sol, vars = (0,22:24), color=:black, style=:solid, linewidth=2, xaxis="Time (t)",yaxis="Torque [Nm]",           label=[L"\tau_r(t)" nothing nothing])
    plot!(      sol, vars = (0,25:27), color=:blue,  style=:solid, linewidth=2, xaxis="Time (t)",yaxis="Torque [Nm]",           label=[L"\tau(t)" nothing nothing])
    plot(pqqr, pwwr, pttr, layout=(3,1), size=(1000,750))
    gui()
    if savePlot
        savefig(joinpath(AerialVehicleControl.CONT_LIB_DOCS, "$(namePlot)_FOF_states.png"))
    end
end

if 2 in showPlot
    using Plots
    plotattd = plot(sol, vars = (0,30), color=:black, linewidth=2, xaxis="Time (t)",yaxis="Attitude error",    label=L"\Psi(\hat{R}(t), R(t))")
    plot!(          sol, vars = (0,31), color=:blue,  linewidth=2, xaxis="Time (t)",yaxis="Attitude error",    label=L"\Gamma(R_r(t), R(t))")
    plotlyap    = plot(sol, vars = (0,29), color=:black, linewidth=2, xaxis="Time (t)",yaxis="Lyapunov function", label=L"\mathcal{V}(t)")
    plotlyaplog = plot(sol.t, log10.(abs.(sol[29,:])), color=:black, linewidth=2, xaxis="Time (t)",yaxis="Lyapunov function", label=L"log_{10}(\mathcal{V}(t))")
    plot(plotattd, plotlyap, plotlyaplog, layout=(3,1), size=(1000,500))
    gui()
    if savePlot
        savefig(joinpath(AerialVehicleControl.CONT_LIB_DOCS, "$(namePlot)_FOF_errors.png"))
    end
end


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
#
# Contributed by Andrei Fomiga, Stefan Karpinski, Viral B. Shah, Jeff
# Bezanson, smallnamespaces, and Adam Beckmeyer.
#
# This implementation was specifically written so that for-loops have
# integer-literal ranges as their iterable. This makes the number of
# loop iterations available as a compile-time constant, so llvm can
# unroll them. To convince llvm to do so, this script should be run
# with the evironment variable:
# JULIA_LLVM_ARGS='-unroll-threshold=500'.

using Printf, Base.Cartesian

const SOLAR_MASS = 4 * pi * pi
const DAYS_PER_YEAR = 365.24
# Precalculate the pairs of bodies that must be compared so that it
# doesn't have to be done each loop.
const PAIRS = Tuple((i, j) for i=1:4 for j=i+1:5)

# Use a struct instead of mutable struct since a struct can be stored
# inline in an array avoiding the overhead of following a pointer
struct Body
    x::NTuple{3,Float64}
    v::NTuple{3,Float64}
    m::Float64
end

function init_sun(bodies)
    p = (0.0, 0.0, 0.0)
    for b in bodies
        p = p .- b.v .* b.m
    end
    Body((0.0, 0.0, 0.0), p ./ SOLAR_MASS, SOLAR_MASS)
end

# Advance all bodies in the system by one timestep of 0.01. This
# function always uses a timestep of 0.01 and assumes that there are
# exactly 5 bodies in the system.
@inline function advance!(bodies)
    # In a system with 5 bodies, there are 10 unique pairs of
    # bodies. Δx holds the difference in position between these pairs
    # of bodies.
    Δx = @ntuple 10 k-> @inbounds bodies[PAIRS[k][1]].x .- bodies[PAIRS[k][2]].x
    dsq = @ntuple 10 k-> sum(Δx[k] .* Δx[k])

    # When called with @fastmath, 1 / sqrt(x::Float32) will use
    # SSE single-precision rsqrt approximation followed by a single
    # iteration of the Newton-Raphson method. This, followed by an
    # additional double-precision Newton-Raphson iteration gives
    # sufficient precision for this problem and is significantly
    # faster than double-precision division and sqrt.
    rd = @ntuple 10 k-> @fastmath Float64(1 / sqrt(Float32(dsq[k])))
    # This is a Newton-Raphson iteration.
    rd = @ntuple 10 k-> 1.5rd[k] - 0.5dsq[k] * rd[k] * (rd[k] * rd[k])

    # Alternatively 0.01rd[k] / dsq[k] may be faster. This is what
    # other fast implementations do, but on my machine the 2
    # multiplications are faster than a single division.
    mag = @ntuple 10 k-> 0.01rd[k] * (rd[k] * rd[k])

    # Update the velocities of each body using the precalculated mag.
    k = 1
    @inbounds for i=1:5
        bi = bodies[i]
        # For body i, since velocity is the only part of the object
        # changing on each iteration, it's more efficient to update
        # outside the vector.
        iv = bi.v
        for j=i+1:5
            iv = iv .- Δx[k] .* (mag[k] * bodies[j].m)
            bodies[j] = Body(bodies[j].x,
                             bodies[j].v .+ Δx[k] .* (mag[k] * bi.m),
                             bodies[j].m)
            k += 1
        end
        # Advance body positions using the updated velocities.
        bodies[i] = Body(bi.x .+ (iv .* 0.01), iv, bi.m)
    end
end

# Total energy of the system
function energy(bodies)
    e = 0.0
    # Kinetic energy of bodies
    @inbounds for b in bodies
        e += 0.5b.m * sum(b.v .* b.v)
    end

    # Potential energy between body i and body j
    @inbounds for (i, j) in PAIRS
        Δx = bodies[i].x .- bodies[j].x
        e -= bodies[i].m * bodies[j].m / √sum(Δx .* Δx)
    end
    e
end

# Mutate bodies array according to symplectic integrator in advance!
# for n iterations.
nbody!(bodies, n) = for i=1:n
    advance!(bodies)
end

# Doing the allocation of the Vector{Body} as a global constant
# instead of within the nbody! function speeds up inference
# considerably. Inference takes less than 60% of the time it would
# otherwise for an overall speedup of 2%-3%.
const bodies = [
    # Jupiter
    Body(( 4.84143144246472090e+0,                # x
          -1.16032004402742839e+0,                # y
          -1.03622044471123109e-1),               # z
         ( 1.66007664274403694e-3DAYS_PER_YEAR,   # vx
           7.69901118419740425e-3DAYS_PER_YEAR,   # vy
          -6.90460016972063023e-5DAYS_PER_YEAR),  # vz
           9.54791938424326609e-4SOLAR_MASS)      # mass
    # Saturn
    Body(( 8.34336671824457987e+0,
           4.12479856412430479e+0,
          -4.03523417114321381e-1),
         (-2.76742510726862411e-3DAYS_PER_YEAR,
           4.99852801234917238e-3DAYS_PER_YEAR,
           2.30417297573763929e-5DAYS_PER_YEAR),
           2.85885980666130812e-4SOLAR_MASS)
    # Uranus
    Body(( 1.28943695621391310e+1,
          -1.51111514016986312e+1,
          -2.23307578892655734e-1),
         ( 2.96460137564761618e-3DAYS_PER_YEAR,
           2.37847173959480950e-3DAYS_PER_YEAR,
          -2.96589568540237556e-5DAYS_PER_YEAR),
           4.36624404335156298e-5SOLAR_MASS)
    # Neptune
    Body(( 1.53796971148509165e+1,
          -2.59193146099879641e+1,
           1.79258772950371181e-1),
         ( 2.68067772490389322e-3DAYS_PER_YEAR,
           1.62824170038242295e-3DAYS_PER_YEAR,
          -9.51592254519715870e-5DAYS_PER_YEAR),
           5.15138902046611451e-5SOLAR_MASS)
]
push!(bodies, init_sun(bodies))

function real_main()
    if !isinteractive()
        @printf("%.9f\n", energy(bodies))
        nbody!(bodies, parse(Int, ARGS[1]))
        @printf("%.9f\n", energy(bodies))
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end


text:
function create_transformation(statistics)
    q005 = zeros(size(statistics,2))
    q995 = similar(q005)
    q50 = similar(q005)
    for i = 1:size(statistics,2)
        q005[i] = quantile(statistics[:,i],0.005)
        q50[i] = quantile(statistics[:,i],0.5)
        q995[i] = quantile(statistics[:,i],0.995)
    end
    info = [q005 q50 q995]
    writedlm("transformation_info", info)
    return info
end


text:
#Demo 12: Eye diagrams, shifting, multiple channels, etc
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
vch_vs_t = cons(:a, labels = set(yaxis="Centered @ ch#", xaxis="Time (s)"))
vch_vs_ui = cons(:a, labels = set(yaxis="Centered @ ch#", xaxis="Symbol Time (UI)"))

#==Input data
===============================================================================#
nchannels = 3
osr = 20 #samples per bit
nsamples = 127
skewmax = 100e-12 #+/- maximum skew a channel can exhibit
tbitVec = 1e-9*[1, 1/3, 1/5] #Bit periods to "simulate"
tau = 1e-9/8 #Circuit time constant
nbit_Π = 5 #Π-pulse length, in number of bits

#==Computations
===============================================================================#
#Generate a different PRBS sequenc for each channel:
seqVec = []
for seed in [11, 35, 2]
	push!(seqVec, 1.0*prbs(reglen=7, seed=seed, nsamples=nsamples))
end

#Generate parameter sweeps:
sweeplist = PSweep[
	PSweep("tbit", tbitVec)
]

#Generate data (like a simulator would do):
#-------------------------------------------------------------------------------
Π = fill(DataHR, sweeplist) do tbit
	tΠ = DataF1(0:(tbit/osr):(nbit_Π*tbit))
	return pulse(tΠ, Pole(1/tau,:rad), tpw=tbit)
end

seq = DataHR[] #Broadcastable version of prbs sequence - include undefined @ start
pat = DataHR[]
tbit = parameter(Π, "tbit")
for i in 1:nchannels
	curseq = fill(DataHR{DataF1}, sweeplist) do tbit
		undefdata = rand(0:1.0, 5) #Each test/channel gets different undefined data
		_seq = vcat(undefdata, zeros(10), seqVec[i])
		return DataF1(collect(1:length(_seq)), _seq)
	end
	push!(seq, curseq)
	Δ = rand(-skewmax:1e-12:skewmax) #Each channel is skewed differently
	push!(pat, pattern(curseq, xshift(Π, Δ), tbit=tbit))
end

#Compute derived data:
#-------------------------------------------------------------------------------
tskipundef = 10*tbit #Starting point for skipping over undefined data
#Reference channel #1 for starting point:
tstart = xcross1(pat[1]-0.5, xstart=tskipundef, allow=CrossType(:rise))
#NOTE: tstart is not really necessary. tskipundef would be sufficient.

#Show that something happened
@show tskipundef
@show tstart

#Pattern suitable for eye diagram (position start time @ t=0 & t-normalized)
eyepat = DataHR[]
for i in 1:nchannels
	curpat = xscale(xshift(pat[i], tbit/4-tstart), 1/tbit)
	push!(eyepat, curpat)
end

#Compute undefined data (normalized by tbit for readability)
#-------------------------------------------------------------------------------
undefdata = DataHR[]
for i in 1:nchannels
	curpat = xscale(clip(pat[i], xmax=20*tbit), 1/tbit)
	push!(undefdata, curpat)
end

#==Helper functions
===============================================================================#
#Scale/shift data to center by channel number:
_scaleshift(data, chnum) = (chnum-1)+(data-0.5)*.8

#TODO: supply eyeparam() as a FoldedAxis or Plot constructor???
eyeparam(tbit; teye=1.5*tbit, tstart=0) = cons(:a,
	xfolded = set(tbit, xstart=tstart, xmax=teye),
	xaxis = set(min=0, max=teye), #Force limits on exact data range.
)

#==Generate plot
===============================================================================#
p1 = cons(:plot, vch_vs_ui, title="Initial Undefined Pattern", legend=false)
for i in 1:nchannels
	push!(p1, cons(:wfrm, _scaleshift(undefdata[i], i), label="pat[$i]"))
end
p2 = cons(:plot, vch_vs_ui, eyeparam(1.0, teye=1.5), title="Eye", legend=false)
for i in 1:nchannels
	push!(p2, cons(:wfrm, _scaleshift(eyepat[i], i), label="pat[$i]"))
end
p3 = cons(:plot, vch_vs_t, title="Pattern", legend=false)
for i in 1:nchannels
	push!(p3, cons(:wfrm, _scaleshift(pat[i], i), label="pat[$i]"))
end

pcoll = push!(cons(:plotcoll, title="Normalized Eye Diagrams (shifted around channel #)"), p1, p2, p3)
	pcoll.ncolumns = 2

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
# ---
# title: 628. Maximum Product of Three Numbers
# id: problem628
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Math
# link: <https://leetcode.com/problems/maximum-product-of-three-numbers/description/>
# hidden: true
# ---
#
# Given an integer array `nums`, _find three numbers whose product is maximum
# and return the maximum product_.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3]
#     Output: 6
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,3,4]
#     Output: 24
#
#
# **Example 3:**
#
#
#
#     Input: nums = [-1,-2,-3]
#     Output: -6
#
#
#
#
# **Constraints:**
#
#   * `3 <= nums.length <= 104`
#   * `-1000 <= nums[i] <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
abstract type Edible end
eat(::Edible) = "Yum!"

mutable struct FoodBox{T<:Edible}
    food::Vector{T}
end

struct Carrot <: Edible
    variety::AbstractString
end

struct Brick
    volume::Float64
end

c = Carrot("Baby")
b = Brick(125.0)
eat(c)
eat(b)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#---------------------------------------------------------
# # Tracer Transport Operators
#---------------------------------------------------------

#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`tracer_transport_operators.ipynb`](@__NBVIEWER_ROOT_URL__/examples/generated/tracer_transport_operators.ipynb)

# To model marine biogeochemical processes on a global scale we need to be able to account for the movement of chemical constituents both horizontally and vertically.
# We do this with a **tracer transport operator**.
# When this operator acts on a tracer field it produces the advective-diffusive divergence of the tracer.

#---------------------------------------------------------
# ## Discretization
#---------------------------------------------------------

# In order to represent the transport operator on a computer one needs to discretize the tracer concentration field and the operator.
# Once discretized the tracer field is represented as a vector and the operator is represented as a sparse matrix.

#md # !!! note
#md #     A sparse matrix behaves the same way as a regular matrix.
#md #     The only difference is that in a sparse matrix the majority of the entries are zeros.
#md #     These zeros are not stored explicitly to save computer memory making it possible to deal with fairly high resolution ocean models.
#nb # > **Note**
#nb # > A sparse matrix behaves the same way as a regular matrix.
#nb # > The only difference is that in a sparse matrix the majority of the entries are zeros.
#nb # > These zeros are not stored explicitly to save computer memory making it possible to deal with fairly high resolution ocean models.

# Mathematically, the discretization converts an expression with partial derivatives into a matrix vector product.
# For the ocean circulation, we do the following conversion
#
# $$\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \cdot \nabla \right] C \longrightarrow \mathbf{T} \, \boldsymbol{C}$$
#
# where $C(\boldsymbol{r})$ is a tracer concentration at location $\boldsymbol{r}$.
# (We often omit the $\boldsymbol{r}$ dependency in equations for brevity.)
# The $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \cdot \nabla \right] C$ term is the flux divergence of the tracer due to the marine currents and turbulent eddies.
# ($\boldsymbol{u}$ is the 3D current velocity and $\mathbf{K}$ the diffusivity matrix.)
# The matrix $\mathbf{T}$ is the flux divergence transport matrix and $\boldsymbol{C}$ is the tracer concentration vector.
#
# One can go a long way towards understanding what a tracer transport operator is by playing with a simple model with only a few boxes, which is the goal of this example.

# The simple box model we consider is embeded in a 2×2×2 "shoebox".
# It has 5 *wet* boxes and 3 *dry* boxes, as illustrated below:

#md # ```@raw html
#md # <img src="https://user-images.githubusercontent.com/4486578/58314610-3b130b80-7e53-11e9-9fe8-9527cdcca2d0.png" width =800>
#md # ```
#nb # <img src="https://user-images.githubusercontent.com/4486578/58314610-3b130b80-7e53-11e9-9fe8-9527cdcca2d0.png" width =800>

# The circulation consists of
# - a meridional overturning circulation flowing in a cycle through boxes 1 → 2 → 6 → 5 → 1 (shown in the "meridional section 1" panel above)
# - a zonal current in a reentrant cycling through boxes 1 → 3 → 1 (shown in the "layer 1" panel above)
# - vertical mixing representing deep convection between boxes 2 ↔ 6 (not shown)

#---------------------------------------------------------
# ## The model grid and the transport matrix
#---------------------------------------------------------

# Like for any models using AIBECS, we start by telling Julia just that:

using AIBECS

# We then load the shoebox model via

wet3D, grd, T = Primeau_2x2x2.load() ;

# where we have loaded 3 objects, `wet3D`, `grd`, and `T`.

# `wet3D` is 3D array representing the 3D ocean with `true` for "wet" boxes and `false` for "dry" boxes.
# Let's have a look at its contents:

wet3D

# It's a 2×2×2 `BitArray`, i.e., an array of bit elements (the `true` and `false` entries).
# You can check that it matches our "shoebox" model.
# (Well, except for the orientation, for which northwards in the box model is downwards in the array.)

# We can find all the wet boxes simply via

findall(wet3D)

# These are the 3D indices of the wet boxes, `(i,j,k)`, called the "cartesian" indices.
# If you want the "linear" indices, i.e., the numbers as shown in the image of the shoebox model, you can simply transform `wet3D` into a vector, via

iwet = findall(vec(wet3D))

# We can also check that we indeed have the expected number of wet boxes via

nb = length(iwet)

# Now let's look at the grid, `grd`:

grd

# It's an `OceanGrid` object, of size 2×2×2, as expected.
# This object is defined in the [OceanGrids](https://github.com/briochemc/OceanGrids.jl) package, on which AIBECS depends.
# There are many ways to look inside the grid, one of which is to iterate over it:

[println(box) for box in grd] ;

# shows some details about all the boxes of the model, one at a time.
# The `grd` object also contains other information about the grid, like the 3D depths of the boxes:

grd.depth_3D

# or the 3D latitudes:

grd.lat_3D

# where you can check that northwards = downwards in the array.
# You may notice these come with units!
# This helps ensure that degrees of latitude are not confused with meters for example.

#md # !!! note
#md #     Julia comes with [Unitful](https://github.com/PainterQubits/Unitful.jl), a package for using units, which AIBECS uses.
#nb # > **Note**
#nb # > Julia comes with [Unitful](https://github.com/PainterQubits/Unitful.jl), a package for using units, which AIBECS uses.

# Finally, let's have a look at the transport matrix `T`, which represents $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \cdot \nabla \right]$, i.e., the flux divergence operator for dissolved tracers:

T

# It's a sparse matrix (a `SparseMatrixCSC` object in Julia), which only stores those entries that are not zero.
# We can display the full matrix via

Matrix(T)

# to check its structure out.
# Note that all the diagonal terms are positive, which is the sign that the transport matrix acts as a divergence, which would be positive for a box containing all the tracer.
# For example, if there was only some tracer in box 2, then that tracer should "diverge" away from that box, resulting in the positiva value in that box:

j = 2                            # index where we put some tracer
x = 1.0 * [i == j for i in 1:nb] # vector of 0's except for index j
x, (T * x)[j]

# shows that the flux divergence `T * x` is positive where we injected the tracer.

#---------------------------------------------------------
# ## Radiocarbon
#---------------------------------------------------------

# Radiocarbon, ¹⁴C, is produced by cosmic rays in the lower stratosphere and upper troposphere.
# It quickly reacts with oxygen to produce ¹⁴CO₂, which is then mixed throughout the troposphere and enters the ocean through air–sea gas exchange.
# Because the halflife of radiocarbon is only 5730 years a significant amount of decay can occur before the dissolved inorganic radiocarbon (DI¹⁴C) can mix uniformally throughout the ocean.
# As such the ¹⁴C serves as a tracer label for water that was recently in contact with the atmosphere.

# ### Tracer Equation

# Mathematically, the ¹⁴C tracer concentration, denoted $R$ (for Radiocarbon), satisfies the following tracer equation:
#
# $$\frac{\partial R}{\partial t} + \nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \cdot \nabla \right] R = \Lambda(R_\mathsf{atm} - R) - R / \tau,$$
#
# where $\Lambda(R_\mathsf{atm} - R)$ represents the air–sea exchanges and $R / \tau$ the radioactive decay rate.
# ($\tau$ is the radioactive decay timescale.)

# The discretized tracer is thus given by
#
# $$\frac{\partial \boldsymbol{R}}{\partial t} + \mathbf{T} \, \boldsymbol{R} = \mathbf{\Lambda}(R_\mathsf{atm} - \boldsymbol{R}) - \boldsymbol{R} / \tau.$$

# ### Translation to AIBECS Code

# We will perform an idealized radiocarbon simulation in our model and use the ocean circulation defined earlier using AIBECS.
# In this model we prescribe the atmospheric concentration, $R_\mathsf{atm}$, to be simply equal to 1.
# (We do not specify its unit or its specific value because it is not important for determining the age of a water parcel — only the decay rate does.)

# To use AIBECS, one must put the equations into the generic form of
#
# $$\frac{\partial \boldsymbol{x}}{\partial t} + \mathbf{T}(\boldsymbol{p}) \, \boldsymbol{x} = \boldsymbol{G}(\boldsymbol{x}, \boldsymbol{p}),$$
#
# where $\boldsymbol{x}$ is the state vector, $\boldsymbol{p}$ is the vector of model parameters, $\mathbf{T}(\boldsymbol{p})$ is the transport operator, and $\boldsymbol{G}(\boldsymbol{x}, \boldsymbol{p})$ is the local sources minus sinks.

# In our radiocarbon-model context, with $\boldsymbol{x} = \boldsymbol{R}$, we have that
#
# $$\boldsymbol{G}(\boldsymbol{x}, \boldsymbol{p}) = \mathbf{\Lambda}(R_\mathsf{atm} - \boldsymbol{x}) - \boldsymbol{x} / \tau.$$
#
# Hence, we must create `T(p)` and `G(x,p)` to give AIBECS the means to simulate the tracer distribution and/or its evolution in time.

# #### Sources and Sinks

# The local sources and sinks are thus simply given by

function G(x, p)
    τ, Ratm = p.τ, p.Ratm
    return Λ(Ratm .- x, p) - x / τ
end

# where `τ` is the decay rate timescale and `Ratm` is the atmospheric concentration of radiocarbon.

# We must define the air–sea exchange rate, `Λ(x,p)`, which requires us to define which boxes are located at the surface first.
# This is done, e.g., via

surface_boxes = grd.depth_3D[iwet] .== grd.depth[1]

# The air–sea exchange rate is then given by

function Λ(x, p)
    λ, h = p.λ, p.h
    return λ / h * surface_boxes .* x
end

# where `λ` is the piston velocity and `h` is the height of the top layer of the model grid.

# #### Parameters

# For the air–sea gas exchange, we use a constant piston velocity $\lambda$ of 50m / 10years, which will happen in the top layer, of height given by, well, the height of the top layer.
# And for the radioactive decay we use a timescale $\tau$ of 5730/log(2) years.
# We define these as parameters using the dedicated API from the AIBECS:

t = empty_parameter_table()                   # initialize an empty table of parameters
add_parameter!(t, :τ, 5730u"yr"/log(2)) # radioactive decay e-folding timescale
add_parameter!(t, :λ, 50u"m" / 10u"yr") # piston velocity
add_parameter!(t, :h, grd.δdepth[1])    # height of top layer
add_parameter!(t, :Ratm, 1.0u"mol/m^3") # atmospheric concentration
t

# shows the parameters that you just created.

# We now generate a new object to contain all these parameters via

initialize_Parameters_type(t, "C14_shoebox_parameters") # creates the type for parameters
p = C14_shoebox_parameters()                            # creates the parameters object

# #### Generate the state function and its Jacobian

# The last step for the setup is for AIBECS to create $\boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p}) = \boldsymbol{G}(\boldsymbol{x}, \boldsymbol{p}) - \mathbf{T}(\boldsymbol{p}) \, \boldsymbol{x}$, which defines the rate of change of the state, $\boldsymbol{x}$.
# This is done via

F, ∇ₓF = state_function_and_Jacobian(p -> T, G) # generates the state function (and its Jacobian!)
x = zeros(nb)
F(x,p)

# Note here that AIBECS has automatically created `∇ₓF`, i.e., $\nabla_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p})$, which is the Jacobian of the system.
# This Jacobian will be useful in the simulations below.

# That's it!
# Your model is entirely setup and ready to be used for simulations!

#---------------------------------------------------------
# ## Time stepping
#---------------------------------------------------------

# One way to see how the tracer evolves with time is to time step it.
# Here we will use the [Crank-Nicolson](https://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method) scheme embedded in AIBECS.
#
# A single step is performed via, e.g.,

δt = ustrip(1.0u"yr" |> u"s")
AIBECS.crank_nicolson_step(x, p, δt, F, ∇ₓF)

# We can write a function to run all the time steps and save them into a `x_hist` object, via

function time_steps(x₀, Δt, n, F, ∇ₓF)
    x_hist = [x₀]
    δt = Δt / n
    for i in 1:n
        push!(x_hist, AIBECS.crank_nicolson_step(last(x_hist), p, δt, F, ∇ₓF))
    end
    return reduce(hcat, x_hist), 0:δt:Δt
end

#md # !!! note
#md #     We store all the history of `x` in `x_hist`.
#md #     At each step, the `push!` function adds a new `x` to `x_hist`.
#md #     Technically, `x_hist` is a vector of vectors, so at the end, we horizontally concatenate it via `reduce(hcat, x_hist)` to rearrange it into a 2D array.
#nb # > **Note**
#nb # > We store all the history of `x` in `x_hist`.
#nb # > At each step, the `push!` function adds a new `x` to `x_hist`.
#nb # > Technically, `x_hist` is a vector of vectors, so at the end, we horizontally concatenate it via `reduce(hcat, x_hist)` to rearrange it into a 2D array.

# Now let's simulate the evolution of radiocarbon for 7500 years, starting from a concentration of 1 everywhere, via

Δt = ustrip(7500u"yr" |> u"s") # 7500 years in seconds
x₀ = ones(5)             #
x_hist, t_hist = time_steps(x₀, Δt, 1000, F, ∇ₓF) # runs the simulation

# This should take a few seconds to run.

#md # !!! note
#md #     AIBECS also has other time-stepping methods available, including Euler-froward (fast but unstable), Euler-backwards (like Crank-Nicolson, slow but stable), and Crank-Nicolson-Leapfrog (Fast-ish and more stable than Euler forward).
#nb # > **Note**
#nb # > AIBECS also has other time-stepping methods available, including Euler-froward (fast but unstable), Euler-backwards (like Crank-Nicolson, slow but stable), and Crank-Nicolson-Leapfrog (Fast-ish and more stable than Euler forward).

# Once it's done, we can plot the evolution of radiocarbon through time via

using PyPlot
clf()
C14age_hist = -log.(x_hist) * ustrip(p.τ * u"s" |> u"yr")
plot(t_hist * ustrip(1u"s" |> u"yr"), C14age_hist')
xlabel("simulation time (years)")
ylabel("¹⁴C age (years)")
legend("box " .* string.(iwet))
title("Simulation of the evolution of ¹⁴C age with Crank-Nicolson time steps")
gcf()

# The box model took more than 4000 years to spin up to equilibrium.
# For a box model that's no big deal because it is not computationally expensive to run the model, but for a big circulation model waiting for the model to spinup is painfully long.
# We therefore want a better way to find the equilibrium solution.

#---------------------------------------------------------
# ## Solving Directly for the Steady State
#---------------------------------------------------------

# With the AIBECS, you can create the steady state problem and solve it in just 2 commands:

prob = SteadyStateProblem(F, ∇ₓF, x, p)
x_final = solve(prob, CTKAlg()).u

# Converting radiocarbon into years gives the following values

C14age_final = -log.(x_final) * p.τ * u"s" .|> u"yr"
println.("box ", iwet, ": ", C14age_final);

# These are exactly the limit that the Radiocarbon age reaches after about 4000 years of simulation!
# So what happened there?

# Well, AIBECS used a version of [Newton's method](https://en.wikipedia.org/wiki/Newton%27s_method) to solve for $\boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p}) = 0$.
# Simply put, Newton's method iterates on the state via the recursion relation
#
# $$\boldsymbol{x}_{n+1} = \boldsymbol{x}_n - \nabla_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p})^{-1} \boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p}).$$

# Here, our radiocarbon model only requires a single iteration.
# This is because the state function, $\boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p})$, is linear.


text:
include("DoME.jl");
include("usefulFunctions.jl");

# Uncomment the desired line to perform the experiments and obtain the results shown in Table 5
datasetName = "1027_ESL";                      Strategy = Strategy2; MinimumReductionMSE = 1e-5; MaximumNodes = 100;
# datasetName = "1028_SWD";                      Strategy = Strategy1; MinimumReductionMSE = 1e-5; MaximumNodes = 60;
# datasetName = "1029_LEV";                      Strategy = Strategy3; MinimumReductionMSE = 1e-4; MaximumNodes = 45;
# datasetName = "1030_ERA";                      Strategy = Strategy2; MinimumReductionMSE = 1e-7; MaximumNodes = 100;
# datasetName = "1089_USCrime";                  Strategy = Strategy3; MinimumReductionMSE = 1e-3; MaximumNodes = 20;
# datasetName = "1096_FacultySalaries";          Strategy = Strategy1; MinimumReductionMSE = 1e-7; MaximumNodes = 155;
# datasetName = "192_vineyard";                  Strategy = Strategy4; MinimumReductionMSE = 1e-3; MaximumNodes = 15;
# datasetName = "195_auto_price";                Strategy = Strategy1; MinimumReductionMSE = 1e-7; MaximumNodes = 65;
# datasetName = "207_autoPrice";                 Strategy = Strategy2; MinimumReductionMSE = 1e-4; MaximumNodes = 160;
# datasetName = "210_cloud";                     Strategy = Strategy4; MinimumReductionMSE = 1e-7; MaximumNodes = 55;
# datasetName = "228_elusage";                   Strategy = Strategy1; MinimumReductionMSE = 1e-3; MaximumNodes = 15;
# datasetName = "230_machine_cpu";               Strategy = Strategy4; MinimumReductionMSE = 1e-7; MaximumNodes = 70;
# datasetName = "485_analcatdata_vehicle";       Strategy = Strategy4; MinimumReductionMSE = 1e-5; MaximumNodes = 145;
# datasetName = "519_vinnie";                    Strategy = Strategy3; MinimumReductionMSE = 1e-4; MaximumNodes = 5;
# datasetName = "522_pm10";                      Strategy = Strategy2; MinimumReductionMSE = 1e-4; MaximumNodes = 30;
# datasetName = "523_analcatdata_neavote";       Strategy = Strategy3; MinimumReductionMSE = 1e-2; MaximumNodes = 10;
# datasetName = "527_analcatdata_election2000";  Strategy = Strategy3; MinimumReductionMSE = 1e-5; MaximumNodes = 100;
# datasetName = "542_pollution";                 Strategy = Strategy4; MinimumReductionMSE = 1e-4; MaximumNodes = 30;
# datasetName = "547_no2";                       Strategy = Strategy1; MinimumReductionMSE = 1e-4; MaximumNodes = 50;
# datasetName = "556_analcatdata_apnea2";        Strategy = Strategy4; MinimumReductionMSE = 1e-4; MaximumNodes = 115;
# datasetName = "557_analcatdata_apnea1";        Strategy = Strategy1; MinimumReductionMSE = 1e-6; MaximumNodes = 125;
# datasetName = "561_cpu";                       Strategy = Strategy3; MinimumReductionMSE = 1e-7; MaximumNodes = 140;
# datasetName = "659_sleuth_ex1714";             Strategy = Strategy3; MinimumReductionMSE = 1e-2; MaximumNodes = 40;
# datasetName = "663_rabe_266";                  Strategy = Strategy1; MinimumReductionMSE = 1e-7; MaximumNodes = 185;
# datasetName = "665_sleuth_case2002";           Strategy = Strategy3; MinimumReductionMSE = 1e-6; MaximumNodes = 30;
# datasetName = "666_rmftsa_ladata";             Strategy = Strategy1; MinimumReductionMSE = 1e-3; MaximumNodes = 70;
# datasetName = "678_visualizing_environmental"; Strategy = Strategy3; MinimumReductionMSE = 1e-5; MaximumNodes = 30;
# datasetName = "687_sleuth_ex1605";             Strategy = Strategy2; MinimumReductionMSE = 1e-7; MaximumNodes = 30;
# datasetName = "690_visualizing_galaxy";        Strategy = Strategy1; MinimumReductionMSE = 1e-7; MaximumNodes = 75;
# datasetName = "695_chatfield_4";               Strategy = Strategy3; MinimumReductionMSE = 1e-6; MaximumNodes = 30;
# datasetName = "706_sleuth_case1202";           Strategy = Strategy3; MinimumReductionMSE = 1e-3; MaximumNodes = 70;
# datasetName = "712_chscase_geyser1";           Strategy = Strategy1; MinimumReductionMSE = 1e-3; MaximumNodes = 45;

# Load the dataset
(inputs, targets) = loadDataset(datasetName);

numFolds = 10;
# Create the same cross-validation indices used in the rest of experiments
indicesKFold = crossvalidationIndices(length(targets), numFolds);

testValues = Array{Float64,1}(undef,numFolds);

println("Dataset \"", datasetNameFromFile(datasetName), "\" with strategy ", Strategy, ", Min. MSE reduction ", MinimumReductionMSE, " and Maximum num. nodes ", MaximumNodes);

for numFold = 1:numFolds

    trainingInputs = inputs[indicesKFold.!=numFold,:];
    trainingTargets = targets[indicesKFold.!=numFold];

    (_, _, _, bestTree) = dome(trainingInputs, trainingTargets;
        minimumReductionMSE = MinimumReductionMSE ,
        maximumNodes = MaximumNodes ,
        strategy = Strategy ,
        showText = false
    );

    # Get this expression as a string with vector operations
    expr = vectorString(bestTree);
    # Convert the variable names "X" to "inputs"
    expr = replace(expr, "X" => "inputs");
    # Evaluate this expression
    outputs = eval(Meta.parse(expr))
    # Calculate MSE in the test samples
    testMSE = mean(((targets .- outputs).^2)[indicesKFold.==numFold]);

    println("   Finished fold $numFold/$numFolds, MSE in test: $testMSE");

    testValues[numFold] = testMSE;

end;

println("Median test MSE: ", median(testValues));


text:
set_cursor_col(n::Int = 0) = "\e[$(n)G"
clear_to_eol() = "\e[0J"
set_cursor_line(n::Int = 1) = "\e[$(n)F"


text:
using GLMakie.AbstractPlotting: backend_display, getscreen

@testset "unit tests" begin
    @testset "Window handling" begin
        AbstractPlotting.inline!(false)
        screen = GLMakie.global_gl_screen((100, 100), false)
        @test isopen(screen)
        scene = scatter(1:4);
        screen2 = display(scene)
        @test screen === screen2
        @test getscreen(scene) === screen
        close(screen)

        # assure we correctly close screen and remove it from plot
        @test getscreen(scene) === nothing
        @test !events(scene).window_open[]
        @test isempty(events(scene).window_open.listeners)
    end

    @testset "Pick a plot element or plot elements inside a rectangle" begin
        N = 100000
        scene = scatter(1:N, 1:N)
        xlims!((99990,100000))
        ylims!((99990,100000))
        screen = display(scene)
        yield()
        # test for pick a single data point (with idx > 65535)
        point_px = AbstractPlotting.project(scene, Point(100000,100000))
        plot,idx = pick(scene, point_px)
        @test idx == 100000

        # test for pick a rectangle of data points (also with some indices > 65535)
        rect = FRect2D(99990.5,99990.5,8,8)
        origin_px = AbstractPlotting.project(scene, Point(origin(rect)))
        tip_px = AbstractPlotting.project(scene, Point(origin(rect) .+ widths(rect)))
        rect_px = IRect2D(round.(origin_px), round.(tip_px .- origin_px))
        #! there is no pick(::Scene,::IRect2D)
        plot_idx = pick(screen, rect_px)

        # objects returned in plot_idx should be either grid lines (i.e. LineSegments) or Scatter points
        @test all(pi-> pi[1] isa Union{LineSegments,Scatter}, plot_idx)
        # scatter points should have indices equal to those in 99991:99998
        scatter_plot_idx = filter(pi -> pi[1] isa Scatter, plot_idx)
        @test Set(last.(scatter_plot_idx)) == Set(99991:99998)
    end
end


text:
"""
    mutable struct RTEMatpowerGenerator <: AbstractNodeLabel

Mutable structure descendant of AbstractNodeLabel

# Fields
- `power_min::Complex128`
- `power_max::Complex128`
- `cost_degree::Int64`
- `cost_coeffs::Array{Float64}`
- `gen_status::Bool`: generator on/off

"""
mutable struct RTEMatpowerGenerator <: AbstractNodeLabel
  power_min::Complex #Smin = Pmin + i Qmin
  power_max::Complex #Smax = Pmax + i Qmax
  active_power_input::Float64
  gen_status::Bool
end

"""
    create_vars!(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String) where T <: RTEMatpowerGenerator

Create power variables n for generator `elemid` at `bus` in `bus_vars` if `elem_formulation == :NewVar`\n
Return nothing

"""
function create_vars!(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String) where T <: RTEMatpowerGenerator
    bus_vars[elemid] = Variable(variable_name("Sgen", bus, elemid, scenario), Complex)
    #bus_vars["lambda_$(sign)"] = Variable("lambda_$(sign)", Real)
  return
end

"""
[cstrname, polynom, lb, ub] = Snodal(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}) where T <: RTEMatpowerGenerator

Return the polynomial contribution in power of generator `elemid` at `bus`(name, value, lower bound, upper bound). Will be used to construct power balance constraints in polynomial problem.\n
If `elem_formulation == :NbMinVar`, return generator bounds `["UNIT", Polynomial() ,Smin, Smax]`\n
If `elem_formulation == :NewVar`, return -Sgen `["UNIT", -Sgen, 0, 0]`\n
Return no contribution `["", Polynomial(), 0, 0]` otherwise
"""
function Snodal(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}) where T <: RTEMatpowerGenerator
  cstrname = "UNIT"
  return [cstrname, -bus_vars[elemid], 0, 0]

end

"""
    constraint(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String, OPFpbs::OPFProblems) where T <: RTEMatpowerGenerator

Return all the constraints defined by generator `elemid` at `bus`. Will be used to construct constraints in polynomial problem.\n
If `elem_formulation == :NewVar`, return generator bounds : "Genbounds" => Smin <= Sgen <= Smax\n
Return empty dictionary otherwise
"""
function constraint(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String, OPFpbs::OPFProblems) where T <: RTEMatpowerGenerator
  #bounds on the reactive power
    Qmin = imag(element.power_min)
    Qmax = imag(element.power_max)
    reactive_power = 0.5*(-im)*(bus_vars[elemid] - conj(bus_vars[elemid]))
    Pmin = real(element.power_min)
    Pmax = real(element.power_max)
    Pinput = element.active_power_input
    if !element.gen_status
      Pmin = Pmax = 0
      Qmin = Qmax = 0
    end
    active_power = 0.5*(bus_vars[elemid]+conj(bus_vars[elemid]))
    return SortedDict{String, Constraint}( "Reactivepower" => Qmin << reactive_power << Qmax)

end

"""
    cost(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_elems_vars::SortedDict{String, Variable}, Snode::Polynomial, lb, ub) where T <: RTEMatpowerGenerator

Return the polynomial contribution in objective generator `elemid` at `bus`.
"""
function cost(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_elems_var::SortedDict{String, Variable}, Snode::Polynomial, lb, ub) where T <: RTEMatpowerGenerator
  Sgen = bus_elems_var[elemid]
  Pgen = 0.5*(Sgen+conj(Sgen))
  if !element.gen_status
    Pgen = 0
  end
  return Pgen
end


text:
using DifferentialEquations
using Plots
gr()

function lif(u,p,t);
    gL, EL, C, Vth, I = p
    (-gL*(u-EL)+I)/C
end

function thr(u,t,integrator)
    integrator.u > integrator.p[4]
end

function reset!(integrator)
    integrator.u = integrator.p[2]
end

threshold = DiscreteCallback(thr,reset!)
current_step= PresetTimeCallback([2,15],integrator -> integrator.p[5] += 210.0)
cb = CallbackSet(current_step,threshold)

u0 = -75
tspan = (0.0, 40.0)
# p = (gL, EL, C, Vth, I)
p = [10.0, -75.0, 5.0, -55.0, 0]

prob = ODEProblem(lif, u0, tspan, p, callback=cb)

sol = solve(prob)

plot(sol)

#Izhikevichch Model
using DifferentialEquations
using Plots

function izh!(du,u,p,t);
    a, b, c, d, I = p

    du[1] = 0.04*u[1]^2+5*u[1]+140-u[2]+I
    du[2] = a*(b*u[1]-u[2])
end

function thr(u,t,integrator)
    integrator.u[1] >= 30
end

function reset!(integrator)
    integrator.u[1] = integrator.p[3]
    integrator.u[2] += integrator.p[4]
end

threshold = DiscreteCallback(thr,reset!)
current_step= PresetTimeCallback(50,integrator -> integrator.p[5] += 10)
cb = CallbackSet(current_step,threshold)

p = [0.02, 0.2, -50, 2, 0]
u0 = [-65, p[2]*-65]
tspan = (0.0, 300)

prob = ODEProblem(izh!, u0, tspan, p, callback=cb)

sol = solve(prob);
plot(sol, vars=1)

plot(sol, vars=2)

p = [0.02, 0.2, -65, 8, 0]
u0 = [-65, p[2]*-65]
tspan = (0.0, 300)

prob = ODEProblem(izh!, u0, tspan, p, callback=cb)
sol = solve(prob);
plot(sol, vars=1)

using DifferentialEquations
using Plots

# Potassium ion-channel rate functions
alpha_n(v) = (0.02 * (v - 25.0)) / (1.0 - exp((-1.0 * (v - 25.0)) / 9.0))
beta_n(v) = (-0.002 * (v - 25.0)) / (1.0 - exp((v - 25.0) / 9.0))

# Sodium ion-channel rate functions
alpha_m(v) = (0.182*(v + 35.0)) / (1.0 - exp((-1.0 * (v + 35.0)) / 9.0))
beta_m(v) = (-0.124 * (v + 35.0)) / (1.0 - exp((v + 35.0) / 9.0))

alpha_h(v) = 0.25 * exp((-1.0 * (v + 90.0)) / 12.0)
beta_h(v) = (0.25 * exp((v + 62.0) / 6.0)) / exp((v + 90.0) / 12.0)

function HH!(du,u,p,t);
    gK, gNa, gL, EK, ENa, EL, C, I = p
    v, n, m, h = u

    du[1] = (-(gK * (n^4.0) * (v - EK)) - (gNa * (m ^ 3.0) * h * (v - ENa)) - (gL * (v - EL)) + I) / C
    du[2] = (alpha_n(v) * (1.0 - n)) - (beta_n(v) * n)
    du[3] = (alpha_m(v) * (1.0 - m)) - (beta_m(v) * m)
    du[4] = (alpha_h(v) * (1.0 - h)) - (beta_h(v) * h)
end

current_step= PresetTimeCallback(100,integrator -> integrator.p[8] += 1)

# n, m & h steady-states
n_inf(v) = alpha_n(v) / (alpha_n(v) + beta_n(v))
m_inf(v) = alpha_m(v) / (alpha_m(v) + beta_m(v))
h_inf(v) = alpha_h(v) / (alpha_h(v) + beta_h(v))

p = [35.0, 40.0, 0.3, -77.0, 55.0, -65.0, 1, 0]
u0 = [-60, n_inf(-60), m_inf(-60), h_inf(-60)]
tspan = (0.0, 1000)

prob = ODEProblem(HH!, u0, tspan, p, callback=current_step)

sol = solve(prob);
plot(sol, vars=1)

plot(sol, vars=[2,3,4], tspan=(105.0,130.0))

function gSyn(max_gsyn, tau, tf, t);
    if t-tf >= 0
        return max_gsyn * exp(-(t-tf)/tau)
    else
        return 0.0
    end
end
function HH!(du,u,p,t);
    gK, gNa, gL, EK, ENa, EL, C, I, max_gSyn, ESyn, tau, tf = p
    v, n, m, h = u

    ISyn = gSyn(max_gSyn, tau, tf, t) * (v - ESyn)

    du[1] = (-(gK * (n^4.0) * (v - EK)) - (gNa * (m ^ 3.0) * h * (v - ENa)) - (gL * (v - EL)) + I - ISyn) / C
    du[2] = (alpha_n(v) * (1.0 - n)) - (beta_n(v) * n)
    du[3] = (alpha_m(v) * (1.0 - m)) - (beta_m(v) * m)
    du[4] = (alpha_h(v) * (1.0 - h)) - (beta_h(v) * h)
end

p = [35.0, 40.0, 0.3, -77.0, 55.0, -65.0, 1, 0, 0.008, 0, 20, 100]
tspan = (0.0, 200)
prob = ODEProblem(HH!, u0, tspan, p)
sol = solve(prob);
plot(sol, vars=1)

p = [35.0, 40.0, 0.3, -77.0, 55.0, -65.0, 1, 0, 0.01, 0, 20, 100]
tspan = (0.0, 200)
prob = ODEProblem(HH!, u0, tspan, p)
sol = solve(prob);
plot!(sol, vars=1)

function HH!(du,u,p,t);
    gK, gNa, gL, EK, ENa, EL, C, I, tau, tau_u, tau_R, u0, gmax, Esyn  = p
    v, n, m, h, u, R, gsyn = u

    du[1] = ((gK * (n^4.0) * (EK - v)) + (gNa * (m ^ 3.0) * h * (ENa - v)) + (gL * (EL - v)) + I + gsyn * (Esyn - v)) / C
    du[2] = (alpha_n(v) * (1.0 - n)) - (beta_n(v) * n)
    du[3] = (alpha_m(v) * (1.0 - m)) - (beta_m(v) * m)
    du[4] = (alpha_h(v) * (1.0 - h)) - (beta_h(v) * h)

    # Synaptic variables
    du[5] = -(u/tau_u)
    du[6] = (1-R)/tau_R
    du[7] = -(gsyn/tau)
end

function epsp!(integrator);
    integrator.u[5] += integrator.p[12] * (1 - integrator.u[5])
    integrator.u[7] += integrator.p[13] * integrator.u[5] * integrator.u[6]
    integrator.u[6] -= integrator.u[5] * integrator.u[6]

end

epsp_ts= PresetTimeCallback(100:100:500, epsp!)

p = [35.0, 40.0, 0.3, -77.0, 55.0, -65.0, 1, 0, 30, 1000, 50, 0.5, 0.005, 0]
u0 = [-60, n_inf(-60), m_inf(-60), h_inf(-60), 0.0, 1.0, 0.0]
tspan = (0.0, 700)
prob = ODEProblem(HH!, u0, tspan, p, callback=epsp_ts)
sol = solve(prob);
plot(sol, vars=1)

plot(sol, vars=7)

plot(sol, vars=[5,6])

epsp_ts= PresetTimeCallback(100:1000:5100, epsp!)

p = [35.0, 40.0, 0.3, -77.0, 55.0, -65.0, 1, 0, 30, 500, 50, 0.5, 0.005, 0]
u0 = [-60, n_inf(-60), m_inf(-60), h_inf(-60), 0.0, 1.0, 0.0]
tspan = (0.0, 5300)
prob = ODEProblem(HH!, u0, tspan, p, callback=epsp_ts)
sol = solve(prob);
plot(sol, vars=7)

plot(sol, vars=[5,6])

epsp_ts= PresetTimeCallback(100:100:500, epsp!)

p = [35.0, 40.0, 0.3, -77.0, 55.0, -65.0, 1, 0, 30, 100, 1000, 0.5, 0.005, 0]
u0 = [-60, n_inf(-60), m_inf(-60), h_inf(-60), 0.0, 1.0, 0.0]
tspan = (0.0, 700)
prob = ODEProblem(HH!, u0, tspan, p, callback=epsp_ts)
sol = solve(prob);
plot(sol, vars=7)

plot(sol, vars=[5,6])

using SciMLTutorials
SciMLTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
"""
    monte_carlo_integration(f::Function, a::Real, b::Real, n::Int)

Monte carlo integration is a very easy and scalable way to do multidimentional integrals.
However, only single variable integrals are considered.

# Arguments
- `f`: the function to integrate. (at the momment only single variable is suported)
- `a`: start in the integration limits.
- `b`: endin the integration limits.
- `N`: Number of points to sample. For most simple functions, 1000 to 10,000 should be okay.

# Examples
```julia
julia> monte_carlo_integration(x -> 3*x^2, 0, 1, 100000) # integrate a polynomial
1.0000037602209

julia> monte_carlo_integration(x -> sin(x), 0, pi, 1000) # integrate the sin function
2.0018927826323756
```

# References
- https://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/monte-carlo-methods-in-practice/monte-carlo-integration
- https://kingaa.github.io/sbied/pfilter/monteCarlo.html

# Contributors
- [AugustoCL](https://github.com/AugustoCL)
- [Ved Mahajan](https://github.com/Ved-Mahajan)
"""
function monte_carlo_integration(f::Function, a::Real, b::Real, n::Int)
    Δₓ = ((b - a) / n)

    Σ = 0.0
    for _ ∈ 1:n
        # generate uniform(a, b) using uniform(0, 1)
        Xᵢ = a + (b-a)*rand()
        Σ += f(Xᵢ)
    end

    return Δₓ * Σ
end

max_stars_repo_name: frankschmitt/Julia

max_issues_repo_name: frankschmitt/Julia

max_forks_repo_name: frankschmitt/Julia


text:
using AWS
using AWS.S3

include("config.jl")

## env=AWSEnv(; id=id, key=key, dbg=dbg)

env = AWSEnv(timeout=60.0)

bkt = "1-test-temp-bkt"

println("List all buckets")
resp=S3.list_all_buckets(env)
println(resp.obj)

println("Create a bucket")
acl=S3.S3_ACL()
acl.acl="public-read"
resp=S3.create_bkt(env, bkt, acl=acl)
println(resp)

po = S3.PutObjectOptions()
po.acl = acl

println("Add a file")
resp = S3.put_object(env, bkt, "file1", "Hello ", options=po)
println(resp)

println("Add another file")
resp = S3.put_object(env, bkt, "file2", "World ")
println(resp)

println("Add another file")
resp = S3.put_object(env, bkt, "file3", "Hello World ")
println(resp)

println("List bucket")
resp = S3.get_bkt(env, bkt)
println(resp)

println("Get file 1")
resp = S3.get_object(env, bkt, "file1")
println(resp)

println("Get file 2")
resp = S3.get_object(env, bkt, "file2")
println(resp)

println("Sleep for 10 secs ....")
sleep(10)

println("Delete file 1")
resp = S3.del_object(env, bkt, "file1")
println(resp)

println("Delete file 2 using the multi api")
resp = S3.del_object_multi(env, bkt, S3.DeleteObjectsType([S3.ObjectType("file2"), S3.ObjectType("file3")]))
println(resp)

println("Delete Bucket")
resp = S3.del_bkt(env, bkt)
println(resp)


text:
module BulkAndTails

  using Distributions, ForwardDiff, Ipopt, Roots, Random
  export BulkAndTailsDist, fitbats, fit_bats_mle_covariates, fitbats_covariates, batspdf, batscdf, batsquantile, batslogpdf, batslogcdf, batsrand

  include("BulkAndTailsDist_struct.jl")
  include("BulkAndTailsDist_functions.jl")
  include("fitting.jl")
  include("ipopt_interface.jl")

end # module


text:
#! /usr/bin/env julia

function print_args(;args...)
    # do some formatting
    str = ""
    for (key, val) in args
        str *= string(key, " " , val, " ")
    end
    println(str)
end

function arg_collector(;args...)

    # do some operation on args
    println(args)
    println(typeof(args))
    print_args(;args...)
end

arg_collector(x=5, y=6, z=7)

max_stars_repo_name: sjkelly/JuliaPlayground

max_issues_repo_name: sjkelly/JuliaPlayground

max_forks_repo_name: sjkelly/JuliaPlayground