text:
using HighDimPDE
using Random
using Test
using Statistics
using Revise

#relative error l2
function rel_error_l2(u, uanal)
    if abs(uanal) >= 10 * eps(eltype(uanal))
        sqrt((u - uanal)^2 / u^2)
    else # overflow
        abs(u-uanal)
    end
end

# https://en.wikipedia.org/wiki/Heat_equation#Fundamental_solutions
@testset "MLP - heat equation - single thread" begin # solving at one unique point
    for d in [1, 10, 50]
        println("test for d = ", d)
        tspan = (0e0, 5e-1)
        dt = 5e-1  # time step
        μ(x, p, t) = 0e0 # advection coefficients
        σ(x, p, t) = 1e-1 # diffusion coefficients

        u_anal(x, t) = sum(x.^2) + d * σ(0., 0., 0.)^2 * t
        g(x) = sum(x.^2)

        # d = 10
        x0 = fill(2e0,d)

        alg = MLP()

        f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0e0 .* v_y #TODO: this fix is not nice

        # defining the problem
        prob = PIDEProblem(g, f, μ, σ, tspan, x = x0)
        # solving
        xs,ts,sol = solve(prob, alg, multithreading = false)
        u1 = sol[end]
        u1_anal = u_anal(x0, tspan[end])
        e_l2 = rel_error_l2(u1, u1_anal)
        println("rel_error_l2 = ", e_l2, "\n")
        @test e_l2 < 0.1
    end
end

# https://en.wikipedia.org/wiki/Heat_equation#Fundamental_solutions
@testset "MLP - heat equation - multi thread" begin # solving at one unique point
    for d in [1, 10, 50]
        println("test for d = ", d)
        tspan = (0e0, 5e-1)

        μ(x, p, t) = 0e0 # advection coefficients
        σ(x, p, t) = 1e-1 # diffusion coefficients

        u_anal(x, t) = sum(x.^2) + d * σ(0., 0., 0.)^2 * t
        g(x) = sum(x.^2)

        # d = 10
        x0 = fill(2e0,d)

        alg = MLP()

        f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0e0 .* v_y #TODO: this fix is not nice

        # defining the problem
        prob = PIDEProblem(g, f, μ, σ, tspan, x = x0)
        # solving
        xs,ts,sol = solve(prob, alg, multithreading = true)
        u1 = sol[end]
        u1_anal = u_anal(x0, tspan[end])
        e_l2 = rel_error_l2(u1, u1_anal)
        println("rel_error_l2 = ", e_l2, "\n")
        @test e_l2 < 0.1
    end
end

@testset "MLP - heat equation - single thread - Neumann BC" begin # solving at one unique point
    tspan = (0f0, 5f-1)
    μ(x, p, t) = 0f0 # advection coefficients
    σ(x, p, t) = 1f-1 #1f-1 # diffusion coefficients
    for d in [1,2,5]
        u1s = []
        for _ in 1:2
            neumann_bc = (fill(-5e-1, d), fill(5e-1, d))
            g(x) = sum(x.^2)
            # d = 10
            x0 = fill(3e-1,d)
            alg = MLP()
            f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = 0e0 .* v_y #TODO: this fix is not nice

            # defining the problem
            prob = PIDEProblem(g, f, μ, σ, tspan, x = x0, neumann_bc = neumann_bc)
            # solving
            xs,ts,sol = solve(prob, alg)
            push!(u1s, sol[end])
            println("d = $d, u1 = $(sol[end])")

        end
        e_l2 = mean(rel_error_l2.(u1s[1], u1s[2]))
        println("rel_error_l2 = ", e_l2, "\n")
        @test e_l2 < 0.1
    end
end

@testset "MLP - exponential - interval" begin
    for d in [1, 3, 10]
        println("test for d = ", d)
        tspan = (0f0, 5f-1)

        μ(x, p, t) = 0f0 # advection coefficients
        σ(x, p, t) = 0f0 # diffusion coefficients
        r = 1f-1
        x = fill(0f0, d)

        u_anal(x, t) = (g(x) * exp(r * t))
        g(x) = sum(x.^2) .+ 2f0

        alg = MLP()

        f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = r * v_y #TODO: this fix is not nice

        # defining the problem
        prob = PIDEProblem(g, f, μ, σ, tspan, x=x)
        # solving
        xs,ts,sol = solve(prob, alg)
        u1 = sol[end]
        u1_anal = u_anal.(xs, tspan[end])
        e_l2 = mean(rel_error_l2.(u1, u1_anal))
        println("rel_error_l2 = ", e_l2, "\n")
        @test e_l2 < 0.1
    end
end

@testset "MLP - allen cahn" begin

    tspan = (3f-1, 6f-1)

    μ(x, p, t) = 0f0 # advection coefficients
    σ(x, p, t) = 1f0 # diffusion coefficients

    for d in [10]

        alg = MLP()

        X0 = fill(0f0,d)  # initial point
        g(X) =  1f0 ./ (2f0 .+ 4f-1 * sum(X.^2))   # initial condition
        a(u) = u - u^3
        f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = - a.(v_y)

        # defining the problem
        prob = PIDEProblem(g, f, μ, σ, tspan, x = X0 )
        # solving
        xs,ts,sol = solve(prob, alg,)
        u1 = sol[end]
        # value coming from \cite{Beck2017a}
        e_l2 = rel_error_l2(u1, 0.30879)
        @test e_l2 < 0.5 # this is quite high as a relative error.
        println("d = $d, rel_error_l2 = $e_l2")
    end
end

@testset "MLP - allen cahn local, Neumann BC" begin

    tspan = (0f0, 5f-1)

    μ(x, p, t) = 0f0 # advection coefficients
    σ(x, p, t) = 1f0 # diffusion coefficients

    for d in [1,2,5]
        neumann = (fill(-5f-1, d), fill(5f-1, d))

        alg = MLP()

        X0 = fill(0f0,d)  # initial point
        g(X) = exp.(-0.25f0 * sum(X.^2))   # initial condition
        a(u) = u - u^3
        f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = a.(v_y) # nonlocal nonlinear part of the

        # defining the problem
        prob = PIDEProblem(g, f, μ, σ, tspan, x = X0, neumann = neumann )
        # solving
        xs,ts,sol = solve(prob, alg, )
        u1 = sol[end]
        @test !isnan(u1)
        println("d = $d, u1 = $(u1)")
    end
end

@testset "MLP algorithm - Black-Scholes Equation with Default Risk" begin
    tspan = (0f0, 1f0)

    μ(x, p, t) = 0.02f0 * x # advection coefficients
    σ(x, p, t) = 0.2f0 * x # diffusion coefficients

    d = 20

    alg = MLP()

    X0 = fill(100f0,d)  # initial point
    g(X) =  minimum(X) # initial condition

    δ = 2.0f0/3f0
    R = 0.02f0

    vh = 50.0f0
    vl = 70.0f0
    γh = 0.2f0
    γl = 0.02f0

    Q(u) = (u .< vh) .* γh .+ (u .>= vl) .*  γl .+ ( (u .>= vh) .* (u .< vl)) .* (((γh - γl) / (vh - vl)) * (u .- vh) .+ γh)

    µc = 0.02f0
    σc = 0.2f0

    f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = -(1f0 - δ) * Q.(v_y) .* v_y .- R * v_y

    # defining the problem
    prob = PIDEProblem(g, f, μ, σ, tspan, x = X0 )
    # solving
    xs,ts,sol = solve(prob, alg, )

    u1 = sol[end]

    analytical_ans = 60.781
    error_l2 = rel_error_l2(u1, analytical_ans)

    @test error_l2 < 0.1

    println("d = $d, error_l2 = $(error_l2)")

