text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# This file depends on base/regex.jl and base/show.jl
# and needs to follow these in base/sysimg.jl

"""
    Module to implement replace method with multiple pairs of patterns-replacements
"""
module ReplaceImpl

export replace_code
import Base: Callable, SubstitutionString

#special type to wrap a function f, which has a method f(::Char)
struct CharMethod{T<:Function} <: Function
    f::T
end
(cm::CharMethod)(x::Union{Char,AbstractString}) = cm.f(x)

# types of input patterns
const TIF = Union{Char,AbstractString,Regex,Callable,Tuple{Vararg{Char}},
                  AbstractVector{Char},Set{Char}}
# types of input replacements
#const TIS = Union{Char,AbstractString,Callable,SubstitutionString}
const TIS = Any

# types input patterns are converted to
const TF = Union{Regex,String,Callable}
# types output patterns are converted to
const TS = Union{Char,String,Callable,SubstitutionString,CharMethod}

# extract pattern from a Pair
pattern(p::Pair{<:TF}) = p.first
# extract replacement from a Pair
repl(p::Pair{<:TF,CharMethod{T}}) where T<:Function = p.second.f
repl(p::Pair{<:TF,<:TS}) = p.second

# type of function, if wrapped in CharMethod
charmethod(x::Type{CharMethod{T}}) where T<:Function = T
charmethod(x) = x

# code generation for all input patterns
function gen_pattern_all(ptypes::NTuple{N,DataType}) where N
    Expr(:block, [gen_pattern(i, ptypes[i]) for i in 1:N]...)
end

# generate code to find next substring complying Regex or identical to String
function gen_pattern(pri::Int, ::Type{<:Pair{Regex}})
    :(  let j = $(getpr(1, pri)), k = $(getpr(2, pri)), r;
        if j < MAX
            if j < pos
                r = findnext($(getpa(pri)), str, pos)
                j, k = first(r), last(r)
                j = ifelse(j == 0, MAX, j)
                $(setpr(1, pri, :j)); $(setpr(2, pri, :k))
            end
            $(gen_min_max(pri))
        end
    end
    )
end

function gen_pattern(pri::Int, ::Type{<:Pair{String}})
    :(  let j = $(getpr(1, pri)), k = $(getpr(2, pri));
        if j < MAX
            if j < pos
                k = endof($(getpa(pri)))
                if k > 0
                    j = Base._searchindex(cstr, $(getpavec(pri)), pos)
                    if j == 0
                        j = k = MAX
                    else
                        k += j - 1
                    end
                else
                    j = pos; k = pos - 1
                end
                $(setpr(1, pri, :j)); $(setpr(2, pri, :k))
            end
            $(gen_min_max(pri))
        end
    end
    )
end

# generate code to find next character honoring predicate function
function gen_pattern(pri::Int, ::Type{<:Pair{<:Function}})
    :(  let j = $(getpr(1, pri)), k;
        if j < MAX && pos != maxk
            if j < pos
                r = findnext($(getpa(pri)), str, pos)
                j = r isa Nothing ? MAX : r
                $(setpr(1, pri, :j)); $(setpr(2, pri, :j))
            end
            k = j
            $(gen_min_max(pri))
        end
    end
    )
end

# generate code to calculate minj, maxk, and minp
function gen_min_max(pri::Int)
    :( (j < minj || j == minj && k > maxk) && ((minj, maxk, minp) = (j, k, $pri)) )
end

# generate code to select one of several options, depending on runtime variable `minp`.
function extree(n::Int, m::Int, ex::Vector)
    if n < m
        k = (n + m) ÷ 2
        #Expr(:if, :(minp <= $k), extree(n, k, ex), extree(k+1, m, ex))
        :( minp <= $k ? $(extree(n, k, ex)) : $(extree(k+1, m, ex)) )
    elseif n == m
        ex[n]
    else
        nothing
    end
end

# generate code to access variables
getpr(k::Int, pri::Int) = Symbol("pr", k, "_", pri)
setpr(k::Int, pri::Int, j) = Expr(:(=), getpr(k, pri), j)
setprall(k::Int, n::Int) = [setpr(k, i, 0) for i = 1:n]

getpa(pri::Int) = Symbol("pa_", pri)
setpa(pri::Int, T::Type, j) = Expr(:(=), Expr(:(::), getpa(pri), T), j)
function setpaall(pt::NTuple{N,DataType}) where N
    [setpa(i, pt[i].types[1], :(pattern(pat_repls[$i]))) for i = 1:N]
end

getpavec(pri::Int) = Symbol("pavec_", pri)
setpavec(pri::Int, T::Type, j) = Expr(:(=), Expr(:(::), getpavec(pri), T), j)
function setpavecall(pt::NTuple{N,DataType}) where N
    ex = Expr[]
    for i = 1:N
        if pt[i].types[1] == String
            push!(ex, setpavec(i, Vector{UInt8},
                               :(unsafe_wrap(Vector{UInt8}, $(getpa(i))))))
        end
    end
    ex
end

getre(pri::Int) = Symbol("re_", pri)
setre(pri::Int, T::Type, j) = Expr(:(=), Expr(:(::), getre(pri), T), j)
function setreall(pt::NTuple{N,DataType}) where N
    [setre(i, charmethod(pt[i].types[2]), :(repl(pat_repls[$i]))) for i = 1:N]
end

function gen_repl_all(ptypes::NTuple{N,DataType}) where N
    extree(1, N, [gen_repl(i, ptypes[i]) for i in 1:N])
end

# code generation depending on pattern_type
function gen_repl(pri::Int, ::Type{<:Pair{<:Any, Char}})
    :( write(out, $(getre(pri))))
end

function gen_repl(pri::Int, ::Type{<:Pair{<:Any, String}})
    :( print(out, $(getre(pri))))
end

# call the output methods for (Char,)
function gen_repl(pri::Int, pt::Type{<:Pair{<:Callable, <:CharMethod}})
    :( write(out, $(getre(pri))(getindex(str, minj))))
end

function gen_repl(pri::Int, pt::Type{<:Pair{<:Any, <:CharMethod}})
    :( minj == maxk ?
        write(out, $(getre(pri))(getindex(str, minj))) :
        write(out, $(getre(pri))(getindex(str, minj:maxk))) )
end

function gen_repl(pri::Int, pt::Type{<:Pair{<:Callable, <:Callable}})
    :( write(out, $(getre(pri))(getindex(str, minj:minj))))
end

function gen_repl(pri::Int, pt::Type{<:Pair{<:Any, <:Callable}})
    :( write(out, $(getre(pri))(getindex(str, minj:maxk))))
end

function gen_repl(pri::Int, ::Type{<:Pair{Regex, <:SubstitutionString}})
    :(Base._replace(out, $(getre(pri)), str, minj:maxk, $(getpa(pri))))
end

# all arguments are types of the real arguments
function replace_gen_impl(str::Type{String}, count::Type{Int}, pat_repls...)
    n = length(pat_repls)
    gen_input = gen_pattern_all(pat_repls)
    gen_output = gen_repl_all(pat_repls)

    quote
        begin
            # The following local variables occur once per pair
            $(setpaall(pat_repls)...)       # patterns
            $(setpavecall(pat_repls)...)    # UInt8-vector forms of string pattern
            $(setreall(pat_repls)...)       # replacements
            $(setprall(1, n)...)            # next start of match
            $(setprall(2, n)...)            # next end of match
            cstr = unsafe_wrap(Vector{UInt8}, str)  # vector form of input string
            MAX = typemax(Int)
            eos = endof(str)
            pos = 1  # start position for next findnext
            wpos = 1 # start of write out in original string
            # buffer size hint greater than original size - empirical value
            out = IOBuffer(Base.StringVector(eos*12÷10), true, true)
            out.size = 0
            ctr = count
            while ctr > 0

                minj = MAX  # first of next found pattern (MAX if nothing found)
                maxk = -1   # end of next found pattern
                minp = 0    # type of repl for found pattern
                $gen_input  # calculates minj, maxk, minp
                minj == MAX && break    # no further pattern found
                if wpos == 1 || wpos <= maxk
                    unsafe_write(out, pointer(str, wpos), UInt(minj-wpos))
                    $gen_output
                end
                if maxk < minj   # empty pattern found
                    wpos = minj  # write this position
                    minj > eos && break
                    pos = nextind(str, minj) # but skip searching next pattern
                else
                    pos = wpos = nextind(str, maxk) # write and search from after patt
                end
                ctr -= 1
            end
            write(out, SubString(str, wpos))
            String(take!(out))
        end
    end
end

@generated function replace_gen(str::String, count::Int, pat_repls::Pair...)
    replace_gen_impl(str, count, pat_repls...)
end

# consolidate patterns
# collections of chars are transformed to Char if collection has only one member
# strings of length 1 are transformed to the character
# single char x is transformed to `equalto(x)`
# abstract strings are converted to strings
# character collection pat is transformed to `occursin(pat)`
# `occursin(character collection)` is transformed like occursin().x
# `occursin(empty string)` is ignored
# Regex and Callables are unchanged
# other patterns are stringified (TODO should be rejected!)
#
function consol_pattern(pat::Union{Tuple{Vararg{Char}},AbstractVector{Char},Set{Char}})
    if length(pat) > 1
        occursin(pat)
    elseif length(pat) == 1
        equalto(next(pat, start(pat))[1])
        # works also for `Set{Char}` as opposed to `first`
    else
        isempty # case of empty collection - `isempty(::Char) == false`
    end
end
function consol_pattern(pat::Base.OccursIn{<:Union{Tuple{Vararg{Char}},
                                                   AbstractVector{Char},Set{Char}}})
    consol_pattern(pat.x)
end
function consol_pattern(pat::Base.OccursIn{<:AbstractString})
    isempty(pat.x) ? isempty : pat
end
consol_pattern(pat::Char) = equalto(pat)
consol_pattern(pat::AbstractString) = length(pat) == 1 ? equalto(pat[1]) : string(pat)
consol_pattern(pat::Union{Regex,Callable}) = pat
consol_pattern(pat) = string(pat)

# consolidate replacements
# SubstitutionString, Char, and functions are unchanged
# strings of length 1 are transformed to the character
# abstract strings are converted to strings
# other replacements are stringified (TODO should be rejected?)
#
consol_repl(repl::Function) = hasmethod(repl, (Char,)) ? CharMethod(repl) : repl
consol_repl(repl::Union{Char,SubstitutionString,Callable}) = repl
consol_repl(repl::AbstractString) = length(repl) == 1 ? repl[1] : string(repl)
consol_repl(repl) = string(repl)

# check if pattern is equal to replacement (no-op replacement)
is_modify(p::Pair{Base.EqualTo{Char},Char}) = p.first.x != p.second
is_modify(p::Pair{typeof(isempty)}) = false
is_modify(p::Pair) = p.first != p.second

# replace collection of Char with function that checks occurrence in this collection
function pat_repl_pair(p::Pair{<:TIF,<:TIS})::Pair{<:TF,<:TS}
    pat::TF  = consol_pattern(p.first)
    repl::TS = consol_repl(p.second)
    repl isa SubstitutionString && !(pat isa Regex) &&
        error("substitution string requires regex")
    repl = funtostring(pat, repl)
    pat === p.first && repl === p.second ? p : pat => repl
end

funtostring(p::Pair) = Pair(p.first, funtostring(p.first, p.second))

funtostring(pat::Base.EqualTo{Char}, fun::Callable) = stringfun(fun, pat.x)
funtostring(pat::Union{Char,String}, fun::Callable) = stringfun(fun, pat)
funtostring(::Any, y) = y

function stringfun(fun::Callable, s::Union{Char,String})
    if hasmethod(fun, (Char,)) && length(s) == 1
        fun(s[1])
    else
        fun(string(s))
    end
end

function check_args(str::String, count::Int, pat_repls::Pair...)
    count < 0 && throw(DomainError(count, "`count` must be non-negative."))
    count == 0 || length(pat_repls) == 0
end

function consolidate_args(pat_repls::Pair...)
    pairs = collect(Pair, pat_repl_pair.(pat_repls))
    pairs = expand_pairs(pairs)
    pairs = funtostring.(pairs)
    filter!(is_modify, pairs)
    condense_pairs(pairs)
end

"""
    replace_code(string, pattern-pairs)

Return `Expr` of the body of the generated function
"""
function replace_code(str::String, pat_repls::Pair...; count::Int=typemax(Int))
    check_args(str, count, pat_repls...) && return :str
    pat_repls2 = consolidate_args(pat_repls...)
    length(pat_repls2) == 0 && return :str
    ReplaceImpl.replace_gen_impl(typeof(str), Int, typeof.(pat_repls2)...)
end

include("condense.jl")

end # module ReplaceImpl

import Base.replace

# cover the multiple pairs case
function replace(str::AbstractString, pat_repls::Pair...; count::Integer=typemax(Int))
    count = Int(clamp(count, typemin(Int), typemax(Int)))
    ReplaceImpl.check_args(str, count, pat_repls...) && return str
    pat_repls2 = ReplaceImpl.consolidate_args(pat_repls...)
    length(pat_repls2) == 0 && return str
    ReplaceImpl.replace_gen(String(str), count, pat_repls2...)
end


text:
"""
    UpdateCamera(camera::RayCamera3D)

Return new camera with updated parameter.
"""
function Binding.UpdateCamera(camera::RayCamera3D)
    new_camera_ref = Ref(camera)
    UpdateCamera(new_camera_ref)

    return new_camera_ref[]
end

"""
    UpdateCamera!(camera::RayCamera3D)

Update camera position for selected mode
"""
function UpdateCamera!(camera::RayCamera3D)
    camera_ptr = convert(Ptr{RayCamera3D}, pointer_from_objref(camera))
    UpdateCamera(camera_ptr)
    return camera
end

"""
    GetDroppedFiles()

Return a list of dropped file paths.
"""
function Binding.GetDroppedFiles()
    count = Ref{Cint}(0)
    fptrs = GetDroppedFiles(count)
    fcstr = Base.unsafe_wrap(Vector{Cstring}, fptrs, count[])
    return map(Base.unsafe_string, fcstr)
end

