text:
include(joinpath(AS.topdirectory, "test", "EpemaData.jl"))

wing = EpemaData.wing
rotor = EpemaData.rotor
setup = EpemaData.setup
results = EpemaData.results

Vinf = setup[:"Vinf"]
vinfs = [Vinf]
plotstepi = [1]
alphas = [setup[:"wing_aoa"]]
wing_b = wing[:"span"]
wing_TR = wing["chord_tip"] / wing[:"chord_root"]
wing_AR = wing["span"]^2 / wing[:"area"]
wing_θroot = 0.0
wing_θtip = 0.0
wing_le_sweep = 0.0
wing_ϕ = 0.0
rotor_omegas = [AS.get_omega(Vinf=setup[:"Vinf"], J=setup[:"J"], D=rotor[:"diameter"])]
nblades = [rotor[:"Nblades"]]
rhub = [rotor[:"radius_hub"] / 2]
rtip = [rotor[:"radius"]]
chord_data = rotor["r/R vs chord/R"]
twist_data = rotor["r/R vs twist"]
q = 1/2 * 1.0 * Vinf^2

radii = [chord_data[:,1] * rtip[1]]
chords = [chord_data[:,2] * rtip[1]]
twists = [FM.linear(twist_data[:,1]*rtip[1], twist_data[:,2], radii[1])]
airfoilcontour = joinpath(AS.topdirectory, "data", "airfoil", "contours", "e212-il.csv")
airfoilcontours = [fill(airfoilcontour, length(radii[1]))]
airfoilname = "eppler212"
airfoilnames = [fill(airfoilname, length(radii[1]))]
index = [1]
rotor_positions = [[0.0; rotor[:"y"]; 0.0]]
rotor_orientations = [[-1.0; 0.0; 0.0]] # positive x downstream
spin_directions = [true]
polar_directory=joinpath(AS.topdirectory, "data","airfoil","polars","20210618")

args = AS.lift_distribution_template(vinfs, plotstepi, alphas, wing_b, wing_TR, wing_AR, wing_θroot, wing_θtip, wing_le_sweep, wing_ϕ,
                        Res_list = [fill([EpemaData.setup[:"Re"]], length(radii))],
                        surfacenames = ["epema wing"], Vinf=Vinf)

outs = AS.runsimulation!(args...)

aircraft = args[1]
parameters = args[2]

span_plot = aircraft.wingsystem.lifting_line_rs[1][2,2:end] / (EpemaData.wing[:"span"]/2)
cls_plot = parameters.cfs[1][3,:]
wing_chords = FM.linear(epema_chords[:,1], epema_chords[:,2], range(epema_chords[1,1], stop=epema_chords[3,1], length=length(span_plot)))
cls_plot = cls_plot .* wing_chords / EpemaData.wing[:"mac"]

epema_chords = EpemaData.wing[:"chords"]
epema_cl = EpemaData.results[:"lift distribution VLM rotors off"][:,2]
epema_x = EpemaData.results[:"lift distribution VLM rotors off"][:,1]
epema_cl_interpolated = FM.linear(epema_x, epema_cl, span_plot)
epema_cl_interpolated = [[0 0 cl] for cl in epema_cl_interpolated]

gammas = [panel.gamma for panel in aircraft.wingsystem.system.properties[1]]
velocities = [panel.velocity for panel in aircraft.wingsystem.system.properties[1]]

cls_new = 2 * gammas' ./ wing_chords .* wing_chords / EpemaData.wing[:"mac"]

panels = aircraft.wingsystem.system.surfaces[1]
Δs = VL.top_vector.(panels)
Vfreestream = VL.freestream_velocity(aircraft.wingsystem.system.freestream[1])

# calculating gammas from epema data
# Vs = [Vfreestream + vel for vel in velocities]

# gammas_epema = zeros(length(wing_chords))
# for i in 1:length(wing_chords)
#     test_gamma = epema_cl_interpolated[i]' * Vinf^2 .* wing_chords[i] ./ (2 * cross(Vs[i], Δs[i]))[3]
#     gammas_epema[i] = test_gamma[3]
# end

figure()
plot(span_plot, cls_plot, label="BEM+VLM")
plot(EpemaData.results[:"lift distribution VLM rotors off"][:,1], epema_cl, label="Epema VLM")
plot(span_plot, cls_new, label=L"2 * Γ / V_{\infty}")
xlabel(L"2y/b")
ylabel(L"normalized\ c_l (c/c_{mac})")
legend()

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
using AstroBase
using Distributions
using JPLEphemeris
using Random

Random.seed!(430)

function runner_ep(times)
    pos = zeros(3)
    vel = zeros(3)
    state!(pos, vel, spk, times[1], luna, mercury)
    total = 0.0
    for t in times
        time = @elapsed state!(pos, vel, spk, t, luna, mercury)
        total += time
    end
    print_time(total / n)
end

function print_time(time)
    if time < 1e-6
        println("Time: $(time*1e9) ns")
    elseif time < 1e-3
        println("Time: $(time*1e6) μs")
    elseif time < 1.0
        println("Time: $(time*1e3) ms")
    else
        println("Time: $time s")
    end
end

kernel = joinpath(dirname(pathof(JPLEphemeris)), "..",  "deps", "de430.bsp")
spk = SPK(joinpath(kernel))
seg = spk.segments[0][3]
n = 1_000_000

first = seg.firstdate
last = seg.lastdate
linear = collect(range(first, stop=last, length=n))
d = Truncated(Normal(), first, last)
random = rand(d, n)

linear_ep = TDBEpoch.(linear, origin=:julian)
random_ep = TDBEpoch.(random, origin=:julian)

runner_ep(linear_ep)
runner_ep(random_ep)

ext: jl


text:
include("test_2d.jl")
include("test_manual.jl")


text:
## low-level API

const nvtxDomainHandle_t = Ptr{Cvoid}

struct Domain
    handle::nvtxDomainHandle_t

    function Domain(name::String)
        info("Creating domain $name")
        handle = ccall((:nvtxDomainCreateA, libnvtx), nvtxDomainHandle_t, (Cstring,), name)
        new(handle)
    end
end

Base.unsafe_convert(::Type{nvtxDomainHandle_t}, dom::Domain) = dom.handle

unsafe_destroy!(dom::Domain) =
    ccall((:nvtxDomainDestroy, libnvtx), Nothing, (nvtxDomainHandle_t,), dom)

## high-level API

function domain(f::Function, name::String)
    dom = Domain(name)
    f(dom)
    unsafe_destroy!(dom)
end

# TODO: support for ex mark/range, with domain argument


text:
# simple methods to process audio signals

import DSP.resample, DSP.arraysplit
export mono, resample, duration, play, pitchshift, speedup, slowdown, zero_crossing_rate

"""
    mono(audio)

convert a multichannel audio to mono
"""
function mono(audio::SampleBuf{T, 2}) where {T <: AbstractFloat}
    SampleBuf{T, 1}(
        vec(monoch(audio).data),
        audio.samplerate
    )
end

# special code for fixed-point samples, to avoid overflow
function mono(audio::SampleBuf{T, 2}) where {T <: Fixed}
    nchannels = nchannels(audio)
    if nchannels == 1
        SampleBuf{T, 1}(vec(audio.data), audio.samplerate)
    elseif nchannels == 2
        nsamples = nframes(audio)
        buffer = Array{T}(undef, nsamples)
        for i = 1:nsamples
            @inbounds a = audio.data[i, 1].i
            @inbounds b = audio.data[i, 2].i
            m = (a >> 1) + (b >> 1) + (a & b & 1)
            @inbounds buffer[i] = T(m, 0)
        end
        SampleBuf{T, 1}(buffer, audio.samplerate)
    else
        SampleBuf{T, 1}(
            map(T, vec(mean(map(Float32, audio.data)))),
            audio.samplerate
        )
    end
end

