text:
"""
    readsiggen(file_path::String[, T::Type=Float64])
Read the '*.config' file in 'file_path' for SigGen and returns a dictionary of all parameters.
Non-existing parameteres are set to 0.
...
# Arguments
- `file_path::String`: file path for the SigGen config file.
- `T::Type=Float64`: type of the parameters in the output dictionary.
...
"""
function readsiggen(file_path::String; T::Type=Float64)
    #detector_name = split(basename(file_path), ".config")[1]
    #In order to use the plotting recipe:
    detector_name = "Public Inverted Coax"

    config = Dict("name"       => detector_name,
        "verbosity_level"      => 0,
        "xtal_length"          => 0,
        "xtal_radius"          => 0,
        "top_bullet_radius"    => 0, # not implemented in later conversions
        "bottom_bullet_radius" => 0, # not implemented in later conversions
        "pc_length"            => 0,
        "pc_radius"            => 0,
        "bulletize_PC"         => 0, # not implemented in later conversions
        "wrap_around_radius"   => 0,
        "ditch_depth"          => 0,
        "ditch_thickness"      => 0,
        "bottom_taper_length"  => 0,
        "hole_length"          => 0,
        "hole_radius"          => 0,
        "outer_taper_length"   => 0,
        "inner_taper_length"   => 0,
        "outer_taper_width"    => 0,
        "inner_taper_width"    => 0,
        "taper_angle"          => 0,
        "taper_length"         => 0,
        "Li_thickness"         => 0,
        "energy"               => 0,
        "xtal_grid"            => 0,
        "impurity_z0"          => 0,
        "impurity_gradient"    => 0,
        "impurity_surface"     => 0,
        "xtal_HV"              => 0,
        "max_iterations"       => 0,
        "write_field"          => 0,
        "write_WP"             => 0,
        "drift_name"           => "",
        "field_name"           => "",
        "wp_name"              => "",
        "xtal_temp"            => 0,
        "preamp_tau"           => 0,
        "time_steps_calc"      => 0,
        "step_time_calc"       => 0,
        "step_time_out"        => 0,
        "charge_cloud_size"    => 0,
        "use_diffusion"        => 0,
        "surface_drift_vel_factor" => 0);

    file = open(file_path) do f
        temp          = readlines(f)
        config_file   = Dict()
        for line in temp
            if line  != "" && string(line[1]) != "#"
                line  = strip(split(line, "#")[1])
                name  = ""
                value = ""
                if length(split(line, "\t")) > 1
                    name  = split(line, "\t")[1]
                    value = strip(split(line, "\t")[2])
                else
                    name  = split(line, " ")[1]
                    value = strip(split(line, " ")[end])
                end

                config_file[name] = value

            end
        end
        return config_file
    end

    config_keys = keys(config)
    file_keys   = keys(file)

    for key in file_keys
        if key in config_keys
            if key == "wp_name" || key == "drift_name" || key == "field_name"
                config[key] = string(file[key])
            elseif key == "taper_angle" && parse(Float64, file["taper_angle"]) != 0
                config["outer_taper_width"] = tan(2*pi*parse(T, file[key])/360)*parse(T, file["outer_taper_length"])
                config[key]                 = parse(T, file[key])

            else
                config[key] = parse(T, file[key])
            end
        end

    end
    return config
end

"""
    siggentodict(config::Dict[, units::Dict, detector_name::String])
Converts the dictionary containing the parameters from a SigGen config file
to a SSD config dictionary. This dictionary can be saved as a JSON file using
the JSON package and 'JSON.print(file, config, 4)'.
The 'detector_name' is set to "Public Inverted Coax" by default to inherit
the colour scheme.
...
# Arguments
- `config::Dict`: dictionary containing SigGen parameters (output of readsiggen()).
- `units::Dict`: units used in SigGen file (set to 'mm', 'deg', 'V' and 'K').
- `detector_name::String`: name of the detector.
...
"""
function siggentodict(config::Dict;
        units::Dict = Dict("length" => "mm","angle" => "deg","potential" => "V","temperature" => "K"))

    #>----------------Settings START-------------------------------------------<

    #>----------------Setting the z and r limits-------------------------------<
    z_limit      = ceil(config["xtal_length"]; sigdigits=1);
    if z_limit - config["xtal_length"] <= 0
        z_limit += 10
    end
    r_limit      = ceil(config["xtal_radius"]; sigdigits=1);
    if r_limit - config["xtal_radius"] <= 0
        r_limit += 10
    end

    #>----------------Setting the Axes & Grid dictionaries---------------------<
    axes         = Dict("r"           => Dict("to"         => r_limit,
                                              "boundaries" => "inf"),
                        "phi"         => Dict("from"       => 0,
                                              "to"         => 0,
                                              "boundaries" => "periodic"),
                        "z"           => Dict("from"       => -10,
                                              "to"         => z_limit,
                                              "boundaries" => Dict("left" => "inf", "right" => "inf"))
                        )

    grid         = Dict("coordinates" => "cylindrical",
                        "axes"        => axes);

    #>----------------Settings END---------------------------------------------<
    #<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><
    #>----------------Volume START---------------------------------------------<

    #>----------------Main Volume -> Initial Cylinder--------------------------<
    # r and L are given by the input file. Cylindrical symmetry is set to default
    main_vol     = Dict("name"        => "Initial Cylinder",
                        "type"        => "tube",
                        "r"           => Dict("from" => 0.0, "to" => config["xtal_radius"]),
                        "phi"         => Dict("from" => 0.0, "to" => 360.0),
                        "h"           => config["xtal_length"]);

    #>----------------Upper Cone / Taper on top--------------------------------<
    # If not existing, everything will be set to 0
    upper_cone   = Dict("name"        => "Upper Cone",
                        "type"        => "cone",
                        "r"           => Dict("bottom" => Dict("from" => config["xtal_radius"],
                                                               "to"   => ceil(config["xtal_radius"] + 1.0)),
                                              "top"    => Dict("from" => config["xtal_radius"] -
                                                                            config["outer_taper_width"],
                                                               "to"   => ceil(config["xtal_radius"] + 1.0))),
                        "phi"         => Dict("from"   => 0.0, "to" => 360.0),
                        "h"           => config["outer_taper_length"],
                        "translate"   => Dict("z"      => config["xtal_length"] -
                                                                config["outer_taper_length"]));

    #>----------------Lower Cone / Taper on bottom-----------------------------<
    # 45-degree taper at bottom of ORTEC-type crystal
    lower_cone   = Dict("name"        => "Lower Cone",
                        "type"        => "cone",
                        "r"           => Dict("top"    => Dict("from" => config["xtal_radius"],
                                                               "to"   => ceil(config["xtal_radius"] + 1.0)),
                                              "bottom" => Dict("from" => config["xtal_radius"] -
                                                                            config["taper_length"],
                                                               "to"   => ceil(config["xtal_radius"] + 1.0))),
                        "phi"         => Dict("from"   => 0.0, "to" => 360.0),
                        "h"           => config["taper_length"],
                        "translate"   => Dict("z"      => 0.0));

    #>----------------Ditch around the bottom contact--------------------------<
    # If not existing, everything will be set to 0
    ditch        = Dict("name"        => "Ditch",
                        "type"        => "tube",
                        "r"           => Dict("from" => config["wrap_around_radius"] -
                                                            config["ditch_thickness"],
                                              "to"   => config["wrap_around_radius"]),
                        "phi"         => Dict("from" => 0.0, "to" => 360.0),
                        "h"           => config["ditch_depth"]);

    #>----------------Borehole on bottom as contact----------------------------<
    # If not existing, everything will be set to 0
    borehole     = Dict("name"        => "Borehole",
                        "type"        => "tube",
                        "r"           => Dict("from" => 0.0, "to" => config["pc_radius"]),
                        "phi"         => Dict("from" => 0.0, "to" => 360.0),
                        "h"           => config["pc_length"]);

    #>----------------Hole on top of the detector------------------------------<
    # If not existing, everything will be set to 0
    hole = Dict()
    if config["inner_taper_width"] != 0 && config["inner_taper_length"] != 0
        hole         = Dict("name"        => "Inner Cone",
                            "type"        => "cone",
                            "r"           => Dict("top"    => Dict("from" => 0,
                                                                   "to"   => config["hole_radius"] +
                                                                                config["inner_taper_width"]),
                                                  "bottom" => Dict("from" => 0,
                                                                   "to"   => config["hole_radius"])),
                            "phi"         => Dict("from"   => 0.0, "to" => 360.0),
                            "h"           => config["inner_taper_length"],
                            "translate"   => Dict("z"      => config["xtal_length"] -
                                                                config["inner_taper_length"]));
    else
        hole         = Dict("name"        => "Top Hole",
                            "type"        => "tube",
                            "r"           => Dict("from" => 0.0, "to" => config["hole_radius"]),
                            "phi"         => Dict("from" => 0.0, "to" => 360.0),
                            "h"           => config["hole_length"],
                            "translate"   => Dict("z"    => config["xtal_length"] -
                                                            config["hole_length"]));
    end

    #>----------------Subtract volumes-----------------------------------------<
    geometry = main_vol;
    if config["outer_taper_length"] != 0 || config["outer_taper_width"] != 0
        geometry = Dict("type" => "difference", "parts" => [geometry, upper_cone])
    end
    if config["inner_taper_length"] != 0 || config["inner_taper_width"] != 0 || config["hole_radius"] != 0 || config["hole_length"] != 0
        geometry = Dict("type" => "difference", "parts" => [geometry, hole])
    end
    if config["taper_length"] != 0
        geometry = Dict("type" => "difference", "parts" => [geometry, lower_cone])
    end
    if config["ditch_thickness"] != 0
        geometry = Dict("type" => "difference", "parts" => [geometry, ditch])
    end
    if config["pc_radius"] != 0 && config["pc_length"] >= 0.1 # below this, it is the thickness of the layer
        geometry = Dict("type" => "difference", "parts" => [geometry, borehole])
    end

    #>----------------Volume END-----------------------------------------------<
    #<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><
    #>----------------Contacts START-------------------------------------------<

    geometry_1   = Any[] # gathers all parts of contact_1

    #>----------------Area/volume of point contact-----------------------------<
    #>----------------Contact 1------------------------------------------------<
    # Depending on whether the contact is a borehole or not
    if config["pc_radius"] != 0 && config["pc_length"] >= 5 # below this, it is the thickness of the layer
        borehole_wall = Dict(    "type" => "tube",
             "r"         => Dict("from" => config["pc_radius"],
                                   "to" => config["pc_radius"]),
             "phi"       => Dict("from" => 0.0,
                                   "to" => 360.0),
             "h"         => config["pc_length"])
        borehole_top = Dict(     "type" => "tube",
             "r"         => Dict("from" => 0.0,
                                   "to" => config["pc_radius"]),
             "phi"       => Dict("from" => 0.0,
                                   "to" => 360.0),
             "h"         => 0.0,
             "translate" => Dict("z"    => config["pc_length"]))
        borehole_around = Dict(  "type" => "tube",
             "r"         => Dict("from" => config["pc_radius"],
                                   "to" => config["wrap_around_radius"] - config["ditch_thickness"]),
             "phi"       => Dict("from" => 0.0,
                                   "to" => 360.0),
             "h"         => 0.0,
             "translate" => Dict("z"    => 0.0))
        push!(geometry_1, borehole_wall)
        push!(geometry_1, borehole_top)
        push!(geometry_1, borehole_around)

    elseif config["pc_radius"] != 0 && config["pc_length"] < 5
        pc_vol = Dict(           "type" => "tube",
             "r"         => Dict("from" => 0.0,
                                   "to" => config["pc_radius"]),
             "phi"       => Dict("from" => 0.0,
                                   "to" => 360.0),
             "h"         => config["pc_length"])
        push!(geometry_1, pc_vol)
    end

    #>----------------Define contact 1-----------------------------------------<

    contact_1 = Dict("type"        => "contact",
                     "material"    => "HPGe",
                     "channel"     => 1,
                     "potential"   => 0.0,
                     "geometry"    => Dict("type"  => "union",
                                           "parts" => geometry_1));

    #>----------------Contact 2------------------------------------------------<
    geometry_2 = Any[]

    if config["wrap_around_radius"] != 0.0
        bottom = Dict("type" => "tube",
                      "r"    => Dict("from" => config["wrap_around_radius"],
                                       "to" => config["xtal_radius"] - config["taper_length"]),
                      "phi"  => Dict("from" => 0.0,
                                       "to" => 360.0),
                      "h"    => 0.0)
        push!(geometry_2, bottom)
    end

    wall   = Dict("type" => "tube",
                  "r"    => Dict("from" => config["xtal_radius"],
                                 "to"   => config["xtal_radius"]),
                  "phi"  => Dict("from" => 0.0,
                                   "to" => 360.0),
                  "h"    => config["xtal_length"] - config["taper_length"] - config["outer_taper_length"],
                  "translate" => Dict("z" => config["taper_length"]))
    push!(geometry_2, wall)

    top    = Dict("type" => "tube",
                  "r"    => Dict("from" => config["hole_radius"] + config["inner_taper_width"],
                                   "to" => config["xtal_radius"] - config["outer_taper_width"]),
                  "phi"  => Dict("from" => 0.0,
                                   "to" => 360.0),
                  "h"    => 0.0,
                  "translate" => Dict("z" => config["xtal_length"]))
    push!(geometry_2, top)
    if config["hole_radius"] != 0 && config["hole_length"] != 0 && config["inner_taper_width"] == 0 && config["inner_taper_length"] == 0
        hole   = Dict("type" => "tube",
                      "r"    => Dict("from" => config["hole_radius"],
                                       "to" => config["hole_radius"]),
                      "phi"  => Dict("from" => 0.0,
                                       "to" => 360.0),
                      "h"    => config["hole_length"],
                      "translate" => Dict("z" => config["xtal_length"] - config["hole_length"]))
        inner_hole_bottom    =  Dict("type" => "tube",
                      "r"    => Dict("from" => 0.0,
                                       "to" => config["hole_radius"]),
                      "phi"  => Dict("from" => 0.0,
                                       "to" => 360.0),
                      "h"    => 0.0,
                      "translate" => Dict("z" => config["xtal_length"] - config["hole_length"]))
        push!(geometry_2, hole)
        push!(geometry_2, inner_hole_bottom)
    end
    if config["hole_radius"] != 0 && config["hole_length"] != 0 && config["inner_taper_width"] != 0 && config["inner_taper_length"] != 0
        inner_taper = Dict("type" => "cone",
                      "r"         => Dict("top"    => Dict("from" => config["hole_radius"] +
                                                                        config["inner_taper_width"],
                                                           "to"   => config["hole_radius"] +
                                                                        config["inner_taper_width"]),
                                          "bottom" => Dict("from" => config["hole_radius"],
                                                           "to"   => config["hole_radius"])),
                      "phi"       => Dict("from"   => 0.0, "to" => 360.0),
                      "h"         => config["inner_taper_length"],
                      "translate" => Dict("z"      => config["xtal_length"] -
                                                        config["inner_taper_length"]));
        inner_hole_bottom   =  Dict("type" => "tube",
                      "r"    => Dict("from" => 0.0,
                                       "to" => config["hole_radius"]),
                      "phi"  => Dict("from" => 0.0,
                                       "to" => 360.0),
                      "h"    => 0.0,
                      "translate" => Dict("z" => config["xtal_length"] - config["hole_length"]))
        push!(geometry_2, inner_taper)
        push!(geometry_2, inner_hole_bottom)
    end
    if config["outer_taper_width"] != 0 && config["outer_taper_length"] != 0
        outer_taper = Dict("type" => "cone",
                      "r"         => Dict("top"    => Dict("from" => config["xtal_radius"] -
                                                                        config["outer_taper_width"],
                                                           "to"   => config["xtal_radius"] -
                                                                        config["outer_taper_width"]),
                                          "bottom" => Dict("from" => config["xtal_radius"],
                                                           "to"   => config["xtal_radius"])),
                      "phi"       => Dict("from"   => 0.0, "to" => 360.0),
                      "h"         => config["outer_taper_length"],
                      "translate" => Dict("z"      => config["xtal_length"] -
                                                        config["outer_taper_length"]));
        push!(geometry_2, outer_taper)
    end
    if config["taper_length"] != 0
        taper = Dict("type"       => "cone",
                      "r"         => Dict("top"    => Dict("from" => config["xtal_radius"],
                                                           "to"   => config["xtal_radius"]),
                                          "bottom" => Dict("from" => config["xtal_radius"] -
                                                                        config["taper_length"],
                                                           "to"   => config["xtal_radius"] -
                                                                        config["taper_length"])),
                      "phi"       => Dict("from"   => 0.0, "to" => 360.0),
                      "h"         => config["taper_length"]);
        push!(geometry_2, taper)
    end

    #>----------------Define contact 2-----------------------------------------<

    contact_2 = Dict("type"        => "contact",
                     "material"    => "HPGe",
                     "channel"     => 2,
                     "potential"   => config["xtal_HV"],
                     "geometry"    => Dict("type"  => "union",
                                           "parts" => geometry_2));

    #>----------------Contacts END---------------------------------------------<
    #<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><
    #>----------------Final configuration--------------------------------------<

    impurity_z0            = config["impurity_z0"]*1e10
    impurity_gradient      = config["impurity_gradient"]*1e10
    if units["length"]     == "mm"
        impurity_z0        /= 1e3
        impurity_gradient  /= 1e4
    elseif units["length"] == "m"
        impurity_z0        *= 1e6
        impurity_gradient  *= 1e8
    end

    objects = Any[Dict("type"        => "semiconductor",
                    "material"    => "HPGe",
                    "bulk_type"   => "p",
                    "temperature" => config["xtal_temp"],
                    "charge_density_model" => Dict(
                                "name"     => "linear",
                                "r"        => Dict("init"     => 0.0,
                                                   "gradient" => 0.0),
                                "phi"      => Dict("init"     => 0.0,
                                                   "gradient" => 0.0),
                                "z"        => Dict("init"     => impurity_z0,
                                                   "gradient" => impurity_gradient)
                            ),
                    "geometry"    => geometry
                    ),
                    contact_1,
                    contact_2];

    conf = Dict();
    conf = Dict("name"    => config["name"],
                "units"   => units,
                "grid"    => grid,
                "medium"  => "vacuum",
                "objects" => objects);

    return conf