"""
    RayFileData(filename::AbstractString)

Read the file data. It will be auto-unloaded when being garbage collected.
"""
struct RayFileData
    data::Vector{UInt8}

    function RayFileData(filename::AbstractString)
        bytes = Ref{Cuint}(0)
        dptr = LoadFileData(filename, bytes)
        fdata = Base.unsafe_wrap(Vector{UInt8}, dptr, bytes[])
        finalizer(fdata) do x
            ptr = pointer(x)
            UnloadFileData(ptr)
        end

        return new(fdata)
    end
end

Base.length(fdata::RayFileData) = length(fdata.data)
Base.pointer(fdata::RayFileData) = pointer(fdata.data)

for func in :(
    IsKeyDown, IsKeyPressed, IsKeyReleased, IsKeyUp
).args
    @eval Binding.$func(k::KeyboardKey) = $func(convert(Integer, k))
end

ext: jl


text:
using FizzBuzz
using Base.Test

const fizzbuzzer = FizzBuzz.fizzbuzzer
const d = Dict(3 => "fizz", 5 => "buzz", 15 => "fizz buzz")

@testset "FizzBuzz" begin

    @testset "Test <:Integer returns $(d[v])" for v = [3, 5, 15]
        @test fizzbuzzer(v) == d[v]
    end

    @testset "Test <:AbstractFloat returns $(d[v])" for v = [3., 5., 15.]
        @test fizzbuzzer(v) == d[v]
    end

    @testset "Test <:Rational returns $(d[v])" for v = [3//1, 5//1, 15//1]
        @test fizzbuzzer(v) == d[v]
    end

    @testset "Test $v returns $v" for v = (1, 1., 1//1)
        @test fizzbuzzer(v) == v
    end

    @testset "Test UnitRange{<:Integer}" begin
        io = IOBuffer()
        fizzbuzz(1:15, io)
        @test String(take!(io)) == "1\n2\nfizz\n4\nbuzz\nfizz\n7\n8\nfizz\nbuzz\n11\nfizz\n13\n14\nfizz buzz\n"
        close(io)
    end

    @testset "Test UnitRange{<:AbstractFloat}" begin
        io = IOBuffer()
        fizzbuzz(1.:15., io)
        @test String(take!(io)) == "1.0\n2.0\nfizz\n4.0\nbuzz\nfizz\n7.0\n8.0\nfizz\nbuzz\n11.0\nfizz\n13.0\n14.0\nfizz buzz\n"
        close(io)
    end

    @testset "Test Array{<:Integer}" begin
        io = IOBuffer()
        fizzbuzz(collect(1:15), io)
        @test String(take!(io)) == "1\n2\nfizz\n4\nbuzz\nfizz\n7\n8\nfizz\nbuzz\n11\nfizz\n13\n14\nfizz buzz\n"
    end

    @testset "Test Array{<:AbstractFloat}" begin
        io = IOBuffer()
        fizzbuzz(collect(1.:15.), io)
        @test String(take!(io)) == "1.0\n2.0\nfizz\n4.0\nbuzz\nfizz\n7.0\n8.0\nfizz\nbuzz\n11.0\nfizz\n13.0\n14.0\nfizz buzz\n"
    end

    @testset "Test not divisible by 3 or 5 returns number" begin
        for v = 1:100
            while true
                r = rand(1:100_000, 1)[1]
                if r % 3 > 0 && r % 5 > 0
                    @test fizzbuzzer(r) == r
                    break
                end
            end
        end
    end

end


text:
export azimuth, final_azimuth, spherical_angle, spherical_excess

"""
    azimuth(point₁::Point, point₂::Point)

    output: angle [deg]

Return the initial `azimuth` [deg], measured clockwise from the north direction,
of the great circle line between the positions `point₁` and `point₂` [deg] on a
unit sphere.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function azimuth(point₁::Point, point₂::Point)
    Δpos = point₂ - point₁
    return mod(atand(sind(Δpos.λ)*cosd(point₂.ϕ), cosd(point₁.ϕ)*sind(point₂.ϕ) -
    sind(point₁.ϕ)*cosd(point₂.ϕ)*cosd(Δpos.λ)), 360.0)
end

"""
    azimuth(line::Line)

    output: angle [deg]

Return the initial `azimuth` [deg], measured clockwise from the north direction,
of the line on a unit sphere.
"""
azimuth(arc::Arc) = azimuth(arc.point₁, arc.point₂)

"""
    final_azimuth(point₁::Point, point₂::Point)

    output: angle [deg]

Return the final `azimuth` [deg], measured clockwise from the north direction,
of the great circle line between the positions `point₁` and `point₂` [deg] on a
unit sphere.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function final_azimuth(point₁::Point, point₂::Point)
    return mod(azimuth(point₂, point₁) + 180.0, 360.0)
end

"""
    final_azimuth(line::Line)

    output: angle [deg]

Return the `final_azimuth` [deg] of the line on a unit sphere.
"""
final_azimuth(arc::Arc) = final_azimuth(arc.point₁, arc.point₂)

"""
    spherical_angle(angular_distance₁₂::Float64, angular_distance₂₃::Float64, angular_distance₁₃::Float64)

    output: angle [deg]

Return the 'spherical_angle' [deg] of angle₁₂₃ given the angular distances [deg] between point₁ and point₂,
point₂ and point₃, and point₁ and point₃.

Source: en.wikipedia.org/wiki/Spherical_trigonometry
"""
function spherical_angle(angular_distance₁₂::Float64, angular_distance₂₃::Float64, angular_distance₁₃::Float64,
    tolerance::Float64=tolerance_deg)
    a = deg2rad(angular_distance₁₃)
    b = deg2rad(angular_distance₂₃)
    c = deg2rad(angular_distance₁₂)
    if b < deg2rad(tolerance) || c < deg2rad(tolerance)
        return acos(1 - deg2rad(a)^2/2)
    end
    cosA = (cos(a) - cos(b)*cos(c)) / (sin(b)*sin(c))
    return rad2deg(acos(cosA))
end

"""
    spherical_angle(point₁::Point, point₂::Point, point₃::Point)

    output: angle [deg]

Return the 'spherical_angle' [deg] of angle₁₂₃ given the triangle point₁ - point₂ - point₃

Source: en.wikipedia.org/wiki/Spherical_trigonometry
"""
spherical_angle(point₁::Point, point₂::Point, point₃::Point, tolerance::Float64=tolerance_deg) =
spherical_angle(angular_distance(point₁, point₂), angular_distance(point₂, point₃), angular_distance(point₁, point₃))

"""
    spherical_excess(angular_distance₁₂::Float64, angular_distance₂₃::Float64, angular_distance₁₃::Float64)

    output: angle [deg]

Return the 'spherical_excess' [deg] of the spherical triangle 123 given the angular distances between the
points point₁, point₂, and point₃

Source: mathworld.wolfram.com/SphericalExcess.html
"""
function spherical_excess(angular_distance₁₂::Float64, angular_distance₂₃::Float64,
    angular_distance₁₃::Float64, tolerance::Float64=tolerance_deg)
    A = spherical_angle(angular_distance₁₂, angular_distance₂₃, angular_distance₁₃)
    B = spherical_angle(angular_distance₁₃, angular_distance₁₂, angular_distance₂₃)
    C = spherical_angle(angular_distance₂₃, angular_distance₁₃, angular_distance₁₂)
    return A+B+C-180.0
end

"""
    spherical_excess(point₁::Point, point₂::Point, point₃::Point)

    output: angle [deg]

Return the 'spherical_excess' [deg] of the spherical triangle 123 given the points point₁ - point₂ - point₃

Source: mathworld.wolfram.com/SphericalExcess.html
"""
spherical_excess(point₁::Point, point₂::Point, point₃::Point) =
spherical_excess(angular_distance(point₁, point₂), angular_distance(point₂, point₃),
angular_distance(point₁, point₃))


text:
# in file 'myscript.jl'
push!( LOAD_PATH, "./" )
import MyPrecompiledModule
println( "Hello from the script. The arguments passed into it were $ARGS" )
MyPrecompiledModule.exportedfun()

max_stars_repo_name: timm/misc

max_issues_repo_name: timm/misc

max_forks_repo_name: timm/misc


text:
@inline nvalues(v::CS.Variable) = v.last_ptr-v.first_ptr+1

"""
    value(v::CS.Variable)

Get the value of the variable if it is fixed. Otherwise one of the possible values is returned.
Can be used if the status is :Solved as then all variables are fixed.
"""
@inline value(v::CS.Variable) = v.values[v.last_ptr]

@inline values(v::CS.Variable) = v.values[v.first_ptr:v.last_ptr]

"""
    values(m::Model, v::VariableRef)

Return all possible values for the variable. (Only one if solved to optimality)
"""
function values(m::Model, v::VariableRef)
    com = backend(m).optimizer.model.inner
    return values(com.search_space[v.index.value])
end

@inline view_values(v::CS.Variable) = @views v.values[v.first_ptr:v.last_ptr]

function issetto(v::CS.Variable, x::Int)
    if !isfixed(v)
        return false
    else
        return x == value(v)
    end
end

function has(v::CS.Variable, x::Int)
    if x > v.max || x < v.min
        return false
    end
    ind = v.indices[x+v.offset]
    return v.first_ptr <= ind <= v.last_ptr
end

function rm!(com::CS.CoM, v::CS.Variable, x::Int; in_remove_several=false, changes=true)
    if !in_remove_several
        # after removing nothing would be possible
        len_vals = nvalues(v)
        if len_vals == 1
            com.bt_infeasible[v.idx] += 1
            return false
        elseif len_vals == 2
            possible = values(v)
            left_over = possible[1] == x ? possible[2] : possible[1]
            if !fulfills_constraints(com, v.idx, left_over)
                com.bt_infeasible[v.idx] += 1
                return false
            end
        end
    end

    ind = v.indices[x+v.offset]
    v.indices[x+v.offset], v.indices[v.values[v.last_ptr]+v.offset] = v.indices[v.values[v.last_ptr]+v.offset], v.indices[x+v.offset]
    v.values[ind], v.values[v.last_ptr] = v.values[v.last_ptr], v.values[ind]
    v.last_ptr -= 1
    if !in_remove_several
        vals = values(v)
        if length(vals) > 0
            if x == v.min
                v.min = minimum(vals)
            end
            if x == v.max
                v.max = maximum(vals)
            end
        end
        changes && push!(v.changes[com.c_backtrack_idx], (:rm, x, 0, 1))
    end
    return true
end

function fix!(com::CS.CoM, v::CS.Variable, x::Int; changes=true)
    if !fulfills_constraints(com, v.idx, x)
        com.bt_infeasible[v.idx] += 1
        return false
    end
    ind = v.indices[x+v.offset]
    pr_below = ind-v.first_ptr
    pr_above = v.last_ptr-ind
    changes && push!(v.changes[com.c_backtrack_idx], (:fix, x, v.last_ptr, 0))
    v.last_ptr = ind
    v.first_ptr = ind
    v.min = x
    v.max = x
    return true
end

function isfixed(v::CS.Variable)
    return v.last_ptr == v.first_ptr
end

function remove_below!(com::CS.CoM, var::CS.Variable, val::Int; changes=true)
    vals = values(var)
    still_possible = filter(v -> v >= val, vals)
    if length(still_possible) == 0
        com.bt_infeasible[var.idx] += 1
        return false
    elseif length(still_possible) == 1
        if !fulfills_constraints(com, var.idx, still_possible[1])
            com.bt_infeasible[var.idx] += 1
            return false
        end
    end

    nremoved = 0
    for v in vals
        if v < val
            rm!(com, var, v; in_remove_several = true)
            nremoved += 1
        end
    end
    if nremoved > 0 && feasible(var)
        var.min = minimum(values(var))
        changes && push!(var.changes[com.c_backtrack_idx], (:remove_below, val, 0, nremoved))
    end
    return true
end

function remove_above!(com::CS.CoM, var::CS.Variable, val::Int; changes=true)
    vals = values(var)
    still_possible = filter(v -> v <= val, vals)
    if length(still_possible) == 0
        com.bt_infeasible[var.idx] += 1
        return false
    elseif length(still_possible) == 1
        if !fulfills_constraints(com, var.idx, still_possible[1])
            com.bt_infeasible[var.idx] += 1
            return false
        end
    end

    nremoved = 0
    for v in vals
        if v > val
            rm!(com, var, v; in_remove_several = true)
            nremoved += 1
        end
    end
    if nremoved > 0 && feasible(var)
        var.max = maximum(values(var))
        changes && push!(var.changes[com.c_backtrack_idx], (:remove_above, val, 0, nremoved))
    end
    return true
end

function feasible(var::CS.Variable)
    return var.last_ptr >= var.first_ptr
end


text:
module DiffSolving

export EulerMethod, InstableAlgorithm, EulerCromer, Verlet, VelocityVerlet

using LaTeXStrings, Plots, ProgressBars, Statistics, StatsBase
cd(dirname(@__FILE__))
#First Order ODE_Comparison between Euler method's and Numerical solution

function EulerMethod(func::Function, xᵢₙᵢₜ::Vector{Float64}, tₘᵢₙ::Float64, tₘₐₓ::Float64, h::Float64)
    t_range=collect(range(tₘᵢₙ, tₘₐₓ, step=h))
    Q=zeros(length(t_range), length(xᵢₙᵢₜ))
    Q[1, :]=xᵢₙᵢₜ

    for (i,t) in enumerate(t_range[1:end-1])
        Q[i+1, :] = Q[i, :] + func(t,Q[i, :])*h
    end
    return t_range, Q
end

function InstableAlgorithm(func::Function, xᵢₙᵢₜ::Float64, tₘᵢₙ::Float64, tₘₐₓ::Float64, h::Float64 )
    t_range=collect(range(tₘᵢₙ, tₘₐₓ, step=h))
    Q=zeros(length(t_range))

    Q[1]=xᵢₙᵢₜ
    Q[2]=xᵢₙᵢₜ + h*func(tₘᵢₙ, xᵢₙᵢₜ )
    for i in 2:length(t_range)-1
        Q[i+1] = Q[i-1] + 2*h*func(t_range[i], Q[i])
    end
    return t_range, Q
end

#Second Order ODE
#Simple Harmonic Oscillator

function EulerCromer(func::Function, xᵢₙᵢₜ::Float64, vᵢₙᵢₜ::Float64, tₘᵢₙ::Float64, tₘₐₓ::Float64, h::Float64)
    t=collect(range(tₘᵢₙ, tₘₐₓ, step=h))
    x=zeros(length(t))
    v=zeros(length(t))

    v[1]= vᵢₙᵢₜ
    x[1]= xᵢₙᵢₜ
    for i in 1:length(t)-1
        v[i+1] = v[i] + h * func(t[i], x[i])
        x[i+1] = x[i] + h * v[i+1]
    end

    return t, x, v
end

function Verlet(func::Function, xᵢₙᵢₜ::Float64, vᵢₙᵢₜ::Float64, tₘᵢₙ::Float64, tₘₐₓ::Float64, h::Float64)
    t=collect(range(tₘᵢₙ, tₘₐₓ, step=h))
    x=zeros(length(t))
    v=zeros(length(t))

    v[1]= vᵢₙᵢₜ
    x[1]= xᵢₙᵢₜ
    x[2]= xᵢₙᵢₜ + vᵢₙᵢₜ * h + func(tₘᵢₙ, xᵢₙᵢₜ ) * h^2 *1/2

    for i in 2: length(t)-1
        x[i+1]= 2 * x[i] - x[i-1] + func(t[i], x[i]) * h^2
        v[i+1]= (x[i+1] - x[i])/h
    end
    return t, x, v
end

function VelocityVerlet(func::Function, xᵢₙᵢₜ::Float64, vᵢₙᵢₜ::Float64, tₘᵢₙ::Float64, tₘₐₓ::Float64, h::Float64)
    t= collect(range(tₘᵢₙ, tₘₐₓ, step=h))
    x=zeros(length(t))
    v=zeros(length(t))

    v[1]= vᵢₙᵢₜ
    x[1]= xᵢₙᵢₜ

    for i in 1:length(t)-1
        x[i+1] = x[i] + v[i] * h + func(t[i], x[i]) * h^2 * 1/2
        v[i+1] = v[i] + (func(t[i], x[i]) + func(t[i+1], x[i+1])) * h * 1/2
    end

    return t, x, v
end

end


text:
"""

    calib_wrapper(param, st_hydro, prec, epot, q_obs, q_sim)

Wrapper function required for calibrating hydrological routing model.

"""

function calib_wrapper(param, st_hydro, prec, epot, q_obs, q_sim, states_sim, warmup, force_states)

    # Assign parameter values

    assign_param(st_hydro, param)

    # Initilize model states

    init_states(st_hydro)

    # Run model

    ntimes = size(prec, 2)

    for itime in 1:ntimes

        get_input(st_hydro, prec, epot, itime)

        run_timestep(st_hydro)

        q_sim[itime] = st_hydro.q_sim

        if ~isempty(states_sim)

            states_hydro = get_states(st_hydro)

            states_sim[:, itime] = vcat(states_hydro)

        end

    end

    # Compute penatly

    if force_states == true && ~isempty(states_sim)

        penalty = zeros(size(states_sim, 1))

        for istate = 1:size(states_sim, 1)

            slope = cov(states_sim[istate, warmup:end], warmup:ntimes) / var(warmup:ntimes)

            penalty[istate] = 2 ./ (1  + exp(-(10 * slope).^2)) - 1

        end

        penalty = maximum(penalty)

        # penalty = mean(abs(slope))

    else

        penalty = 0.0

    end

    # Compute performance measure

    perf_measure = 1.0 - nse(q_sim[warmup:end], q_obs[warmup:end]) + penalty

    return perf_measure

end

"""

    calib_wrapper(param, st_snow, st_hydro, tair, prec, epot, q_obs)

Wrapper function required for calibrating snow and hydrological routing model.

"""

function calib_wrapper(param, st_snow, st_hydro, tair, prec, epot,
                       q_obs, q_sim, states_sim, warmup, force_states)

    # Assign parameter values

    assign_param(st_snow, param[1:length(st_snow.param)])
    assign_param(st_hydro, param[length(st_snow.param)+1:end])

    # Initilize model states

    init_states(st_snow)
    init_states(st_hydro)

    # Run model

    ntimes = size(prec, 2)

    for itime in 1:ntimes

        get_input(st_snow, prec, tair, itime)

        run_timestep(st_snow)

        get_input(st_snow, st_hydro, epot, itime)

        run_timestep(st_hydro)

        q_sim[itime] = st_hydro.q_sim

        states_snow = get_states(st_snow)
        states_hydro = get_states(st_hydro)

        states_sim[:, itime] = vcat(states_snow, states_hydro)

    end

    # Compute penatly

    if force_states == true

        penalty = zeros(size(states_sim, 1))

        for istate = 1:size(states_sim, 1)

            slope = cov(states_sim[istate, warmup:end], warmup:ntimes) / var(warmup:ntimes)

            penalty[istate] = 2 ./ (1  + exp(-(10 * slope).^2)) - 1

        end

        penalty = maximum(penalty)

        # penalty = mean(abs(slope))

    else

        penalty = 0.0

    end

    # Compute performance measure

    perf_measure = 1.0 - nse(q_sim[warmup:end], q_obs[warmup:end]) + penalty

    return perf_measure

end

"""

    run_model_calib(st_hydro::Hydro, prec, epot, q_obs)

Run calibration of hydrological routing model, for example HBV.

"""

function run_model_calib(st_hydro::Hydro, prec, epot, q_obs;
                         verbose = :silent, warmup = 3*365, force_states = false)

    # Get parameter range

    param_range = get_param_range(st_hydro)

    # Allocate output array

    states_hydro = Vann.get_states(st_hydro)

    ntimes = size(prec, 2)
    nstates = length(states_hydro)

    q_sim = zeros(ntimes)

    states_sim = zeros(nstates, ntimes)

    # Run calibration

    calib_wrapper_tmp(param) = calib_wrapper(param, st_hydro, prec, epot, q_obs,
                                             q_sim, states_sim, warmup, force_states)

    res = bboptimize(calib_wrapper_tmp; SearchRange = param_range, TraceMode = verbose)

    param_hydro = best_candidate(res)

    return param_hydro

end

"""

    run_model_calib(st_snow::Snow, st_hydro::Hydro, tair, prec, epot, q_obs)

Run calibration of snow and hydrological routing model.

"""

function run_model_calib(st_snow::Snow, st_hydro::Hydro, tair, prec, epot, q_obs;
                         verbose = :silent, warmup = 3*365, force_states = false)

    # Get parameter range

    param_range_snow  = get_param_range(st_snow)
    param_range_hydro = get_param_range(st_hydro)

    param_range = vcat(param_range_snow, param_range_hydro)

    # Allocate output array

    states_snow = Vann.get_states(st_snow)
    states_hydro = Vann.get_states(st_hydro)

    ntimes = size(prec, 2)
    nstates = length(vcat(states_snow, states_hydro))

    q_sim = zeros(ntimes)
    states_sim = zeros(nstates, ntimes)

    # Run calibration

    calib_wrapper_tmp(param) = calib_wrapper(param, st_snow, st_hydro, tair, prec, epot,
                                             q_obs, q_sim, states_sim, warmup, force_states)

    res = bboptimize(calib_wrapper_tmp; SearchRange = param_range, TraceMode = verbose)

    # Extract parameters for snow and hydrological routing model

    param = best_candidate(res)

    param_snow = param[1:length(st_snow.param)]
    param_hydro = param[length(st_snow.param)+1:end]

    return param_snow, param_hydro

end

"""

    kge(q_sim, q_obs)

Compute modified Kling-Gupta efficiency

"""

function kge(q_sim, q_obs)

    # ikeep = q_obs .!= -999.

    if all(isnan, q_sim) || all(isnan, q_obs)

        kge = NaN

    else

        ikeep = !isnan(q_obs)

        q_sim = q_sim[ikeep]
        q_obs = q_obs[ikeep]

        r = cor(q_sim, q_obs)

        beta = mean(q_sim) / mean(q_obs)

        gamma = (std(q_sim) / mean(q_sim)) / (std(q_obs) / mean(q_obs))

        kge = 1 - sqrt( (r-1)^2 + (beta-1)^2 + (gamma-1)^2 )

    end

    return kge

end

"""

    nse(q_sim, q_obs)

Compute Nash-Sutcliffe efficiency

"""

function nse(q_sim, q_obs)

    # ikeep = q_obs .!= -999.

    ikeep = !isnan(q_obs)

    q_sim = q_sim[ikeep]
    q_obs = q_obs[ikeep]

    ns = 1 - sum((q_sim - q_obs).^2) / sum((q_obs - mean(q_obs)).^2)

end

#=
### Snow and hydrological model

# Run model and compute performance measure

function calib_wrapper(param::Vector, grad::Vector, st_snow, st_hydro, date, tair, prec, epot, q_obs)

  # Assign parameter values

  assign_param(st_snow, param[1:length(st_snow.param)]);
  assign_param(st_hydro, param[length(st_snow.param)+1:end]);

  # Initilize model states

  init_states(st_snow);
  init_states(st_hydro);

  # Run model

  ntimes = size(prec, 2);

  q_sim = zeros(Float64, ntimes);

  for itime in 1:ntimes

    get_input(st_snow, prec, tair, date, itime);

    run_timestep(st_snow);

    get_input(st_snow, st_hydro, epot, itime);

    run_timestep(st_hydro);

    q_sim[itime] = st_hydro.q_sim

  end

  return(kge(q_sim, q_obs));

end

# Run calibration

function run_model_calib(st_snow::Snow, st_hydro::Hydro, date, tair, prec, epot, q_obs)

  # Get parameter range

  param_range_snow  = get_param_range(st_snow);
  param_range_hydro = get_param_range(st_hydro);

  param_range = vcat(param_range_snow, param_range_hydro);

  # Lower and upper bounds for parameters

  param_lower = [param_range[i][1] for i in eachindex(param_range)];
  param_upper = [param_range[i][2] for i in eachindex(param_range)];

  # Starting point

  param_start = (param_lower + param_upper) / 2;

  # Wrapper function

  calib_wrapper_tmp(param::Vector, grad::Vector) = calib_wrapper(param::Vector, grad::Vector, st_snow, st_hydro, date, tair, prec, epot, q_obs);

  # Perform global optimization

  opt_global = Opt(:GN_ESCH, length(param_start));

  max_objective!(opt_global, calib_wrapper_tmp);

  maxeval!(opt_global, 5000);

  lower_bounds!(opt_global, param_lower);
  upper_bounds!(opt_global, param_upper);

  (min_func, best_global, ret_nlopt) = optimize(opt_global, param_start);

  # Check that optimal value is inside bounds

  for iparam = 1:length(best_global)

    if best_global[iparam] < param_lower[iparam] + 0.2*(param_upper[iparam] - param_lower[iparam])
      best_global[iparam] = param_lower[iparam] + 0.2*(param_upper[iparam] - param_lower[iparam]);
    end

    if best_global[iparam] > param_upper[iparam] - 0.2*(param_upper[iparam] - param_lower[iparam]);
      best_global[iparam] = param_upper[iparam] - 0.2*(param_upper[iparam] - param_lower[iparam]);
    end

  end

  # Perform local optimization

  opt_local = Opt(:LN_NELDERMEAD, length(param_start));

  max_objective!(opt_local, calib_wrapper_tmp);

  lower_bounds!(opt_local, param_lower);
  upper_bounds!(opt_local, param_upper);

  (min_func, best_local, ret_nlopt) = optimize(opt_local, best_global);

  # Return best parameter values

  return(best_local)

end

### Hydrological model

# Run model and compute performance measure

function calib_wrapper(param::Vector, grad::Vector, st_hydro, prec, epot, q_obs)

  # Assign parameter values

  assign_param(st_hydro, param);

  # Initilize model states

  init_states(st_hydro);

  # Run model

  ntimes = size(prec, 2);

  q_sim = zeros(Float64, ntimes);

  for itime in 1:ntimes

    get_input(st_hydro, prec, epot, itime);

    q_sim[itime] = run_timestep(st_hydro);

  end

  return(kge(q_sim, q_obs));

end

# Run calibration

function run_model_calib(st_hydro::Hydro, prec, epot, q_obs)

  # Get parameter range

  param_range = get_param_range(st_hydro);

  # Lower and upper bounds for parameters

  param_lower = [param_range[i][1] for i in eachindex(param_range)];
  param_upper = [param_range[i][2] for i in eachindex(param_range)];

  # Starting point

  param_start = (param_lower + param_upper) / 2;

  # Wrapper function

  calib_wrapper_tmp(param::Vector, grad::Vector) = calib_wrapper(param::Vector, grad::Vector, st_hydro, prec, epot, q_obs);

  # Perform global optimization

  opt_global = Opt(:GN_ESCH, length(param_start));

  max_objective!(opt_global, calib_wrapper_tmp);

  maxeval!(opt_global, 5000);

  lower_bounds!(opt_global, param_lower);
  upper_bounds!(opt_global, param_upper);

  (min_func, best_global, ret_nlopt) = optimize(opt_global, param_start);

  # Check that optimal value is inside bounds

  for iparam = 1:length(best_global)

    if best_global[iparam] < param_lower[iparam] + 0.2*(param_upper[iparam] - param_lower[iparam])
      best_global[iparam] = param_lower[iparam] + 0.2*(param_upper[iparam] - param_lower[iparam]);
    end

    if best_global[iparam] > param_upper[iparam] - 0.2*(param_upper[iparam] - param_lower[iparam]);
      best_global[iparam] = param_upper[iparam] - 0.2*(param_upper[iparam] - param_lower[iparam]);
    end

  end

  # Perform local optimization

  opt_local = Opt(:LN_NELDERMEAD, length(param_start));

  max_objective!(opt_local, calib_wrapper_tmp);

  lower_bounds!(opt_local, param_lower);
  upper_bounds!(opt_local, param_upper);

  (min_func, best_local, ret_nlopt) = optimize(opt_local, best_global);

  # Return best parameter values

  return(best_local)

end

=#

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
module LibSerialPort

using Libdl
using libserialport_jll

export
    # Enum types
    SPReturn,
    SPMode,
    SPEvent,
    SPBuffer,
    SPParity,
    SPrts,
    SPcts,
    SPdtr,
    SPdsr,
    SPXonXoff,
    SPFlowControl,
    SPSignal,
    SPTransport,

    # Enum values
    SP_OK,
    SP_ERR_ARG,
    SP_ERR_FAIL,
    SP_ERR_MEM,
    SP_ERR_SUPP,

    SP_MODE_READ,
    SP_MODE_WRITE,
    SP_MODE_READ_WRITE,

    SP_EVENT_RX_READY,
    SP_EVENT_TX_READY,
    SP_EVENT_ERROR,

    SP_BUF_INPUT,
    SP_BUF_OUTPUT,
    SP_BUF_BOTH,

    SP_PARITY_INVALID,
    SP_PARITY_NONE,
    SP_PARITY_ODD,
    SP_PARITY_EVEN,
    SP_PARITY_MARK,
    SP_PARITY_SPACE,

    SP_RTS_INVALID,
    SP_RTS_OFF,
    SP_RTS_ON,
    SP_RTS_FLOW_CONTROL,

    SP_CTS_INVALID,
    SP_CTS_IGNORE,
    SP_CTS_FLOW_CONTROL,

    SP_DTR_INVALID,
    SP_DTR_OFF,
    SP_DTR_ON,
    SP_DTR_FLOW_CONTROL,

    SP_DSR_INVALID,
    SP_DSR_IGNORE,
    SP_DSR_FLOW_CONTROL,

    SP_XONXOFF_INVALID,
    SP_XONXOFF_DISABLED,
    SP_XONXOFF_IN,
    SP_XONXOFF_OUT,
    SP_XONXOFF_INOUT,

    SP_FLOWCONTROL_NONE,
    SP_FLOWCONTROL_XONXOFF,
    SP_FLOWCONTROL_RTSCTS,
    SP_FLOWCONTROL_DTRDSR,

    SP_SIG_CTS,
    SP_SIG_DSR,
    SP_SIG_DCD,
    SP_SIG_RI,

    SP_TRANSPORT_NATIVE,
    SP_TRANSPORT_USB,
    SP_TRANSPORT_BLUETOOTH,

    # Functions from libserialport API
    sp_get_port_by_name,
    sp_free_port,
    sp_list_ports,
    sp_copy_port,
    sp_free_port_list,
    sp_open,
    sp_close,
    sp_get_port_name,
    sp_get_port_description,
    sp_get_port_transport,
    sp_get_port_usb_bus_address,
    sp_get_port_usb_vid_pid,
    sp_get_port_usb_manufacturer,
    sp_get_port_usb_product,
    sp_get_port_usb_serial,
    sp_get_port_bluetooth_address,
    sp_get_port_handle,
    sp_new_config,
    sp_free_config,
    sp_get_config,
    sp_set_config,
    sp_set_baudrate,
    sp_get_config_baudrate,
    sp_set_config_baudrate,
    sp_set_bits,
    sp_get_config_bits,
    sp_set_config_bits,
    sp_set_parity,
    sp_get_config_parity,
    sp_set_config_parity,
    sp_set_stopbits,
    sp_get_config_stopbits,
    sp_set_config_stopbits,
    sp_set_rts,
    sp_get_config_rts,
    sp_set_config_rts,
    sp_set_cts,
    sp_get_config_cts,
    sp_set_config_cts,
    sp_set_dtr,
    sp_get_config_dtr,
    sp_set_config_dtr,
    sp_set_dsr,
    sp_get_config_dsr,
    sp_set_config_dsr,
    sp_set_xon_xoff,
    sp_get_config_xon_xoff,
    sp_set_config_xon_xoff,
    sp_set_config_flowcontrol,
    sp_set_flowcontrol,
    sp_blocking_read,
    sp_blocking_read_next,
    sp_nonblocking_read,
    sp_blocking_write,
    sp_nonblocking_write,
    sp_input_waiting,
    sp_output_waiting,
    sp_flush,
    sp_drain,
    sp_new_event_set,
    sp_add_port_events,
    sp_wait,
    sp_free_event_set,
    sp_get_signals,
    sp_start_break,
    sp_end_break,
    sp_last_error_code,
    sp_last_error_message,
    sp_get_major_package_version,
    sp_get_minor_package_version,
    sp_get_micro_package_version,
    sp_get_package_version_string,
    sp_get_current_lib_version,
    sp_get_revision_lib_version,
    sp_get_age_lib_version,
    sp_get_lib_version_string,

    # Wrapper object for high-level API
    SerialPort,

    # Functions from high-level API
    list_ports,
    get_port_list,
    print_port_metadata,
    print_port_settings,
    set_speed,
    set_frame,
    set_flow_control,
    seteof,
    reseteof

include("wrap.jl")
include("high-level-api.jl")
include("utils.jl")

end # LibSerialPort


text:
function part1(input)
    nothing
end

function part2(input)
    x = parse_input(input)

    mask = 1 << 11
    start, stop = 1, length(x)
    last0 = partialsort(x, mask, start, stop)
    if last0 - start > stop - last0 - 1
        start2, stop2 = last0 + 1, stop
        stop = last0
    else
        start2, stop2 = start, last0
        start = last0 + 1
    end

    while start < stop
        mask >>= 1
        last0 = partialsort(x, mask, start, stop)
        if last0 - start > stop - last0 - 1
            stop = last0
        else
            start = last0 + 1
        end
    end
    oxygen = x[start]

    start, stop = start2, stop2
    mask = 1 << 10
    while start < stop
        last0 = partialsort(x, mask, start, stop)
        if last0 - start > stop - last0 - 1
            start = last0 + 1
        else
            stop = last0
        end
        mask >>= 1
    end
    CO2 = x[start]

    return oxygen * CO2
end

function parse_input(input)
    x = Int[]
    sizehint!(x, 1000)
    n = 0
    for c in read(input)
        if c == UInt8('\n')
            push!(x, n)
            n = 0
        else
            n <<= 1
            n |= c == UInt8('1')
        end
    end
    return x
end

function partialsort(x, mask, start, stop)
    last0 = start - 1
    while start < stop
        while start < stop && (x[start] & mask) == 0
            last0 = start
            start += 1
        end
        while start < stop && (x[stop] & mask) != 0
            stop -= 1
        end
        if start < stop
            x[start], x[stop] =  x[stop], x[start]
            last0 = start
            start += 1
            stop -= 1
        end
        if start >= stop
            if (x[stop] & mask) == 0
                last0 = stop
            end
            return last0
        end
    end
end


text:
include("testfunctions.jl")

module opti
using LinearAlgebra

# include("line_search.jl")
include("opti_1d.jl")

function bracketing_line_search(fcn, x0, step_direction; max_iter = 100, alpha = 1.0, x_tol = 1e-2)

    obj_fcn(alpha) = fcn( x0 + step_direction * alpha )
    x_lo = 0.0
    x_hi = alpha
    bracketing_interval = bracket_min(obj_fcn, x_lo, x_hi, max_iter = max_iter)

    stat = golden_section_search(
        obj_fcn,
        bracketing_interval.x_lo,
        bracketing_interval.x_hi,
        max_iter = max_iter,
        x_tol = x_tol)

    x_sol = x0 + stat.x_sol * step_direction
    return (x_sol, stat.y_sol)
end

function basis(n_dim::Int, ii::Int)
    e_i = zeros(n_dim)
    e_i[ii] = 1.0
    return e_i
end

struct DirectMethodToleranceOptions
    f_abs_tol::Float64
    f_rel_tol::Float64
end
function DirectMethodToleranceOptions()
    f_abs_tol = 1e-10
    f_rel_tol = 1e-10
    return DirectMethodToleranceOptions(f_abs_tol, f_rel_tol)
end

function powell(
    fcn,
    x0;
    max_iter = 100,
    step_size_tol = 1e-10,
    line_search_x_tol = 1e-2,
    reset_iter = 5 * length(x0),
    tol_options = DirectMethodToleranceOptions()
    )

    x = copy(x0)
    mach_eps = tol_options.f_abs_tol * 1e-5

    n_dim = length(x0)
    U = Matrix( one(eltype(x0)) * I , n_dim, n_dim )

    break_flag = false

    f_prev = fcn(x0)
    iter = 0
    for outer iter = 1:max_iter

        for ii = 1:n_dim
            (x, f_act) = bracketing_line_search(
                fcn,
                x,
                U[:, ii];
                max_iter = max_iter,
                x_tol = line_search_x_tol)

                # stopping criteria for function value
            if abs(f_act - f_prev) <= tol_options.f_abs_tol || abs(f_act - f_prev) / ( abs(f_prev) == 0 ? mach_eps : abs(f_prev)) < tol_options.f_rel_tol
                break_flag = true
                break
            end

            println(" x = $(x) ")
        end

        if break_flag
            break
        end

        for ii = 2:n_dim
            U[:, ii-1] = U[:, ii]
        end
        d = x - U[:,n_dim]
        U[:,n_dim] = d
        (x_hat, f_act) = bracketing_line_search(
                fcn,
                x,
                U[:, n_dim];
                max_iter = max_iter,
                x_tol = line_search_x_tol)
        println(" x_hat = $(x_hat) ")
        step = x_hat - x
        if norm(step) <= step_size_tol
            break
        end

        if iter % reset_iter == 0
            U = Matrix( one(eltype(x0)) * I , n_dim, n_dim )
        end

    end

    return (x, iter)

end

function exponential_annealing_schedule(iter, start_value = 1.0, gamma = 1.0)
    return gamma * start_value^iter
end

function fast_annealing_schedule(iter, start_value = 1.0)
    return start_value / iter
end

function logarithmic_annealing_schedule(iter, gamma = 1.0)
    return log( gamma * iter )
end

function simulated_annealing(
    fcn,
    x0,
    rand_fcn;
    max_iter = 100,
    annealing_schedule = exponential_annealing_schedule
    )

    # gamma = 10
    # annealing_schedule(iter) = gamma * iter

    n_dim = length(x0)
    x_path = zeros(n_dim, max_iter)

    f_old = fcn(x0)
    x_best = x0
    f_best = f_old

    x_old = copy(x0)
    iter = 0
    for outer iter = 1:max_iter

        x_new = x_old + rand_fcn(x_old, iter)
        x_path[:, iter] = x_new

        f_new = fcn(x_new)
        df = f_new - f_old

        if df <= 0
            x_old = x_new
            f_old = f_new
        elseif rand() < exp( -df / annealing_schedule(iter) )
            x_old = x_new
            f_old = f_new
        end

        if f_new <= f_best
            x_best = x_new
            f_best = f_new
        end
    end

    return (x_best, f_best, iter, x_path)

end

function corana_update!(v, a, c, n_crit_cycles)
    for ii = 1:length(v)

        if a[ii] > 0.6 * n_crit_cycles
            v[ii] *= (1 + c[ii] / 0.4 * ( a[ii]/n_crit_cycles - 0.6 )  )
        elseif a[ii] < 0.4 * n_crit_cycles
            v[ii] /= (1 + c[ii] / 0.4 * (0.4 - a[ii]/n_crit_cycles))
        end

    end
    return v
end

function adaptive_simulated_annealing(
    fcn,
    x0,
    v,
    temperature,
    eps;
    n_crit_cycles = 20,
    neps = 4,
    n_crit_resets = max( 100, 5 * length(x0) ),
    gamma = 0.85,
    c = fill(2, length(x0)),
    max_iter = 100
    )

    x_best = copy(x0)
    x_old  = copy(x0)

    f_old  = fcn(x0)
    f_best = f_old

    y_arr = []
    n_dim = length(x0)
    uniform_random() = 2 * (rand() - 0.5)

    a = zeros(n_dim)
    counts_cycles = 0
    counts_resets = 0

    basis_vectors = Matrix{eltype(x0)}( one(eltype(x0)) * I, n_dim, n_dim )

    iter = 0
    for outer iter = 1:max_iter

        for ii = 1:n_dim
            x_new = x_old + basis_vectors[:,ii] * uniform_random() * v[ii]
            f_new = fcn(x_new)

            df = f_new - f_old
            if df < 0 || rand() < exp( df / temperature )
                x_old = x_new
                f_old = f_new
                a[ii] += 1
                if f_new < f_best
                    x_best = x_new
                    f_best = f_new
                end
            end
        end

        counts_cycles += 1
        # counts_cycles >= n_crit_cycles || continue
        if counts_cycles < n_crit_cycles
            continue
        end

        counts_cycles = 0
        corana_update!(v, a, c, n_crit_cycles)
        fill!(a, 0)
        counts_resets += 1
        if counts_resets < n_crit_resets
            continue
        end

        temperature *= gamma
        counts_resets = 0
        push!(y_arr, f_old)

        Ly = length(y_arr)
        cond = Ly > neps && y_arr[end] - y_best <= eps &&
         (all(abs(y_arr[end] - y_arr[end-u])) <= eps for u in 1:neps)

        if cond
            x_old = x_best
            f_old = f_best
        else
            break
        end

    end

    return (x_best, f_best)
end

end # opti

module otest
using ..optitestfun
using ..opti
using PyPlot
PyPlot.pygui(true)

fcn = optitestfun.rosenbrock # optitestfun.quadratic
x0 = [-1.8, -2.0] # [10, 8.0]

rand_fcn(x, iter) = 1 / log(iter+1) * randn( length(x) )

# (x_best, y_best, iter, x_path) = opti.simulated_annealing(
#     fcn,
#     x0,
#     rand_fcn;
#     max_iter = 300,
#     annealing_schedule = opti.fast_annealing_schedule) # opti.exponential_annealing_schedule)

v = [1.0, 1.0]
temp = 1.0
eps = 1e-1
(x_best, f_best) = opti.adaptive_simulated_annealing(
    fcn,
    x0,
    v,
    temp,
    eps;
    max_iter = 1000)

# PyPlot.figure()
# PyPlot.grid()
# PyPlot.plot(
#     x_path[1,:], x_path[2,:], linestyle = :none, marker = :.)
# PyPlot.plot( x_best[1], x_best[2], linestyle = :none, marker = :., color = :r, markersize = 15)

function test_powell()
    # fcn = optitestfun.quadratic
    # x0 = [10, 8.0]

    fcn(x) = optitestfun.rosenbrock(x, 2)
    x0 = [-1.8, -2.0]

    (x_sol, iter) = opti.powell( fcn, x0, step_size_tol = 1e-2, line_search_x_tol = 1e-3 )

    # xs = collect( range(-1.0, 10, 30) )
    # ys = collect( range(-1.0, 10, 30) )

    # zs = zeros( length(xs), length(ys) )
    # ii = 0
    # for xx in xs
    #     for yy in ys
    #         ii += 1
    #         zs[ii] = fcn([xx, yy])
    #     end
    # end

    # PyPlot.figure()
    # PyPlot.grid()
    # PyPlot.contour(xs, ys, zs, [1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 0.2, 0.5, 1, 2, 5, 10, 20, 30, 50, 100])
end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using Dispersal, Test

@testset "chain growth and allee" begin
    init = [1.0 4.0 7.0;
            2.0 5.0 8.0;
            3.0 6.0 9.0]

    output = ArrayOutput(init; tspan=1:3)

    growthRule = ExactExponentialGrowth(intrinsicrate=-log(2.0))
    alleeRule = AlleeExtinction(minfounders = 2.0)
    ruleset = Ruleset(Chain(growthRule, alleeRule))
    sim!(output, ruleset)
    @test output[1] == init
    @test output[2] == [0.0  2.0  3.5;
                        0.0  2.5  4.0;
                        0.0  3.0  4.5]
    @test output[3] == [0.0  0.0  0.0;
                        0.0  0.0  2.0;
                        0.0  0.0  2.25]
end

@testset "chain discrete growth and survival" begin
    init = [1.0 4.0 7.0;
            2.0 5.0 8.0;
            3.0 6.0 9.0]

    suit = [10.0 10.0 10.0;
            5.0 5.0 5.0;
            2.0 2.0 2.0]

    output = ArrayOutput(init; tspan=1:3,aux=(suit=suit,))

    growthRule = DiscreteGrowth(intrinsicrate=2.0)
    survivalRule = SurvLogLogisticMap(layerkey=Val(:suit), LC50=5.0)
    ruleset = Ruleset(Chain(growthRule, survivalRule))
    sim!(output, ruleset)
    @test output[1] == init
    @test output[2] ≈ [0.965357  3.86143  6.7575;
                        2.0       5.0      8.0;
                        3.13735   6.2747   9.41204] atol=1e-4
    @test output[3] ≈ [0.931913  3.72765  6.52339;
                        2.0       5.0      8.0;
                        3.28098   6.56197  9.84295] atol=1e-4
end

@testset "chain SimData growth and survival" begin
    init = [1.0 4.0 7.0;
            2.0 5.0 8.0;
            3.0 6.0 9.0]

    suit = [10.0 10.0 10.0;
            5.0 5.0 5.0;
            2.0 2.0 2.0]

    output = ArrayOutput(init; tspan=1:3,aux=(suit=suit,))

    growthRule = DiscreteGrowth(intrinsicrate=2.0)
    survivalRule = SurvLogLogisticMap(layerkey=Val(:suit), LC50=5.0)
    ruleset = Ruleset(Chain(growthRule, survivalRule))
    sim!(output, ruleset)
    @test output[1] == init
    @test output[2] ≈ [0.965357  3.86143  6.7575;
                        2.0       5.0      8.0;
                        3.13735   6.2747   9.41204] atol=1e-4
    @test output[3] ≈ [0.931913  3.72765  6.52339;
                        2.0       5.0      8.0;
                        3.28098   6.56197  9.84295] atol=1e-4
end


text:
using SpecialFunctions: besselix, besselkx
using ArbNumerics: besseli, besselk, ArbFloat

dImx(m, γ, ρ) = γ*besselix(m + 1, γ*ρ) + m*besselix(m, γ*ρ)/ρ
dKmx(m, γ, ρ) = -γ*besselkx(m + 1, γ*ρ) + m*besselkx(m, γ*ρ)/ρ

dIm(m, γ, ρ) = γ*besseli(m + 1, γ*ρ) + m*besseli(m, γ*ρ)/ρ
dKm(m, γ, ρ) = -γ*besselk(m + 1, γ*ρ) + m*besselk(m, γ*ρ)/ρ

function _get_αβ_start(m, γ, D, ρ, n, ind)

    α = D[ind - 1] * (n[ind - 1])^2 * dKm(m, γ[ind - 1], ρ[ind]) * besseli(m, γ[ind]*ρ[ind])
    α -= D[ind] * (n[ind])^2 * besselk(m, γ[ind - 1]*ρ[ind]) * dIm(m, γ[ind], ρ[ind])

    β = D[ind - 1] * (n[ind - 1])^2 * dIm(m, γ[ind - 1], ρ[ind]) * besseli(m, γ[ind]*ρ[ind])
    β -= D[ind] * (n[ind])^2 * besseli(m, γ[ind - 1]*ρ[ind]) * dIm(m, γ[ind], ρ[ind])

    return α, β
end

##reduced form
function _get_αβ_startx(m, γ, D, ρ, n, ind)

    α = D[ind - 1] * (n[ind - 1])^2 * dKmx(m, γ[ind - 1], ρ[ind]) * besselix(m, γ[ind]*ρ[ind])
    α -= D[ind] * (n[ind])^2 * besselkx(m, γ[ind - 1]*ρ[ind]) * dImx(m, γ[ind], ρ[ind])
    α *= exp(ρ[ind]*(γ[ind] - γ[ind - 1]))

    β = D[ind - 1] * (n[ind - 1])^2 * dImx(m, γ[ind - 1], ρ[ind]) * besselix(m, γ[ind]*ρ[ind])
    β -= D[ind] * (n[ind])^2 * besselix(m, γ[ind - 1]*ρ[ind]) * dImx(m, γ[ind], ρ[ind])
    β *= exp(ρ[ind]*(γ[ind] - γ[ind - 1]))

    return α, β
end

function _get_αn_βn(m, αnp1, βnp1, γ, D, ρ, n, ind)

    αn = D[ind] * (n[ind])^2 * dKm(m, γ[ind], ρ[ind + 1])
    αn *= (αnp1*besseli(m, γ[ind + 1]*ρ[ind + 1]) - βnp1*besselk(m, γ[ind + 1]*ρ[ind + 1]))
    αn_aux = D[ind + 1] * (n[ind + 1])^2 * besselk(m, γ[ind]*ρ[ind + 1])
    αn_aux *= (αnp1*dIm(m, γ[ind + 1], ρ[ind + 1]) - βnp1*dKm(m, γ[ind + 1], ρ[ind + 1]))
    αn -= αn_aux

    βn = D[ind] * (n[ind])^2 * dIm(m, γ[ind], ρ[ind + 1])
    βn *= (αnp1*besseli(m, γ[ind + 1]*ρ[ind + 1]) - βnp1*besselk(m, γ[ind + 1]*ρ[ind + 1]))
    βn_aux = D[ind + 1] * (n[ind + 1])^2 * besseli(m, γ[ind]*ρ[ind + 1])
    βn_aux *= (αnp1*dIm(m, γ[ind + 1], ρ[ind + 1]) - βnp1*dKm(m, γ[ind + 1], ρ[ind + 1]))
    βn -= βn_aux

    return αn, βn
end
##reduced form
function _get_αn_βnx(m, αnp1, βnp1, γ, D, ρ, n, ind)

    αn = D[ind] * (n[ind])^2 * dKmx(m, γ[ind], ρ[ind + 1])
    αn *= (αnp1*besselix(m, γ[ind + 1]*ρ[ind + 1]) - βnp1*exp(-2*ρ[ind + 1]*γ[ind + 1])*besselkx(m, γ[ind + 1]*ρ[ind + 1]))
    αn_aux = D[ind + 1] * (n[ind + 1])^2 * besselkx(m, γ[ind]*ρ[ind + 1])
    αn_aux *= (αnp1*dImx(m, γ[ind + 1], ρ[ind + 1]) - βnp1*exp(-2*ρ[ind + 1]*γ[ind + 1])*dKmx(m, γ[ind + 1], ρ[ind + 1]))
    αn -= αn_aux
    αn *= exp(ρ[ind + 1]*(γ[ind + 1] - γ[ind]))

    βn = D[ind] * (n[ind])^2 * dImx(m, γ[ind], ρ[ind + 1])
    βn *= (αnp1*besselix(m, γ[ind + 1]*ρ[ind + 1]) - βnp1*exp(-2*ρ[ind + 1]*γ[ind + 1])*besselkx(m, γ[ind + 1]*ρ[ind + 1]))
    βn_aux = D[ind + 1] * (n[ind + 1])^2 * besselix(m, γ[ind]*ρ[ind + 1])
    βn_aux *= (αnp1*dImx(m, γ[ind + 1], ρ[ind + 1]) - βnp1*exp(-2*ρ[ind + 1]*γ[ind + 1])*dKmx(m, γ[ind + 1], ρ[ind + 1]))
    βn -= βn_aux
    βn *= exp(ρ[ind + 1]*(γ[ind + 1] + γ[ind]))

    return αn, βn
end

function _get_α1_β1(m, γ, D, ρ, n)

    ind = length(γ)
    α, β = _get_αβ_start(m, γ, D, ρ, n, ind)
    ind -= 1
    for i in ind:-1:2
        println(α, β)
        α, β = _get_αn_βn(m, α, β, γ, D, ρ, n, i)
    end

    return α, β
end
##reduced form
function _get_α1_β1x(m, γ, D, ρ, n)

    ind = length(γ)
    α, β = _get_αβ_startx(m, γ, D, ρ, n, ind)
    ind -= 1
    for i in ind:-1:2
        println(α, β)
        α, β = _get_αn_βnx(m, α, β, γ, D, ρ, n, i)
    end

    return α, β
end

function _get_A1_B1(m, γ, D, ρ, n, ρ0)

    α1, β1 = _get_α1_β1(m, γ, D, ρ, n)

    A1 = besselk(m, γ[1]*ρ[1])
    A1 *= (β1*besselk(m, γ[1]*ρ0) - α1*besseli(m, γ[1]*ρ0))
    A1 /= D[1]*(α1*besseli(m, γ[1]*ρ[1]) - β1*besselk(m, γ[1]*ρ[1]))

    B1 = β1/D[1]
    B1 /= (α1*besseli(m, γ[1]*ρ[1]) - β1*besselk(m, γ[1]*ρ[1]))
    B1 *= (besseli(m, γ[1]*ρ0)*besselk(m, γ[1]*ρ[1]) - besselk(m, γ[1]*ρ0)*besseli(m, γ[1]*ρ[1]))

    return A1, B1
end

function _integrand_modbessel(m, A1, B1, γ, ρ, ϕ)

    out = A1*besseli(m, γ[1]*ρ[1]) + B1*besselk(m, γ[1]*ρ[1])
    out *= cos(m*ϕ)

    return out
end

##reduced form (calculates A1 and B1 together)
function _integrand_modbesselx(m, γ, D, ρ, n, ρ0, ϕ)

    α1, β1 = _get_α1_β1x(m, γ, D, ρ, n)

    A1 = (β1*besselkx(m, γ[1]*ρ0)*exp(-2*γ[1]*ρ0) - α1*besselix(m, γ[1]*ρ0))
    A1 /= (α1*besselix(m, γ[1]*ρ[1]) - β1*besselkx(m, γ[1]*ρ[1])*exp(-2*γ[1]*ρ[1]))
    A1 *= besselkx(m, γ[1]*ρ[1])*besselix(m, γ[1]*ρ[1])/D[1]
    A1 *= exp(γ[1]*(ρ0 - ρ[1]))

    B1 = (besselix(m, γ[1]*ρ0)*besselkx(m, γ[1]*ρ[1])*exp(γ[1]*(ρ0 - ρ[1])) - besselkx(m, γ[1]*ρ0)*besselix(m, γ[1]*ρ[1])*exp(γ[1]*(-ρ0 + ρ[1])))
    B1 /= (α1*besselix(m, γ[1]*ρ[1]) - β1*besselkx(m, γ[1]*ρ[1])*exp(-2*γ[1]*ρ[1]))
    B1 *=  besselkx(m, γ[1]*ρ[1])*β1/D[1]
    B1 *= exp(-2*γ[1]*ρ[1])

    return A1*cos(m*ϕ)
end
function sum_modbessel(k, μa, D, ρ, n, ρ0, ϕ)

    a = 0.0
    A1 = 0.0
    B1 = 0.0
    γ = @. sqrt(μa/D + k^2)

    bessel_index = -100:100
    for m in bessel_index
        A1, B1 = _get_A1_B1(m, γ, D, ρ, n, ρ0)
        a += _integrand_modbessel(m, A1, B1, γ, ρ, ϕ)
    end

    return a
end

##reduced form
function sum_modbesselx(k, μa, D, ρ, n, ρ0, ϕ)

    a = 0.0
    A1 = 0.0
    B1 = 0.0
    γ = @. sqrt(μa/D + k^2)

    bessel_index = 0:200
    for m in bessel_index
        a += (2 - (m == 0 ? 1 : 0)) * _integrand_modbesselx(m, γ, D, ρ, n, ρ0, ϕ)
    end

    return a
end
## switch to arbitrary precision
function sum_modbessel(k, μa, D, ρ, n, ρ0, ϕ)
    a = 0.0
    A1 = 0.0
    B1 = 0.0
    γ = @. sqrt(μa/D + k^2)
    γ1 = @. sqrt(ArbFloat.(μa)/D + k^2)
    if k < 10.0
        bessel_index = 0:500
    else
        bessel_index = 0:500
    end

    for m in bessel_index

        if besselix(m,minimum(γ)*minimum([ρ; ρ0])) < floatmin(eltype(μa))*50000 # besselkx(m,maximum(γ)*maximum([ρ; ρ0])) > typemax(eltype(μa))/2 #besselix(m, minimum(γ)*minimum([ρ; ρ0])) == 0
            #println(m); println("breaking")
            #break
            A1, B1 = _get_A1_B1(m, γ1, D, ρ, n, ρ0)
            #println(besseli(m, γ1[2]*ρ[2]))
            a += convert(eltype(μa), (2 - (m == 0 ? 1 : 0)) * _integrand_modbessel(m, A1, B1, γ1, ρ, ϕ))
            #println(a)
        else
            besselix(m,minimum(γ)*minimum([ρ; ρ0]))
            a += (2 - (m == 0 ? 1 : 0)) * _integrand_modbesselx(m, γ, D, ρ, n, ρ0, ϕ)
           # println(a)
        end

    end

    return a
end

function _greens_homogenous(μa, D, ρ, n, ρ0, ϕ, z)

    g = quadgk(k -> cos(k*z)*sum_modbessel(k, μa, D, ρ, n, ρ0, ϕ), 0, 400; maxevals=10^3)

    return g[1]/(2*π^2)
end

function _green_Nradialcyl(μa, D, ρ, n, ρ0, ϕ, z)

    r = sqrt(ρ[1]^2 + ρ0^2 - 2*ρ[1]*ρ0*cosh(ϕ) + z^2)

    g = exp(-r*sqrt(μa[1]/D[1]))/(4*π*D[1]*r)
    println(g)
    gh = _greens_homogenous(μa, D, ρ, n, ρ0, ϕ, z)

    return g + gh
end

function _greens_homogenous(μa, D, ρ, n, ρ0, ϕ, z)

    N = 190
    x, w = gausslegendre(N)
    ub = 150

    green = zeros(BigFloat, N)

    Threads.@threads for m in eachindex(x)
        green[m] = cos(x[m]*ub/2 +ub/2*z)*sum_modbessel(x[m]*ub/2 +ub/2, μa, D, ρ, n, ρ0, ϕ)*w[m]*ub
    end

    return sum(green)/(2*pi^2)
end

max_issues_repo_name: heltonmc/LightPropagation


text:
# basic test of multiplying big and small together should give small

using ApproxManifoldProducts
using Test

@testset "multiply big and small" begin

N = 100

X1 = randn(3,N)
X2 = 0.1.*randn(3,N)

p = manikde!(X1, (:Euclid, :Euclid, :Circular))
q = manikde!(X2, (:Euclid, :Euclid, :Circular))

pq = manifoldProduct([p;q],(:Euclid,:Euclid,:Circular))

X12 = getPoints(pq)

@test 0.7*N < sum(abs.(X12[1,:]) .< 0.3)
@test 0.7*N < sum(abs.(X12[2,:]) .< 0.3)
@test 0.7*N < sum(abs.(X12[3,:]) .< 0.3)

end


text:
# Broadcasting over pipes

add1(x) = x + 1
mul2(x) = x * 2
add1mul2v(xs) = xs |> add1v |> mul2v

add1v(x) = x .+ 1
mul2v(x) = x .* 2
add1mul2(xs) = xs .|> add1 .|> mul2

using BenchmarkTools
xs = collect(1:10000);
@btime add1mul2v($xs);
@btime add1mul2($xs);

#=
julia> @btime add1mul2v($xs);
  16.128 μs (4 allocations: 156.41 KiB)

julia> @btime add1mul2($xs);
  8.417 μs (2 allocations: 78.20 KiB)
=#


text:
#using Gadfly
using PyPlot
f = open("time2")
a = Float64[]
for ln in eachline(f)
    n = parse(Float64, ln)
    push!(a, n)
end
plot(a)
show()

max_stars_repo_name: frostRed/Exercises

max_issues_repo_name: frostRed/Exercises

max_forks_repo_name: frostRed/Exercises


text:
# Parameters:
#	N::Int
#	- the problem size is N
#	SNF::Array{Foat64, 1}
#	- SNF[i] is the value associated with the Permutation that permutation_index() maps to NZL[i]
#	NZL::Array{Int, 1}
#	- NZL[i] the set of NonZeroLocations for the sparse function over Sn
#	- NZL must be in increasing order
#	YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1}
#	- output1 from yor()
#	PT::Array{Array{Array{Int, 1}, 1}, 1}
#	- output2 from yor()
# Return Values:
#	FFT::Array{Float64, 2}
#	- FFT is the Fast Fourier Transform of SNF
function sn_fft_sp(N::Int, SNF::Array{Float64, 1}, NZL::Array{Int, 1}, YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1}, PT::Array{Array{Array{Int, 1}, 1}, 1})
	np = nprocs()
	if np == 1 || N < 10
		return compute_fft_sp(N, SNF, YOR, PT, NZL, 1, Counter(1))
	else
		sFFT = Array(Array{Array{Float64, 2}, 1}, N)
		RR_YOR = Array(RemoteRef, np)
		RR_PT = Array(RemoteRef, np)
		for p = 1:np
			if p != myid()
				RR_YOR[p] = RemoteRef(p)
				put!(RR_YOR[p], YOR)
				RR_PT[p] = RemoteRef(p)
				put!(RR_PT[p], PT)
			end
		end
		BS = factorial(N - 1)
		sSNF = Array(Array{Float64, 1}, N)
		sNZL = Array(Array{Int, 1}, N)
		SM = zeros(Int, N)
		si = 1
		ei = 1
		for n = 1:N
			if si > length(NZL)
				sSNF[n] = Array(Float64, 0)
				sNZL[n] = Array(Int, 0)
			else
				ni = NZL[ei]
				lb = (n - 1) * BS + 1
				ub = lb + BS - 1
				while lb <= ni && ni <= ub
					SM[n] = 1
					ei += 1
					if ei > length(NZL)
						break
					end
					ni = NZL[ei]
				end
				sSNF[n] = SNF[si:(ei - 1)]
				sNZL[n] = NZL[si:(ei - 1)]
				si = ei
			end
		end
		i = 1
		nextidx() = (idx = i; i += 1; idx)
		@sync begin
			for p = 1:np
				if p != myid()
					@async begin
						while true
							idx = nextidx()
							if idx > N
								break
							end
							if SM[idx] == 1
								lb = (idx - 1) * BS + 1
								sFFT[idx] = remotecall_fetch(p, compute_fft_sp_remote, N - 1, sSNF[idx], RR_YOR[p], RR_PT[p], sNZL[idx], lb, Counter(1))
							else
								sFFT[idx] = Array(Array{Float64, 2}, 0)
							end
						end
					end
				end
			end
		end
		YORn = YOR[N]
		NP = length(YORn)
		PTn = PT[N]
		RR_sFFT = Array(RemoteRef, np)
		for p = 1:np
			if p != myid()
				RR_sFFT[p] = RemoteRef(p)
				put!(RR_sFFT[p], sFFT)
			end
		end
		FFT = Array(Array{Float64, 2}, NP)
		i = 1
		nextidx() = (idx = i; i += 1; idx)
		@sync begin
			for p = 1:np
				if p != myid()
					@async begin
						while true
							idx = nextidx()
							if idx > NP
								break
							end
							FFT[idx] = remotecall_fetch(p, fc_sp_remote, N, idx, RR_YOR[p], RR_PT[p], RR_sFFT[p], SM)
						end
					end
				end
			end
		end
		return FFT
	end
end

# Parameters:
#	N::Int
#	- the problem size is N
#	SNF::Array{Foat64, 1}
#	- SNF[i] is the value associated with the Permutation that permutation_index() maps to NZL[i]
#	NZL::Array{Int, 1}
#	- NZL[i] the set of NonZeroLocations for the sparse function over Sn
#	- NZL must be in increasing order
#	YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1}
#	- output1 from yor()
#	PT::Array{Array{Array{Int, 1}, 1}, 1}
#	- output2 from yor()
#	LB::Int
#	- the LowerBound on the value of permutation_index() for Permutations in this subgroup
#	C::Counter
#	- a wrapper for an Int that is used to coordinate the recursive access to the elements of SNF
# Return Values:
#	FFT::Array{Float64, 2}
#	- FFT is the Fast Fourier Transform of SNF
function compute_fft_sp(N::Int, SNF::Array{Float64, 1}, YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1}, PT::Array{Array{Array{Int, 1}, 1}, 1}, NZL::Array{Int, 1}, LB::Int, C::Counter)
	if N == 1
		sFFT = Array(Array{Float64, 2}, 1)
		sFFTi = Array(Float64, 1, 1)
		sFFTi[1, 1] = SNF[C.N]
		sFFT[1] = sFFTi
		C.N += 1
		return sFFT
	end

	sFFT = Array(Array{Array{Float64, 2}, 1}, N)
	SM = zeros(Int, N)
	lbi = factorial(N - 1) - 1
	for n = 1:N
		if C.N > length(NZL)
			sFFT[n] = Array(Array{Float64, 2}, 0)
		else
			index = NZL[C.N]
			ub = LB + lbi
			if LB <= index && index <= ub
				SM[n] = 1
				sFFT[n] = compute_fft_sp(N - 1, SNF, YOR, PT, NZL, LB, C)
			else
				sFFT[n] = Array(Array{Float64, 2}, 0)
			end
			LB = ub + 1
		end
	end
	FFT = combine_sfft_sp(N, YOR[N], PT[N], sFFT, SM)
	return FFT