"""resample audio with a different sample rate"""
function resample(audio::SampleBuf{T, 2}, samplerate::Real) where T
    rate = samplerate / audio.samplerate
    filter = DSP.resample_filter(rate, 512, 1.0, 140)
    SampleBuf{T, 2}(
        mapslices(audio.data, dims=1) do data
            DSP.resample(data, rate, filter)
        end,
        samplerate
    )
end

"""resample audio with a different sample rate"""
function resample(audio::SampleBuf{T, 1}, samplerate::Real) where {T, F}
    SampleBuf{T, 1}(
        DSP.resample(audio.data, samplerate / audio.samplerate),
        samplerate
    )
end

"""returns the duration of given audio, in seconds"""
function duration(audio::SampleBuf)
    nframes(audio) / samplerate(audio)
end

"""
    play(audio)

play the audio on local computer using PortAudio
"""
function play(audio::SampleBuf{Float32})
    # import PortAudio on-demand
    @eval import PortAudio
    nchannels = nchannels(audio)
    stream = PortAudio.PortAudioStream(2, nchannels)
    try
        write(stream, audio)
    finally
        close(stream)
    end
end
play(audio::SampleBuf{T}) where T = play(map(Float32, audio))

""""""
function pitchshift(audio::SampleBuf{T, N}, semitones::Real = 8) where {T, N}
    rate = 2.0 ^ (semitones / 12.0)
    shifted = resample(slowdown(audio, rate), audio.samplerate / rate)
    SampleBuf{T, N}(
        shifted.data,
        audio.samplerate
    )
end

""""""
function speedup(audio::SampleBuf, speed::Real, windowsize::Int = 1024, hopsize::Int = windowsize >> 2; kwargs...)
    S = stft(audio, windowsize, hopsize; kwargs...)
    S = phase_vocoder(S, speed, hopsize)
    istft(S, audio.samplerate, windowsize, hopsize; kwargs...)
end

""""""
function slowdown(audio::SampleBuf, ratio::Real, windowsize::Int = 1024, hopsize::Int = windowsize >> 2; kwargs...)
    speedup(audio, 1.0 / ratio, windowsize, hopsize; kwargs...)
end

""""""
function zero_crossing_rate(audio::SampleBuf{T, 1}, framesize::Int = 1024, hopsize::Int = framesize >> 2) where T
    nframes = MusicProcessing.nframes(length(audio.data), framesize, hopsize)
    result = Array{Float32}(undef,nframes)

    offset = 0
    for i = 1:nframes
        result[i] = zero_crossings(audio.data, framesize, offset) / framesize
        offset += hopsize
    end

    result
end


text:
module BuildTools

export install, build

using Logging
using PackageCompiler
using Pkg.TOML
using Pkg.PlatformEngines
using ..Comonicon
using ..Comonicon.Parse
using ..Comonicon.CodeGen
using ..Comonicon.PATH
using ..Comonicon.Types
using ..Comonicon.Tools

const COMONICON_URL = "https://github.com/Roger-luo/Comonicon.jl"
# Comonicon.toml
# name = "ion"

# [install]
# bin = "~/.julia/bin"
# completion=true
# quiet=false
# compile="min"
# optimize=2

# [sysimg]
# path="deps/lib"
# incremental=false
# filter_stdlibs=true
# cpu_target="native"

# [download]
# host="github.com"
# user="Roger-luo"

const DEFAULT_INSTALL_CONFIG = Dict(
    "bin" => PATH.default_julia_bin(),
    "completion" => true,
    "quiet" => false,
    "compile" => nothing,
    "optimize" => 2,
)

const DEFAULT_SYSIMG_CONFIG = Dict(
    "path" => "deps/lib",
    "incremental" => true,
    "filter_stdlibs" => false,
    "cpu_target" => "native",
)

const COMONICON_TOML = ["Comonicon.toml", "JuliaComonicon.toml"]

function build(mod, sysimg = true; incremental = true, filter_stdlibs = false, kwargs...)
    configs = read_configs(mod; incremental = incremental, filter_stdlibs = filter_stdlibs, kwargs...)

    validate_toml(mod, configs)
    configs = merge_defaults(mod, configs)

    if sysimg && !haskey(configs, "sysimg")
        configs["sysimg"] = DEFAULT_SYSIMG_CONFIG
    elseif !sysimg
        delete!(configs, "sysimg")
    end
    # do not download in manual mode
    delete!(configs, "download")
    return install(mod, configs)
end

function install(mod; kwargs...)
    configs = read_configs(mod; kwargs...)
    configs = merge_defaults(mod, configs)
    validate_toml(mod, configs)
    return install(mod, configs)
end

function read_configs(mod; kwargs...)
    configs = read_toml(mod)
    for (k, v) in kwargs
        if k == :name
            configs["name"] = v
        end

        if k in [:bin, :completion, :quiet, :compile, :optimize]
            install_configs = get!(configs, "install", Dict{String,Any}())
            install_configs[string(k)] = v
        end

        if k in [:path, :incremental, :filter_stdlibs, :cpu_target]
            sysimg_configs = get!(configs, "sysimg", Dict{String,Any}())
            sysimg_configs[string(k)] = v
        end

        if k in [:host, :repo, :user]
            download_config = get!(configs, "download", Dict{String,Any}())
            download_config[string(k)] = v
        end
    end
    return configs
end

function read_toml(mod)
    path = nothing
    for file in COMONICON_TOML
        _path = PATH.project(mod, file)
        if ispath(_path)
            path = _path
            break
        end
    end

    path === nothing && return Dict()
    configs = TOML.parsefile(path)
    return configs
end

function merge_defaults(mod, configs)
    if haskey(configs, "sysimg")
        configs["sysimg"] = merge(DEFAULT_SYSIMG_CONFIG, configs["sysimg"])
    end

    configs["install"] = merge(DEFAULT_INSTALL_CONFIG, configs["install"])
    if !haskey(configs, "name")
        configs["name"] = Parse.default_name(mod)
    end
    return configs
end

function validate_toml(mod, configs)
    haskey(configs, "name") || error("missing key \"name\" in Comonicon.toml or kwargs")
    got = configs["name"]
    exp = Types.cmd_name(mod.CASTED_COMMANDS["main"])
    got == exp || error("field name must be the same with entry name, expect $exp got $got")

    _check(configs["install"], "compile") do x
        x in [nothing, "min", "no", "all", "yes"]
    end

    for key in ["completion", "quiet"]
        _check(configs["install"], key) do x
            x isa Bool
        end
    end

    _check(configs["install"], "optimize") do x
        x isa Int
    end

    haskey(configs, "sysimg") || return

    for key in ["incremental", "filter_stdlibs"]
        _check(configs["sysimg"], key) do x
            x isa Bool
        end
    end

    haskey(configs, "download") || return

    for key in ["host", "user", "repo"]
        _check(configs["download"], key) do x
            x isa String
        end
    end

    return
end

function _check(f, configs, key)
    haskey(configs, key) || error("missing key $key in Comonicon.toml or kwargs")
    got = configs[key]
    f(got) || error("invalid value $got for field \"$key\" in Comonicon.toml or kwargs")
    return
end

function install(mod::Module, configs::Dict)
    if configs["install"]["quiet"]
        logger = NullLogger()
    else
        logger = ConsoleLogger()
    end

    with_logger(logger) do
        if haskey(configs, "sysimg")
            install_sysimg(mod, configs)
        end

        # if the system image is required by the developer,
        # when system image installation
        # errors, the CLI will not be installed.
        # User will need to use mod.build([sysimg=false]) to install it
        # manually, with an option to install without system image
        # or build it locally.

        # do not install script while building tarball
        if !("sysimg" in ARGS)
            install_script(mod, configs)
        end
    end
    return
end