end


text:
using Garamon
using Documenter

makedocs(;
    modules=[Garamon],
    authors="Jean-François Baffier",
    repo="https://github.com/JuliaConstraints/Garamon.jl/blob/{commit}{path}#L{line}",
    sitename="Garamon.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", nothing) == "true",
        canonical="https://JuliaConstraints.github.io/Garamon.jl",
        assets = ["assets/favicon.ico"; "assets/github_buttons.js"; "assets/custom.css"],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaConstraints/Garamon.jl",
    devbranch="main",
)

ext: jl


text:
using ClangCompiler
using Test

@testset "boot" begin
    @test_nowarn ClangCompiler.print_julia_version()
end


text:
module consoles

export Console, console, err_console, console_height, console_width

const STDOUT = stdout
const STDERR = stderr

# --------------------------------- controls --------------------------------- #

"""
Get cursor position
"""
cursor_position() = cursor_position(stdout)
cursor_position(io::IO) = print(io, "\e[6n")

"""
Move cursor up one line
"""
up() = up(stdout)
up(io::IO) = print(io, "\e[A")

"""
Move cursor to the beginning of the previous line
"""
beginning_previous_line() = beginning_previous_line(stdout)
beginning_previous_line(io::IO) = print(io, "\e[F")

"""
Move cursor down one line
"""
down() = down(stdout)
down(io::IO) = print(io, "\e[B")

"""
Clear terminal.
"""
clear() = clear(stdout)
clear(io::IO) = print(io, "\x1b[2J")

"""
Hide cursor
"""
hide_cursor() = hide_cursor(stdout)
hide_cursor(io::IO) = print(io, "\x1b[?25l")

"""
Show cursor
"""
show_cursor() = show_cursor(stdout)
show_cursor(io::IO) = print(io, "\x1b[?25h")

"""
Print a new line.
"""
line(; i=1) = line(stdout; i=i)
line(io::IO; i=1) = print(io, "\n"^i)

"""
Erase last line in console.
"""
erase_line() = erase_line(stdout)
erase_line(io::IO) = print(io, "\e[2K\r")

# ---------------------------------------------------------------------------- #
#                                    CONSOLE                                   #
# ---------------------------------------------------------------------------- #
"""
    Console

The `Console` object stores information about the dimensions of the output(::IO)
where objects will be printed
"""
struct Console
    io::IO
    width::Int
    height::Int
end

Console(io::IO) = Console(io, displaysize(io)[2], displaysize(io)[1])
Console() = Console(stdout)

console = Console(stdout)
err_console = Console(stderr)

"""
    console_height()

Get the current console height.
"""
console_height() = displaysize(stdout)[1]
console_height(io::IO) = displaysize(io)[1]

"""
    console_width()

Get the current console width.
"""
console_width() = displaysize(stdout)[2]
console_width(io::IO) = displaysize(io)[2]

end


text:
module TestDoctest

using ConcurrentUtils
using Documenter

test() = doctest(ConcurrentUtils; manual = false)

end  # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

julia_version = v"1.6.0"        # Must probably be the same as in libcxxwrap_julia

name = "openPMD_api"
version = v"0.14.2"