end

function compute_fft_sp_remote(N::Int, SNF::Array{Float64, 1}, RR_YOR::RemoteRef, RR_PT::RemoteRef, NZL::Array{Int, 1}, LB::Int, C::Counter)
	YOR = fetch(RR_YOR)
	PT = fetch(RR_PT)
	FFT = compute_fft_sp(N, SNF, YOR, PT, NZL, LB, C)
	return FFT
end

# Parameters:
#	N::Int
#	- the size of the FFT being calculated is N
#	YORn::Array{Array{SparseMatrixCSC, 1}, 1}, 1}
#	- Young's Orthogonal Representations for the Partitions of N
#	PTn::Array{Array{Int, 1}, 1}
#	- the decomposition indices for the Partitions of N
#	sFFT::Array{Array{Array{Float64, 2}, 1}, 1}
#	- the array of N FFT's of size N-1 that is used to compute this FFT of size N
#	SM::Array{Int, 1}
#	- SM[n] == 0 if sFFT[n] corresponded to a homogenous subgroup with a value of 0
#	- SM[n] == 1 otherwise
# Return Values:
#	FFT::Array{Array{Float64, 2}, 1}
#	- FFT is the Fast Fourier Transform for this group
function combine_sfft_sp(N::Int, YORn::Array{Array{SparseMatrixCSC, 1}, 1}, PTn::Array{Array{Int, 1}, 1}, sFFT::Array{Array{Array{Float64, 2}, 1}, 1}, SM::Array{Int, 1})
	NP = length(YORn)
	FFT = Array(Array{Float64, 2}, NP)
	for p = 1:NP
		YORnp = YORn[p]
		PTnp = PTn[p]
		FC = fc_sp(N, YORnp, PTnp, sFFT, SM)
		FFT[p] = FC
	end
	return FFT