function install_script(mod::Module, configs::Dict)
    if haskey(configs, "sysimg") && isfile(sysimage_path(mod, configs))
        sysimg_path = sysimage_path(mod, configs)
    else
        sysimg_path = nothing
    end

    install_configs = configs["install"]
    name = configs["name"]
    bin = install_configs["bin"]

    shadow = joinpath(bin, name * ".jl")
    if install_configs["compile"] === nothing
        compile = nothing
    else
        compile = Symbol(install_configs["compile"])
    end

    shell_script = cmd_script(
        mod,
        shadow;
        sysimg = sysimg_path,
        compile = compile,
        optimize = install_configs["optimize"],
    )

    file = joinpath(bin, name)

    if !ispath(bin)
        @info "cannot find Julia bin folder creating .julia/bin"
        mkpath(bin)
    end

    # generate contents
    @info "generating $shadow"
    open(shadow, "w+") do f
        println(f, "#= generated by Comonicon for $name =# using $mod; exit($mod.command_main())")
    end

    @info "generating $file"
    open(file, "w+") do f
        println(f, shell_script)
    end

    if install_configs["completion"]
        install_completion(mod, joinpath(dirname(bin), "completions"))
    end

    chmod(file, 0o777)
    return
end

function install_sysimg(mod::Module, configs::Dict)
    # if sysimg will be downloaded on user side
    # when we build the sysimg, a tarball should
    # be generated.
    if haskey(configs, "download")
        # we create a system image tarball
        # via an argument sysimg
        if "sysimg" in ARGS
            build_sysimg(mod, configs)
            create_tarball(mod, configs["name"])
        else
            download_sysimg(mod, configs)
        end
    else # manually triggered
        build_sysimg(mod, configs)
        if "tarball" in ARGS
            create_tarball(mod, configs["name"])
        end
    end
    return
end

function sysimage_path(mod, configs)
    return PATH.project(mod, configs["sysimg"]["path"], PATH.sysimg(configs["name"]))
end

function download_sysimg(mod::Module, configs::Dict)
    sysimg_configs = configs["sysimg"]
    name = configs["name"]
    os = osname()
    tarball_name = "$name-$(VERSION)-$os-$(Sys.ARCH).tar.gz"

    url = sysimg_url(mod, configs)
    tarball = joinpath(PATH.project(mod, "deps", tarball_name))
    PlatformEngines.probe_platform_engines!()

    try
        download(url, tarball)
        unpack(tarball, PATH.project(mod, "deps"))
    catch e
        @warn "fail to download $url, building the system image locally"
        # force incremental build
        sysimg_configs["incremental"] = true
        sysimg_configs["filter_stdlibs"] = false
        build_sysimg(mod, configs)
    end

    if ispath(tarball)
        rm(tarball)
    end
    return
end

function sysimg_url(mod, configs)
    name = configs["name"]
    host = configs["download"]["host"]
    if host == "github.com"
        url =
            "https://github.com/" *
            configs["download"]["user"] *
            "/" *
            configs["download"]["repo"] *
            "/releases/download/"
    else
        error("host $host is not supported, please open an issue at $COMONICON_URL")
    end

    tarball = tarball_name(name)
    url *= "v$(Comonicon.get_version(mod))/$tarball"
    return url
end

function build_sysimg(mod::Module, configs::Dict)
    sysimg_configs = configs["sysimg"]
    lib_path = PATH.project(mod, sysimg_configs["path"])
    if !ispath(lib_path)
        @info "creating library path: $lib_path"
        mkpath(lib_path)
    end

    # install precompile script
    precompile_jl = PATH.project(mod, "deps", "precompile.jl")
    @info "generating precompile execution file: $precompile_jl"
    open(precompile_jl, "w+") do f
        print(f, precompile_script(mod))
    end

    project = PATH.project(mod)
    incremental = sysimg_configs["incremental"]
    filter_stdlibs = sysimg_configs["filter_stdlibs"]
    cpu_target = sysimg_configs["cpu_target"]
    image_path = sysimage_path(mod, configs)

    @info "compile under project: $project"
    @info "incremental: $incremental"
    @info "filter stdlibs: $filter_stdlibs"
    @info "system image path: $image_path"

    create_sysimage(
        nameof(mod);
        sysimage_path = image_path,
        incremental = incremental,
        project = project,
        precompile_execution_file = [precompile_jl, PATH.project(mod, "test", "runtests.jl")],
        cpu_target = cpu_target,
        filter_stdlibs = filter_stdlibs,
    )

    return
end

function find_project(path)
    file = joinpath(path, "JuliaProject.toml")
    isfile(file) && return file
    file = joinpath(path, "Project.toml")
    isfile(file) && return file
    return
end

function create_tarball(mod::Module, name)
    version = get_version(mod)
    tarball = tarball_name(name)
    @info "creating tarball $tarball"
    cd(PATH.project(mod, "deps")) do
        run(`tar -czvf $tarball lib`)
    end
    return
end

function tarball_name(name)
    return "$name-$VERSION-$(osname())-$(Sys.ARCH).tar.gz"
end

"""
    osname()

Return the name of OS, will be used in building tarball.
"""
function osname()
    return Sys.isapple() ? "darwin" :
           Sys.islinux() ? "linux" :
           error("unsupported OS, please open an issue to request support at $COMONICON_URL")
end

"""
    cmd_script(mod, shadow; kwargs...)

Generates a shell script that can be use as the entry of
`mod.command_main`.

# Arguments

- `mod`: a module that contains the commands and the entry.
- `shadow`: location of a Julia script that calls the actual `mod.command_main`.

# Keywords

- `exename`: The julia executable name, default is [`PATH.default_exename`](@ref).
- `sysimg`: System image to use, default is `nothing`.
- `project`: the project path of the CLI.
- `compile`: julia compile level, can be [:yes, :no, :all, :min]
- `optimize`: julia optimization level, default is 2.
"""
function cmd_script(
    mod::Module,
    shadow::String;
    project::String = PATH.project(mod),
    exename::String = PATH.default_exename(),
    sysimg = nothing,
    compile = nothing,
    optimize = 2,
)

    head = "#!/bin/sh\n"
    if (project !== nothing) && ispath(project)
        head *= "JULIA_PROJECT=$project "
    end
    head *= exename
    script = String[head]

    if sysimg !== nothing
        push!(script, "-J$sysimg")
    end

    if compile in [:yes, :no, :all, :min]
        push!(script, "--compile=$compile")
    end

    push!(script, "-O$optimize")
    push!(script, "--startup-file=no")
    push!(script, "-- $shadow \$@")

    return join(script, " \\\n    ")
end

"""
    precompile_script(mod)

Generates a script to execute as `precompile_execution_file` for all the commands.
"""
function precompile_script(mod::Module)
    script = "using $mod;\n$mod.command_main([\"-h\"]);\n"

    if isdefined(mod, :CASTED_COMMANDS)
        for (name, cmd) in mod.CASTED_COMMANDS
            if name != "main" # skip main command
                script *= "$mod.command_main([$(precompile_script(mod, cmd))]);\n"
            end
        end
    end
    return script
end

function precompile_script(mod::Module, cmd::EntryCommand)
    return precompile_script(mod, cmd.root)
end

function precompile_script(mod::Module, cmd::LeafCommand)
    return "\"$(cmd_name(cmd))\", \"-h\""
end

function precompile_script(mod::Module, cmd::NodeCommand)
    return join(map(x -> "\"$(cmd_name(cmd))\", " * precompile_script(mod, x), cmd.subcmds))
end

Base.write(x::EntryCommand) = write(cachefile(), x)

"""
    write([io], cmd::EntryCommand)

Write the generated CLI script into a Julia script file. Default is the [`cachefile`](@ref).
"""
function Base.write(io::IO, x::EntryCommand)
    println(io, "#= generated by Comonicon =#")
    println(io, prettify(codegen(x)))
    println(io, "command_main()")
end

"""
    detect_shell()

Detect shell type via `SHELL` environment variable.
"""
function detect_shell()
    haskey(ENV, "SHELL") || error("cannot find available shell command")
    return basename(ENV["SHELL"])
end