# Collection of sources required to complete build
sources = [
    # ArchiveSource("https://github.com/openPMD/openPMD-api/archive/refs/tags/0.13.4.tar.gz",
    #               "46c013be5cda670f21969675ce839315d4f5ada0406a6546a91ec3441402cf5e"),
    # We temporarily use a feature branch instead of a released
    # version because the Julia bindings are not release yet
    ArchiveSource("https://github.com/eschnett/openPMD-api/archive/0df99cc11cfe7555e6d099d62b4016412b650f34.tar.gz",
                  "421538de3a8cba2171a1f2deb5774a340896aa4579e77178a9df27b634fb980d"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd openPMD-api-*
mkdir build
cd build
mpiopts=
if [[ "$target" == *-apple-* ]]; then
    # MPICH's pkgconfig file "mpich.pc" lists these options:
    #     Libs:     -framework OpenCL -Wl,-flat_namespace -Wl,-commons,use_dylibs -L${libdir} -lmpi -lpmpi -lm    -lpthread
    #     Cflags:   -I${includedir}
    # cmake doesn't know how to handle the "-framework OpenCL" option
    # and wants to use "-framework" as a stand-alone option. This fails,
    # and cmake concludes that MPI is not available.
    mpiopts="-DMPI_C_ADDITIONAL_INCLUDE_DIRS='' -DMPI_C_LIBRARIES='-Wl,-flat_namespace;-Wl,-commons,use_dylibs;-lmpi;-lpmpi' -DMPI_CXX_ADDITIONAL_INCLUDE_DIRS='' -DMPI_CXX_LIBRARIES='-Wl,-flat_namespace;-Wl,-commons,use_dylibs;-lmpi;-lpmpi'"
elif [[ "$target" == x86_64-w64-mingw32 ]]; then
    mpiopts="-DMPI_GUESS_LIBRARY_NAME=MSMPI -DMPI_C_LIBRARIES=msmpi64 -DMPI_CXX_LIBRARIES=msmpi64"
elif [[ "$target" == *-mingw* ]]; then
    mpiopts="-DMPI_GUESS_LIBRARY_NAME=MSMPI"
fi
testopts=
if [[ "$target" == *-apple-* ]]; then
    testopts="-DBUILD_TESTING=OFF"
fi
# -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN%.*}_gcc.cmake
cmake \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_FIND_ROOT_PATH=$prefix \
    -DJulia_PREFIX=$prefix \
    -DopenPMD_USE_MPI=ON \
    -DMPI_HOME=$prefix \
    ${mpiopts} \
    ${testopts} \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    ..
cmake --build . --config RelWithDebInfo --parallel $nproc
cmake --build . --config RelWithDebInfo --parallel $nproc --target install
install_license ../COPYING*
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# platforms = supported_platforms()
# Use only platforms where libcxxwrap_julia is supported.
include("../../L/libjulia/common.jl")
platforms = libjulia_platforms(julia_version)
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libopenPMD", :libopenPMD),
    LibraryProduct("libopenPMD_jl", :libopenPMD_jl),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version)),
    Dependency(PackageSpec(name="ADIOS2_jll")),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    # We would need a parallel version of HDF5
    # Dependency(PackageSpec(name="HDF5_jll")),
    Dependency(PackageSpec(name="MPICH_jll")),
    Dependency(PackageSpec(name="MicrosoftMPI_jll")),
    Dependency(PackageSpec(name="libcxxwrap_julia_jll")),
]

# See <https://github.com/JuliaPackaging/Yggdrasil/blob/master/Q/Qt5Base/build_tarballs.jl> for building on macOS 10.14

# Build the tarballs, and possibly a `build.jl` as well.
# We need C++14, which requires at least GCC 5.
# GCC 5 reports incompatible signatures for `posix_memalign` on linux/musl, fixed on GCC 6
# GCC 5 has a bug regarding `std::to_string` on freebsd, fixed on GCC 6
# macos encounters an ICE in GCC 6; switching to GCC 7 instead
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version = v"7")

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
using ClipData
using Test
using InteractiveUtils: clipboard
using Tables

@testset "cliptable" begin
    """
    a b
    1 2
    3 4
    """ |> clipboard

    @test Tables.columntable(cliptable()) == (a = [1, 3], b = [2, 4])

    """
    a\tb
    1\t2
    3\t4
    """ |> clipboard

    @test Tables.columntable(cliptable()) == (a = [1, 3], b = [2, 4])

    t = (a = [1, 3], b = [2, 4])
    cliptable(t)
    @test clipboard() == "a\tb\n1\t2\n3\t4"
end

@testset "cliparray" begin
    """
    1 2
    3 4
    """ |> clipboard

    @test cliparray() == [1 2; 3 4]

    """
    1\t2
    3\t4
    """ |> clipboard

    @test cliparray() == [1 2; 3 4]

    """
    1
    2
    3
    4
    """ |> clipboard

    @test cliparray() == [1, 2, 3, 4]

    X = [1 2; 3 4]
    cliparray(X)
    @test clipboard() == "1\t2\n3\t4"

    x = [1, 2, 3, 4]
    cliparray(x)
    @test clipboard() == "1\n2\n3\n4"

    """
    1 2 3 4
    """ |> clipboard

    @test cliparray() == [1, 2, 3, 4]
end

@testset "mwetable" begin
    """
    a b
    1 2
    3 4
    """ |> clipboard

    s = mwetable(; returnstring=true)
    s_correct =
"""
df = \"\"\"
a,b
1,2
3,4
\"\"\" |> IOBuffer |> CSV.File"""

    @test s == s_correct

    t = (a = [1, 3], b = [2, 4])
    s = mwetable(t; returnstring = true)

    @test s == s_correct
end

@testset "mwearray" begin
    """
    1 2
    3 4
    """ |> clipboard

    s = mwearray(; returnstring=true)
    s_correct =
"""
X = \"\"\"
1,2
3,4
\"\"\" |> IOBuffer |> CSV.File |> Tables.matrix"""

    @test s == s_correct

    t = [1 2; 3 4]
    s = mwearray(t; returnstring = true)

    @test s == s_correct

    """
    1
    2
    3
    4
    """ |> clipboard

    s = mwearray(; returnstring=true)

    s_correct =
"""
x = \"\"\"
1
2
3
4
\"\"\" |> IOBuffer |> CSV.File |> Tables.matrix |> vec"""

    @test s == s_correct

    x = [1, 2, 3, 4]

    s = mwearray(x; returnstring=true)
    @test s == s_correct
end

@testset "@mwetable" begin
    # Can't think of a way to test. Just check
    # for errors.

    mytable = (a = [1, 2], b = [3, 4])

    @mwetable mytable
end

@testset "@mwearray" begin
    # Can't think of a way to test. Just check
    # for errors.

    myarray = [1 2; 3 4]

    @mwearray myarray

    myvector = [1, 2, 3, 4]

    @mwearray myvector
end

max_forks_repo_name: pdeffebach/DataClip


text:
export tod2map, combinemaps!

@doc raw"""
    tod2map{T,O}(pixidx, tod::Array{T}; nside=128) :: (map, hits)

Create a binned map for a TOD and return a tuple containing the map
itself and the hit map.
"""
function tod2map(pixidx, tod::Array{T}; nside=128, ordering=Healpix.Ring) where {T}
    @assert length(pixidx) == length(tod)

    binnedmap = Map{T,ordering}(nside)
    hitmap = Map{Int,ordering}(nside)

    @inbounds for i in eachindex(pixidx)
        if !ismissing(tod[i])
            binnedmap.pixels[pixidx[i]] += tod[i]
            hitmap.pixels[pixidx[i]] += 1
        end
    end

    @inbounds for i in eachindex(binnedmap.pixels)
        if hitmap.pixels[i] > 0
            binnedmap.pixels[i] /= hitmap.pixels[i]
        end
    end

    (binnedmap, hitmap)
end

@doc raw"""
    combinemaps{T, O, H}(destmap::Map{T, O}, desthitmap::Map{H, O}, othermap::Map{T, O}, otherhitmap::Map{H, O})

Sum "othermap" to "destmap", assuming that both maps have been
produced by binning TODs. The parameters `desthitmap` and
`otherhitmap` are the two hit maps. At the end of the call, `destmap`
and `desthitmap` are updated.
"""
function combinemaps!(destmap::Map{T,O},
                      desthitmap::Map{H,O},
                      othermap::Map{T,O},
                      otherhitmap::Map{H,O}) where {T,O,H}
    (conformables(destmap, othermap)
     && conformables(desthitmap, otherhitmap)
     && conformables(destmap, desthitmap)) || throw(DomainError())

    for i in eachindex(destmap.pixels)
        (pix1, pix2) = (destmap.pixels[i], othermap.pixels[i])
        (hit1, hit2) = (desthitmap.pixels[i], otherhitmap.pixels[i])
        hitsum = hit1 + hit2
        if hitsum > 0
            destmap.pixels[i] = (pix1 * hit1 + pix2 * hit2) / hitsum
        end
        desthitmap.pixels[i] = hitsum
    end
end


text:
function applyheartics!(system::CVSystem,old=Dict("a"=>0),restart="no")
    if restart == "no"
        # volumes based on end-diastole
        system.heart.lv.V[1] = 120*cm3Tom3;
        # system.heart.rv.V[1] = 125*cm3Tom3;
        system.heart.rv.V[1] = 125*cm3Tom3;
        # system.heart.la.V[1] = 36*cm3Tom3;
        system.heart.la.V[1] = 36.5*cm3Tom3;
        system.heart.ra.V[1] = 32.7*cm3Tom3;
        # system.heart.ra.V[1] = 35.2*cm3Tom3;

        # aortic valve state
        system.heart.av.zeta[1] = 0;
    elseif restart == "yes"
        lv = old["lv"];
        rv = old["rv"];
        la = old["la"];
        ra = old["ra"];
        av = old["av"];
        system.heart.lv.V[1] = lv["V"][end];
        system.heart.rv.V[1] = rv["V"][end];
        system.heart.la.V[1] = la["V"][end];
        system.heart.ra.V[1] = ra["V"][end];
        system.heart.av.zeta[1] = av["zeta"][end];
    end

    # associated pressures
    CVModule.elastancemodel!(system,0);
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
export BrickletIO4Configuration
struct BrickletIO4Configuration
    direction_mask::Integer
    value_mask::Integer
end

export BrickletIO4Monoflop
struct BrickletIO4Monoflop
    value::Integer
    time::Integer
    time_remaining::Integer
end

export BrickletIO4EdgeCountConfig
struct BrickletIO4EdgeCountConfig
    edge_type::Integer
    debounce::Integer
end