end

# Parameters:
#	N::Int
#	- the Fourier Component that is being calculated corresponds to a Partition, P, of N
#	YORnp::Array{SparseMatrixCSC, 1}
#	- Young's Orthogonal Representations for P
#	PTnp::Array{Int, 1}
#	- the indices for the Partitions that P decomposes into
#	sFFT::Array{Array{Array{Float64, 2}, 1}, 1}
#	- the array of N FFT's of size N-1 that is used to compute this FFT of size N
#	SM::Array{Int, 1}
#	- SM[n] == 0 if sFFT[n] corresponded to a homogenous subgroup with a value of 0
#	- SM[n] == 1 otherwise
# Return Values:
#	FC::Array{Float64, 2}
#	- FC is the Fourier Coefficient corresponding to the Partition P
function fc_sp(N::Int, YORnp::Array{SparseMatrixCSC, 1}, PTnp::Array{Int, 1}, sFFT::Array{Array{Array{Float64, 2}, 1}, 1}, SM::Array{Int, 1})
	Dim = size(YORnp[1], 1)
	if SM[N] == 1
		FC = dsm(Dim, sFFT[N], PTnp)
	else
		FC = zeros(Float64, Dim, Dim)
	end
	CCM = eye(Dim)
	for n = (N - 1):-1:1
		CCM = YORnp[n] * CCM
		if SM[n] == 1
			DSM = dsm(Dim, sFFT[n], PTnp)
			FC_n = CCM * DSM
			FC += FC_n
		end
	end
	return FC