"""
    install_completion(m::Module[, path::String=PATH.default_julia_fpath()])

Install completion script at `path`. Default path is [`PATH.default_julia_fpath()`](@ref).
"""
function install_completion(m::Module, path::String = PATH.default_julia_fpath())
    isdefined(m, :CASTED_COMMANDS) || error("cannot find Comonicon CLI entry")
    haskey(m.CASTED_COMMANDS, "main") || error("cannot find Comonicon CLI entry")

    main = m.CASTED_COMMANDS["main"]
    shell = detect_shell()
    shell === nothing && return

    if shell == "zsh"
        ctx = CodeGen.ZSHCompletionCtx()
    else
        @warn "$shell completion is not supported yet"
        return
    end

    script = CodeGen.codegen(ctx, main)

    if !ispath(path)
        mkpath(path)
    end

    write(joinpath(path, "_" * cmd_name(main)), script)
    return
end

function contain_comonicon_path(rcfile, env = ENV)
    if !haskey(env, "PATH")
        _contain_path(rcfile) && return true
        return false
    end

    for each in split(env["PATH"], ":")
        each == PATH.default_julia_bin() && return true
    end
    return false
end

function contain_comonicon_fpath(rcfile, env = ENV)
    if !haskey(env, "FPATH")
        _contain_fpath(rcfile) && return true
        return false
    end

    for each in split(env["FPATH"], ":")
        each == PATH.default_julia_fpath() && return true
    end
    return false
end

function _contain_path(rcfile)
    for line in readlines(rcfile)
        if strip(line) == "export PATH=\"\$HOME/.julia/bin:\$PATH\"" ||
           strip(line) == "export PATH=\"$(PATH.default_julia_bin()):\$PATH\""
            return true
        end
    end
    return false
end

function _contain_fpath(rcfile)
    for line in readlines(rcfile)
        if strip(line) == "export FPATH=\$HOME/.julia/completions:\$FPATH" ||
           strip(line) == "export FPATH=\"$(PATH.default_julia_fpath()):\$FPATH\""
            return true
        end
    end
    return false
end

function install_env_path(; yes::Bool = false)
    shell = detect_shell()

    config_file = ""
    if shell == "zsh"
        config_file = joinpath((haskey(ENV, "ZDOTDIR") ? ENV["ZDOTDIR"] : homedir()), ".zshrc")
    elseif shell == "bash"
        config_file = joinpath(homedir(), ".bashrc")
    else
        @warn "auto installation for $shell is not supported, please open an issue under Comonicon.jl"
    end

    write_path(joinpath(homedir(), config_file), yes)
end

"""
    write_path(rcfile[, yes=false])

Write `PATH` and `FPATH` to current shell's rc files (.zshrc, .bashrc)
if they do not exists.
"""
function write_path(rcfile, yes::Bool = false, env = ENV)
    isempty(rcfile) && return

    script = []
    msg = "cannot detect $(PATH.default_julia_bin()) in PATH, do you want to add it in PATH?"

    if !contain_comonicon_path(rcfile, env) && Tools.prompt(msg, yes)
        push!(
            script,
            """
# generated by Comonicon
# Julia bin PATH
export PATH="$(PATH.default_julia_bin()):\$PATH"
""",
        )
        @info "adding PATH to $rcfile"
    end

    msg = "cannot detect $(PATH.default_julia_fpath()) in FPATH, do you want to add it in FPATH?"
    if !contain_comonicon_fpath(rcfile, env) && Tools.prompt(msg, yes)
        push!(
            script,
            """
# generated by Comonicon
# Julia autocompletion PATH
export FPATH="$(PATH.default_julia_fpath()):\$FPATH"
autoload -Uz compinit && compinit
""",
        )
        @info "adding FPATH to $rcfile"
    end

    # exit if nothing to add
    isempty(script) && return
    # NOTE: we don't create the file if not exists
    open(rcfile, "a") do io
        write(io, "\n" * join(script, "\n"))
    end
    @info "open a new terminal, or source $rcfile to enable the new PATH."
    return
end

end # BuildTools


text:
import GeometryBasics, GeometryTypes
using ..TopOptProblems: getdh

"""
map Ferrite cell type to VTKDataTypes cell type
"""
const ferrite_to_vtk = Dict( Triangle => 5,
                            QuadraticTriangle => 22,
                            Quadrilateral => 9,
                            QuadraticQuadrilateral => 23,
                            Tetrahedron => 10,
                            QuadraticTetrahedron => 24,
                            Hexahedron => 12,
                            QuadraticHexahedron => 25
                          )

"""
Converting a Ferrite grid to a VTKUnstructuredData from [VTKDataTypes](https://github.com/mohamed82008/VTKDataTypes.jl).
"""
function VTKDataTypes.VTKUnstructuredData(grid::Ferrite.Grid{dim, <:Ferrite.Cell{dim,N,M}, T}) where {dim, N, M, T}
    celltype = ferrite_to_vtk[eltype(grid.cells)]
    celltypes = [celltype for i in 1:length(grid.cells)]
    connectivity = copy(reinterpret(NTuple{N, Int}, grid.cells))
    node_coords = copy(reshape(reinterpret(Float64, grid.nodes), dim, length(grid.nodes)))
    return VTKUnstructuredData(node_coords, celltypes, connectivity)
end
function VTKDataTypes.VTKUnstructuredData(problem::AbstractTopOptProblem)
    return VTKUnstructuredData(getdh(problem).grid)
end
VTKDataTypes.GLMesh(grid::Ferrite.Grid; kwargs...) = GLMesh(VTKUnstructuredData(grid); kwargs...)
VTKDataTypes.GLMesh(problem::AbstractTopOptProblem; kwargs...) = GLMesh(VTKUnstructuredData(problem); kwargs...)

```
workaround taken from https://github.com/JuliaPlots/Makie.jl/issues/647
TODO: should directly convert VTKDataTypes.VTKUnstructuredData to GeometryBasics.Mesh
Do not want to spend more time on this now...
```
function GeometryBasics.Mesh(glmesh::GeometryTypes.GLNormalVertexcolorMesh)
    newverts = reinterpret(GeometryBasics.Point{3, Float32}, glmesh.vertices)
    newfaces = reinterpret(GeometryBasics.NgonFace{3, GeometryBasics.OffsetInteger{-1, UInt32}}, glmesh.faces)
    newnormals = reinterpret(GeometryBasics.Vec{3, Float32}, glmesh.normals)
    return GeometryBasics.Mesh(GeometryBasics.meta(newverts; normals = newnormals, color = glmesh.color), newfaces)
end

"""
Get mesh of the topopt problem with a given topology indicator vector
"""
function GeometryBasics.Mesh(problem::AbstractTopOptProblem, topology::Array{T,1}; kwargs...) where {T}
    mesh = VTKUnstructuredData(problem)
    topology = round.(topology)
    inds = findall(isequal(0), topology)
    deleteat!(mesh.cell_connectivity, inds)
    deleteat!(mesh.cell_types, inds)
    topology = topology[setdiff(1:length(topology), inds)]
    mesh.cell_data["topology"] = topology
    glmesh = GLMesh(mesh, color = "topology"; kwargs...)
    return GeometryBasics.Mesh(glmesh)
end


text:
function sumsqrtn_fixed(n)
r = 0.0
for i = 1:n
r = r + sqrt(i)
end
return r
end

max_stars_repo_path: Module 2/Chapter03/sumsqrtn_fixed(n).jl

max_issues_repo_path: Module 2/Chapter03/sumsqrtn_fixed(n).jl

max_forks_repo_path: Module 2/Chapter03/sumsqrtn_fixed(n).jl