export BrickletIO4Identity
struct BrickletIO4Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletIO4
"""
4-channel digital input/output
"""
mutable struct BrickletIO4 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletIO4(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_io4")
        deviceInternal = package.BrickletIO4(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_value
"""
    $(SIGNATURES)

Sets the output value (high or low) with a bitmask (4bit). A 1 in the bitmask
means high and a 0 in the bitmask means low.

For example: The value 3 or 0b0011 will turn the pins 0-1 high and the
pins 2-3 low.

All running monoflop timers will be aborted if this function is called.

.. note::
 This function does nothing for pins that are configured as input.
 Pull-up resistors can be switched on with :func:`Set Configuration`.
"""
function set_value(device::BrickletIO4, value_mask)
    device.deviceInternal.set_value(value_mask)
end

export get_value
"""
    $(SIGNATURES)

Returns a bitmask of the values that are currently measured.
This function works if the pin is configured to input
as well as if it is configured to output.
"""
function get_value(device::BrickletIO4)
    return device.deviceInternal.get_value()
end

export set_configuration
"""
    $(SIGNATURES)

Configures the value and direction of the specified pins. Possible directions
are 'i' and 'o' for input and output.

If the direction is configured as output, the value is either high or low
(set as *true* or *false*).

If the direction is configured as input, the value is either pull-up or
default (set as *true* or *false*).

For example:

* (15, 'i', true) or (0b1111, 'i', true) will set all pins of as input pull-up.
* (8, 'i', false) or (0b1000, 'i', false) will set pin 3 of as input default (floating if nothing is connected).
* (3, 'o', false) or (0b0011, 'o', false) will set pins 0 and 1 as output low.
* (4, 'o', true) or (0b0100, 'o', true) will set pin 2 of as output high.

Running monoflop timers for the specified pins will be aborted if this
function is called.
"""
function set_configuration(device::BrickletIO4, selection_mask, direction, value)
    device.deviceInternal.set_configuration(selection_mask, direction, value)
end

export get_configuration
"""
    $(SIGNATURES)

Returns a value bitmask and a direction bitmask. A 1 in the direction bitmask
means input and a 0 in the bitmask means output.

For example: A return value of (3, 5) or (0b0011, 0b0101) for direction and
value means that:

* pin 0 is configured as input pull-up,
* pin 1 is configured as input default,
* pin 2 is configured as output high and
* pin 3 is are configured as output low.
"""
function get_configuration(device::BrickletIO4)
    return device.deviceInternal.get_configuration()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the debounce period of the :cb:`Interrupt` callback.

For example: If you set this value to 100, you will get the interrupt
maximal every 100ms. This is necessary if something that bounces is
connected to the IO-4 Bricklet, such as a button.
"""
function set_debounce_period(device::BrickletIO4, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletIO4)
    return device.deviceInternal.get_debounce_period()
end

export set_interrupt
"""
    $(SIGNATURES)

Sets the pins on which an interrupt is activated with a bitmask.
Interrupts are triggered on changes of the voltage level of the pin,
i.e. changes from high to low and low to high.

For example: An interrupt bitmask of 10 or 0b1010 will enable the interrupt for
pins 1 and 3.

The interrupt is delivered with the :cb:`Interrupt` callback.
"""
function set_interrupt(device::BrickletIO4, interrupt_mask)
    device.deviceInternal.set_interrupt(interrupt_mask)
end

export get_interrupt
"""
    $(SIGNATURES)

Returns the interrupt bitmask as set by :func:`Set Interrupt`.
"""
function get_interrupt(device::BrickletIO4)
    return device.deviceInternal.get_interrupt()
end

export set_monoflop
"""
    $(SIGNATURES)

Configures a monoflop of the pins specified by the first parameter as 4 bit
long bitmask. The specified pins must be configured for output. Non-output
pins will be ignored.

The second parameter is a bitmask with the desired value of the specified
output pins. A 1 in the bitmask means high and a 0 in the bitmask means low.

The third parameter indicates the time that the pins should hold
the value.

If this function is called with the parameters (9, 1, 1500) or
(0b1001, 0b0001, 1500): Pin 0 will get high and pin 3 will get low. In 1.5s pin
0 will get low and pin 3 will get high again.

A monoflop can be used as a fail-safe mechanism. For example: Lets assume you
have a RS485 bus and an IO-4 Bricklet connected to one of the slave
stacks. You can now call this function every second, with a time parameter
of two seconds and pin 0 set to high. Pin 0 will be high all the time. If now
the RS485 connection is lost, then pin 0 will get low in at most two seconds.
"""
function set_monoflop(device::BrickletIO4, selection_mask, value_mask, time)
    device.deviceInternal.set_monoflop(selection_mask, value_mask, time)
end

export get_monoflop
"""
    $(SIGNATURES)

Returns (for the given pin) the current value and the time as set by
:func:`Set Monoflop` as well as the remaining time until the value flips.

If the timer is not running currently, the remaining time will be returned
as 0.
"""
function get_monoflop(device::BrickletIO4, pin)
    return device.deviceInternal.get_monoflop(pin)
end

export set_selected_values
"""
    $(SIGNATURES)

Sets the output value (high or low) with a bitmask, according to
the selection mask. The bitmask is 4 bit long, *true* refers to high
and *false* refers to low.

For example: The parameters (9, 4) or (0b0110, 0b0100) will turn
pin 1 low and pin 2 high, pin 0 and 3 will remain untouched.

Running monoflop timers for the selected pins will be aborted if this
function is called.

.. note::
 This function does nothing for pins that are configured as input.
 Pull-up resistors can be switched on with :func:`Set Configuration`.
"""
function set_selected_values(device::BrickletIO4, selection_mask, value_mask)
    device.deviceInternal.set_selected_values(selection_mask, value_mask)
end

export get_edge_count
"""
    $(SIGNATURES)

Returns the current value of the edge counter for the selected pin. You can
configure the edges that are counted with :func:`Set Edge Count Config`.

If you set the reset counter to *true*, the count is set back to 0
directly after it is read.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function get_edge_count(device::BrickletIO4, pin, reset_counter)
    return device.deviceInternal.get_edge_count(pin, reset_counter)
end

export set_edge_count_config
"""
    $(SIGNATURES)

Configures the edge counter for the selected pins.

The edge type parameter configures if rising edges, falling edges or
both are counted if the pin is configured for input. Possible edge types are:

* 0 = rising (default)
* 1 = falling
* 2 = both

Configuring an edge counter resets its value to 0.

If you don't know what any of this means, just leave it at default. The
default configuration is very likely OK for you.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function set_edge_count_config(device::BrickletIO4, selection_mask, edge_type, debounce)
    device.deviceInternal.set_edge_count_config(selection_mask, edge_type, debounce)
end

export get_edge_count_config
"""
    $(SIGNATURES)

Returns the edge type and debounce time for the selected pin as set by
:func:`Set Edge Count Config`.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function get_edge_count_config(device::BrickletIO4, pin)
    return device.deviceInternal.get_edge_count_config(pin)
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletIO4)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletIO4, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
@inline function cti(cone::Cone, i)
	return i + cone.offs
end

#e vector generation

@unroll function make_e!(cones::Tuple{Vararg{Cone}}, r)
	@unroll for c in cones
		make_e!(c, r)
	end
end

function make_e!(cone::POC{dim}, r) where dim
	for i=cti(cone, 1):cti(cone, dim)
		r[i] = 1.0
	end
end

function make_e!(cone::SOC{dim}, r) where dim
	r[cti(cone,1)] = 1.0
	for i=cti(cone, 2):cti(cone, dim)
		r[i] = 0.0
	end
end

function make_e(cone::Union{SOC,POC})
	res = zeros(cone.offs + conedim(cone))
	make_e!(cone, res)
	return @view res[cti(cone, 1):cti(cone, conedim(cone))]
end

function make_e(cones::Tuple{Vararg{Cone}})
	tvpres = zeros(cones[end].offs + conedim(cones[end]))
	make_e!(cones, tvpres)
	return tvpres
end

# vector products

function vprod(cone::POC{dim}, u, v) where dim
	res = zeros(dim+cone.offs)
	vprod!(cone, res, u, v)
	return res
end

function vprod(cone::SOC{dim}, u, v) where dim
	res = zeros(dim+cone.offs)
	vprod!(cone, res, u, v)
	return res
end

@unroll function vprod(cones::Tuple{Vararg{Cone}}, u, v)
	tvpres = zeros(cones[end].offs + conedim(cones[end]))
	vprod!(cones, tvpres, u, v)
	return tvpres
end

function vprod!(cone::POC{dim}, t, u, v) where dim
	for i in cti(cone, 1):cti(cone, dim)
		t[i] = u[i]*v[i]
	end
end

function vprod!(cone::SOC{dim}, t, u, v) where dim
	i1 = cti(cone, 1)
	t[i1] = 0.0
	for i = cti(cone, 1):cti(cone, dim)
		t[i1] += u[i] * v[i]
	end
	iu = u[i1]
	iv = v[i1]
	for i in cti(cone, 2):cti(cone, dim)
		t[i] = iu * v[i] + iv * u[i]
	end
end

@unroll function vprod!(cones::Tuple{Vararg{Cone}}, t, u, v)
	@unroll for cone in cones
		vprod!(cone, t, u, v)
	end
end

# inverse vector products

function iprod(cone::Cone{dim}, lam::Vector{Float64}, v::Vector{Float64}) where dim
	tres = zeros(dim+cone.offs)
	iprod!(cone, tres, lam, v)
	return tres
end

@unroll function iprod(cones::Tuple{Vararg{Cone}}, u, v)
	tres = zeros(conedim(cones[end])+cones[end].offs)
	@unroll for c in cones
		iprod!(c, tres, u, v)
	end
	return tres
end

function iprod!(cone::POC{dim}, t::Vector{Float64}, lam::Vector{Float64}, v::Vector{Float64}) where dim
	for i=cti(cone, 1):cti(cone, dim)
		t[i] = v[i]/lam[i]
	end
end

function iprod!(cone::SOC{dim}, t, lam::Vector{Float64}, v::Vector{Float64}) where dim
	i1 = cti(cone, 1)
	l1 = lam[i1]
	a = l1^2
	for i=cti(cone, 2):cti(cone, dim)
		a -= lam[i] * lam[i]
	end
	for i=i1:cti(cone, dim)
		t[i] = 0
	end
	t[i1] += v[i1] * l1/a
	for j=cti(cone, 2):cti(cone, dim)
		t[i1] -= v[j] * lam[j]/a
	end
	for i=cti(cone, 2):cti(cone, dim)
		t[i] -= v[i1] * lam[i]/a
		for j=cti(cone, 2):cti(cone, dim)
			t[i] += v[j] * ((i == j ? a : 0.0) + lam[i]*lam[j])/(l1*a)
		end
	end
end

@unroll function iprod!(cones::Tuple{Vararg{Cone}}, t, u, v)
	@unroll for cone in cones
		iprod!(cone, t, u, v)
	end
end

# per-cone gt defns

function cgt(c::POC{dim}, x, dx) where dim
	for i in cti(c, 1):cti(c, dim)
		if x[i] + dx[i] < 0
			return false
		end
	end
	return true
end

function cgt(c::SOC{dim}, x, dx) where dim
	tot = 0.0
	for i in cti(c, 2):cti(c, dim)
		val = x[i] + dx[i]
		tot += val * val
	end
	return sqrt(tot) <= x[cti(c, 1)] + dx[cti(c, 1)]
end

@unroll function cgt(cs::Tuple{Vararg{Cone}}, x, dx)
	@unroll for c in cs
		if !cgt(c, x, dx)
			return false
		end
	end
	return true
end

# degree

function deg(c::POC{d}) where d
	return d
end

function deg(c::SOC)
	return 1
end

@unroll function deg(cs::Tuple{Vararg{Cone}})
	dg = 0
	@unroll for c in cs
		dg += deg(c)
	end
	return dg
end


text:
# Base case to test the distribution network model.
#
# Runs a single-period problem without storage and flexible loads.
# Uses the original CIGRE network with the following changes:
# - added generator cost: all have equal prices, HV grid exchanges cost twice;
# - a fixed 1.0 tap ratio is assigned to the transformer of branch 17;
# - added some candidate branches.

## Import packages and choose a solver

import PowerModels; const _PM = PowerModels
import FlexPlan; const _FP = FlexPlan
import Cbc
optimizer = _FP.optimizer_with_attributes(Cbc.Optimizer, "logLevel"=>0)

## Settings

bus_1_voltage_fixed        = false
gen_7_power_reactive_fixed = false
branch_1_impedance_zeroed  = false
branch_12_rating_zeroed    = false
branch_16_tap_ratio_fixed  = false
load_bus_1_big             = false
load_bus_12_big            = false
load_bus_13_big            = false
load_bus_14_big            = false

result_file = "test/data/output_files/test_dist_nw_model_result.txt"

## Load and preprocess data

data_file = "test/data/cigre_mv_eu/cigre_mv_eu_unit_test.m"
data = _FP.parse_file(data_file)
_FP.add_dimension!(data, :hour, 1)
_FP.add_dimension!(data, :year, 1)
data = _FP.make_multinetwork(data)

## Apply changes for testing purpose

branch = data["nw"]["1"]["branch"]
bus    = data["nw"]["1"]["bus"]
gen    = data["nw"]["1"]["gen"]
load   = data["nw"]["1"]["load"]
if bus_1_voltage_fixed
    bus["1"]["vmin"] = 1.0
    bus["1"]["vmax"] = 1.0
end
if gen_7_power_reactive_fixed
    gen["7"]["qmin"] = 0.0
    gen["7"]["qmax"] = 0.0
end
if branch_1_impedance_zeroed
    branch["1"]["br_r"] = 0.0
    branch["1"]["br_x"] = 0.0
end
if branch_12_rating_zeroed
    branch["12"]["rate_a"] = 0.0
    branch["12"]["rate_b"] = 0.0
    branch["12"]["rate_c"] = 0.0
end
if branch_16_tap_ratio_fixed
    branch["16"]["tm_min"] = 0.95
    branch["16"]["tap"]    = 0.95
    branch["16"]["tm_max"] = 0.95
end
if load_bus_1_big
    load["1"]["pd"] += 10.0
end
if load_bus_12_big
    load["11"]["pd"] += 10.0
end
if load_bus_13_big
    load["12"]["pd"] += 8.0
end
if load_bus_14_big
    load["13"]["pd"] += 1.0
end

## Solve problem

result = _FP.strg_tnep(data, _FP.BFARadPowerModel, optimizer)
# Two-step alternative
#pm = _PM.instantiate_model(data, _FP.BFARadPowerModel, _FP.post_strg_tnep; ref_extensions=[_FP.ref_add_gen!, _FP.ref_add_storage!, _FP.ref_add_ne_storage!, _PM.ref_add_on_off_va_bounds!, _FP.ref_add_ne_branch_allbranches!, _FP.ref_add_frb_branch!, _FP.ref_add_oltc_branch!])
#result = _PM.optimize_model!(pm; optimizer, solution_processors=[_PM.sol_data_model!])
@assert result["termination_status"] ∈ (_PM.OPTIMAL, _PM.LOCALLY_SOLVED) "$(result["optimizer"]) termination status: $(result["termination_status"])"

## Write result

_PM.print_summary(result["solution"]["nw"]["1"])
mkpath(dirname(result_file))
open(result_file, "w") do io
    _PM.summary(io, result["solution"]["nw"]["1"])
end

## Perform other unit tests

gen = result["solution"]["nw"]["1"]["gen"]

pg = sum(g["pg"] for g in values(gen))
pd = sum(l["pd"] for l in values(load))
@assert pg ≈ pd "generated active power ($pg) does not match demand ($pd)"

qg = sum(g["qg"] for g in values(gen))
qd = sum(l["qd"] for l in values(load))
@assert qg ≈ qd "generated reactive power ($qg) does not match demand ($qd)"


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -6.019960436071491 <= q <= 11.696817216382863)

                     add_NL_constraint(m, :(sigmoid(0.9042935565441375 + 0.31236426273908524*sigmoid(0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4])) + -0.4128166127064836*sigmoid(0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4])) + -0.2476169775927386*sigmoid(0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))) + sigmoid(0.8600467319105785 + -0.4006275013725942*sigmoid(0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4])) + 0.22106887260345687*sigmoid(0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4])) + 0.7210213948433983*sigmoid(0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))) + sigmoid(0.982756228863964 + 0.34792310275142935*sigmoid(0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4])) + 0.7924515917531365*sigmoid(0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4])) + -0.9399735455074207*sigmoid(0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using IntervalSets
using Base.Test

@testset "IntervalSets" begin
    io = IOBuffer()
    @test ordered(2, 1) == (1, 2)
    @test ordered(1, 2) == (1, 2)
    @test ordered(Float16(1), 2) == (1, 2)

    @testset "Closed Sets" begin
        @test_throws ArgumentError :a .. "b"
        I = 0..3
        print(io, I)
        @test takebuf_string(io) == "0..3"
        J = 3..2
        K = 5..4
        L = 3 ± 2
        M = ClosedInterval(2, 5.0)
        print(io, M)
        @test takebuf_string(io) == "2.0..5.0"
        N = ClosedInterval(UInt8(255), 300)
        O = CartesianIndex(1, 2, 3, 4) ± 2
        @test O == (-1..3, 0..4, 1..5, 2..6)

        @test eltype(I) == Int
        @test eltype(M) == Float64
        @test convert(ClosedInterval{Float64}, I) === 0.0..3.0
        @test !(convert(ClosedInterval{Float64}, I) === 0..3)
        @test ClosedInterval{Float64}(1,3) === 1.0..3.0
        @test ClosedInterval(0.5..2.5) === 0.5..2.5
        @test ClosedInterval{Int}(1.0..3.0) === 1..3

        @test !isempty(I)
        @test isempty(J)
        @test J == K
        @test I != K
        @test I == I
        @test J == J
        @test L == L
        @test isequal(I, I)
        @test isequal(J, K)

        @test typeof(M.left) == typeof(M.right) && typeof(M.left) == Float64
        @test typeof(N.left) == typeof(N.right) && typeof(N.left) == Int

        @test maximum(I) === 3
        @test minimum(I) === 0

        @test 2 in I
        @test 1..2 in 0.5..2.5

        @test I ∪ L == ClosedInterval(0, 5)
        @test I ∩ L == ClosedInterval(1, 3)
        @test isempty(J ∩ K)
        @test isempty((2..5) ∩ (7..10))
        @test isempty((1..10) ∩ (7..2))
        A = Float16(1.1)..Float16(1.234)
        B = Float16(1.235)..Float16(1.3)
        C = Float16(1.236)..Float16(1.3)
        D = Float16(1.1)..Float16(1.236)
        @test A ∪ B == Float16(1.1)..Float16(1.3)
        @test D ∪ B == Float16(1.1)..Float16(1.3)
        @test D ∪ C == Float16(1.1)..Float16(1.3)
        @test_throws(ArgumentError, A ∪ C)

        @test J ⊆ L
        @test (L ⊆ J) == false
        @test K ⊆ I
        @test ClosedInterval(1, 2) ⊆ I
        @test I ⊆ I
        @test (ClosedInterval(7, 9) ⊆ I) == false
        @test I ⊇ I
        @test I ⊇ ClosedInterval(1, 2)

        @test hash(1..3) == hash(1.0..3.0)
    end
end


text:
# The structure of the fractional derivative:
# (1) Type definition.
# (2) Checks function to ensure the input parameters are correct and handle some special cases, for example, 0 end point case.
# (3) Multiple dispatch for fracdiff of different algorithms, the first is used to compute the value of a specific point, the second is used to compute the fractional derivative of a Vector.
# (4) Fractional derivative macros for convenient computing.

using QuadGK, SpecialFunctions
using LinearAlgebra, InvertedIndices

"""
Base type of fractional differentiation algorithms, in FractionalCalculus.jl, all of the fractional derivative algorithms belong to ```FracDiffAlg```
"""
abstract type FracDiffAlg end

"""
Caputo sense fractional derivative algorithms, please refer to [Caputo derivative](https://en.wikipedia.org/wiki/Fractional_calculus#Caputo_fractional_derivative) for more details.
"""
abstract type Caputo <: FracDiffAlg end

"""
Grünwald–Letnikov sense fractional derivative algorithms, please refer to [Grünwald–Letnikov derivative](https://en.wikipedia.org/wiki/Gr%C3%BCnwald%E2%80%93Letnikov_derivative) for more details
"""
abstract type GL <: FracDiffAlg end

"""
Riemann-Liouville sense fractional derivative algorithms, please refer to [Riemann-Liouville derivative](https://en.wikipedia.org/wiki/Fractional_calculus#Riemann%E2%80%93Liouville_fractional_derivative)
"""
abstract type RLDiff <: FracDiffAlg end

"""
Note *Caputo Direct* algorithms belong to direct computing, precise are ensured, but maybe cause more memory allocation and take more compilation time.

Using the direct mathematic expression:

```math
^CD_t^α=\\frac{1}{\\Gamma(n-α)}\\int_0^t\\frac{f^{(n)}(τ)}{(t-τ)^{α+1-n}}
```

As for the derivative inside the integral, we use the **Complex Step Differentiation** to obtain a more accurate value.
"""
struct Caputo_Direct <: Caputo end

"""
With first derivative known, we can direct use the derivative to compute a more accurate result.
"""
struct Caputo_Direct_First_Diff_Known <: Caputo end

"""
With first and second derivative known, we can direct use derivative and second derivative to compute a more accurate result.
"""
struct Caputo_Direct_First_Second_Diff_Known <: Caputo end

"""
Grunwald Letnikov direct compute method to obtain fractional derivative, precision are guaranteed but cause more memory allocation and compilation time.
"""
struct GL_Direct <: GL end

"""
@article{LI20113352,
title = {Numerical approaches to fractional calculus and fractional ordinary differential equation},
author = {Changpin Li and An Chen and Junjie Ye},
}
"""

"""
Using piecewise linear interpolation function to approximate input function and combining Caputo derivative then implement summation.
"""
struct Caputo_Piecewise <: Caputo end

"""
@article{10.1093/imanum/17.3.479,
author = {DIETTELM, KAI},
title = "{Generalized compound quadrature formulae for finite-part integrals}",
journal = {IMA Journal of Numerical Analysis},
doi = {10.1093/imanum/17.3.479},
}
"""

"""
Diethelm's method for computingn Caputo sense fractional derivative using [product trapezoidal rule](https://en.wikipedia.org/wiki/Trapezoidal_rule) and [Richardsin extrapolation](https://en.wikipedia.org/wiki/Richardson_extrapolation)
"""
struct Caputo_Diethelm <: Caputo end

"""
@book{oldham_spanier_1984,
title={The fractional calculus: Theory and applications of differentiation and integration to arbitrary order},
author={Oldham, Keith B. and Spanier, Jerome},
year={1984}}
"""

"""
Using Grünwald–Letnikov multiplication-addition-multiplication-addition··· method to compute fractional derivative in Grünwald Letnikov sense.
"""
struct GL_Multiplicative_Additive <: GL end

"""
Using Lagrange three points interpolation method to compute fractional derivative in Grünwald Letnikov sense.
"""
struct GL_Lagrange_Three_Point_Interp <: GL end

#RLDiff_LinearInterp maybe??
"""
Using Linear interpolation to approximate fractional derivative in Riemann Liouville  fractional derivative sense.
"""
struct RLDiff_Approx <: RLDiff end
struct RLDiff_Approx_2 <: RLDiff end
struct RLDiff_Approx_3 <: RLDiff end
struct RLDiff_Approx_Any <: RLDiff end

"""
@article{2009,
title={Matrix approach to discrete fractional calculus II: Partial fractional differential equations},
DOI={10.1016/j.jcp.2009.01.014},
author={Podlubny, Igor and Chechkin, Aleksei and Skovranek, Tomas and Chen, YangQuan and Vinagre Jara, Blas M.},
}

"""

"""
Using [Triangular Strip Matrix](https://en.wikipedia.org/wiki/Triangle_strip) to discrete the derivative.
"""
struct RLDiff_Matrix <: RLDiff end

"""
Using **Grünwald–Letnikov finite difference method** to compute Grünwald–Letnikov sense fractional derivative.
"""
struct GL_Finite_Difference <: GL end

################################################################
###                    Type defination done                  ###
################################################################

"""
Check if the format of nargins is correct.
"""
function checks(f::Union{Function, Number}, α::Float64, start_point::Float64, end_point::Float64)
    α % 1 != 0 ? nothing : error("α must be a decimal number")
    if isa(end_point, Number)
        start_point < end_point ? nothing : DomainError("Domain error! Start point must smaller than end point")
    elseif isa(end_point, AbstractArray)
        start_point < minimum(end_point) ? nothing : DomainError("Vector domain error! Start point must smaller than end point")
    else
        ErrorException("Please input correct point you want to compute")
    end

    end_point == 0 ? 0 : nothing
end

function checks(f::Number, α::Float64, start_point::Float64, end_point::Float64)
    f = f == 0 ? 0 : f/sqrt(pi*end_point)

end

"""
# FracDiffAlg

    fracdiff(f::Function, α, start_point, end_point, FracDiffAlg())

### Example

```julia-repl
julia> fracdiff(Function, Order, Start_Point, End_Point, AlgType)
```

Compute the α-order fractional derivative of f from start point to end point with specific algorithm.
"""
function fracdiff(f, α, end_point, h)
    return fracdiff(f, α, end_point, h, RLDiff_Approx())
end

"""
# Caputo sense fractional derivative.

    fracdiff(f::Function, α, start_point, end_point, step_size ::Caputo)

### Example:

```julia-repl
julia> fracdiff(x->x^5, 0.5, 0, 2.5, 0.0001, Caputo_Direct())
```

Returns a tuple (**result**, **error**), which means the value of this derivative is 141.59714979764541, and the error estimate is 1.1532243848672914e-6.

When the input end points is an array, **fracdiff** will compute

Refer to [Caputo derivative](https://en.wikipedia.org/wiki/Fractional_calculus#Caputo_fractional_derivative)
"""
function fracdiff(f::Union{Function, Number}, α::Float64, start_point, end_point, h::Float64, ::Caputo_Direct)
    checks(f, α, start_point, end_point)

    #Using complex step differentiation to calculate the first order differentiation
    g(τ) = imag(f(τ+1*im*h) ./ h) ./ ((end_point-τ) .^α)
    result = quadgk(g, start_point, end_point) ./gamma(1-α)
    return result
end

function fracdiff(f::Union{Number, Function}, α::Float64, start_point, end_point::AbstractArray, h::Float64, ::Caputo_Direct)::Vector
    ResultArray = Float64[]
    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracdiff(f, α, start_point, value, h, Caputo_Direct()))
    end
    return ResultArray