end

function fc_sp_remote(N::Int, p::Int, RR_YOR::RemoteRef, RR_PT::RemoteRef, RR_sFFT::RemoteRef, SM::Array{Int, 1})
	YOR = fetch(RR_YOR)
	PT = fetch(RR_PT)
	sFFT = 	fetch(RR_sFFT)
	FC = fc_sp(N, YOR[N][p], PT[N][p], sFFT, SM)
	return FC
end


text:
using DamagedShearBand; const DSB = DamagedShearBand
using Tensors
using StaticArrays
using Test

# Data :
S, σ₃ = 5, -1e6
r = DSB.Rheology()
D = r.D₀
σᵢⱼ = SymmetricTensor{2,3}([S*σ₃ 0 0 ; 0 σ₃ 0 ; 0 0 r.ν*(S+1)*σ₃])

@testset "get_plane_strain_oop_stress" begin
    atol=1e-12
    σᵢⱼ_cor = DSB.set_plane_strain_oop_stress(σᵢⱼ, r, D ; abstol=atol, maxiter=100, σoop_guess=nothing)
    ϵoop = DSB.compute_ϵ_oop(r,D,σᵢⱼ_cor)
    @test abs(ϵoop)<=atol
end

@testset "get_damage_onset" begin
    Sc = DSB.get_damage_onset(r,σ₃)

    # KI at damage onset
    σᵢⱼ_onset = DSB.build_principal_stress_tensor(r,Sc,σ₃)
    KI_onset = DSB.compute_KI(r, σᵢⱼ_onset, r.D₀)
    @test KI_onset < 1e-10

    # KI after damage onset
    σᵢⱼ_up = DSB.build_principal_stress_tensor(r,Sc+0.1,σ₃)
    KI_up = DSB.compute_KI(r, σᵢⱼ_up, r.D₀)
    @test KI_up > 0

    # KI before damage onset
    σᵢⱼ_down = DSB.build_principal_stress_tensor(r,Sc-0.1,σ₃)
    KI_down = DSB.compute_KI(r, σᵢⱼ_down, r.D₀)
    @test KI_down < 0