text:
@doc raw"""
`SPSDE`: Stratonovich Split Partitioned Stochastic Differential Equation

Defines a partitioned stochastic differential initial value problem
```math
\begin{align*}
\dq (t) &=   v(t, q(t)) \, dt + B(t, q(t)) \circ dW , & q(t_{0}) &= q_{0} ,
\dp (t) &= [ f1(t, q(t)) + f2(t, q(t)) ] \, dt + [ G1(t, q(t)) + G2(t, q(t)) ] \circ dW , & p(t_{0}) &= p_{0}
\end{align*}
```
with the drift vector fields ``v`` and ``fi``, diffusion matrices ``B`` and ``Gi``,
initial conditions ``q_{0}`` and ``p_{0}``, the dynamical variables ``(q,p)`` taking
values in ``\mathbb{R}^{d} \times \mathbb{R}^{d}``, and the m-dimensional Wiener process W

### Fields

* `d`:  dimension of dynamical variable ``q`` and the vector fields ``vi``
* `m`:  dimension of the Wiener process
* `n`:  number of initial conditions
* `ns`: number of sample paths
* `v` :  function computing the drift vector field for the position variable q
* `f1`:  function computing the drift vector field for the momentum variable p
* `f2`:  function computing the drift vector field for the momentum variable p
* `B` :  function computing the d x m diffusion matrix for the position variable q
* `G1`:  function computing the d x m diffusion matrix for the momentum variable p
* `G2`:  function computing the d x m diffusion matrix for the momentum variable p
* `t₀`: initial time
* `q₀`: initial condition for dynamical variable ``q`` (may be a random variable itself)
* `p₀`: initial condition for dynamical variable ``p`` (may be a random variable itself)

The functions `v`, `f`, 'B' and `G`, providing the drift vector fields and diffusion matrices, take four arguments,
`v(t, q, p, v)`, `f(t, q, p, f)`, `B(t, q, p,  B)` and `G(t, q, p, G)`, where `t` is the current time, `(q, p)` is the
current solution vector, and `v`, `f`, 'B' and `G` are the variables which hold the result
of evaluating the vector fields ``v``, ``f`` and the matrices ``B``, ``G`` on `t` and `(q,p)`.

### Example

```julia
    function v(λ, t, q, v)
        v[1] = λ*q[1]
        v[2] = λ*q[2]
    end

    function B(μ, t, q, B)
        B[1] = μ*q[1]
        B[2] = μ*q[2]
    end

    t₀ = 0.
    q₀ = [1., 1.]
    λ  = 2.
    μ  = 1.

    v_sde = (t, q, v) -> v(λ, t, q, v)
    B_sde = (t, q, B) -> B(μ, t, q, B)

    sde = SDE(v_sde, B_sde, t₀, q₀)
```
"""
struct SPSDE{dType <: Number, tType <: Number, vType <: Function, f1Type <: Function, f2Type <: Function, BType <: Function, G1Type <: Function, G2Type <: Function, N} <: Equation{dType, tType}
    d::Int
    m::Int
    n::Int
    ns::Int
    v ::vType
    f1::f1Type
    f2::f2Type
    B ::BType
    G1::G1Type
    G2::G2Type
    t₀::tType
    q₀::Array{dType, N}           #Initial condition: N=1 - single deterministic, N=2 - single random or multiple deterministic, N=3 - multiple deterministic
    p₀::Array{dType, N}
    periodicity::Vector{dType}

    function SPSDE{dType,tType,vType,f1Type,f2Type,BType,G1Type,G2Type,N}(d, m, n, ns, v, f1, f2, B, G1, G2, t₀, q₀, p₀; periodicity=[]) where {dType <: Number, tType <: Number, vType <: Function, f1Type <: Function, f2Type <: Function, BType <: Function, G1Type <: Function, G2Type <: Function, N}

        @assert dType == eltype(q₀) == eltype(p₀)
        @assert tType == typeof(t₀)
        @assert ndims(q₀) == ndims(p₀) == N
        @assert d == size(q₀,1) == size(p₀,1)

        if ns==1 && n==1
            # single sample path and single initial condition, therefore N=1
            @assert N == 1
        elseif ns>1 && n==1
            # single initial condition, but may be random (N=2) or deterministic (N=1)
            @assert N ∈ (1,2)
            if N==2
                @assert ns == size(q₀,2)
            end
        elseif ns==1 && n>1
            # multiple  deterministic initial conditions, so N=2
            @assert N == 2
            @assert n == size(q₀,2)
        elseif ns>1 && n>1
            # either multiple random initial conditions (N=3),
            # or multiple deterministic initial conditions (N=2)
            @assert N ∈ (2,3)

            if N==2
                @assert n == size(q₀,2)
            else
                @assert ns == size(q₀,2)
                @assert n == size(q₀,3)
            end
        end

        if !(length(periodicity) == d)
            periodicity = zeros(dType, d)
        end

        new(d, m, n, ns, v, f1, f2, B, G1, G2, t₀, q₀, p₀, periodicity)
    end
end

function SPSDE(m::Int, ns::Int, v::VT, f1::F1T, f2::F2T, B::BT, G1::G1T, G2::G2T, t₀::TT, q₀::DenseArray{DT,1}, p₀::DenseArray{DT,1}; periodicity=[]) where {DT,TT,VT,F1T,F2T,BT,G1T,G2T}
    # A 1D array q₀ contains a single deterministic initial condition, so n=1, but we still need to specify
    # the number of sample paths ns
    @assert size(q₀) == size(p₀)
    SPSDE{DT, TT, VT, F1T, F2T, BT, G1T, G2T, 1}(size(q₀, 1), m, 1, ns, v, f1, f2, B, G1, G2, t₀, q₀, p₀, periodicity=periodicity)
end

# A 2-dimensional matrix q0 can represent a single random initial condition with ns>1 and n=1,
# or a set of deterministic initial conditions with n>1 (for which we can have both ns=1 and ns>1)
# The function below assumes q0 to represent a single random initial condition (n=1, ns=size(q₀, 2))
function SPSDE(m::Int, v::VT, f1::F1T, f2::F2T, B::BT, G1::G1T, G2::G2T, t₀::TT, q₀::DenseArray{DT,2}, p₀::DenseArray{DT,2}; periodicity=[]) where {DT,TT,VT,F1T,F2T,BT,G1T,G2T}
    SPSDE{DT, TT, VT, F1T, F2T, BT, G1T, G2T, 2}(size(q₀, 1), m, 1, size(q₀, 2), v, f1, f2, B, G1, G2, t₀, q₀, p₀, periodicity=periodicity)
end

# On the other hand, the function below assumes q₀ represents multiple deterministic initial conditions
# (n=size(q₀, 2)), but these initial conditions may be run an arbitrary number ns of sample paths, so ns has to be explicitly specified
function SPSDE(m::Int, ns::Int, v::VT, f1::F1T, f2::F2T, B::BT, G1::G1T, G2::G2T, t₀::TT, q₀::DenseArray{DT,2}, p₀::DenseArray{DT,2}; periodicity=[]) where {DT,TT,VT,F1T,F2T,BT,G1T,G2T}
    SPSDE{DT, TT, VT, F1T, F2T, BT, G1T, G2T, 2}(size(q₀, 1), m, size(q₀, 2), ns, v, f1, f2, B, G1, G2, t₀, q₀, p₀, periodicity=periodicity)
end

function SPSDE(m::Int, v::VT, f1::F1T, f2::F2T, B::BT, G1::G1T, G2::G2T, t₀::TT, q₀::DenseArray{DT,3}, p₀::DenseArray{DT,3}; periodicity=[]) where {DT,TT,VT,F1T,F2T,BT,G1T,G2T}
    @assert size(q₀) == size(p₀)
    SPSDE{DT, TT, VT, F1T, F2T, BT, G1T, G2T, 3}(size(q₀, 1), m, size(q₀,3), size(q₀,2), v, f1, f2, B, G1, G2, t₀, q₀, p₀, periodicity=periodicity)
end

function SPSDE(m::Int, ns::Int, v::VT, f1::F1T, f2::F2T, B::BT, G1::G1T, G2::G2T, q₀::DenseArray{DT,1}, p₀::DenseArray{DT,1}; periodicity=[]) where {DT,VT,F1T,F2T,BT,G1T,G2T}
    SPSDE(m, ns, v, f1, f2, B, G1, G2, zero(DT), q₀, p₀, periodicity=periodicity)
end