end

"""
# Grünwald–Letnikov sense fractional dervivative.

    fracdiff(f, α, start_point, end_point, GL_Direct())

### Example:

```julia-repl
julia> fracdiff(x->x^5, 0, 0.5, GL_Direct())
```

> Please note Grunwald-Letnikov sense fracdiff only support 0 < α < 1.

Please refer to [Grünwald–Letnikov derivative](https://en.wikipedia.org/wiki/Gr%C3%BCnwald%E2%80%93Letnikov_derivative) for more details.
"""
function fracdiff(f::Union{Function, Number}, α::Float64, start_point, end_point, ::GL_Direct)
    checks(f, α, start_point, end_point)

    g(τ) = (f(end_point)-f(τ)) ./ (end_point - τ) .^ (1+α)
    result = f(end_point)/(gamma(1-α) .* (end_point - start_point) .^ α) .+ quadgk(g, start_point, end_point) .* α ./ gamma(1-α)
    return result
end

function fracdiff(f::Union{Number, Function}, α::Float64, start_point, end_point::AbstractArray, ::GL_Direct)::Vector
    ResultArray = Float64[]
    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracdiff(f, α, start_point, value, GL_Direct()))
    end
    return ResultArray
end

"""
# Caputo sense fractional derivative with first derivative known.

    fracdiff(fd, α, start_point, end_point, Caputo_Direct_First_Diff_Known())

If the first order derivative of a function is already known, we can use this method to compute the fractional order derivative more precisely.

The inout function should be the first order derivative of a function

### Example

```julia-repl
julia> fracdiff(x->5*x^4, 0.5, 0, 2.5, Caputo_Direct_First_Diff_Known())
```
Return the semi-derivative of ``f(x)=x^5`` at ``x=2.5``.

Compared with **Caputo_Direct** method, this method don't need to specify step size, more precision are guaranteed.
"""
function fracdiff(fd::Function, α, start_point, end_point, ::Caputo_Direct_First_Diff_Known)

    g(τ) = (end_point-τ) .^ (-α) * fd(τ)
    result = quadgk(g, start_point, end_point) ./gamma(1-α)
    return result