end

@testset "band_coords" begin
    σ0 = DSB.band_coords(σᵢⱼ,0) # plane is vertical,
    @test σ0[1,1] == σᵢⱼ[2,2]
    σm90 = DSB.band_coords(σᵢⱼ,-90) # plane is horizontal,
    @test all(σm90 .== σᵢⱼ)
    σ90 = DSB.band_coords(σᵢⱼ,90) # plane is horizontal,
    @test all(σ90 .== σᵢⱼ)
    σ90 = DSB.band_coords(σᵢⱼ,90+360) # plane is horizontal,
    @test all(σ90 .== σᵢⱼ)
    σ_rotated = DSB.band_coords(σᵢⱼ,30)
    σ_back = DSB.principal_coords(σ_rotated,30)
    @test σ_back ≈ σᵢⱼ
end

# @testset "solving using given ϵ̇11. D=D₀=0 test vs elastic rheology" begin
#     tspan = (0.0,1.0)
#     Δt = 0.1
#     ϵ̇11 = -1e-5
#     r = DSB.Rheology(D₀=0.0) # object containing elastic moduli and damage parameters.
#     D_i = r.D₀
#     σ₃ = -1e6
#     S = 1

#     n_iter_saved = 5000
#     n_iter_printed = 10
#     op = DSB.OutputParams(save_frequency = n_iter_saved/(tspan[2]-tspan[1]),
#                         save_period = nothing,
#                         print_frequency = n_iter_printed/(tspan[2]-tspan[1]),
#                         print_period = nothing)