# Assumes q0 represents a single random initial condition (n=1, ns=size(q₀, 2))
function SPSDE(m::Int, v::VT, f1::F1T, f2::F2T, B::BT, G1::G1T, G2::G2T, q₀::DenseArray{DT,2}, p₀::DenseArray{DT,2}; periodicity=[]) where {DT,VT,F1T,F2T,BT,G1T,G2T}
    SPSDE(m, v, f1, f2, B, G1, G2, zero(DT), q₀, p₀, periodicity=periodicity)
end

# Assumes q₀ represents multiple deterministic initial conditions (n=size(q₀, 2))
function SPSDE(m::Int, ns::Int, v::VT, f1::F1T, f2::F2T, B::BT, G1::G1T, G2::G2T, q₀::DenseArray{DT,2}, p₀::DenseArray{DT,2}; periodicity=[]) where {DT,VT,F1T,F2T,BT,G1T,G2T}
    SPSDE(m, ns, v, f1, f2, B, G1, G2, zero(DT), q₀, p₀, periodicity=periodicity)
end

function SPSDE(m::Int, v::VT, f1::F1T, f2::F2T, B::BT, G1::G1T, G2::G2T, q₀::DenseArray{DT,3}, p₀::DenseArray{DT,3}; periodicity=[]) where {DT,VT,F1T,F2T,BT,G1T,G2T}
    SPSDE(m, v, f1, f2, B, G1, G2, zero(DT), q₀, p₀, periodicity=periodicity)
end

Base.hash(sde::SPSDE, h::UInt) = hash(sde.d, hash(sde.m, hash(sde.n, hash(sde.ns, hash(sde.v, hash(sde.f1, hash(sde.f2, hash(sde.B, hash(sde.G1, hash(sde.G2, hash(sde.t₀, hash(sde.q₀, hash(sde.p₀, hash(sde.periodicity, h))))))))))))))

Base.:(==)(sde1::SPSDE, sde2::SPSDE) = (
                                sde1.d == sde2.d
                             && sde1.m == sde2.m
                             && sde1.n == sde2.n
                             && sde1.ns == sde2.ns
                             && sde1.v  == sde2.v
                             && sde1.f1 == sde2.f1
                             && sde1.f2 == sde2.f2
                             && sde1.B  == sde2.B
                             && sde1.G1 == sde2.G1
                             && sde1.G2 == sde2.G2
                             && sde1.t₀ == sde2.t₀
                             && sde1.q₀ == sde2.q₀
                             && sde1.p₀ == sde2.p₀
                             && sde1.periodicity == sde2.periodicity)

function Base.similar(sde::SPSDE{DT,TT,VT,F1T,F2T,BT,G1T,G2T}, q₀::DenseArray{DT,1}, p₀::DenseArray{DT,1}, ns::Int) where {DT, TT, VT, F1T, F2T, BT, G1T, G2T}
    similar(sde, sde.t₀, q₀, p₀, ns)
end

# Assumes q0 represents a single random initial condition (n=1, ns=size(q₀, 2))
function Base.similar(sde::SPSDE{DT,TT,VT,F1T,F2T,BT,G1T,G2T}, q₀::DenseArray{DT,2}, p₀::DenseArray{DT,2}) where {DT, TT, VT, F1T, F2T, BT, G1T, G2T}
    similar(sde, sde.t₀, q₀, p₀)
end

# Assumes q₀ represents multiple deterministic initial conditions (n=size(q₀, 2))
function Base.similar(sde::SPSDE{DT,TT,VT,F1T,F2T,BT,G1T,G2T}, q₀::DenseArray{DT,2}, p₀::DenseArray{DT,2}, ns::Int) where {DT, TT, VT, F1T, F2T, BT, G1T, G2T}
    similar(sde, sde.t₀, q₀, p₀, ns)
end

function Base.similar(sde::SPSDE{DT,TT,VT,F1T,F2T,BT,G1T,G2T}, q₀::DenseArray{DT,3}, p₀::DenseArray{DT,3}) where {DT, TT, VT, F1T, F2T, BT, G1T, G2T}
    similar(sde, sde.t₀, q₀, p₀)
end

function Base.similar(sde::SPSDE{DT,TT,VT,F1T,F2T,BT,G1T,G2T}, t₀::TT, q₀::DenseArray{DT,1}, p₀::DenseArray{DT,1}, ns::Int) where {DT, TT, VT, F1T, F2T, BT, G1T, G2T}
    @assert size(q₀) == size(p₀)
    @assert sde.d == size(q₀,1)
    SPSDE(sde.m, ns, sde.v, sde.f1, sde.f2, sde.B, sde.G1, sde.G2, t₀, q₀, p₀, periodicity=sde.periodicity)
end

# Assumes q0 represents a single random initial condition (n=1, ns=size(q₀, 2))
function Base.similar(sde::SPSDE{DT,TT,VT,F1T,F2T,BT,G1T,G2T}, t₀::TT, q₀::DenseArray{DT,2}, p₀::DenseArray{DT,2}) where {DT, TT, VT, F1T, F2T, BT, G1T, G2T}
    @assert size(q₀) == size(p₀)
    @assert sde.d == size(q₀,1)
    SPSDE(sde.m, sde.v, sde.f1, sde.f2, sde.B, sde.G1, sde.G2, t₀, q₀, p₀, periodicity=sde.periodicity)
end

# Assumes q₀ represents multiple deterministic initial conditions (n=size(q₀, 2))
function Base.similar(sde::SPSDE{DT,TT,VT,F1T,F2T,BT,G1T,G2T}, t₀::TT, q₀::DenseArray{DT,2}, p₀::DenseArray{DT,2}, ns::Int) where {DT, TT, VT, F1T, F2T, BT, G1T, G2T}
    @assert size(q₀) == size(p₀)
    @assert sde.d == size(q₀,1)
    SPSDE(sde.m, ns, sde.v, sde.f1, sde.f2, sde.B, sde.G1, sde.G2, t₀, q₀, p₀, periodicity=sde.periodicity)
end

function Base.similar(sde::SPSDE{DT,TT,VT,F1T,F2T,BT,G1T,G2T}, t₀::TT, q₀::DenseArray{DT,3}, p₀::DenseArray{DT,3}) where {DT, TT, VT, F1T, F2T, BT, G1T, G2T}
    @assert size(q₀) == size(p₀)
    @assert sde.d == size(q₀,1)
    PSDE(sde.m, sde.v, sde.f1, sde.f2, sde.B, sde.G1, sde.G2, t₀, q₀, p₀, periodicity=sde.periodicity)
end

Base.ndims(sde::SPSDE) = sde.d


text:
if !("." in LOAD_PATH) # for easier local testing
    push!(LOAD_PATH, ".")
end

include("TemperatureProfiles.jl")
include("relations.jl")


text:
using QuantumControlBase
using QuantumPropagators
using Documenter

DocMeta.setdocmeta!(QuantumControlBase, :DocTestSetup, :(using QuantumControlBase); recursive=true)