end

function fracdiff(f::Function, α, start_point, end_point::AbstractArray, step_size, ::Caputo_Direct_First_Diff_Known)::Vector
    ResultArray = Float64[]
    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracdiff(f, α, start_point, value, step_size, Caputo_Direct_First_Diff_Known()))
    end
    return ResultArray
end

"""
# Caputo sense fractional derivative with first and second derivative known.

    fracdiff(fd1, fd2, α, start_point, end_point, Caputo_Direct_First_Second_Diff_Known)

If the first and second order derivative of a function is already known, we can use this method to compute the fractional order derivative more precisely.

### Example

```julia-repl
julia> fracdiff(x->5*x^4, x->20*x^3, 0.5, 0, 2.5, Caputo_Direct_First_Second_Diff_known())
```
Return the semi-derivative of ``f(x)=x^5`` at ``x=2.5``.

Compared with **Caputo_Direct** method, this method don't need to specify step size, more precision are guaranteed.
"""
function fracdiff(fd1::Function, fd2, α, start_point, end_point, ::Caputo_Direct_First_Second_Diff_Known)

    temp1 = fd1(start_point) .* (end_point-start_point) .^ (1-α) ./ gamma(2-α)
    g(τ) = fd2(τ) .* ((end_point-τ) .^ (1-α))
    temp2 = quadgk(g, start_point, end_point) ./ gamma(2-α)
    result = temp1 .+ temp2
    return result
end

function fracdiff(f::Function, α::Float64, start_point, end_point::AbstractArray, h, ::Caputo_Direct_First_Second_Diff_Known)::Vector
    ResultArray = Float64[]
    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracdiff(f, α, start_point, value, h, Caputo_Direct_First_Second_Diff_Known()))
    end
    return ResultArray
end

"""
# Caputo sense Piecewise algorithm

    fracdiff(f, α, end_point, h, Caputo_Piecewise())

Using the **piecewise algorithm** to obtain the fractional derivative at a specific point.

### Example

```julia-repl
julia> fracdiff(x->x^5, 0.5, 2.5, 0.001, Caputo_Piecewise())
```

Return the fractional derivative of ``f(x)=x^5`` at point ``x=2.5``.

"""
function fracdiff(f::Union{Function, Number}, α::Float64, end_point, h::Float64, ::Caputo_Piecewise)::Float64
    checks(f, α, 0, end_point)

    m = floor(α)+1

    summation = 0
    n = Int64(end_point/h)

    for i in range(0, n, step=1)
        summation +=W₅(i, n, m, α)*first_order(f, i*h, h)
    end

    result=summation*h^(m-α)/gamma(m-α+2)
    return result
end
function W₅(i, n, m, α)
    if i==0
        return n^(m-α)*(m-α+1-n)+(n-1)^(m-α+1)
    elseif i==n
        return 1
    else
        (n-i-1)^(m-α+1)+(n-i+1)^(m-α+1)-2*(n-i)^(m-α+1)
    end
end
function first_order(f, point, h)
    return (f(point+h)-f(point-h))/(2*h)
end

function fracdiff(f::Union{Number, Function}, α, end_point::AbstractArray, h, ::Caputo_Piecewise)::Vector
    ResultArray = Float64[]
    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracdiff(f, α, value, h, Caputo_Piecewise()))
    end
    return ResultArray
end

#TODO: Use the improved alg!! This algorithm is not accurate
#This algorithm is not good, still more to do
"""
# Grünwald Letnikov sense derivative approximation

    fracdiff(f, α, end_point, h, GL_Multiplicative_Additive())

Grünwald–Letnikov multiplication-addition-multiplication-addition··· method to approximate fractional derivative.

### Example

```julia-repl
julia> fracdiff(x->x, 0.5, 1, 0.007, GL_Multiplicative_Additive())
1.127403405642918
```

!!! danger "Inaccurate"
    The `GL_Multiplicative_Additive` method is not accruate, please use it at your own risk.
"""
function fracdiff(f::Union{Function, Number}, α, end_point, h, ::GL_Multiplicative_Additive)::Float64
    checks(f, α, 0, end_point)

    summation = 0
    n = Int64(end_point/h)

    for i in range(0, n-1, step=1)
        summation+=gamma(i-α)/gamma(i+1)*f(end_point-i*end_point/n)
    end
    result=summation*end_point^(-α)*n^α/gamma(-α)
    return result
end

function fracdiff(f::Union{Number, Function}, α::Float64, end_point::AbstractArray, h, ::GL_Multiplicative_Additive)::Vector
    ResultArray = Float64[]
    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracdiff(f, α, value, h, GL_Multiplicative_Additive()))
    end
    return ResultArray
end

#TODO: This algorithm is same with the above one, not accurate!!!
#This algorithm is not good, still more to do
"""
# Grünwald Letnikov sense derivative approximation

    fracdiff(f, α, end_point, h, GL_Lagrange_Three_Point_Interp())

Using Lagrange three poitns interpolation to approximate the fractional derivative.

### Example

```julia-repl
julia> fracdiff(x->x, 0.5, 1, 0.007, GL_Lagrange_Three_Point_Interp())
1.1283963376811044
```

!!! danger "Inaccurate"
    The `GL_Multiplicative_Additive` method is not accruate, please use it at your own risk.
"""
function fracdiff(f::Union{Function, Number}, α::Float64, end_point, h, ::GL_Lagrange_Three_Point_Interp)::Float64
    checks(f, α, 0, end_point)

    n = end_point/h
    summation=0

    for i in range(0, n-1, step=1)
        summation += gamma(i-α)/gamma(i+1)*(f(end_point-i*end_point/n)+1/4*α*(f(end_point-(i-1)*end_point/n)-f(end_point-(i+1)*end_point/n))+1/8*α^2*(f(end_point-(i-1)*end_point/n)-2*f(end_point-i*end_point/n)+f(end_point-(i+1)*end_point/n)))
    end

    result = summation*end_point^(-α)*n^α/gamma(-α)
    return result
end

function fracdiff(f::Union{Number, Function}, α::Float64, end_point::AbstractArray, h, ::GL_Lagrange_Three_Point_Interp)::Vector
    ResultArray = Float64[]
    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracdiff(f, α, value, h, GL_Lagrange_Three_Point_Interp()))
    end
    return ResultArray
end