#     sp = DSB.SolverParams(newton_abstol = 1e-12,
#                         newton_maxiter = 100,
#                         time_maxiter = nothing,
#                         e₀ = (D=1e-1, σ=100.0, ϵ=1e-3))

#     p = DSB.Params(op,sp)

#     σᵢⱼ_i = DSB.build_principal_stress_tensor(r,S,σ₃,D_i ; abstol=1e-15) # takes care of the plane strain constraint
#     ϵᵢⱼ_i = DSB.compute_ϵij(r,D_i,σᵢⱼ_i)
#     t_vec, σᵢⱼ_vec, ϵᵢⱼ_vec, D_vec = DSB.time_integration(r,p,σᵢⱼ_i,ϵᵢⱼ_i,D_i,ϵ̇11,Δt,tspan)

#     Ce = DSB.get_elastic_stiffness_tensor(r)
#     σᵢⱼ_vec_elast = Ref(Ce) .⊡ ϵᵢⱼ_vec

#     @test all([σᵢⱼ[1,1] for σᵢⱼ in σᵢⱼ_vec] .≈ [σᵢⱼ[1,1] for σᵢⱼ in σᵢⱼ_vec_elast])
#     @test all([σᵢⱼ[2,2] for σᵢⱼ in σᵢⱼ_vec] .≈ [σᵢⱼ[2,2] for σᵢⱼ in σᵢⱼ_vec_elast])
#     @test all([σᵢⱼ[3,3] for σᵢⱼ in σᵢⱼ_vec] .≈ [σᵢⱼ[3,3] for σᵢⱼ in σᵢⱼ_vec_elast])
#     @test all([σᵢⱼ[1,2] for σᵢⱼ in σᵢⱼ_vec] .≈ [σᵢⱼ[1,2] for σᵢⱼ in σᵢⱼ_vec_elast])
#     @test all([σᵢⱼ[1,3] for σᵢⱼ in σᵢⱼ_vec] .≈ [σᵢⱼ[1,3] for σᵢⱼ in σᵢⱼ_vec_elast])
#     @test all([σᵢⱼ[2,3] for σᵢⱼ in σᵢⱼ_vec] .≈ [σᵢⱼ[2,3] for σᵢⱼ in σᵢⱼ_vec_elast])

#     # check tangent modulus along xx (σxx/ϵxx)
#     Emod = DSB.E_from_Gν(r.G,r.ν)
#     tan_mod_elast_model = (σᵢⱼ_vec_elast[end][1,1]-σᵢⱼ_vec_elast[1][1,1]) / (ϵᵢⱼ_vec[end][1,1]-ϵᵢⱼ_vec[1][1,1])
#     tan_mod_Dat0_model = (σᵢⱼ_vec[end][1,1]-σᵢⱼ_vec_elast[1][1,1]) / (ϵᵢⱼ_vec[end][1,1]-ϵᵢⱼ_vec[1][1,1])
#     tan_mod_theoretical = Emod/(1-r.ν^2)
#     @test abs(tan_mod_Dat0_model-tan_mod_theoretical) < 1e-2
#     @test abs(tan_mod_elast_model-tan_mod_theoretical) < 1e-2
# end

@testset "from principal stress rate to band orientation" begin
    θ = 60.0
    Ṡ = 1e-2
    Δt = 1
    σᵢⱼ_principal = DSB.build_principal_stress_tensor(r,S,σ₃,D ; abstol=1e-15)
    σᵢⱼ1 = DSB.compute_rotated_stress_rate_from_principal_coords(r,Ṡ,σ₃,σᵢⱼ_principal,D,Δt,θ ; damaged_allowed=true)
    σᵢⱼ2 = DSB.compute_rotated_stress_rate_from_band_coords(r,Ṡ,σ₃,DSB.band_coords(σᵢⱼ_principal,θ),D,Δt,θ; damaged_allowed=true)
    @test all(σᵢⱼ1 .≈ σᵢⱼ2)
end

@testset "get_stress_deviation_from_far_field" begin
    σᵢⱼ = DSB.build_principal_stress_tensor(r,S,σ₃,D)
    σᵢⱼ_rotated_90 = DSB.insert_into(σᵢⱼ, (σᵢⱼ[1,1],σᵢⱼ[2,2]), ((2,2),(1,1))) # swap two principals in plane stresses
    σᵢⱼ_rotated_30 = DSB.band_coords(σᵢⱼ,60)
    angle_90 = DSB.get_stress_deviation_from_far_field(σᵢⱼ_rotated_90 ; offdiagtol=1e-5, compression_axis=:x)
    angle_30 = DSB.get_stress_deviation_from_far_field(σᵢⱼ_rotated_30 ; offdiagtol=1e-5, compression_axis=:x)
    @test abs(angle_90) ≈ 90
    @test abs(angle_30) ≈ 30

    # # from simple shear test :
    Dᵢ, _ = DSB.get_KI_minimizer_D_S(r,σ₃)
    σᵢⱼ0 = DSB.build_principal_stress_tensor(r,2,σ₃,Dᵢ ; abstol=1e-15) # σ1 along x
    σᵢⱼ1 = DSB.build_principal_stress_tensor(r,0.5,σ₃,Dᵢ ; abstol=1e-15) # σ1 along y
    σᵢⱼ2p = SymmetricTensor{2,3,Float64}(SA[σᵢⱼ0[1,1], -0.1*σᵢⱼ0[1,1], 0, σᵢⱼ0[2,2], 0, σᵢⱼ0[3,3]]) # σᵢⱼ0 with low positive shear stress
    σᵢⱼ2n = SymmetricTensor{2,3,Float64}(SA[σᵢⱼ0[1,1], 0.1*σᵢⱼ0[1,1], 0, σᵢⱼ0[2,2], 0, σᵢⱼ0[3,3]]) # σᵢⱼ0 with low negative shear stress
    σᵢⱼ3 = DSB.build_principal_stress_tensor(r,1,σ₃,Dᵢ ; abstol=1e-15) # hydrostatic in plane stress
    σ_angle_0 = DSB.get_stress_deviation_from_y(σᵢⱼ0)
    σ_angle_1 = DSB.get_stress_deviation_from_y(σᵢⱼ1)
    σ_angle_2p = DSB.get_stress_deviation_from_y(σᵢⱼ2p)
    σ_angle_2n = DSB.get_stress_deviation_from_y(σᵢⱼ2n)
    σ_angle_3 = DSB.get_stress_deviation_from_y(σᵢⱼ3 ; shear_mode=:simple)
    @test σ_angle_0 == 90 # σ1 along x
    @test σ_angle_1 == 0 # σ1 along y
    @test 0 < σ_angle_2p < 90 # positive trigonometric angle from σyy
    @test -90 < σ_angle_2n < 0 # negative trigonometric angle from σyy => because the angle is expressed in [-π/2 ; π/2]
    @test σ_angle_3 == 45 # set to 45° for hydrostatic conditions in simple shear test because it goes to this value as soon as shear terms become non zero
end

@testset "get_KI_minimizer_D" begin
    r = DSB.Rheology(D₀=0.25)
    Dc = DSB.get_KI_minimizer_D(r,S,σ₃)
    σᵢⱼ_Dc = DSB.build_principal_stress_tensor(r,S,σ₃,Dc)
    KI_Dc = DSB.compute_KI(r, σᵢⱼ_Dc, Dc)

    # KI after damage onset
    KI_up = DSB.KI_from_external_load(r, S, σ₃, Dc+0.01)
    @test KI_up > KI_Dc

    # KI before damage onset
    KI_down = DSB.KI_from_external_load(r, S, σ₃, Dc-0.01)
    @test KI_down > KI_Dc
end

@testset "get_D_and_S_that_minimize_KI" begin
    σ₃ = -1e6
    r = DSB.Rheology(D₀=0.175)
    Dc, Sc = DSB.get_KI_minimizer_D_S(r,σ₃)
    σᵢⱼ_min = DSB.build_principal_stress_tensor(r,Sc,σ₃,Dc)
    KI_min = DSB.compute_KI(r,σᵢⱼ_min,Dc)
    @test abs(KI_min) < 1e-8

    eps = 1e-6
    S_more = Sc + eps
    S_less = Sc - eps
    D_more = Dc + eps
    D_less = Dc - eps
    KI_S_more = DSB.compute_KI(r,DSB.build_principal_stress_tensor(r,S_more,σ₃,Dc),Dc)
    KI_S_less = DSB.compute_KI(r,DSB.build_principal_stress_tensor(r,S_less,σ₃,Dc),Dc)
    KI_D_more = DSB.compute_KI(r,DSB.build_principal_stress_tensor(r,Sc,σ₃,D_more),D_more)
    KI_D_less = DSB.compute_KI(r,DSB.build_principal_stress_tensor(r,Sc,σ₃,D_less),D_less)

    @test abs(KI_S_more) > abs(KI_min)
    @test abs(KI_S_less) > abs(KI_min)
    @test abs(KI_D_more) > abs(KI_min)
    @test abs(KI_S_less) > abs(KI_min)
end

@testset "initialize stresses and strains from D_in and D_out" begin
    # parameters initialization
    σ₃ = -1e6
    S_i = 1.0
    Dⁱ = 0.5
    Dᵒ = 0.4
    θ = 60.0
    op = DSB.OutputParams()
    sp = DSB.SolverParams()
    p = DSB.Params(op,sp)
    r = DSB.Rheology()

    # tested function
    σᵒᵢⱼ, σⁱᵢⱼ, ϵᵒᵢⱼ, ϵⁱᵢⱼ = DSB.initialize_state_var_D(r,p,S_i,σ₃,Dⁱ,Dᵒ,θ ; coords=:band)

    # tests on stress continuity
    @test σᵒᵢⱼ[1,2] ≈ σⁱᵢⱼ[1,2]
    @test σᵒᵢⱼ[2,2] ≈ σⁱᵢⱼ[2,2]

    # test on strain compatibility
    @test ϵᵒᵢⱼ[1,1] ≈ ϵⁱᵢⱼ[1,1]

    #test that epsilon oop is zero
    @test ϵᵒᵢⱼ[3,3] ≈ 0
    @test ϵⁱᵢⱼ[3,3] ≈ 0

    # test that functional form of epsilon applied to sigma gives plane strain
    ϵ_out = DSB.compute_ϵij(r,Dᵒ,σᵒᵢⱼ)
    ϵ_in = DSB.compute_ϵij(r,Dⁱ,σⁱᵢⱼ)
    abstol = 1e-15
    @test abs(ϵ_out[3,3]) < abstol*abs(ϵ_out[1,1])
    @test abs(ϵ_in[3,3]) < abstol*abs(ϵ_in[1,1])
    @test all(i -> i < abstol, abs.(ϵ_out .- ϵᵒᵢⱼ)) == true
    @test all(i -> i < abstol, abs.(ϵ_in .- ϵⁱᵢⱼ)) == true
end

@testset "comparison between two residual forms" begin
    # input :
    r = DSB.Rheology(
      G = 28e9,
      ν = 0.25,
      μ = 0.6, # Friction coef
      β = 0.1, # Correction factor
      K₁c = 1.74e6, # Critical stress intensity factor (Pa.m^(1/2))
      a = 0.1, # Initial flaw size (m)
      ψ = 0.5*atand(1/0.6),# crack angle to the principal stress (radians)
      D₀ = 0.55,# Initial flaw density
      n = 5.5, # Stress corrosion index
      l̇₀ = 7e-5, # Ref. crack growth rate (m/s)
      H = 50e3, # Activation enthalpy (J/mol)
      A = 5.71 # Preexponential factor (m/s)
    )

    σ₃ = -1e6
    du = [-17.060700132915812, -4.572215062877375, -12.380368686287465, 1.530946116594802e-8, 8.045129698897287e-7]
    u = [-1.8458631082604998e6, 878593.6772462169, -1.942972357187292e6, 0.00024964206656874537, 0.9899999999999999]
    σ₁₁, σ₁₂, σₒₒₚ, ϵ₂₂, D = u
    σ̇₁₁, σ̇₁₂, σ̇ₒₒₚ, ϵ̇₂₂, Ḋ = du
    σᵢⱼ = SymmetricTensor{2,3}(SA[σ₁₁ σ₁₂ 0 ; σ₁₂ σ₃ 0 ; 0 0 σₒₒₚ])
    Ḋ = du[5]
    du_nl = SA[σ̇₁₁, σ̇₁₂, σ̇ₒₒₚ, ϵ̇₂₂*r.G]

    σ̇ᵢⱼ = SymmetricTensor{2,3}([ σ̇₁₁    σ̇₁₂   0.0  ;
                                 σ̇₁₂    0.0   0.0  ;
                                 0.0    0.0   σ̇ₒₒₚ ])

    𝕀 = SymmetricTensor{2,3}(DSB.δ) # Second order identity tensor

    # stress at previous timestep
    σ = (1/3)*(tr(σᵢⱼ))
    sᵢⱼ = dev(σᵢⱼ)
    τ = DSB.get_τ(sᵢⱼ)

    # stress derivatives
    σ̇ = (1/3)*(tr(σ̇ᵢⱼ))
    τ̇ = sᵢⱼ ⊡ σ̇ᵢⱼ / (2*τ)

    # damage constants
    c1, c2, c3 = DSB.compute_c1c2c3(r,D)
    A, B = DSB.compute_AB(r,c1,c2,c3)
    A1, B1 = DSB.compute_A1B1(r,A,B)

    dc1dD = DSB.compute_dc1dD(r,D)
    dc2dD = DSB.compute_dc2dD(r,D)
    dc3dD = DSB.compute_dc3dD(r,D)
    dA1dD = DSB.compute_dA1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)
    dB1dD = DSB.compute_dB1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)

    la1 = DSB.λ₁(A1,B1,σ,τ)
    la2 = DSB.λ₂(r,A1,B1,σ,τ)
    la3 = DSB.λ₃(A1,B1)
    dla1dD  = DSB.dλ₁dD(A1,B1,dA1dD,dB1dD,σ,τ)
    dla2dD  = DSB.dλ₂dD(A1,B1,dA1dD,dB1dD,σ,τ)
    dla3dD  = DSB.dλ₃dD(A1,B1,dA1dD,dB1dD)
    dla1dσ  = DSB.dλ₁dσ(A1,B1,τ)
    dla1dτ  = DSB.dλ₁dτ(A1,B1,σ,τ)

    # RES JAO
    res1 = SA[ la1*σ̇₁₁ - la2*σ̇ + la3*τ̇ + Ḋ * (dla1dD*σᵢⱼ[1,1] - dla2dD*σ + dla3dD*τ) + (dla1dσ*σ̇ + dla1dτ*τ̇)*sᵢⱼ[1,1],
              -2*r.G*ϵ̇₂₂ - la2*σ̇ + la3*τ̇ + Ḋ * (dla1dD*params.σ₃ - dla2dD*σ + dla3dD*τ) + (dla1dσ*σ̇ + dla1dτ*τ̇)*sᵢⱼ[2,2],
              la1*σ̇ₒₒₚ - la2*σ̇ + la3*τ̇ + Ḋ * (dla1dD*σᵢⱼ[3,3] - dla2dD*σ + dla3dD*τ) + (dla1dσ*σ̇ + dla1dτ*τ̇)*sᵢⱼ[3,3],
              -2*r.G*params.ϵ̇₁₂ + la1*σ̇₁₂ + Ḋ*dla1dD*σᵢⱼ[1,2] + (dla1dσ*σ̇ + dla1dτ*τ̇)*sᵢⱼ[1,2] ]

    # RES initial
    t1 = DSB.λ₁(A1,B1,σ,τ)*σ̇ᵢⱼ - ( DSB.λ₂(r,A1,B1,σ,τ)*σ̇ - (1/3)*A1*B1*τ̇ )*𝕀
    t2 = ( DSB.dλ₁dσ(A1,B1,τ)*σ̇ + DSB.dλ₁dτ(A1,B1,σ,τ)*τ̇ )*sᵢⱼ
    t3 = DSB.dλ₁dD(A1,B1,dA1dD,dB1dD,σ,τ)*Ḋ*σᵢⱼ - ( DSB.dλ₂dD(A1,B1,dA1dD,dB1dD,σ,τ)*Ḋ*σ - (1/3)*Ḋ*(dA1dD*B1 + A1*dB1dD)*τ )*𝕀
    ϵ̇ᵢⱼ = t1 + t2 + t3

    res2 = SA[ϵ̇ᵢⱼ[1,1],
              -2*r.G*ϵ̇₂₂ + ϵ̇ᵢⱼ[2,2],
              ϵ̇ᵢⱼ[3,3],
              -2*r.G*params.ϵ̇₁₂ + ϵ̇ᵢⱼ[1,2]]

    @test maximum(abs.(res1 .- res2)) <=1e-12