end

###################################################
########### NON LOCAL #############################
###################################################

@testset "MLP - replicator mutator" begin
    tspan = (0e0,5e-1)
    σ_sampling = 1e0
    μ(x, p, t) = 0e0 # advection coefficients
    σ(x, p, t) = 1e-1 #1f-1 # diffusion coefficients
    ss0 = 1e-2#std g0

    # Analytic sol
    function _SS(x, t, p)
        d = length(x)
        MM = σ(x, p, t) * ones(d)
        SSt = MM .* ((MM .* sinh.(MM *t) .+ ss0 .* cosh.( MM * t)) ./ (MM .* cosh.(MM * t ) .+ ss0 .* sinh.(MM * t)))
        return SSt
    end

    function uanal(x, t, p)
            d = length(x)
            return (2*π)^(-d/2) * prod(_SS(x, t, p) .^(-1/2)) * exp(-0.5 *sum(x .^2 ./ _SS(x, t, p)) )
    end

    sols = []
    xs = []
    for d in [1, 2, 5]

        x = fill(0e0, d)

        alg = MLP(M=4, K=10, L = 4, mc_sample = NormalSampling(σ_sampling) )

        g(x) = (2*π)^(-d/2) * ss0^(- d * 5e-1) * exp.(-5e-1 *sum(x .^2e0 / ss0)) # initial condition
        m(x) = - 5e-1 * sum(x.^2)

        f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = max(0.0, v_y) * (m(y) - max(0.0, v_z) * m(z) * (2.0 * π)^(d/2) * σ_sampling^d * exp(0.5 * sum(z.^2) / σ_sampling^2)) # nonlocal nonlinear part of the

        # defining the problem
        prob = PIDEProblem(g, f, μ, σ, tspan,x=x)
        # solving
        xgrid,ts,sol = solve(prob, alg,)
        u1 = sol[end]
        u1_anal = uanal(x,tspan[2],Dict())
        e_l2 = mean(rel_error_l2.(u1, u1_anal))
        println("rel_error_l2 = ", e_l2, "\n")
        @test e_l2 < 0.1
        push!(sols, sol[end])
    end
end

@testset "MLP- allen cahn non local Neumann BC" begin

    tspan = (0f0, 5f-1)

    μ(x, p, t) = 0f0 # advection coefficients
    σ(x, p, t) = 1f-1 #1f-1 # diffusion coefficients

    for d in [1,2,5]
        u1s = []
        for _ in 1:2
            neumann = (fill(-5f-1, d), fill(5f-1, d))

            alg = MLP(M=4, K=10, L = 4, mc_sample = UniformSampling(neumann...) )

            x = fill(0f0,d)  # initial point
            g(X) = exp.(-0.25f0 * sum(X.^2))   # initial condition
            a(u) = u - u^3
            f(y,z,v_y,v_z,∇v_y,∇v_z, p, t) = a.(v_y) .- a.(v_z) #.* Float32(π^(d/2)) * σ_sampling^d .* exp.(sum(z.^2, dims = 1) / σ_sampling^2) # nonlocal nonlinear part of the

            # defining the problem
            prob = PIDEProblem(g, f, μ, σ, tspan, x = x, neumann = neumann)
            # solving
            xs,ts,sol = solve(prob, alg)
            push!(u1s, sol[end])
            println("d = $d, u1 = $(sol[end])")

        end
        e_l2 = mean(rel_error_l2.(u1s[1], u1s[2]))
        println("rel_error_l2 = ", e_l2, "\n")
        @test e_l2 < 0.1
    end
end


text:
struct Gencost
    "cost model, 1 = piecewise linear, 2 = polynomial "
    MODEL::Vector{Int64}
    "startup cost in US dollars"
    STARTUP::Vector{Float64}
    "shutdown cost in US dollars"
    SHUTDOWN::Vector{Float64}
    "number breakpoints in piecewise linear cost function, or number of coefficients in polynomial cost function"
    NCOST::Vector{Int64}
    "parameters defining total cost function begin in this col
        (MODEL = 1) : p0, f0, p1, f1, ..., pn, fn
             where p0 < p1 < ... < pn and the cost f(p) is defined
             by the coordinates (p0,f0), (p1,f1), ..., (pn,fn) of
             the end/break-points of the piecewise linear cost
        (MODEL = 2) : cn, ..., c1, c0
             n+1 coefficients of an n-th order polynomial cost fcn,
             starting with highest order, where cost is
             f(p) = cn*p^n + ... + c1*p + c0"
    COST::Vector{Vector{Float64}}
end

function extract_gencost(mpc)
    gencost = mpc["gencost"]
    ng = size(gencost, 1)
    coef = Vector{Vector{Float64}}([gencost[i,:][5:end] for i in 1:ng])
    return Gencost([gencost[:, i] for i in 1:4]..., coef)
end

"""
Shortcut for returning gencost info for one of the cases listed by
`casenames()`.
"""
function extract_gencost(cname::AbstractString)
    mpc = loadcase(cname; describe=false)
    extract_gencost(mpc)
end

ext: jl


text:
using Parameters

# """
# Newton-Krylov based saddle search method
# """
# @with_kw type NK
#    tol::Float64 = 1e-5
#    maxnumdE::Int = 1000
#    len::Float64 = 1e-7
#    precon = I
#    precon_prep! = (P, x) -> P
#    verbose::Int = 1
#    krylovinit::Symbol = :resrot  # allow res, rand, rot, resrot
#    maxstep::Float64 = Inf
#    eigatol::Float64 = 1e-1
#    eigrtol::Float64 = 1e-1
# end