"""
# Riemann Liouville sense derivative approximation

    fracdiff(f, α, end_point, h, RLDiff_Approx())

Using approximation to obtain fractional derivative value.

### Example

```julia-repl
julia> fracdiff(x->x^5, 0.5, 2.5, 0.0001, RLDiff_Approx())
```

!!! warning
    The RLDiff_Approx algorithm only support for 0 < α < 1.

"""
function fracdiff(f::Union{Number, Function}, α, end_point, h, ::RLDiff_Approx)::Float64
    checks(f, α, 0, end_point)

    summation = 0
    n = end_point/h

    @fastmath @inbounds @simd for i in range(0, n-1, step=1)
        summation+=(f(end_point-i*end_point/n)-f(end_point-(i+1)*end_point/n))*((i+1)^(1-α)-i^(1-α))
    end

    result=((1-α)*f(0)/n^α+summation)*end_point^(-α)*n^α/gamma(2-α)
    return result
end

function fracdiff(f::Union{Number, Function}, α, end_point::AbstractArray, h, ::RLDiff_Approx)::Vector
    ResultArray = Float64[]
    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracdiff(f, α, value, h, RLDiff_Approx()))
    end
    return ResultArray
end
#TODO: Add arbitrary order based on Riemann Liouville approximation
function fracdiff(f::Union{Number, Function}, α, end_point, h, ::RLDiff_Approx_2)::Float64

    #The fractional derivative of number is relating with the end_point.
    if typeof(f) <: Number
        if f == 0
            return 0
        else
            return f/sqrt(pi*end_point)
        end
    end

    if end_point == 0
        return 0
    end

    summation = 0
    n = end_point/h

    for i in range(0, n-1, step=1)
        summation+=(f(end_point-(i-1)*h) - 2*f(end_point-i*h) + f(end_point-(i+1)*h))*((i+1)^(2-α)-i^(2-α))
    end

    result=((1-α)*(2-α)*f(0)/n^α + (2-α)*(f(h)-f(0))/n^(α-1) + summation)*end_point^(-α)*n^α/gamma(3-α)
    return result
end

function fracdiff(f::Union{Number, Function}, α, end_point, h, ::RLDiff_Approx_3)::Float64

    #The fractional derivative of number is relating with the end_point.
    if typeof(f) <: Number
        if f == 0
            return 0
        else
            return f/sqrt(pi*end_point)
        end
    end

    if end_point == 0
        return 0
    end

    summation = 0
    n = end_point/h

    for i in range(0, n-1, step=1)
        summation+=(f(end_point-i*end_point/n)-f(end_point-(i+1)*end_point/n))*((i+1)^(1-α)-i^(1-α))
    end

    result=((1-α)*f(0)/n^α+summation)*end_point^(-α)*n^α/gamma(2-α)
    return result
end

function fracdiff(f, α, end_point, ::RLDiff_Approx_Any)
    leftsummation = 0
    rightsummation = 0
    diffgen = 0

    n=end_point/h

    #TODO:
    for i in range(0, n-1, step=1)
        rightsummation+=1
    end

end

"""
https://en.wikipedia.org/wiki/Numerical_differentiation#Higher_derivatives
"""
function diffgen(f, k, n, h, end_point)

    summation=0

    for i in range(0, k, step=1)
        summation+=(-1)^(k+n)*binomial(n, k)*f(end_point)
    end
end

"""
# Grünwald Letnikov sense derivative approximation

    fracdiff(f::Union{Function, Number}, α::AbstractArray, end_point, h, ::GL_Finite_Difference)::Vector

Use finite difference method to obtain Grünwald Letnikov sense fractional derivative.

### Example

```julia-repl
julia> fracdiff(x->x, 0.5, 1, 0.01, GL_Finite_Difference())
1.1269695801851276
```

!!! danger "Inaccurate"
    `GL_Finite_Difference` method is not accruate, please use it at your own risk.
"""
function fracdiff(f::Union{Number, Function}, α::Float64, end_point::Real, h, ::GL_Finite_Difference)::Float64
    checks(f, α, 0, end_point)

    n = end_point/h
    result=0

    @fastmath @simd for i in range(0, n, step=1)
        result+=(-1)^i*gamma(α+1)/(gamma(i+1)*gamma(α-i+1))*f(end_point-i*h)
    end

    result1=result*h^(-α)
    return result1
end

function fracdiff(f::Union{Function, Number}, α::AbstractArray, end_point, h, ::GL_Finite_Difference)::Vector
    ResultArray = Float64[]

    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracint(f, α, value, h, GL_Finite_Difference()))
    end
    return ResultArray
end

"""
# Caputo sense Diethelm computation

    fracdiff(f, α, end_point, h, Caputo_Diethelm())

Using quadrature weights(derived from product trapezoidal rule) to approximate the derivative.

### Example

```julia-repl
julia> fracdiff(x->x, 0.5, 1, 0.007, Caputo_Diethelm())
1.128378318687192
```

!!!info "Scope"
    0 < α < 1
"""
function fracdiff(f::Union{Function, Number}, α::Float64, end_point, h, ::Caputo_Diethelm)
    checks(f, α, 0, end_point)

    N = end_point/h

    result=0

    @fastmath @inbounds @simd for i in range(0, N, step=1)
        result+=quadweights(i, N, α)*(f(end_point-i*h)-f(0))
    end
    return result*h^(-α)/gamma(2-α)

end
function quadweights(n, N, α)
    if n == 0
        return 1
    elseif  0 < n < N
        return (n+1)^(1-α)-2*n^(1-α)+(n-1)^(1-α)
    elseif n == N
        return (1-α)*N^(-α)-N^(1-α)+(N-1)^(1-α)
    end
end
function nderivative(f, n, end_point, h)
    temp = 0

    @fastmath @inbounds @simd for k in range(0, n, step=1)
        temp += (-1)^(k+n)*binomial(n, k)*f(end_point+k*h)
    end
    return h^(-n)*temp
end

"""
# Riemann Liouville sense derivative using Triangular Strip Matrix to discrete and compute.

    fracdiff(f, α, end_point, h, RLDiff_Matrix())

Using Triangular Strip Matrix to approximate fractional derivative.

### Example

```julia-repl
julia> fracdiff(x->x^5, 0.5, 2.5, 0.0001, RLInt_Matrix())
```

!!! info
    Triangular Strip Matrix method returns the derivative in the interval ``[0, T]`` in ```Vector```

!!! tip
    With the advancing Triangular Strip Matrix method, you can not only compute fractional derivatives, integer order, higher order derivative is also supported!!
Try to set α as an integer, arbitrary integer of course! I promise you would enjoy it😏
"""
function fracdiff(f, α, end_point, h::Float64, ::RLDiff_Matrix)
    N=Int64(end_point/h+1)
    @views tspan=collect(0:h:end_point)
    return B(N, α, h)*f.(tspan)
end

#Compute the eliminator matrix Sₖ by omiting n-th row
function eliminator(n, row)
    temp = zeros(n, n)+I
    return @views temp[Not(row), :]
end

#Generate elements in Matrix.
function omega(n, p)
    omega = zeros(n+1)

    omega[1]=1
    @fastmath @inbounds @simd for i in range(1, n, step=1)
        omega[i+1]=(1-(p+1)/i)*omega[i]
    end

    return omega

end
function B(N, p, h::Float64)
    result=zeros(N, N)
    temp=omega(N, p)

    @inbounds @simd for i in range(1, N, step=1)
        @views result[i, 1:i]=reverse(temp[1:i])
    end

    return h^(-p)*result
end

## Macros for convenient computing.
"""
    @fracdiff(f, α, point)

Return the α-order derivative of **f** at specific point.

```julia-repl
julia> @fracdiff(x->x, 0.5, 1)
1.1283791670955188
```
"""
macro fracdiff(f::Union{Number, Function}, α::Float64, point)
    return :(fracdiff($f, $α, $point, 0.0001, RLDiff_Approx()))
end

"""
    @semifracdiff(f, point)

Return the semi-derivative of **f** at spedific point.

```julia-repl
julia> @semifracdiff(x->x, 1)
1.1283791670955188
```
"""
macro semifracdiff(f::Union{Number, Function}, point)
    return :(fracdiff($f, 0.5, $point, 0.0001, RLDiff_Approx()))
end


text:
# RRR
# ---
#
# Indexable and compressible bit vectors. The compression is based on the
# abandance of 1s (or 0s) in a block. If 1s or 0s are clustered in regions, RRR
# can achieve high compression.  This file implements two variants of the RRR:
# * RRR: RRR of Claude and Navarro (2008)
#   - The block size is 15 bits and the original block is decoded from a universal table.
# * LargeRRR: RRR of Navarro and Providel (2012)
#   - The block size is 63 bits and the original block is decoded on the fly.
#
# Data Layout
# bits:      ..|...............|...............|... ...|...............|..
# blocks:      |   block j+1   |   block j+2   |  ...  |  block j+ssr  |
# superblocks: |                       superblock                      |
#
# TODO: formatting citations
# Raman, R., Raman, V., & Satti, S. R. (2007).
# Succinct Indexable Dictionaries with Applications to Encoding k-ary Trees, Prefix Sums and Multisets,
# 3(4), 1–25. doi:10.1145/1290672.1290680
#
# Claude, F., & Navarro, G. (2008).
# Practical Rank / Select Queries over Arbitrary Sequences,
# 080019(i), 176–187.
#
# Navarro, G., & Providel, E. (2012). Fast , Small , Simple Rank / Select on Bitmaps, (1).

const superblock_sampling_rate = 16  # blocks per superblock

# sampling rate controls the tradeoff between performance and space
struct SuperBlock
    # partial sum of rank values
    rank::UInt64
    # offset of the starting position of r-index
    offset::UInt64
end

"""
Compressed indexable bit vector.

`RRR` compresses a bit vector in an information-theoretically optimal representation.
This compression is based on local abundance of bits; if 0s or 1s are clustered in a
bit vector, it will be compressed even when the numbers of 0s and 1s are globally equal.

Let `n` be the length of a `RRR`, the asymptotic query times are

* getindex: `O(1)`
* rank: `O(1)`
* select: `O(log n)`

See (Raman et al, 2007, doi:10.1145/1290672.1290680) for more details.
"""
struct RRR <: AbstractIndexableBitVector
    # class values for each block
    # class values are encoded in 4 bits, hence `ks` is a packed array of 4-bit elements
    ks::Vector{UInt8}
    # r-indices (or offset) for each block
    # r-index is a variable-length encoding of position in a class. This array is also packed without delimiters.
    rs::Vector{UInt16}
    # sampled blocks
    superblocks::Vector{SuperBlock}
    # length of bits
    len::Int
end

blocksizeof(::Type{RRR}) = 15
empty_rs(::Type{RRR}) = UInt16[]

Base.copy(rrr::RRR) = RRR(copy(rrr.ks), copy(rrr.rs), copy(rrr.superblocks), rrr.len)

RRR() = RRR(UInt8[], UInt8[], SuperBlock[], 0)
RRR(vec::AbstractVector{Bool}) = make_rrr(RRR, vec)
Base.convert(::Type{RRR}, vec::AbstractVector{Bool}) = RRR(vec)

Base.length(rrr::RRR) = rrr.len

# data size of payload in bytes
Base.sizeof(rrr::RRR) = sizeof(rrr.ks) + sizeof(rrr.rs) + sizeof(rrr.superblocks)