end
# @testset "residual function DiffEq form" begin
#     # parameters initialization
#     σ₃ = -1e6
#     S_i = 1.0
#     Dⁱ = 0.5
#     Dᵒ = 0.4
#     θ = 60.0
#     ϵ̇ⁱξη = 0.0
#     op = DSB.OutputParams()
#     sp = DSB.SolverParams()
#     mp = DSB.Params(op,sp)
#     r = DSB.Rheology()

#     # parameters
#     p = DSB.DiffEqParams(r = r,
#                     mp = mp,
#                     du = fill(ϵ̇ⁱξη,6),
#                     scalars = SA[σ₃, ϵ̇ⁱξη, θ],
#                     allow_Ḋᵒ = true)
#     σᵒᵢⱼ, σⁱᵢⱼ, ϵᵒᵢⱼ, ϵⁱᵢⱼ = DSB.initialize_state_var_D(r,mp,S_i,σ₃,Dⁱ,Dᵒ,θ ; coords=:band)
#     u0 = SA_F64[S_i, σⁱᵢⱼ[1,1], σⁱᵢⱼ[3,3], ϵⁱᵢⱼ[2,2], Dᵒ, Dⁱ]
#     du_nl = SA_F64[0.0, 0.0, 0.0, 0.0]

#     KIᵒ = DSB.compute_KI(r,σᵒᵢⱼ,Dᵒ)
#     KIⁱ = DSB.compute_KI(r,σⁱᵢⱼ,Dⁱ)
#     Ḋᵒ = DSB.compute_subcrit_damage_rate(r, KIᵒ, Dᵒ)
#     Ḋⁱ = DSB.compute_subcrit_damage_rate(r, KIⁱ, Dⁱ)

#     res  = DSB.residual_2_points_2(du_nl, σᵒᵢⱼ, σⁱᵢⱼ, Dⁱ, Ḋⁱ, Dᵒ, Ḋᵒ, p)
#     @test norm(res) ≈ 0
# end


text:
function getEmbeddedIn3D(
            unitcell :: U,
            z :: Real = 0.0
        ) :: Unitcell{Site{LS,3}, B} where {LS,B,S<:AbstractSite{LS,2},U<:AbstractUnitcell{S,B}}

    # build a new list of sites
    sites_new = Site{LS,3}[
        newSite(Site{LS,3}, [point(s)[1], point(s)[2], z], label(s)) for s in sites(unitcell)
    ]

    # create and return a new Unitcell
    return newUnitcell(
        Unitcell{Site{LS,3}, B},
        latticeVectors(unitcell),
        sites_new,
        bonds(unitcell)
    )
end

function getEmbeddedIn3D(
            lattice :: L,
            z :: Real = 0.0
        ) :: Lattice{Site{LS,3}, B, U} where {LS,B,S<:AbstractSite{LS,2},U,L<:AbstractLattice{S,B,U}}

    # build a new list of sites
    sites_new = Site{LS,3}[
        newSite(Site{LS,3}, [point(s)[1], point(s)[2], z], label(s)) for s in sites(lattice)
    ]

    # create and return a new Unitcell
    return newLattice(
        Lattice{Site{LS,3}, B, U},
        latticeVectors(lattice),
        sites_new,
        bonds(lattice),
        unitcell(lattice)
    )
end

export getEmbeddedIn3D


text:
## subsystem initialization

export Shutdown, ismultithreaded

ismultithreaded() = convert(Core.Bool, API.LLVMIsMultithreaded())

for subsystem in [:Core, :TransformUtils, :ScalarOpts, :ObjCARCOpts, :Vectorization,
                  :InstCombine, :IPO, :Instrumentation, :Analysis, :IPA, :CodeGen, :Target]
    jl_fname = Symbol(:Initialize, subsystem)
    api_fname = Symbol(:LLVM, jl_fname)
    @eval begin
        export $jl_fname
        $jl_fname(R::PassRegistry) = API.$api_fname(R)
    end
end

## back-end initialization

export backends

const libllvm_backends = [:AArch64, :AMDGPU, :ARC, :ARM, :AVR, :BPF, :Hexagon, :Lanai,
                          :MSP430, :Mips, :NVPTX, :PowerPC, :RISCV, :Sparc, :SystemZ,
                          :VE, :WebAssembly, :X86, :XCore]
const libllvm_components = [:Target, :TargetInfo, :TargetMC, :AsmPrinter, :AsmParser, :Disassembler]

# discover supported back-ends and their components by looking at available symbols.
# this reimplements LLVM macros and `static inline` functions that are hard to call.
Libdl.dlopen(libllvm) do library
    supported_backends = filter(libllvm_backends) do backend
        initializer = "LLVMInitialize$(backend)Target"
        Libdl.dlsym(library, initializer; throw_error=false) !== nothing
    end
    @eval backends() = $supported_backends

    # generate subsystem initialization routines for every back-end
    supported_components = Dict(component => [] for component in libllvm_components)
    for backend in libllvm_backends
        backend_supported = backend in supported_backends

        for component in libllvm_components
            jl_fname = Symbol(:Initialize, backend, component)
            @eval export $jl_fname

            api_fname = Symbol(:LLVM, jl_fname)
            if backend_supported
                component_supported = Libdl.dlsym(library, api_fname; throw_error=false) !== nothing
                if component_supported
                    push!(supported_components[component], backend)
                    @eval $jl_fname() = ccall(($(QuoteNode(api_fname)),libllvm), Cvoid, ())
                else
                    @eval $jl_fname() = error($"The $backend back-end does not contain a $component component.")
                end
            else
                @eval $jl_fname() = error($"The $backend back-end is not part of your LLVM library.")
            end
        end
    end

    # same, for initializing subsystems for all back-ends at once
    for component in keys(supported_components)
        jl_fname = Symbol(:Initialize, :All, component, :s)
        exprs = Expr[]
        for backend in supported_components[component]
            fname = Symbol(:Initialize, backend, component)
            push!(exprs, :($fname()))
        end

        @eval begin
            export $jl_fname
            function $jl_fname()
                $(exprs...)
                return
            end
        end
    end
end

# same, for the native back-end
for component in [:Target, :AsmPrinter, :AsmParser, :Disassembler]
    jl_fname = Symbol(:Initialize, :Native, component)
    api_fname = Symbol(:LLVM, jl_fname)
    @eval begin
        export $jl_fname
        $jl_fname() = convert(Core.Bool, API.$api_fname()) &&
                      throw(LLVMException($"Could not initialize native $component"))
    end
end


text:
module DeepGP

using LinearAlgebra
using KernelFunctions
using ZygoteRules
using AbstractGPs: AbstractGPs, AbstractGP, GP
using Distributions
using Statistics
using StatsBase
using ParameterHandling
using PDMats

export mean, cov, var
export predict

include("utils.jl")
include("layers.jl")
include("deepgp_base.jl")
include("training.jl")
end # module


text:
using MyExample
using Test

@testset "f_tests" begin
    @test my_f(2,1) == 7
    @test my_f(2,3) == 13
    @test my_f(1,3) == 11
    @test my_f(3,3) == 15
end

@testset "Derivative tests" begin
    @test derivative_of_my_f(2,1) == 2
end


text:
#=――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
    NOTE
    The algorithm used here only works up to vigintillions. Beyond that, it is
    difficult to find consensus on how numbers should be displayed. It is also
    of dubious real-world utility.
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――=#

module Numeric

export english, unenglish

const ONES = [
    "one", "two", "three", "four", "five", "six", "seven", "eight", "nine"]
const TEENS = [
    "eleven", "twelve", "thirteen", "fourteen", "fifteen", "sixteen",
    "seventeen", "eighteen", "nineteen"]
const TENS = [
    "ten", "twenty", "thirty", "forty", "fifty", "sixty", "seventy", "eighty",
    "ninety"]
const POWERS = [
    "thousand", "million", "billion", "trillion", "quadrillion", "quintillion",
    "sextillion", "septillion", "octillion", "nonillion", "decillion",
    "undecillion", "duodecillion", "tredecillion", "quattuordecillion",
    "quindecillion", "sexdecillion", "septendecillion", "octodecillion",
    "novemdecillion", "vigintillion"]

function reverse_word_lookup(ones, teens, tens)
    result = Dict{String, Int}()
    for (i, o) in enumerate(ones)
        result[o] = i
    end
    for (i, o) in enumerate(teens)
        result[o] = 10 + i
    end
    for (i, o) in enumerate(tens)
        result[o] = i * 10
    end
    result
end

function reverse_power_lookup(powers)
    result = Dict{String, Int}()
    for (i, o) in enumerate(powers)
        result[o] = 3i
    end
    result
end

const REVERSE_WORD = reverse_word_lookup(ONES, TEENS, TENS)
const REVERSE_POWER = reverse_power_lookup(POWERS)

"""
Convert ``n`` to English, given that ``0 < n < 100``.
"""
function lth(n::Integer)
    if n < 10
        ONES[n]
    elseif n % 10 == 0
        TENS[n ÷ 10]
    elseif n < 20
        TEENS[n % 10]
    else
        string(TENS[n ÷ 10], '-', ONES[n % 10])
    end
end

"""
Convert ``n`` to English, given that ``0 < n < 1000``.
"""
function ltt(n::Integer)
    if n < 100
        lth(n)
    elseif n % 100 == 0
        string(ONES[n ÷ 100], " hundred")
    else
        string(ONES[n ÷ 100], " hundred ", lth(n % 100))
    end
end

"""
    english(n::Integer)

Convert ``n`` to English, given that ``0 \\le n < 10^{66}``.

```jldoctest
julia> using EnglishText

julia> english(16)
"sixteen"
```
"""
function english(n::Integer)
    if n == 0
        "zero"
    else
        result = String[]
        segments = digits(n, base=1000)
        if segments[1] ≠ 0
            push!(result, ltt(segments[1]))
        end
        for (i, part) in enumerate(segments[2:end])
            if part ≠ 0
                push!(result, POWERS[i], ltt(part))
            end
        end
        join(reverse(result), " ")
    end
end

"""
    unenglish(T <: Integer, data::AbstractString) → T

Convert `data` to an integral type. This function has the guarantee that
`unenglish(Int, english(x)) == x`, modulo any type differences. It is not
guaranteed to work well or throw exceptions on other inputs.

```jldoctest
julia> using EnglishText

julia> unenglish(Int, "sixteen")
16
```
"""
function (unenglish(::Type{T}, data::AbstractString)::T) where T <: Integer
    words = split(data)
    bigpart::T = zero(T)
    smallpart::T = zero(T)
    for word in words
        pieces = split(word, '-')
        for piece in pieces
            piece = lowercase(piece)
            if haskey(REVERSE_WORD, piece)
                smallpart += REVERSE_WORD[piece]
            elseif piece == "hundred"
                smallpart *= 100
            elseif haskey(REVERSE_POWER, piece)
                bigpart += T(10) ^ REVERSE_POWER[piece] * smallpart
                smallpart = zero(T)
            end
        end
    end
    bigpart + smallpart
end

end  # module Numeric

import .Numeric: english, unenglish
export english, unenglish


text:
using HTTP
using JSON

function anki_suspend_card(card_id)
	if length(card_id) == 1
		card_id = "[$card_id]"
	end

	response = HTTP.request(:POST, "http://localhost:8765", [], """{
				\"action\": \"suspend\",
				\"version\": 6,
				\"params\": {
					\"cards\": $card_id
					}
				}""")
	JSON.parse(String(response.body))
end