text:
import QuantumInformation: _mdbit, _mdxit
b = SpinBasis(1//2)
up = spinup(b)

@test try tex(up); true catch false end
@test try tex(up'); true catch false end
@test try tex(im * up); true catch false end
@test try tex((1 + im) * up); true catch false end
@test try tex(-up); true catch false end
@test try tex(im * up); true catch false end
@test try tex(-im * up); true catch false end
@test try tex((-1 - im) * up); true catch false end
@test try tex(qubit("010")); true catch false end
@test try tex(qubit("010"), "ϕ"); true catch false end
@test try tex(qudit("010", 3)); true catch false end
@test try tex(quxit("010", [3,4,3])); true catch false end
@test try tex(hadamard() * up); true catch false end
@test try tex(hadamard() * up, "ρ"); true catch false end
@test try tex(sigmay() * up); true catch false end
@test try tex(sigmay() * hadamard() * up); true catch false end
@test try tex(sigmaz() * hadamard() * up); true catch false end
@test try tex(sigmax()); true catch false end
@test try tex(hadamard(2)); true catch false end
@test try tex((qubit("01") ⊗ basisstate(NLevelBasis(3), 3)) ⊗ qubit("10")'); true catch false end
@test try tex((qubit("01") ⊗ basisstate(NLevelBasis(3), 3)) ⊗ qubit("10")', "ρ"); true catch false end
@test try tex(sparse((qubit("01") ⊗ basisstate(NLevelBasis(3), 3)) ⊗ qubit("10")')); true catch false end
@test try tex(sparse((qubit("01") ⊗ basisstate(NLevelBasis(3), 3)) ⊗ qubit("10")'), "ρ"); true catch false end
@test _mdbit(qubit("0"), "psi") == "| psi \\rangle = | 0 \\rangle"
@test _mdbit(qubit("0100"), "psi") == "| psi \\rangle = | 0100 \\rangle"
@test _mdbit(qubit("0100")', "psi") == "\\langle psi | = \\langle 0100 |"
@test _mdbit(dm(qubit("0100")), "psi") == "psi = | 0100 \\rangle \\langle 0100 |"
@test _mdbit(sparse(dm(qubit("0100"))), "psi") == "psi = | 0100 \\rangle \\langle 0100 |"
@test _mdbit(H(4) * qubit("0100"), "psi") == "| psi \\rangle = 0.25 | 0000 \\rangle+0.25 | 0001 \\rangle+0.25 | 0010 \\rangle+0.25 | 0011 \\rangle-0.25 | 0100 \\rangle-0.25 | 0101 \\rangle-0.25 | 0110 \\rangle-0.25 | 0111 \\rangle+0.25 | 1000 \\rangle+0.25 | 1001 \\rangle+0.25 | 1010 \\rangle+0.25 | 1011 \\rangle-0.25 | 1100 \\rangle-0.25 | 1101 \\rangle-0.25 | 1110 \\rangle-0.25 | 1111 \\rangle"
@test _mdxit(quxit("0", [2]), "psi") == "| psi \\rangle = | 0 \\rangle"
@test _mdxit(quxit("02", [2,3]), "psi") == "| psi \\rangle = | 02 \\rangle"
@test _mdxit(quxit("02", [2,3])', "psi") == "\\langle psi | = \\langle 02 |"
@test _mdxit(dm(quxit("02", [2,3])), "psi") == "psi = | 02 \\rangle \\langle 02 |"
@test _mdxit(sparse(dm(quxit("02", [2,3]))), "psi") == "psi = | 02 \\rangle \\langle 02 |"


text:
"""
Core genie configuration / settings functionality.
"""
module Configuration

import Pkg
import Dates

using Random
using AMQPClient

import VersionCheck

function __init__()
  try
    @async VersionCheck.newversion("Genie", url = "https://genieframework.com/CHANGELOG.html")
  catch
  end
end

"""
    pkginfo(pkg::String)

Returns installed package information for `pkg`
"""
pkginfo(pkg::String) = filter(x -> x.name == pkg && x.is_direct_dep, values(Pkg.dependencies()) |> collect)

import Logging
import Genie
import MbedTLS

export isdev, isprod, istest, env
export Settings, DEV, PROD, TEST

# app environments
const DEV   = "dev"
const PROD  = "prod"
const TEST  = "test"

"""
    isdev()  :: Bool

Set of utility functions that return whether or not the current environment is development, production or testing.

# Examples
```julia
julia> Configuration.isdev()
true

julia> Configuration.isprod()
false
```
"""
isdev() :: Bool  = (Genie.config.app_env == DEV)

"""
    isprod() :: Bool

Set of utility functions that return whether or not the current environment is development, production or testing.

# Examples
```julia
julia> Configuration.isdev()
true

julia> Configuration.isprod()
false
```
"""
isprod():: Bool = (Genie.config.app_env == PROD)

"""
    istest() :: Bool

Set of utility functions that return whether or not the current environment is development, production or testing.

# Examples
```julia
julia> Configuration.isdev()
true

julia> Configuration.isprod()
false
```
"""
istest():: Bool = (Genie.config.app_env == TEST)

"""
    env() :: String

Returns the current Genie environment.

# Examples
```julia
julia> Configuration.env()
"dev"
```
"""
env() :: String = Genie.config.app_env

"""
    buildpath() :: String

Constructs the temp dir where Genie's view files are built.
"""
buildpath() :: String = Base.Filesystem.mktempdir(prefix = "jl_genie_build_")

"""
    mutable struct Settings

App configuration - sets up the app's defaults. Individual options are overwritten in the corresponding environment file.

# Arguments
- `server_port::Int`: the port for running the web server (default 8000)
- `server_host::String`: the host for running the web server (default "127.0.0.1")
- `server_document_root::String`: path to the document root (default "public/")
- `server_handle_static_files::Bool`: if `true`, Genie will also serve static files. In production, it is recommended to serve static files with a web server like Nginx.
- `server_signature::String`: Genie's signature used for tagging the HTTP responses. If empty, it will not be added.
- `app_env::String`: the environment in which the app is running (dev, test, or prod)
- `cors_headers::Dict{String,String}`: default `Access-Control-*` CORS settings
- `cors_allowed_origins::Vector{String}`: allowed origin hosts for CORS settings
- `cache_duration::Int`: cache expiration time in seconds
- `log_level::Logging.LogLevel`: logging severity level
- `log_to_file::Bool`: if true, information will be logged to file besides REPL
- `session_key_name::String`: the name of the session cookie
- `session_storage::Symbol`: the backend adapter for session storage (default File)
- `inflector_irregulars::Vector{Tuple{String,String}}`: additional irregular singular-plural forms to be used by the Inflector
- `run_as_server::Bool`: when true the server thread is launched synchronously to avoid that the script exits
- `websockets_server::Bool`: if true, the websocket server is also started together with the web server
- `html_parser_close_tag::String`: default " /". Can be changed to an empty string "" so the single tags would not be closed.
- `ssl_enabled::Bool`: default false. Server runs over SSL/HTTPS in development.
- `ssl_config::MbedTLS.SSLConfig`: default `nothing`. If not `nothing` and `ssl_enabled`, it will use the config to start the server over HTTPS.
- `webchannels_keepalive_frequency::Int`: default `30000`. Frequency in miliseconds to send keepalive messages to webchannel/websocket to keep the connection alive. Set to `0` to disable keepalive messages.
"""
Base.@kwdef mutable struct Settings
  server_port::Int                                    = (haskey(ENV, "PORT") ? parse(Int, ENV["PORT"]) : 8000) # default port for binding the web server
  server_host::String                                 = haskey(ENV, "HOST") ? ENV["HOST"] : "127.0.0.1"
  server_document_root::String                        = "public"
  server_handle_static_files::Bool                    = true
  server_signature::String                            = "Genie/Julia/$VERSION"

  app_env::String                                     = haskey(ENV, "GENIE_ENV") ? ENV["GENIE_ENV"] : DEV

  cors_headers::Dict{String,String}                   = Dict{String,String}(
                                                          "Access-Control-Allow-Origin"       => "", # ex: "*" or "http://mozilla.org"
                                                          "Access-Control-Expose-Headers"     => "", # ex: "X-My-Custom-Header, X-Another-Custom-Header"
                                                          "Access-Control-Max-Age"            => "86400", # 24 hours
                                                          "Access-Control-Allow-Credentials"  => "", # "true" or "false"
                                                          "Access-Control-Allow-Methods"      => "", # ex: "POST, GET"
                                                          "Access-Control-Allow-Headers"      => "", # ex: "X-PINGOTHER, Content-Type"
                                                        )
  cors_allowed_origins::Vector{String}                = String[]

  cache_duration::Int                                 = 0
  cache_storage::Union{Symbol,Nothing}                = nothing

  log_level::Logging.LogLevel                         = Logging.Debug
  log_to_file::Bool                                   = false
  log_requests::Bool                                  = true

  inflector_irregulars::Vector{Tuple{String,String}}  = Tuple{String,String}[]

  run_as_server::Bool                                 = false

  websockets_server::Bool                             = false
  websockets_port::Int                                = server_port

  initializers_folder::String                         = "initializers"

  path_config::String                                 = "config"
  path_env::String                                    = joinpath(path_config, "env")
  path_app::String                                    = "app"
  path_resources::String                              = joinpath(path_app, "resources")
  path_lib::String                                    = "lib"
  path_helpers::String                                = joinpath(path_app, "helpers")
  path_log::String                                    = "log"
  path_tasks::String                                  = "tasks"
  path_build::String                                  = buildpath()
  path_plugins::String                                = "plugins"
  path_cache::String                                  = "cache"
  path_initializers::String                           = joinpath(path_config, initializers_folder)
  path_db::String                                     = "db"
  path_bin::String                                    = "bin"
  path_src::String                                    = "src"

  webchannels_default_route::String                   = "____"
  webchannels_js_file::String                         = "channels.js"
  webchannels_subscribe_channel::String               = "subscribe"
  webchannels_unsubscribe_channel::String             = "unsubscribe"
  webchannels_autosubscribe::Bool                     = true
  webchannels_eval_command::String                    = ">eval:"
  webchannels_timeout::Int                            = 1_000
  webchannels_keepalive_frequency::Int                = 30_000 # 30 seconds

  webthreads_default_route::String                    = webchannels_default_route
  webthreads_js_file::String                          = "webthreads.js"
  webthreads_pull_route::String                       = "pull"
  webthreads_push_route::String                       = "push"
  webthreads_connection_threshold::Dates.Millisecond  = Dates.Millisecond(60_000) # 1 minute

  html_attributes_replacements::Dict{String,String}   = Dict("v__on!" => "v-on:")
  html_parser_close_tag::String                       = " /"
  html_parser_char_at::String                         = "!!"
  html_parser_char_dot::String                        = "!"
  html_parser_char_column::String                     = "!"
  html_parser_char_dash::String                       = "__"

  ssl_enabled::Bool                                   = false
  ssl_config::Union{MbedTLS.SSLConfig,Nothing}        = nothing

  session_key_name::String                            = "__geniesid"
  session_storage::Union{Symbol,Nothing}              = nothing
  session_options::Dict{String,Any}                   = Dict{String,Any}("Path" => "/", "HttpOnly" => true, "Secure" => ssl_enabled)

  base_path::String                                   = haskey(ENV, "BASEPATH") ? ENV["BASEPATH"] : ""

  features_peerinfo::Bool                             = false

  format_julia_builds::Bool                           = false
  format_html_output::Bool                            = true
  format_html_indentation_string::String              = "  "

  amqp_port::Int                                      = (haskey(ENV, "AMQP_PORT") ? parse(Int, ENV["AMQP_PORT"]) : 5672)
  amqp_username::String                               = haskey(ENV, "AMQP_USERNAME") ? ENV["AMQP_USERNAME"] : "guest"
  amqp_password::String                               = haskey(ENV, "AMQP_PASSWORD") ? ENV["AMQP_PASSWORD"] : "guest"
  amqp_virtualhost::String                            = haskey(ENV, "AMQP_VIRTUALHOST") ? ENV["AMQP_VIRTUALHOST"] : "/"
  amqp_host::String                                   = haskey(ENV, "AMQP_HOST") ? ENV["AMQP_HOST"] : "localhost"
end

end


text:
#=
# 2. Potential flow with an impenetrable surface

This part introduces the treatment of rigid impenetrable surfaces in `GridPotentialFlow`. We can impose the no-penetration constraint by setting the fluid streamfunction equal to that of the surface $\mathfrak{s}_b$, up to a uniform value $s_0$. The discrete no-penetration constraint is thus

$\mathsf{Es}=\mathfrak{s}_{b}-\mathsf{E} \mathsf{s}_{\infty}-\mathfrak{s}_{0},$

where $\mathsf{E}$ is the interpolation operator that interpolates grid data to the surface points, $\mathsf{s}_{\infty}$ is the streamfunction of the uniform flow. For a body translating at velocity $(U,V)$ and rotating at angular velocity $\Omega$, this streamfunction would be

$\mathfrak{s}_{b,k} = U \mathfrak{r}_y - V \mathfrak{r}_x - \frac{1}{2} \Omega (\mathfrak{r}_x^2+\mathfrak{r}_y^2).$

For later shorthand, we will denote the difference between the body motion streamfunction and interpolated uniform flow streamfunction by $\mathfrak{s}'_b \equiv \mathfrak{s}_b - \mathsf{E} \mathsf{s}_{\infty}$. The no-penetration constraint is enforced in the basic potential flow problem with the help of a vector of Lagrange multipliers, $\mathfrak{f}$, on the surface points. The modified potential flow problem is thus

$\mathsf{Ls} + \mathsf{R}\mathfrak{f} = -\mathsf{w},$

where $\mathsf{R}$ is the regularization operator that transfers data from surface data to nodes. $\mathsf{E}$ can be constructed (and we will assume it has) so that it is the transpose of the interpolation operator, $\mathsf{E} = \mathsf{R}^{T}$.

From the previous equation, it is clear (by simple comparison with the vorticity on the right-hand side) that the vector $\mathfrak{f}$ represents the strength of the discrete bound vortex sheet on the surface. Suppose we consider the bound vortex sheet $\gamma(s)$ that emerges from the analogous continuous problem on the undiscretized surface, where $s$ is the arc-length parameter along the surface. At each point $p$, the discrete solution $\mathfrak{f}$ is approximately equal to this continuous solution, multiplied by the length $\delta S_p$ of the small segment surrounding the point:

$\mathfrak{e}_{p}^{T} \mathfrak{f} \approx \gamma(s) \delta S_p.$

Thus, the potential flow problem in the presence of the impenetrable surface is

$\begin{bmatrix}
\mathsf{L} & \mathsf{R} \\
\mathsf{E} &  0
\end{bmatrix} \begin{pmatrix} \mathsf{s} \\ \mathfrak{f} \end{pmatrix} =
\begin{pmatrix} -\mathsf{w} \\ \mathfrak{s}'_b - \mathfrak{s}_0 \end{pmatrix}.$

This problem has the structure of a generic saddle-point problem. `GridPotentialFlow.jl` automatically formulates the no-penetration constraint on the streamfunction and solves the saddle system when it needs to. See the next page for the role of $\mathfrak{s}_0$ and how it will be treated internally.
=#

# We will now consider two examples. In these examples, we use the following grid.
#md # ```@setup 2.-Potential-flow-with-an-impenetrable-surface
#md # using GridPotentialFlow
#md # using Plots
#md # ```
#!md using GridPotentialFlow
#!md using Plots
Δx = 0.02
Lx = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx);

#=
## Vortex near a cylinder

As a basic example, consider now a point vortex near a circular cylinder.
=#

# A body in `GridPotentialFlow.jl` is represented by `PotentialFlowBody`. This is a wrapper around a `Body` from `RigidBodyTools.jl` and can be endowed with a linear velocity, rotational velocity, circulation, and other properties using the appropriate keywords during construction or by mutating its fields using its setter methods. For our circular cylinder, we use the `Circle` shape.
Rc = Lx/4
Δs = 2Δx
body = PotentialFlowBody(Circle(Rc,Δs))

# We choose an initial position $R_v$ of the vortex at $3/2 R_c$ from the origin.
Rv = 3/2*Rc
Γv = 1.0
v = Vortex(Rv,0.0,Γv);

# And we set an equal and opposite circulation about the body
setΓ(body,-Γv)

# We construct the `VortexModel` by using the optional `bodies` keyword. In this case, we are interested in the bound vortex sheet strength $\mathfrak{f}$, instead of just the streamfunction field. For this reason, we use `solve`, which in this case returns a `ConstrainedIBPoissonSolution`.
model = VortexModel(g,vortices=[v],bodies=[body]);
sol = solve(model);

#md # ```@setup 2.-Potential-flow-with-an-impenetrable-surface
#md # # Analytical solution (vortex at Rv and image vortex at Rc^2/Rv)
#md # r1 = sqrt.((body.points.x.-Rv).^2 .+ (body.points.y).^2) # distance from vortex to points on circle
#md # r2 = sqrt.((body.points.x.-Rc^2/Rv).^2 .+ (body.points.y).^2) # distance from image vortex to points on circle
#md # θ1 = π.+atan.((body.points.y)./(body.points.x.-Rv)) # angle from vortex to points on circle
#md # θ2 = atan.((body.points.y),(body.points.x.-Rc^2/Rv)) # angle from image vortex to points on circle
#md # v1x = -v.Γ./(2π*r1).*cos.(θ1.-π/2) # x velocity on circle induced by vortex
#md # v1y = -v.Γ./(2π*r1).*sin.(θ1.-π/2) # y velocity on circle induced by vortex
#md # v2x = v.Γ./(2π*r2).*cos.(θ2.-π/2) # x velocity on circle induced by image vortex
#md # v2y = v.Γ./(2π*r2).*sin.(θ2.-π/2) # y velocity on circle induced by image vortex
#md # V = sqrt.((v1x.+v2x).^2+(v1y.+v2y).^2) # velocity magnitude on circle
#md # γ = -V; # bound vortex sheet strength on circle (velocity on circle is clockwise if positive vortex is to the right of it)
#md # ```

#!md # Analytical solution (vortex at Rv and image vortex at Rc^2/Rv)
#!md r1 = sqrt.((body.points.x.-Rv).^2 .+ (body.points.y).^2) # distance from vortex to points on circle
#!md r2 = sqrt.((body.points.x.-Rc^2/Rv).^2 .+ (body.points.y).^2) # distance from image vortex to points on circle
#!md θ1 = π.+atan.((body.points.y)./(body.points.x.-Rv)) # angle from vortex to points on circle
#!md θ2 = atan.((body.points.y),(body.points.x.-Rc^2/Rv)) # angle from image vortex to points on circle
#!md v1x = -v.Γ./(2π*r1).*cos.(θ1.-π/2) # x velocity on circle induced by vortex
#!md v1y = -v.Γ./(2π*r1).*sin.(θ1.-π/2) # y velocity on circle induced by vortex
#!md v2x = v.Γ./(2π*r2).*cos.(θ2.-π/2) # x velocity on circle induced by image vortex
#!md v2y = v.Γ./(2π*r2).*sin.(θ2.-π/2) # y velocity on circle induced by image vortex
#!md V = sqrt.((v1x.+v2x).^2+(v1y.+v2y).^2) # velocity magnitude on circle
#!md γ = -V; # bound vortex sheet strength on circle (velocity on circle is clockwise if positive vortex is to the right of it)

# We can then easily retrieve the streamfunction and the bound vortex sheet strength from the fields of the solution variable.
plot(sol.ψ,g)
plot!(body,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2)
scatter!([v.x],[v.y],color=:black,markersize=2,xlabel="x",ylabel="y")

# For this example, we can easily calculate the analytical bound vortex sheet strength $\gamma$ and compare it with our numerical solution $f$.
#md # See the notebook in the examples folder for the analytical solution.
plot(sol.f./Δs,label="f/ds",xlabel="body point index")
plot!(γ,label="gamma")

# Finally, we can also check if the vortex will move around the cylinder on a circular path if we advance the system in time, as predicted by analytical potential flow theory. The period of the circular motion can be easily determined by calculating the velocity induced by the image vortex at $R_c^2/R_v$:

Vv = Γv/(2π*(Rv-Rc^2/Rv))
T = 2π*Rv/Vv

# For the time stepping, we again use the fourth-order Runge Kutta scheme from `OrdinaryDiffEq.jl`.
import OrdinaryDiffEq
function rhs(X,model,t)
    setvortexpositions!(model,X)
    Ẋ = vortexvelocities!(model)
    return Ẋ
end
X = getvortexpositions(model)
prob = OrdinaryDiffEq.ODEProblem(rhs,X,(0.0,T),model);
sol = OrdinaryDiffEq.solve(prob,dt=0.1,OrdinaryDiffEq.RK4(),dense=false,adaptive=false);
plot(body,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2)
plot!(map(s->s.u[1],sol.u),map(s->s.v[1],sol.u))
scatter!([v.x],[v.y],color=:black,markersize=2,xlabel="x",ylabel="y")

#jl @testset "Vortex near cylinder" begin
#jl     import OrdinaryDiffEq
#jl
#jl     # analytical solution
#jl     Rv = sqrt(v.x^2+v.y^2)
#jl     Vθ = -v.Γ/(2π*(Rv-Rc^2/Rv))
#jl     Tv = 2π*Rv/abs(Vθ)
#jl
#jl     # numerical solution
#jl     tspan = (0,Tv)
#jl     model = VortexModel(g,bodies=[body],vortices=[v])
#jl
#jl     function rhs(X,model,t)
#jl         setvortexpositions!(model,X)
#jl         Ẋ = vortexvelocities!(model)
#jl         return Ẋ
#jl     end
#jl
#jl     X = getvortexpositions(model)
#jl     prob = OrdinaryDiffEq.ODEProblem(rhs,X,tspan,model);
#jl     sol = OrdinaryDiffEq.solve(prob,dt=0.1,OrdinaryDiffEq.RK4(),dense=false,adaptive=false);
#jl
#jl     @test isapprox(sol.u[end].u[1], 1.5; atol = 1e-2)
#jl     @test isapprox(sol.u[end].v[1], 0.0; atol = 1e-1)
#jl end

#=
## Moving cylinder

In this example we will create a model with a moving body.
=#

# Here we use the keyword `U` to add create a body with a linear velocity.
body = PotentialFlowBody(Circle(Rc,Δs), U=(1.0,0.0))
model = VortexModel(g,bodies=[body]);
sol = solve(model);

#md # ```@setup 2.-Potential-flow-with-an-impenetrable-surface
#md # # Analytical solution
#md # θ = atan.(body.points.y,body.points.x)
#md # γ = 2*sin.(θ);
#md # ```

#!md # # Analytical solution
#!md θ = atan.(body.points.y,body.points.x)
#!md γ = 2*sin.(θ);

plot(sol.ψ,g)
plot!(body,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2,xlabel="x",ylabel="y")

# Again, we can compare the discrete bound vortex sheet strength with the analytical vortex sheet strength.
plot(sol.f./Δs,label="f/ds",xlabel="body point index")
plot!(γ,label="gamma")


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

@with_kw struct MassMatrixAdaptor <: HMCTuningAlgorithm
    target_acceptance::Float64 = 0.8
end

function ahmc_adaptor(
    tuning::MassMatrixAdaptor,
    metric::AdvancedHMC.AbstractMetric,
    integrator::AdvancedHMC.AbstractIntegrator
)
    return AdvancedHMC.MassMatrixAdaptor(metric)
end

@with_kw struct StepSizeAdaptor <: HMCTuningAlgorithm
    target_acceptance::Float64 = 0.8
end

function ahmc_adaptor(
    tuning::StepSizeAdaptor,
    metric::AdvancedHMC.AbstractMetric,
    integrator::AdvancedHMC.AbstractIntegrator
)
    return AdvancedHMC.StepSizeAdaptor(tuning.target_acceptance, integrator)
end

@with_kw struct NaiveHMCTuning <: HMCTuningAlgorithm
    target_acceptance::Float64 = 0.8
end

function ahmc_adaptor(
    tuning::NaiveHMCTuning,
    metric::AdvancedHMC.AbstractMetric,
    integrator::AdvancedHMC.AbstractIntegrator
)
    mma = AdvancedHMC.MassMatrixAdaptor(metric)
    ssa = AdvancedHMC.StepSizeAdaptor(tuning.target_acceptance, integrator)
    return AdvancedHMC.NaiveHMCTuning(mma, ssa)
end

# Uses Stan (also AdvancedHMC) defaults
# (see https://mc-stan.org/docs/2_26/reference-manual/hmc-algorithm-parameters.html):
@with_kw struct StanHMCTuning <: HMCTuningAlgorithm
    "target acceptance rate"
    target_acceptance::Float64 = 0.8

    "width of initial fast adaptation interval"
    initial_bufsize::Int = 75

    "width of final fast adaptation interval"
    term_bufsize::Int = 50

    "initial width of slow adaptation interval"
    window_size::Int = 25
end

function ahmc_adaptor(
    tuning::StanHMCTuning,
    metric::AdvancedHMC.AbstractMetric,
    integrator::AdvancedHMC.AbstractIntegrator,
)
    mma = AdvancedHMC.MassMatrixAdaptor(metric)
    ssa = AdvancedHMC.StepSizeAdaptor(tuning.target_acceptance, integrator)
    stan_adaptor = AdvancedHMC.StanHMCAdaptor(
        mma, ssa,
        init_buffer = Int(tuning.initial_bufsize), term_buffer = Int(tuning.term_bufsize), window_size = Int(tuning.window_size)
    )
    return stan_adaptor
end

mutable struct AHMCTuner{A<:AdvancedHMC.AbstractAdaptor} <: AbstractMCMCTunerInstance
    target_acceptance::Float64
    adaptor::A
end

function (tuning::HMCTuningAlgorithm)(chain::MCMCIterator)
    adaptor = ahmc_adaptor(tuning, chain.hamiltonian.metric, chain.proposal.integrator)
    AHMCTuner(tuning.target_acceptance, adaptor)
end

function tuning_init!(tuner::AHMCTuner, chain::MCMCIterator, max_nsteps::Integer)
    AdvancedHMC.Adaptation.initialize!(tuner.adaptor, Int(max_nsteps - 1))
    nothing
end

tuning_postinit!(tuner::AHMCTuner, chain::MCMCIterator, samples::DensitySampleVector) = nothing

function tuning_reinit!(tuner::AHMCTuner, chain::MCMCIterator, max_nsteps::Integer)
    AdvancedHMC.Adaptation.initialize!(tuner.adaptor, Int(max_nsteps - 1))
    nothing
end

function tuning_update!(tuner::AHMCTuner, chain::MCMCIterator, samples::DensitySampleVector)
    max_log_posterior = maximum(samples.logd)
    accept_ratio = eff_acceptance_ratio(chain)
    if accept_ratio >= 0.9 * tuner.target_acceptance
        chain.info = MCMCIteratorInfo(chain.info, tuned = true)
        @debug "MCMC chain $(chain.info.id) tuned, acceptance ratio = $(Float32(accept_ratio)), integrator = $(chain.proposal.integrator), max. log posterior = $(Float32(max_log_posterior))"
    else
        chain.info = MCMCIteratorInfo(chain.info, tuned = false)
        @debug "MCMC chain $(chain.info.id) *not* tuned, acceptance ratio = $(Float32(accept_ratio)), integrator = $(chain.proposal.integrator), max. log posterior = $(Float32(max_log_posterior))"
    end
    nothing
end

function tuning_finalize!(tuner::AHMCTuner, chain::MCMCIterator)
    adaptor = tuner.adaptor
    AdvancedHMC.finalize!(adaptor)
    chain.hamiltonian = AdvancedHMC.reconstruct(chain.hamiltonian, adaptor)
    chain.proposal = AdvancedHMC.reconstruct(chain.proposal, adaptor)
    nothing
end

tuning_callback(tuner::AHMCTuner) = AHMCTunerCallback(tuner)

struct AHMCTunerCallback{T<:AHMCTuner} <: Function
    tuner::T
end

function (callback::AHMCTunerCallback)(::Val{:mcmc_step}, chain::AHMCIterator)
    adaptor = callback.tuner.adaptor
    tstat = AdvancedHMC.stat(chain.transition)

    AdvancedHMC.adapt!(adaptor, chain.transition.z.θ, tstat.acceptance_rate)
    chain.hamiltonian = AdvancedHMC.reconstruct(chain.hamiltonian, adaptor)
    chain.proposal = AdvancedHMC.reconstruct(chain.proposal, adaptor)
    tstat = merge(tstat, (is_adapt =true,))

    nothing
end


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

name = "Fontconfig"
version = v"2.13.1"

# Collection of sources required to build Fontconfig
sources = [
    "https://www.freedesktop.org/software/fontconfig/release/fontconfig-2.13.1.tar.gz" =>
    "9f0d852b39d75fc655f9f53850eb32555394f36104a044bb2b2fc9e66dbbfa7f",

]

# Bash recipe for building across all platforms
 # weird issue with uuid/uuid.h not being found without C_INCLUDE_PATH
script = raw"""
cd $WORKSPACE/srcdir/fontconfig-2.13.1/
export C_INCLUDE_PATH=$prefix/include:$C_INCLUDE_PATH
./configure --prefix=$prefix --host=$target
make
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libfontconfig", :libfontconfig),
    ExecutableProduct(prefix, "fc-query", :fc_query),
    ExecutableProduct(prefix, "fc-list", :fc_list),
    ExecutableProduct(prefix, "fc-validate", :fc_validate),
    ExecutableProduct(prefix, "fc-cat", :fc_cat),
    ExecutableProduct(prefix, "fc-conflist", :fc_conflist),
    ExecutableProduct(prefix, "fc-match", :fc_match),
    ExecutableProduct(prefix, "fc-scan", :fc_scan),
    ExecutableProduct(prefix, "fc-pattern", :fc_pattern),
    ExecutableProduct(prefix, "fc-cache", :fc_cache)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/bjarthur/ExpatBuilder/releases/download/v2.2.6/build_Expat.v2.2.6.jl",
    "https://github.com/bjarthur/UtilLinuxBuilder/releases/download/v2.32/build_UtilLinux.v2.32.0.jl",
    "https://github.com/JuliaGraphics/FreeTypeBuilder/releases/download/v2.9.1-1/build_FreeType2.v2.9.1.jl",
    "https://github.com/bjarthur/GperfBuilder/releases/download/v3.1.0/build_GperfBuilder.v3.1.0.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: SimonDanisch/FontconfigBuilder

max_issues_repo_name: SimonDanisch/FontconfigBuilder

max_forks_repo_name: SimonDanisch/FontconfigBuilder


text:
using SpikingNeuralNetworks
using Plots
SNN = SpikingNeuralNetworks
SNN.@load_units
using SpikeNetOpt
SNO = SpikeNetOpt
@testset "IZHI_spike_search" begin

    cell_type = "IZHI"
    ngt_spikes = rand(5:15)
    RS = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.02, b = 0.2, c = -65, d = 8))
    current_ = SNO.current_search(cell_type, RS, ngt_spikes)
    RS.I = [current_ * nA]
    SNN.monitor(RS, [:v])
    SNN.monitor(RS, [:fire])
    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms

    SNN.sim!(
        [RS];
        dt = 1 * ms,
        delay = ALLEN_DELAY,
        stimulus_duration = ALLEN_DURATION,
        simulation_duration = ALLEN_DURATION + ALLEN_DELAY + 443ms,
    )

    spikes = SNO.get_spikes(RS)
    spikes = [s * ms for s in spikes]
    nspk = size(spikes)[1]
    @test nspk == ngt_spikes

end


text:
"""
A `ProductGrid` represents the cartesian product of other grids.

`struct ProductGrid{TG,T,N} <: AbstractGrid{T,N}`

Parameters:
- TG is a tuple of (grid) types
- T is the element type of the grid
- N is the dimension of the grid layout
"""
struct ProductGrid{TG,T,N} <: AbstractGrid{T,N}
	grids	::	TG
end

# Generic functions for composite types:
elements(grid::ProductGrid) = grid.grids
element(grid::ProductGrid, j::Int) = grid.grids[j]
element(grid::ProductGrid, range::AbstractRange) = cartesianproduct(grid.grids[range]...)
iscomposite(::ProductGrid) = true

function ProductGrid(grids...)
	TG = typeof(grids)
	T1 = Tuple{map(eltype, grids)...}
	T2 = DomainSets.simplify_product_eltype(T1)
	ProductGrid{typeof(grids),T2,length(grids)}(grids)
end

size(g::ProductGrid) = map(length, g.grids)
size(g::ProductGrid, j::Int) = length(g.grids[j])

support(g::ProductGrid) = cartesianproduct(map(support, elements(g))...)
isperiodic(g::ProductGrid) = reduce(&, map(isperiodic, elements(g)))

getindex(g::ProductGrid{TG,T,N}, I::Vararg{Int,N}) where {TG,T,N} =
	convert(T, map(getindex, g.grids, I))

similargrid(grid::ProductGrid, ::Type{T}, dims...) where T = error()#ProductGrid([similargrid(g, eltype(T), dims[i]) for (i,g) in enumerate(elements(grid))]...)

# Flatten a sequence of elements that may be recursively composite
# For example: a ProductDomain of ProductDomains will yield a list of each of the
# individual domains, like the leafs of a tree structure.
function flatten(::Type{T}, elements::Array, BaseType = Any) where {T}
    flattened = BaseType[]
    for element in elements
        append_flattened!(T, flattened, element)
    end
    flattened
end

flatten(T, elements...) = tuple(flatten(T, [el for el in elements])...)

function append_flattened!(::Type{T}, flattened::Vector, element::T) where {T}
    for el in elements(element)
        append_flattened!(T, flattened, el)
    end
end

function append_flattened!(::Type{T}, flattened::Vector, element) where {T}
    append!(flattened, [element])
end

for (BaseType,TPType) in [ (:AbstractGrid, :ProductGrid)]
    # Override × for grids
    @eval cross(args::$BaseType...) = cartesianproduct(args...)
	@eval ^(arg::$BaseType, n::Int) = cartesianproduct(arg, n)
    # In order to avoid strange nested structures, we flatten the arguments
    @eval cartesianproduct(args::$BaseType...) = $TPType(flatten($TPType, args...)...)
    @eval cartesianproduct(arg::$BaseType, n::Int) = cartesianproduct(ntuple(k->arg,Val(n)))
    # Disallow cartesian products with just one argument
    @eval cartesianproduct(arg::$BaseType) = arg
end


text:
using Images, TestImages, HistogramThresholding

function gaussian(z;P₁,μ₁,σ₁)
    (P₁ / (σ₁*√(2*π))) * exp( -(z-μ₁)^2 / (2*σ₁^2) )
end

function surge(t;A,k)
    A*t*exp(-k*t)
end

# FInd threshold for moon surface image
img = testimage("moonsurface")
edges, counts = build_histogram(img,256)
index = find_threshold(UniThreshold(),counts[1:end], edges)
t = edges[index]

#= FInd threshold for synthetic data
edges = 0:255
counts = collect(255:-1:0)
counts[2] = 1
counts[100] = 10
counts[150] = 0
index = find_threshold(UniThreshold(),counts, edges)
=#

# Making bin after max very low
edges = 0:255
samples = collect(255:-1:0)
samples[2] = 1
samples
t = find_threshold(UniThreshold(),samples,edges)

# Surge Function (also technically a test for no 0 bins as well)
edges = 0:255
samples = map(x->surge(x, A = 20 ,k=0.04),edges)
t = find_threshold(UniThreshold(),samples,edges)

# multiple zero bins

samples = fill(10,256)
samples[5] = 5
samples[100] = 0
samples[120] = 1
samples[150] = 0
t = find_threshold(UniThreshold(),samples,edges)

# Different number of bins

edges = 0:200
samples = collect(200:-1:0)
samples[2] = 1
t = find_threshold(UniThreshold(),samples,edges)

# binarises the image
img_binary = map(img) do val
    if val < t
        return Gray(0)
    else
        return Gray(1)
    end
end
# displays it
img_binary


text:
using CSV
using LinearAlgebra: dot, I
using CairoMakie

const Σ = [1.0  0.1  0.1;
           0.1  1.0  0.1;
           0.1  0.1  1.0]

const v = [1.0/sqrt(3), 1.0/sqrt(3), 1.0/sqrt(3)]
const σ² = 1

const Σ2 = -0.95 * v * v' + I(3)

f(x) = abs2(dot(v, Σ * x)) / (dot(x, Σ * x) + σ² * dot(x, x))

g(x) = abs2(dot(v, Σ2 * x)) / (dot(x, Σ2 * x) + σ² * dot(x, x))

# Helper function to compute extra coordinate
function getZ(x,y)
    if abs2(x) + abs2(y) > 1 - 10*eps()
        return NaN
    end
    return sqrt(1 - abs2(x) - abs2(y))
end

function rewardPositive(x,y)
    return f([x,y,getZ(x,y)])
end

function rewardNegative(x,y)
    return f([x,y,-getZ(x,y)])
end

function rewardPositive2(x,y)
    return g([x,y,getZ(x,y)])
end

function rewardNegative2(x,y)
    return g([x,y,-getZ(x,y)])
end

# Sample Circle
x = -1.1:0.01:1.1
y = -1.1:0.01:1.1

noto_sans = "./resources/NotoSans-Regular.ttf"
noto_sans_bold = "./resources/NotoSans-Bold.ttf"

tickfontsize    = 26
labelfontsize   = 28
basewidth = 2

fig = Figure(font=noto_sans, resolution=(1500,525), figure_padding=40)

ax1 = Axis( fig[1,1],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="x", ylabel = "y",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
            titlesize=labelfontsize,
            ratio=AxisAspect(1))

hm = heatmap!(ax1,x,y,rewardPositive)

ax2 = Axis( fig[1,2],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="x", ylabel = "y",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
            titlesize=labelfontsize,
            ratio=AxisAspect(1))

heatmap!(ax2,x,y,rewardNegative)

ax3 = Axis( fig[1,3],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="x", ylabel = "y",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
            titlesize=labelfontsize,
            ratio=AxisAspect(1))

hm2 = heatmap!(ax3,x,y,rewardPositive2)

ax4 = Axis( fig[1,4],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="x", ylabel = "y",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
            titlesize=labelfontsize,
            ratio=AxisAspect(1))

heatmap!(ax4,x,y,rewardNegative2)

Colorbar(fig[2, 1:2], hm, vertical = false, flipaxis = false, ticklabelsize=tickfontsize)
Colorbar(fig[2, 3:4], hm2, vertical = false, flipaxis = false, ticklabelsize=tickfontsize)

colgap!(fig.layout,15)
colgap!(fig.layout,2,100)

Label(fig[1,1,TopLeft()], "A", font=noto_sans_bold, textsize = 36, halign=:left, valign=:bottom, padding=(25,0,10,0))
Label(fig[1,3,TopLeft()], "B", font=noto_sans_bold, textsize = 36, halign=:left, valign=:bottom, padding=(25,0,10,0))

save("out/OptSurface.pdf",fig)


text:
# f(x,y)=x+y'

FXPYtranspose(x::Float64,y::Float64)=Fxpy(x,y)
FXPYtranspose(x::Float64,y)=Fxpy(x,y')

FXPYtranspose(x,y::Float64)=Fxpy(x,y)

FXPYtranspose(x::Array,y::Array)=Fxpy(x,y')

FXPYtranspose_inplace(value,auxvalue,x::Array,y::Array)=Fxpy_inplace(value,auxvalue,x,y')

function DXPYtranspose(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::Array,y::Array)
    if derivativeIDX==1
        if size(x)==(1,1)
            axpy!(1.0,[sum(grad_c)],grad_n)
        else
            axpy!(1.0,grad_c,grad_n)
        end
    elseif derivativeIDX==2
        if size(y)==(1,1)
            axpy!(1.0,[sum(grad_c)],grad_n)
        else
            axpy!(1.0,grad_c',grad_n)
        end
    end
end

if PROC=="GPU"

FXPYtranspose_inplace(value,auxvalue,x::CudaArray,y::CudaArray)=FXtransposePY_inplace(value,auxvalue,y,x)

function DXPYtranspose(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray,y::CudaArray)
    if derivativeIDX==1
        if size(x)==(1,1)
            axpy!(1.0,sum(grad_c),grad_n)
        else
            axpy!(1.0,grad_c,grad_n)
        end
    elseif derivativeIDX==2
        if size(y)==(1,1)
            axpy!(1.0,sum(grad_c),grad_n)
        else
            tmp=CudaArray(Float64,size(grad_c))
            CUBLAS.geam!('T','N',1.0,grad_c,0.0,tmp,tmp) # tmp=grad_c'
            axpy!(1.0,tmp,grad_n); free(tmp)
        end
    end
end

end

Derivative[FXPYtranspose]=DXPYtranspose
Inplace[FXPYtranspose]=FXPYtranspose_inplace

export FXPYtranspose

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using BinaryProvider, Compat
using CMakeWrapper: cmake_executable
using Compat.Libdl: dlext

function compile(libname, tarball_url, hash; prefix=BinaryProvider.global_prefix, verbose=false)
    # download to tarball_path
    tarball_path = joinpath(prefix, "downloads", "src.tar.gz")
    download_verify(tarball_url, hash, tarball_path; force=true, verbose=verbose)

    # unpack into source_path
    tarball_dir = joinpath(prefix, "downloads", dirname(first(list_tarball_files(tarball_path)))) # e.g. "c-blosc-1.14.3"
    source_path = joinpath(prefix, "downloads", "src")
    verbose && Compat.@info("Unpacking $tarball_path into $source_path")
    rm(tarball_dir, force=true, recursive=true)
    rm(source_path, force=true, recursive=true)
    unpack(tarball_path, dirname(tarball_dir); verbose=verbose)
    mv(tarball_dir, source_path)

    build_dir = joinpath(source_path, "build")
    mkdir(build_dir)
    verbose && Compat.@info("Compiling in $build_dir...")
    cd(build_dir) do
        run(`$cmake_executable -DBUILD_TESTS=Off -DBUILD_BENCHMARKS=Off ..`)
        run(`$cmake_executable --build .`)
        mkpath(libdir(prefix))
        Compat.cp("blosc/libblosc.$dlext", joinpath(libdir(prefix), libname*"."*dlext),
            force=true, follow_symlinks=true)
    end
end


text:
module FastConv

# package code goes here
using Base.Cartesian
import DSP: conv, conv2

include("./utils.jl");

end # module