function Base.getindex(rrr::RRR, i::Integer)
    checkbounds(rrr, i)
    blocksize = blocksizeof(RRR)
    j = div(i - 1, blocksize) + 1
    k, r, _ = jthblock(rrr, j)
    bits = E[K[k+1]+r+1] << (16 - blocksize)
    return bitat(UInt16, bits, rem(i - 1, blocksize) + 1)
end

function rank1(rrr::RRR, i::Int)
    i = clamp(i, 0, rrr.len)
    if i == 0
        return 0
    end
    blocksize = blocksizeof(RRR)
    j, rem = divrem(i - 1, blocksize)
    k, r, rank = jthblock(rrr, j + 1)
    bits = E[K[k+1]+r+1] << (16 - blocksize)
    rank += count_ones(bits & lmask(UInt16, rem + 1))
    return convert(Int, rank)
end

rank1(rrr::RRR, i::Integer) = rank1(rrr, Int(i))

# return the class of j-th block
@inline function classof(rrr::RRR, j::Int)
    ki, rem = divrem(j - 1, 2)
    @inbounds k = rrr.ks[ki+1]
    if rem == 0
        k >>= 4
    else
        k &= 0x0f
    end
    return k
end

# RRR of Navarro and Providel
# the block size is 63 and the block is decoded on the fly.
struct LargeRRR <: AbstractIndexableBitVector
    ks::Vector{UInt8}
    rs::Vector{UInt64}
    superblocks::Vector{SuperBlock}
    len::Int
end

blocksizeof(::Type{LargeRRR}) = 63
empty_rs(::Type{LargeRRR}) = UInt64[]

Base.copy(rrr::LargeRRR) = LargeRRR(
    copy(rrr.ks), copy(rrr.rs), copy(rrr.superblocks), rrr.len)

# 3 elements of ks store 4 classes
# ks:      |........|........|........|  8bits for each
# classes: |......       .... ..      |  6bits for each
#          |      .. ....       ......|

LargeRRR() = LargeRRR(UInt8[], UInt64[], SuperBlock[], 0)
Base.convert(::Type{LargeRRR}, vec::Union{BitVector,Vector{Bool}}) = make_rrr(LargeRRR, vec)

Base.length(rrr::LargeRRR) = rrr.len

function Base.getindex(rrr::LargeRRR, i::Integer)
    checkbounds(rrr, i)
    blocksize = blocksizeof(LargeRRR)
    j = div(i - 1, blocksize) + 1
    k, r, _ = jthblock(rrr, j)
    bits = rindex2bits(r, blocksize, convert(Int, k)) << (64 - blocksize)
    return bitat(UInt64, bits, rem(i - 1, blocksize) + 1)
end

function rank1(rrr::LargeRRR, i::Integer)
    i = clamp(i, 0, rrr.len)
    if i == 0
        return 0
    end
    blocksize = blocksizeof(LargeRRR)
    j, rem = divrem(i - 1, blocksize)
    k, r, rank = jthblock(rrr, j + 1)
    bits = rindex2bits(r, blocksize, convert(Int, k)) << (64 - blocksize)
    rank += count_ones(bits & lmask(UInt64, rem + 1))
    return convert(Int, rank)
end

function classof(rrr::LargeRRR, j::Int)
    ki, rem = divrem(j - 1, 4)
    if rem == 0
        k = rrr.ks[3ki+1] >> 2
    elseif rem == 1
        k1 = (rmask(UInt8, 2) & rrr.ks[3ki+1]) << 4
        k2 = rrr.ks[3ki+2] >> 4
        k = convert(UInt8, k1 + k2)
    elseif rem == 2
        k1 = (rmask(UInt8, 4) & rrr.ks[3ki+2]) << 2
        k2 = rrr.ks[3ki+3] >> 6
        k = convert(UInt8, k1 + k2)
    else  # rem == 3
        k = rrr.ks[3ki+3] & rmask(UInt8, 6)
    end
    return k
end

function make_rrr(::Type{T}, src::AbstractVector{Bool}) where {T<:Union{RRR,LargeRRR}}
    len = length(src)
    if len > typemax(Int)
        error("the bit vector is too large")
    end
    blocksize = blocksizeof(T)
    n_blocks = div(len - 1, blocksize) + 1
    ks = UInt8[]
    rs = empty_rs(T)
    superblocks = SuperBlock[]
    rs_el_bits = bitsof(eltype(rs))
    rank = 0
    n_rembits = 0  # the number of remaining bits in rs
    for i in 1:n_blocks
        # sample superblock
        if i % superblock_sampling_rate == 1
            push!(superblocks, SuperBlock(rank, sizeof(rs) * 8 - n_rembits + 1))
        end

        # bits and class
        bits = read_bits(src, (i - 1) * blocksize + 1, blocksize)
        k = convert(eltype(ks), count_ones(bits))
        @assert 0 ≤ k ≤ blocksize
        rank += k

        # store the class to ks
        if T === RRR
            if isodd(i)
                # use upper 4 bits: here----
                push!(ks, k << 4)
            else
                # use lower 4 bits: ----here
                ks[end] |= k
            end
        elseif T === LargeRRR
            rem = i % 4
            if rem == 1
                push!(ks, 0, 0, 0)
                ks[end-2] |= k << 2
            elseif rem == 2
                ks[end-2] |= k >> 4
                ks[end-1] |= k << 4
            elseif rem == 3
                ks[end-1] |= k >> 2
                ks[end]   |= k << 6
            else  # rem == 4
                ks[end]   |= k
            end
        else
            error()
        end

        # store the r-index to rs
        r = convert(eltype(rs), bits2rindex(bits, blocksize, k))
        @assert 0 ≤ r < Comb[blocksize,k]
        if isempty(rs)
            push!(rs, 0)
            n_rembits += rs_el_bits
        end
        # the number of bits required to store r-index
        n_rbits = nbits(blocksize, k)
        @assert n_rbits ≤ rs_el_bits
        if n_rembits ≥ n_rbits
            # this r-index can be stored in the last element
            rs[end] |= r << (n_rembits - n_rbits)
        else
            # this r-index spans the two last elements
            rs[end] |= r >> (n_rbits - n_rembits)
            push!(rs, 0)
            rs[end] |= r << (rs_el_bits - (n_rbits - n_rembits))
            n_rembits += rs_el_bits
        end
        n_rembits -= n_rbits
    end
    return T(ks, rs, superblocks, len)
end

# Compute the class and r-index of the j-th block, and
# the rank value at the beginning of the j-th block
function jthblock(rrr::Union{RRR,LargeRRR}, j::Int)
    @assert 1 ≤ j
    i = div(j - 1, superblock_sampling_rate)
    superblock = rrr.superblocks[i+1]
    blocksize = blocksizeof(typeof(rrr))
    # read blocks just before the j-th block in the superblock
    rank = convert(Int, superblock.rank)
    offset = convert(Int, superblock.offset)
    lo = i * superblock_sampling_rate + 1
    hi = j - 1
    for j′ in lo:hi
        k = classof(rrr, j′)
        offset += nbits(blocksize, k)
        rank += k
    end
    # read the j-th block
    k = classof(rrr, j)
    r = read_rindex(rrr.rs, offset, nbits(blocksize, k))
    return k, r, rank
end

# return left-aligned bits
function read_bits(src, from, len)
    @assert 1 ≤ len ≤ 64
    bits = zero(UInt64)
    to = from + len - 1
    for i in from:min(to, lastindex(src))
        bits <<= 1
        if src[i]
            bits |= 1
        end
    end
    if to > lastindex(src)
        # align bits to the left
        bits <<= to - lastindex(src)
    end
    return bits
end

# offset is 1-based index
function read_rindex(rs::Vector{T}, offset::Int, len::Int) where {T<:Unsigned}
    w = bitsof(T)
    ri, rem = divrem(offset - 1, w)
    if w - rem ≥ len
        # stored in an element (T = UInt16, w = 16)
        # |      ri+1      |      ri+2      |
        # |................|................|
        # |    xxxxxx      |                |
        #----->|offset
        # |<->|rem
        #      |<-->|len
        @inbounds r = rs[ri+1] & rmask(T, w - rem)
        r >>= w - (rem + len)
    else
        # spans two elements (T = UInt16, w = 16)
        # |      ri+1      |      ri+2      |
        # |................|................|
        # |          xxxxxx|xxxx            |
        #----------->|offset
        # |<------->|rem
        #            |<---- -->|len
        # |<--------------- -->|rem+len
        @inbounds r1 = rs[ri+1] & rmask(T, w - rem)
        r1 <<= rem + len - w
        @inbounds r2 = rs[ri+2] & lmask(T, rem + len - w)
        r2 >>= 2w - (rem + len)
        r = convert(T, r1 + r2)
    end
    return convert(Int, r)
end

# Bit Encoding
# ------------

# Example: t = 4, k = 2
#   r | bits
#   - | ----
#   0 | 0011
#   1 | 0101
#   2 | 0110
#   3 | 1001
#   4 | 1010
#   5 | 1100

# Inverse transformation of bits2rindex (encoding).
# `bits` should be filled from the least significant bit.
# For example, 0b00...001 is the first bits of the class k=1.
function bits2rindex(bits::UInt64, t::Int, k::Integer)
    @assert count_ones(bits) == k
    @assert 0 ≤ k ≤ t ≤ 64
    r = 0
    i = 65 - t
    while k > 0
        if bitat(bits, 65 - t)
            r += Comb[t-1,k]
            k -= 1
        end
        i += 1
        t -= 1
    end
    return r
end

# Inverse transformation of bits2rindex (decoding).
function rindex2bits(r::Int, t::Int, k::Integer)
    @assert 0 ≤ r < Comb[t,k]
    @assert 0 ≤ k ≤ t ≤ 64
    bits = zero(UInt64)
    @inbounds while r > 0
        c = Comb[t-1,k]
        if r ≥ c
            # the first bit is 1
            bits |= 1 << (t - 1)
            r -= c
            k -= 1
        end
        t -= 1
    end
    bits |= typemax(UInt64) >>> (64 - k)
    return bits
end

# Look-up Tables
# --------------

struct CombinationTable
    table::Matrix{Int}
end

@inline function Base.getindex(comb::CombinationTable, t::Int, k::Integer)
    @inbounds c = comb.table[t+1,k+1]
    return c
end

@assert blocksizeof(RRR) ≤ blocksizeof(LargeRRR)

const Comb = CombinationTable([binomial(t, k) for t in 0:blocksizeof(LargeRRR), k in 0:blocksizeof(LargeRRR)])

# enumeration of bit patterns for blocks, sorted by class and r-index
const E, K = let
    t = blocksizeof(RRR) + 1
    bitss = UInt16[]
    offsets = Int[]
    for k in 0:t
        push!(offsets, length(bitss))
        for r in 0:Comb[t,k]-1
            bits = rindex2bits(r, t, k)
            push!(bitss, bits)
        end
    end
    bitss, offsets
end

# Lookup table to know the number of bits to encode class k's r-indices with blocksize t
const NBitsTable = [t ≥ k ? ceil(Int, log2(Comb[t,k])) : 0 for t in 1:blocksizeof(LargeRRR), k in 0:blocksizeof(LargeRRR)]

function nbits(t, k)
    return NBitsTable[t,k+1]
end