makedocs(;
    modules=[QuantumControlBase],
    authors="Michael Goerz <mail@michaelgoerz.net>, Alastair Marshall <alastair@nvision-imaging.com>, and contributors",
    sitename="QuantumControlBase.jl",
    format=Documenter.HTML(;
        prettyurls=true,
        canonical="https://juliaquantumcontrol.github.io/QuantumControlBase.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        hide("QuantumPropagators" => "quantumpropagators.md"),
        "History" => "history.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaQuantumControl/QuantumControlBase.jl",
)


text:
using BenchmarkTools
using InteractiveUtils

function trapz(a::Float64,b::Float64,n::Int64,f::Function)
    h::Float64=(b-a)/n
    sum=0.5*(f(a)+f(b))
    for i=1:n-1
        sum+=f(i*h)
    end
    sum*=h
end

@inline f(x::Float64)=+exp(-x)*x^2

function g(x::Float64)
    h=0.0
    if x<0.5
        h=-exp(-x)
    else
        h= exp(x)
    end
    return h*x^2
end

# Note: replacing f by f1 and g by g1 do not change the computing time:
# f1(x)=x < 1.5 ? +exp(-x)*x^2 :  +exp(-x)*x^2
# g1(x::Float64)=x < 0.5 ? -exp(-x)*x^2 :  exp(x)*x^2

function implicit(t::Float64)
    # implicit = root of  4*sin(x)-exp(x)+t
    # Newton iterations, starting from zero:
    x= 0.0
    F= 4*sin(x)-exp(x)+t
    while abs(F)> 1.e-15
        #
        x-= F/(4*cos(x) - exp(x))
        F=  4*sin(x)-exp(x)+t
        #
    end
    x
end

io = IOContext(stdout, :compact => false)

fw=open("RunningOn"*gethostname(),"w")
# Note: it seems that a loop like:
#
# for F in [f,g,implicit]
#   bench_res = @benchmark trapz(0.,1.,1000,F)
# end
#
# is incompatible with @benchmark... who knows why?

f(0.5)
println("f:")
bench_res = @benchmark trapz(0.,1.,1000,f);
show(io, bench_res)
write(fw,"\nf: "*string(bench_res)*"\n")

g(0.5)
println("\ng:")
bench_res = @benchmark trapz(0.,1.,1000,g);
show(io, bench_res)
write(fw,"g: "*string(bench_res)*"\n")

implicit(0.5)
println("\nimplicit:")
bench_res = @benchmark trapz(0.,1.,1000,implicit)
show(io, bench_res)
write(fw,"implicit: "*string(bench_res)*"\n")

println("\nend.")
close(fw)

max_stars_repo_name: danielcort/benchmarks-python-julia-c

max_issues_repo_name: danielcort/benchmarks-python-julia-c

max_forks_repo_name: danielcort/benchmarks-python-julia-c


text:
using CSV, DataFrames, StatsPlots, StanSample
using Distributions, Statistics, Random

ProjDir = @__DIR__

# Simulate the data

#Random.seed!(123)

countries = ["Guatamala", "Panama", "Kenya", "Ethiopia", "Papua New Guinea"]
brands = [4, 3, 4, 5, 1]
intercepts = [500, 300, 600, 700, 900] # Quantity at cheapest price_point
slopes = [-6, -8, -10, -15, -20] # Assuming constant price elasticity
price_points = [12.5, 12.75, 13.02, 14.20, 15.15] # Selected price points
costs = [4.0, 4.2, 4.6, 4.7, 5.2]
channels = ["retail", "wholesale", "shop"]

df = DataFrame()
for n in 1:1500
	yr = 2019
	wk = rand(1:52)
	c_ind = rand(1:length(countries))
	p_ind = rand(1:length(price_points))
	c = countries[c_ind]							# country
	b = rand(1:brands[c_ind])						# brand
	ch = rand(channels)								# channel
	cst = rand(Normal(costs[c_ind], 0.1))			# cost
	p = price_points[p_ind]							# ppu
	if rand(Bernoulli(0.2))
		p = rand(Normal(p, 0.07))
	end
	p = round(p, digits=2)
	q = intercepts[c_ind] + slopes[c_ind] * p_ind	# quantity sold
	if rand(Bernoulli(0.2))
		q = round(rand(Normal(q, 30)))
	end
	append!(df, DataFrame(
		year=yr, week=wk,
		country=c, brand=b, channel=ch,
		quantity=q, ppu=p)
	)
end

df[!, :ppu_l] = log.(df[:, :ppu])
df[!, :quantity_l] = log.(df[:, :quantity])

# Clean up

df = df[df.ppu_l .> 0.04, :]
df = df[df.quantity_l .> 0.0, :]

CSV.write(joinpath(ProjDir, "data", "data.csv"), df)

bglm = "
data {
  int<lower=0> N;
  vector[N] q; // Outcome
  vector[N] p; // Predictor
}
parameters {
  real<lower=0> a;
  real b;
  real<lower=0> sigma;
}
model {
  vector[N] mu;            		// mu is a vector
  a ~ normal(2.0, 10.0);		// intercept positive
  b ~ normal(-2.0, 10.0);		// demand drops as price increases
  sigma ~ uniform(0, 15);
  mu = a + b * p;
  q ~ normal(mu, sigma);
}
// generated quantities {
//    vector[N] mu_pred;
//    mu_pred = exp(p);
//  }
";

mean_q_l = mean(df[:, :quantity_l])
mean_p_l = mean(df[:, :ppu_l])
data = Dict(
  :N => size(df, 1),
  :q => df[:, :quantity_l] .- mean_q_l,
  :p => df[:, :ppu_l] .- mean_p_l
)

p1 = scatter(df[:, :ppu], df[:, :quantity],
	xlabel="ppu", ylabel="quantity sold", leg=false)
p2 = scatter(data[:p], data[:q],
	xlabel="log(ppu)", ylabel="log(quantity sold)", leg=false)
plot(p1, p2, layout=(2, 1))
savefig(joinpath(ProjDir, "plots", "data.png"))

sm = SampleModel("bglm", bglm)

rc = stan_sample(sm, data=data)

if success(rc)
	chns = read_samples(sm)
	plot(chns)
	savefig(joinpath(ProjDir, "plots", "chains.png"))

	dfs = DataFrame(chns)
	show(chns)

	p = -0.1:0.01:0.15
	p1 = scatter(data[:p], data[:q],
		xlims=(-0.1, 0.15), ylims=(-1, 1),
		xlabel="ppu_l", ylabel="quantity_l sold", leg=false)
	a_mean = mean(dfs[:, :a])
	b_mean = mean(dfs[:, :b])
	mu_q = a_mean .+ b_mean .* p
	plot!(p1, p, mu_q)
	savefig(joinpath(ProjDir, "plots", "regression.png"))

	price_points_l = log.(price_points) .- mean(log.(price_points))
	scatter(price_points_l, a_mean .+ b_mean .* price_points_l,
		xlims = (-0.1, 0.15), ylims= (-0.1, 0.1),
		leg=false, xlabel="log(ppu)", ylabel="log(quantity sold)")
	selected_rows = rand(1:size(dfs, 1), 75)
	for row in eachrow(dfs[selected_rows,:])
		plot!(p, row[:a] .+ row[:b] .* p, color=:gray)
	end
	scatter!(price_points_l, a_mean .+ b_mean .* price_points_l)
	savefig(joinpath(ProjDir, "plots", "sampling.png"))
end


text:
# This is just a bit of code to set up a temporary package environment since not everyone
# will have the same environment set up; this is largely unnecessary when working only on
# your own code.
# See https://bkamins.github.io/julialang/2020/06/28/automatic-project-environments.html
using Pkg: Pkg
# Change to a newly created temporary directory
cd(mktempdir()) do
    Pkg.activate(".")  # activate a new project in the temporary directory
    Pkg.UPDATED_REGISTRY_THIS_SESSION[] = true  # don't update the registry; on Windows it's slow...
    # Add the required packages
    Pkg.add(Pkg.PackageSpec(name="Plots"))
end

# Load required packages
using Plots

function integrate_sde(sde, x0, t1; dt=0.01)
    # Euler-Maruyama for a 1D SDE
    x = x0
    for t in range(0, t1, step=dt)
        (f, g) = sde(x)
        x = x + f*dt + g*sqrt(dt)*randn()
    end
    return x
end

function integrate_many(sde, x0, t1; n = 1, kwargs...)
    # Integrate many sample paths using Euler-Maruyama
    x = zeros(n)
    Threads.@threads for i in eachindex(x)  # use threads to speed up the computation
        x[i] = integrate_sde(sde, x0, t1; kwargs...)
    end
    return x
end

ornstein_uhlenbeck(x, θ, σ) = (-θ*x, σ) # Ornstein-Uhlenbeck equation (this is short hand notation for a function)
cubic(x, θ, σ) = (θ*x - x^3, σ) # SDE with two potential wells

println("Using $(Threads.nthreads()) threads")

x = integrate_many(x->ornstein_uhlenbeck(x, 1.0, 0.1), 0.0, 10.0, n = 1_000_000) # integrate a million sample paths
plt1 = histogram(x)
title!("Ornstein-Uhlenbeck: θ=1.0, σ=0.1")

x = integrate_many(x->cubic(x, 1.0, 0.1), 0.0, 10.0, n = 1_000_000) # integrate a million sample paths
plt2 = histogram(x)
title!("Cubic: θ=1.0, σ=0.1")

x = integrate_many(x->cubic(x, 1.0, 1.5), 0.0, 10.0, n = 1_000_000) # integrate a million sample paths
plt3 = histogram(x)
title!("Cubic: θ=1.0, σ=1.5")

# Use `gui(plt1)` etc to show the different figures

max_stars_repo_name: dawbarton/julia-de-examples

max_issues_repo_name: dawbarton/julia-de-examples

max_forks_repo_name: dawbarton/julia-de-examples


text:
export AvgCutPruningAlgo, AvgCutPruner

"""
    AvgCutPruningAlgo <: AbstractCutPruningAlgo

Removes the cuts with lower trust where the trust is: nused / nwith + bonus
where the cut has been used `nused` times amoung `nwith` optimization done with it.
We say that the cut was used if its dual value is nonzero.
It has a bonus equal to `mycutbonus` if the cut was generated using a trial given by the problem using this cut.
If `nwidth` is zero, `nused/nwith` is replaced by `newcuttrust`.
"""
struct AvgCutPruningAlgo <: AbstractCutPruningAlgo
    # maximum number of cuts
    maxncuts::Int
    newcuttrust::Float64
    mycutbonus::Float64
    function AvgCutPruningAlgo(maxncuts::Int, newcuttrust=3/4, mycutbonus=1/4)
        new(maxncuts, newcuttrust, mycutbonus)
    end
end

mutable struct AvgCutPruner{N, T} <: AbstractCutPruner{N, T}
    # used to generate cuts
    isfun::Bool
    islb::Bool
    lazy_minus::Bool
    A::AbstractMatrix{T}
    b::AbstractVector{T}

    # maximum number of cuts
    maxncuts::Int

    # number of optimization performed
    nwith::Vector{Int}
    # number of times where the cuts have been used
    nused::Vector{Int}
    mycut::Vector{Bool}
    trust::Union{Vector{Float64}, Nothing}
    ids::Vector{Int} # small id means old
    id::Int # current id

    newcuttrust::Float64
    mycutbonus::Float64

    # check redundancy with old cuts when adding new cuts
    excheck::Bool

    # tolerance to check redundancy between two cuts
    TOL_EPS::Float64

    function (::Type{AvgCutPruner{N, T}})(sense::Symbol, maxncuts::Int, newcuttrust=3/4, mycutbonus=1/4, lazy_minus::Bool=false, tol=1e-6, excheck::Bool=false) where {N, T}
        isfun, islb = gettype(sense)
        new{N, T}(isfun, islb, lazy_minus, spzeros(T, 0, N), T[], maxncuts, Int[], Int[], Bool[], nothing, Int[], 0, newcuttrust, mycutbonus, excheck, tol)
    end
end

(::Type{CutPruner{N, T}})(algo::AvgCutPruningAlgo, sense::Symbol, lazy_minus::Bool=false) where {N, T} = AvgCutPruner{N, T}(sense, algo.maxncuts, algo.newcuttrust, algo.mycutbonus, lazy_minus)

# COMPARISON
"""Update cuts relevantness after a solver's call returning dual vector `σρ`."""
function addusage!(man::AvgCutPruner, σρ)
    if ncuts(man) > 0
        man.nwith .+= 1
        # TODO: dry 1e-6 in CutPruner?
        man.nused[σρ .> 1e-6] .+= 1
        man.trust = nothing # need to be recomputed
    end
end

function gettrustof(man::AvgCutPruner, nwith, nused, mycut)
    (nwith == 0 ? man.newcuttrust : nused / nwith) + (mycut ? man.mycutbonus : 0)
end
function initialtrust(man::AvgCutPruner, mycut)
    gettrustof(man, 0, 0, mycut)
end
function hastrust(man::AvgCutPruner)
    man.trust !== nothing
end
function gettrust(man::AvgCutPruner)
    if !hastrust(man)
        trust = man.nused ./ man.nwith
        trust[man.nwith .== 0] .= man.newcuttrust
        trust[man.mycut] .+= man.mycutbonus
        man.trust = trust
    end
    man.trust
end

# CHANGE

function keeponlycuts!(man::AvgCutPruner, K::AbstractVector{Int})
    man.nwith = man.nwith[K]
    man.nused = man.nused[K]
    man.mycut = man.mycut[K]
    _keeponlycuts!(man, K)
end

function replacecuts!(man::AvgCutPruner, K::AbstractVector{Int}, A, b, mycut::AbstractVector{Bool})
    man.nwith[K] .= 0
    man.nused[K] .= 0
    man.mycut[K] .= mycut
    _replacecuts!(man, K, A, b)
    if hastrust(man)
        man.trust[K] .= initialtrusts(man, mycut)
    end
end

function appendcuts!(man::AvgCutPruner, A, b, mycut::AbstractVector{Bool})
    n = length(mycut)
    append!(man.nwith, zeros(n))
    append!(man.nused, zeros(n))
    append!(man.mycut, mycut)
    _appendcuts!(man, A, b)
    if hastrust(man)
        append!(man.trust, initialtrusts(man, mycut))
    end
end


text:
using OneDNN

# stdlib
using Test
using Random
import Statistics

# deps
import CEnum
import Flux
import ProgressMeter
using Zygote

# include("tiled.jl")
include("utils.jl")
include("memory.jl")
include("ops/simple.jl")
include("ops/matmul.jl")
include("ops/batchnorm.jl")
include("ops/innerproduct.jl")
include("ops/pool.jl")
include("ops/concat.jl")
include("ops/convolution.jl")

# include("bf16.jl")


text:
macro declare_file_product(product_name)
    get_path_name = Symbol(string("get_", product_name, "_path"))
    path_name = Symbol(string(product_name, "_path"))
    return esc(quote
        # These will be filled in by init_file_product().
        $(path_name) = ""
        $(product_name) = ""
        function $(get_path_name)()
            return $(path_name)::String
        end
    end)
end

macro init_file_product(product_name, product_path)
    path_name = Symbol(string(product_name, "_path"))
    return esc(quote
        # FileProducts are very simple, and we maintain the `_path` suffix version for consistency
        global $(path_name) = joinpath(artifact_dir, $(product_path))
        global $(product_name) = $(path_name)
    end)
end


text:
using BinDeps

CSIP_VERSION = "0.3.4"
CSIP_URL = "https://github.com/SCIP-Interfaces/CSIP/archive/v$(CSIP_VERSION).zip"
CSIP_LIB = "libcsip"
CSIP_UNPACKED = "CSIP-$(CSIP_VERSION)"

@BinDeps.setup

@assert haskey(ENV, "SCIPOPTDIR") "Environment variable `SCIPOPTDIR` not set!"

function validate_csip(name, handle)
    csip_version = ccall(Libdl.dlsym(handle, :CSIPgetVersion), Cint, ())
    csip_version == 034
end

csipdep = library_dependency(CSIP_LIB, validate=validate_csip)

provides(Sources, Dict(URI(CSIP_URL) => csipdep),
                 unpacked_dir=CSIP_UNPACKED)

provides(SimpleBuild,
    (@build_steps begin
         GetSources(csipdep)
         @build_steps begin
             ChangeDirectory(joinpath(BinDeps.srcdir(csipdep), CSIP_UNPACKED))
             `make`
             `mkdir -p $(libdir(csipdep))`
             `install lib/libcsip.so $(libdir(csipdep))`
         end
    end), csipdep, os = :Unix)

@BinDeps.install Dict(:libcsip => :libcsip)