text:
#   This file is part of Reduce.jl. It is licensed under the MIT license
#   Copyright (C) 2017 Michael Reed

export unfold, mat, treecombine!, irr

"""
REDUCE begin and end marker counter for parsegen
"""
@inline function becount(js,openpar,closepar)
    if typeof(openpar) == String
        sh = split(js,r"[ ]+")
        sum(sh.==openpar)-sum(sh.==closepar)
    else
        count(z->z==openpar,js)-count(z->z==closepar,js)
    end
end

"""
REDUCE begin and end marker match for parsegen
"""
@noinline function bematch(js,sexpr,iter,next,openpar,closepar)
    (h,state) = next
    y = h
    c = becount(js,openpar,closepar)
    flag = c > 0
    nxt = iterate(iter, state)
    while (nxt !== nothing) & flag
        (y,state) = nxt
        c += becount(sexpr[y],openpar,closepar)
        flag = c > 0
        flag && (nxt = iterate(iter, state))
    end
    return (y,state)
end

const prefix = r"(?<!\))(([A-Za-z_][A-Za-z_0-9]*)|([\^+\/-])|([*]{1,2})|(- ))(?=\()"
const parens = r"\(((?>[^\(\)]+)|(?R))*\)"
const braces = r"{((?>[^{}]+)|(?R))*}"
const infix1 = r"^(([\^\+\/])|([*]{1,2})|( -)|( \+)|( [*]{1,2})|( /)|( \^))"
const infix2 = r"(([\^+\/])|([*]{1,2}))$"
const assign = r"^([A-Za-z_ ][A-Za-z_0-9 ]*)(:=)"

@inline function argrfun(mode::Symbol,rfun::Symbol,sep,be=:be)
    if mode == :expr
        :($rfun(fun,$sep;be=$be))
    elseif mode == :args
        :($rfun(fun,$sep,s;be=$be) |> string)
    else
        :($rfun(fun,$sep;be=$be) |> string)
    end
end

function loopshift(js,openpar,closepar,T,sexpr,iter,next)
    ep = Array{T,1}(undef,0)
    c = becount(js,openpar,closepar)
    flag = c ≥ 0
    !flag && (js = join(split(js,closepar)[1:end+c],closepar))
    (h,state) = next
    y = h
    nxt = bematch(js,sexpr,iter,next,openpar,closepar)
    (h,state) = nxt
    push!(ep,js,sexpr[y+1:h]...)
    ep[1] == nothing && popfirst!(ep)
    nxt = iterate(iter, state)
    while (nxt !== nothing) & flag
        (h,state) = nxt
        cQ = c
        js = sexpr[h]
        c += becount(js,openpar,closepar)
        if c ≤ -1
            js = join(split(js,closepar)[1:end+c],closepar)
            flag = false
        end
        y = h
        nxt = bematch(js,sexpr,iter,nxt,openpar,closepar)
        (h,state) = nxt
        epr = vcat(js,sexpr[y+1:h]...)
        epr ≠ nothing && push!(ep,epr...)
        flag && (nxt = iterate(iter, state))
    end
    return ((h,state),ep)
end

for mode ∈ [:expr,:unary,:switch,:args]
    rfun = Symbol(:r,"_",mode)
    modefun = Symbol(:parse,"_",mode)
    argfun = Symbol(:args,"_",mode)
    arty = (mode == :expr) ? :Any : :String
    exec = if mode == :expr
        :(Meta.parse(RSymReplace(js)))
    elseif mode == :unary
        :("$fun($js)" |> RExpr |> rcall)
    elseif mode == :switch
        :(rcall("$js" |> RExpr, fun))
    elseif mode == :args
        :("$fun($(join([js,s...],',')))" |> RExpr |> rcall)
    end
    mode != :args ? (fargs = [:(r::RExpr)]) : (fargs = [:(r::RExpr),Expr(:...,:sr)])
    mode != :args ? (aargs = [:be]) : (aargs = [:be,Expr(:...,:s)])
    quote
        export $modefun
        @noinline function $modefun(fun::String,$(fargs...);be=0)
            nsr = $arty[]
            $(if mode == :args
                quote
                    s = Array{String,1}()
                    for rs ∈ sr
                        push!(s,rs |> RExpr |> string)
                    end
                end
            end)
            sexpr = split(r).str
            iter = 1:length(sexpr)
            next = iterate(iter); #show(sexpr)
            while next !== nothing
                (h,state) = next
                sh = split(sexpr[h],r"[ ]+")
                en = 1
                isempty(replace(sh[en]," "=>"")) && (en = 2); #show(sh[en])
                if occursin(r"^procedure",sh[en])
                    js = join(split(sexpr[h],"procedure")[2:end],"procedure")
                    next = iterate(iter, state)
                    (h,state) = next
                    y = h
                    next = bematch(sexpr[h],sexpr,iter,next,"begin","end")
                    (h,state) = next
                    $(mode != :expr ? :(push!(nsr,String("procedure "*js))) : :(nothing))
                    push!(nsr,$(if mode == :expr
                        :(Expr(:function,Meta.parse(js),$rfun(fun,sexpr[y:h];be=be)))
                    elseif mode == :args
                        :($rfun(fun,sexpr[y:h],s;be=be) |> string)
                    else
                        :($rfun(fun,sexpr[y:h];be=be) |> string)
                    end))
                elseif occursin(r"^begin",sh[en])
                    js = join(split(sexpr[h],"begin")[2:end],"begin")
                    ep = Array{$arty,1}(undef,0)
                    c = becount(js,"begin","end")
                    flag = c ≥ 0
                    !flag && (js = join(split(js,"end")[1:end+c],"end"))
                    y = h
                    next = bematch(js,sexpr,iter,next,"begin","end")
                    (h,state) = next
                    $(mode != :expr ? :(push!(nsr,String("begin "*js))) : :(nothing))
                    push!(ep,$(argrfun(mode,rfun,:(vcat(js,sexpr[y+1:h]...)),:(be+1))))
                    ep[1] == nothing && popfirst!(ep)
                    next = iterate(iter, state)
                    while (next !== nothing) & flag
                        (h,state) = next
                        cQ = c
                        js = sexpr[h]
                        c += becount(js,"begin","end")
                        if c ≤ -1
                            js = join(split(js,"end")[1:end+c],"end")
                            flag = false
                        end
                        y = h
                        next = bematch(js,sexpr,iter,next,"begin","end")
                        (h,state) = next
                        epr = $(argrfun(mode,rfun,:(vcat(js,sexpr[y+1:h]...)),:cQ))
                        epr ≠ nothing && push!(ep,epr)
                        flag && (next = iterate(iter, state))
                    end
                    next = (h,state)
                    push!(nsr,$((mode == :expr) ? :(Expr(:block,ep...)) : Expr(:...,:ep)))
                    $(mode != :expr ? :(push!(nsr,String("end"))) : :(nothing))
                elseif occursin(r"^return",sh[en])
                    js = join(split(sexpr[h],"return")[2:end],"return")
                    y = h
                    next = bematch(js,sexpr,iter,next,"begin","end")
                    (h,state) = next
                    rp = $(argrfun(mode,rfun,:(vcat(js,sexpr[y+1:h]...))))
                    $(mode != :expr ? :(push!(nsr,"return "*rp)) : :(push!(nsr,Expr(:return,rp))))
                elseif occursin(assign,sexpr[h])
                    sp = split(sexpr[h], ":=",limit=2)
                    push!(nsr,$(if mode == :expr
                        :(Expr(:(=),Meta.parse(sp[1]),$rfun(fun,sp[2];be=be)))
                    elseif mode == :args
                        :(String(sp[1]) * ":=" * string($rfun(fun,sp[2] |> String |> RExpr,s;be=be)))
                    else
                        :(String(sp[1]) * ":=" * string($rfun(fun,sp[2] |> String |> RExpr;be=be)))
                    end))
                elseif occursin("for",sh[en])
                    throw(ReduceError("for block parsing not yet supported"))
                elseif occursin(braces,$((mode == :expr) ? :(sexpr[h]) : :("")))
                    $(if mode == :expr; quote
                        ts = sexpr[h]
                        (next,mp) = loopshift(ts,'{','}',String,sexpr,iter,next)
                        (h,state) = next
                        smp = match(braces,join(mp,";\n")).match[2:end-1]
                        ListPrint(ListPrint()+1)
                        args = Array{$arty,1}(undef,0)
                        while smp ≠ ""
                            lsM = match(braces,smp)
                            lsm = lsM ≠ nothing ? lsM.match[2:end-1] : smp
                            lsd = split(smp,braces;limit=2)
                            pre = length(lsd)≠1 ? lsd[1][1:end-1] : ""
                            lsM == nothing && (pre = join([pre,lsm],','))
                            if pre ≠ ""
                                af = $argfun(fun,$(string(mode)),split(pre,','),be)
                                for k ∈ af
                                    k≠nothing && push!(args, k≠[nothing] ? k : Array{Any,1}(undef,0))
                                end
                            end
                            smp = length(lsd)≠1 ? lsd[end][2:end] : ""
                            if lsM ≠ nothing
                                af = $(argrfun(mode,rfun,:(lsM.match)))
                                push!(args, af≠[nothing] ? af : Array{Any,1}(undef,0))
                            end
                        end
                        ListPrint(ListPrint()-1)
                        length(args)==1 && typeof(args[1]) <: Tuple && (args = args[1])
                        push!(nsr,Tuple(args))
                    end; else; :(nothing); end)
                elseif occursin("=>",$((mode == :expr) ? :(sexpr[h]) : :("")))
                    sp = split(sexpr[h],r"=>")
                    stuff = String.(split(sp[2],r"(when)|(and)"))
                    ar = $rfun.(fun,stuff;be=be)
                    push!(nsr,$rfun(fun,sp[1];be=be) => length(ar) ≠ 1 ? ar : ar[1])
                elseif occursin(prefix,$((mode == :expr) ? :(sexpr[h]) : :("")))
                    $(if mode == :expr; quote
                    ts = sexpr[h]
                    (next,mp) = loopshift(ts,'(',')',String,sexpr,iter,next)
                    (h,state) = next
                    smp = replace(join(mp,";\n"),"**"=>'^')
                    qr = IOBuffer()
                    while smp ≠ ""
                        args = Array{$arty,1}(undef,0)
                        if occursin(r"^\s?\(",smp)
                            push!(args,$(argrfun(mode,rfun,:(match(parens,smp).match[2:end-1]))))
                            smp = split(smp,parens;limit=2)[end]
                            print(qr,"($(args...))")
                            if !occursin(prefix,smp)
                                $(if mode == :expr; quote
                                    if occursin(infix1,smp)
                                        print(qr, RSymReplace(match(infix1,smp).match), RSymReplace(split(smp,infix1)[end]))
                                        smp = ""
                                    else
                                        print(qr, RSymReplace(smp))
                                        smp = ""
                                    end; end
                                else
                                    :(print(qr, smp); smp = "")
                                end)
                            elseif occursin(infix1,smp)
                                print(qr, RSymReplace(match(infix1,smp).match))
                                smp = split(smp,infix1)[end]
                            end
                            continue
                        end
                        pf = match(prefix,smp).match |> String
                        sp = split(smp,prefix;limit=2)
                        $(if mode == :expr; quote
                            if occursin(infix2,sp[1])
                                rq = split(sp[1],infix2)[1]
                                if occursin(infix1,rq)
                                    rq = RSymReplace(match(infix1,rq).match) * RSymReplace(split(rq,infix1)[end])
                                else
                                    rq = RSymReplace(rq)
                                end
                                print(qr, rq, RSymReplace(match(infix2,sp[1]).match))
                            elseif occursin(infix1,sp[1])
                                print(qr, RSymReplace(match(infix1,sp[1]).match), RSymReplace(split(sp[1],infix1)[end]))
                            else
                                print(qr, RSymReplace(sp[1]))
                            end; end
                        else
                            :(print(qr, sp[1]))
                        end)
                        smp = split(sp[2],parens;limit=2)[end]
                        lsm = match(parens,sp[2]).match[2:end-1]
                        if pf == "mat"
                            mt = collect((m.match for m=eachmatch(parens,lsm)))
                            for row ∈ mt
                                elm = split(row[2:end-1],',')
                                push!(args,$(if mode == :expr
                                    :(Expr(:row,$argfun(fun,$(string(mode)),elm,be)...))
                                elseif mode == :args
                                    :($argfun(fun,$(string(mode)),elm,be,s) |> string)
                                else
                                    :($argfun(fun,$(string(mode)),elm,be) |> string)
                                end))
                            end
                            print(qr, $(if mode == :expr
                                :(Expr(:vcat,args...) |> string)
                            else
                                :("$pf($(join(args,',')))")
                            end))
                        else
                            ls = split(lsm,',')
                            push!(args,$(if mode == :expr
                                :($argfun(fun,$(string(mode)),ls,be))
                            elseif mode == :args
                                :($argfun(fun,$(string(mode)),ls,be,s) |> string)
                            else
                                :($argfun(fun,$(string(mode)),ls,be) |> string)
                            end)...)
                            rq = "$(RSymReplace(pf))($(join(args,',')))"
                            print(qr, $(if mode == :expr
                                :(((isinfix(pf) && length(args) == 1) ? rq : "($rq)"))
                            else
                                :("$pf($(join(args,',')))")
                            end))
                        end
                        !occursin(prefix,smp) && ($(if mode == :expr; quote
                            if occursin(infix1,smp)
                                print(qr, RSymReplace(match(infix1,smp).match), RSymReplace(split(smp,infix1)[end]))
                                smp = ""
                            else
                                print(qr, RSymReplace(smp))
                                smp = ""
                            end; end
                        else
                            :(print(qr, smp); smp = "")
                        end))
                    end
                    push!(nsr,$((mode == :expr) ? :("("*String(take!(qr))*")" |> Meta.parse |> linefilter!) : :qr))
                    end; else; :(nothing); end)
                elseif occursin("end",sh[en])
                    nothing
                elseif isempty(sh[en])
                    nothing
                elseif occursin("=",$((mode == :expr) ? :(sexpr[h]) : :("")))
                    $(if mode == :expr; quote
                    sp = split(sexpr[h],"=")
                    push!(nsr,$(:(Expr((ListPrint()>0 ? (:(=),) : (:call,:(==)))...,
                            $rfun(fun,sp[1];be=be),$rfun(fun,sp[2];be=be)))))
                    end; end)
                elseif occursin(":",sexpr[h])
                    sp = split(sexpr[h],":")
                    push!(nsr,$(if mode == :expr
                        :(Expr(:(:),$rfun(fun,sp[1];be=be),$rfun(fun,sp[2];be=be)))
                    elseif mode == :args
                        :(($rfun(fun,sp[1],[];be=be)|>string)*":"*($rfun(fun,sp[2],s;be=be)|>string))
                    else
                        :(($rfun(fun,sp[1];be=be)|>string)*":"*($rfun(fun,sp[2];be=be)|>string))
                    end))
                else
                    js=sexpr[h]
                    se=sum(sh.=="end")
                    0<se≤be ? (js=replace(js,"end","")) :
                        (se>be && (js=join(split(js,"end")[1:end-be],"end")))
                    exc = $exec
                    $(if mode ≠ :expr
                        :(if string(exc) == ""
                            c = 1
                            f = SubFail()
                            H = SubHold()
                            while string(exc) == "" && c < f
                                sleep(sqrt(c)*H)
                                exc = $exec
                                c += 1
                            end
                            PipeClogged(string(exc) ≠ "", c, "$fun function")
                            string(exc) == "" && throw(ReduceError(if fun == "//"
                                "If generated code has many calls to $fun, try setting `Reduce.Rational(false)` in code, since rational division is the default; or try `Reduce.Reset()`."
                            else
                                "If generated code has many calls to $fun, try to minimize the number of calls with REDUCE switches and use `Reduce.Reset()` if you'd like to start a new pipe."
                            end))
                        end)
                    else; nothing
                    end)
                    push!(nsr, exc)
                end
                next = iterate(iter, state)
            end
            $(if mode == :expr
                quote
                    u = length(nsr)
                    return u==1 ? nsr[1] : (u==0 ? nothing : Expr(:block,nsr...))
                end
            else
                :(if fun in ["nat","latex"]
                    return nsr |> RExpr |> split |> string
                else
                    return nsr |> RExpr |> split
                end)
            end)
        end
        $(if mode == :args
            quote
                $rfun(fun::String,r::Array{String,1},s;be=0) = $modefun(fun,RExpr(r),s...;be=be)
                $rfun(fun::String,r,s;be=0) = $modefun(fun,r |> String |> RExpr,s...;be=be)
            end
        else
            quote
                $rfun(fun::String,r::Array{String,1};be=0) = $modefun(fun,RExpr(r);be=be)
                $rfun(fun::String,r;be=0) = $modefun(fun,r |> String |> RExpr;be=be)
            end
        end)
        @noinline function $argfun(fun::String,mod::String,ls::Array{SubString{String},1},$(aargs...))
            mode = Symbol(mod)
            args = Array{$arty,1}(undef,0)
            lsi = 1:length(ls)
            nxt = iterate(lsi)
            (lsh,lss) = nxt
            while nxt !== nothing
                (lsh,lss) = nxt
                if occursin(r"^begin",ls[lsh])
                    js = join(split(ls[lsh],"begin")[2:end],"begin")
                    (nxt,ep) = loopshift(js,"begin","end",$arty,ls,lsi,nxt)
                    (lsh,lss) = nxt
                    sep = "begin $(join(ep,',')) end"
                    push!(args,$(argrfun(mode,rfun,:sep)))
                elseif occursin(prefix,ls[lsh])
                    js = ls[lsh]
                    ep = Array{$arty,1}(undef,0)
                    c = count(z->z=='(',js)-count(z->z==')',js)-1
                    flag = c ≥ -1
                    !flag && (js = join(split(js,')')[1:end+c],')'))
                    lsy = lsh
                    nxt = bematch(js,ls,lsi,nxt,'(',')')
                    (lsh,lss) = nxt
                    push!(ep,js,ls[lsy+1:lsh]...)
                    ep[1] == nothing && popfirst!(ep)
                    sep = join(ep,',')
                    push!(args,$(argrfun(mode,rfun,:sep)))
                else
                    push!(args,$(argrfun(mode,rfun,:(ls[lsh]))))
                end
                nxt = iterate(lsi, lss)
            end
            return args
        end
    end |> eval
end

"""
    parsegen(::Symbol,::Symbol)

Parser generator that outputs code to walk and manipulate REDUCE expressions
"""
function parsegen(fun::Symbol,mode::Symbol)
    fune = (fun == ://) ? :/ : (fun == :rlet) ? :let : fun
    mf = Symbol(:parse,"_",mode)
    a = mode != :args ? [:(r::RExpr)] : [:(r::RExpr),Expr(:...,:s)]
    return mode ≠ :expr ? :($fun($(a...);be=0) = $mf($(string(fune)),$(a...);be=0)) :
    :($fun($(a...);be=0) = $mf($(string(fune)),$(a...);be=0) |> treecombine! |> irr)
end

"""
    treecombine!(::Expr)

Recursively simplifies out extra edges from Expr objects
"""
@noinline function treecombine!(e::Expr,redo=[false])
    for i ∈ 1:length(e.args)
        if e.args[i] |> typeof == Expr && e.args[i].head == :call
            if e.head == :call
                s = e.args[i].args[1]
                if s ∈ [://,:/]
                    if e.args[1] == :*
                        d = e.args[i].args[3]
                        e.args[i] = e.args[i].args[2]
                        e.args = [s,deepcopy(e),d]
                        redo[1] = true
                        return treecombine!(e)
                    elseif e.args[1] ∈ [://,:/]
                        if i == 2
                            e.args[3] = e.args[3]*e.args[i].args[3]
                            e.args[2] = deepcopy(e.args[i].args[2])
                        elseif i == 3
                            e.args[2] = e.args[2]*e.args[i].args[3]
                            e.args[3] = deepcopy(e.args[i].args[2])
                        end
                        redo[1] = true
                        return treecombine!(e)
                    end
                elseif s == :-
                    if e.args[1] == :-
                        if length(e.args) == 2
                            if length(e.args[i].args) == 3
                                push!(e.args,deepcopy(e.args[i].args[2]))
                                e.args[2] = deepcopy(e.args[i].args[3])
                                redo[1] = true
                                return treecombine!(e)
                            end
                        #elseif i == 2 && length(e.args) > 2
                        #    push!(e.args[i].args,e.args[3:end]...)
                        #    e.args = deepcopy(e.args[i].args)
                        #    redo[1] = true
                        #    return treecombine!(e)
                        end
                    end
                end
            elseif e.head == :(=) &&
                    e.args[2] |> typeof == Expr && e.args[2].head == :block
                    length(e.args[2].args) == 1 && (e.args[2] = e.args[2].args[1])
            end
            treecombine!(e.args[i],redo)
            d = detectinf(e)
            d ≠ nothing && (e.args[i] = d)
        else
            typeof(e.args[i]) == Expr && treecombine!(e.args[i],redo)
        end
    end
    d = detectinf(e)
    return redo[1] ? treecombine!(e) : d ≠ nothing ? d : e
end
treecombine!(e,redo=[false]) = e
treecombine!(e::T) where T <: Tuple = treecombine!.(e)

function detectinf(e)
    if typeof(e) == Expr && e.head == :call
        if e.args[1] ∈ [://,:/]
            if (e.args[2] == :NaN) | (e.args[3] == :NaN)
                return :NaN
            elseif e.args[2] == :Inf
                if typeof(e.args[3]) <: Number
                    if e.args[3] == :Inf
                        return :NaN
                    else
                        return :Inf
                    end
                elseif (typeof(e.args[3]) == Expr) && (e.args[3].head == :macrocall)
                    if e.args[3].args[1] == Symbol("@big_str") || e.args[3].args[1] == Symbol("@int128_str")
                        return :Inf
                    end
                end
            elseif e.args[3] == :Inf
                if typeof(e.args[2]) <: Number
                    return 0
                elseif (typeof(e.args[2]) == Expr) && (e.args[2].head == :macrocall)
                    if e.args[2].args[1] == Symbol("@big_str") || e.args[2].args[1] == Symbol("@int128_str")
                        return 0
                    end
                end
            end
        elseif e.args[1] ∈ [:*,:+,:-]
            found = false
            for arg ∈ e.args[2:end]
                arg == :NaN && (found = true; return :NaN)
            end
        end
    end
    return nothing
end

@inline irr(expr) = expr ∈ [:ℯ,:π,:(MathConstants.γ),:(MathConstants.φ),:Inf,:NaN] ? eval(expr) : expr

parsegen(:parse,:expr) |> eval

@doc """
    Reduce.parse(r::RExpr)

Parse a Reduce expression into a Julia expression

# Examples
```julia-repl
julia> Reduce.parse(R\"sin(i*x)\")
:(sin(im * x))
```
""" Reduce.parse

function print_args(io::IO,a::Array{Any,1})
    print(io, "(")
    for (i, arg) in enumerate(a)
        show_expr(io, arg)
        i ≠ lastindex(a) ? print(io,",") : print(io,")")
    end
end

@noinline function show_expr(io::IO, expr::Expr) # recursively unparse Julia expression
    if expr.head == :call
        if expr.args[1] == :(:)
            show_expr(io,expr.args[2])
            print(io,":")
            show_expr(io,expr.args[3])
            print(io," ")
        elseif expr.args[1] == :(==)
            show_expr(io,expr.args[2])
            print(io,"=")
            show_expr(io,expr.args[3])
        else
            show_expr(io, expr.args[1])
            print_args(io,expr.args[2:end])
        end
    elseif expr.head == :(=)
        if (typeof(expr.args[1]) == Expr) && (expr.args[1].head == :call) && ListPrint()<1
            show_expr(io,Expr(:function,expr.args[1],expr.args[2]))
        else
            show_expr(io,expr.args[1])
            print(io,ListPrint()>0 ? "=" : ":=")
            show_expr(io,expr.args[2])
        end
    elseif occursin(r"[*\/+-^]=$",string(expr.head))
        if (typeof(expr.args[1]) == Expr) && (expr.args[1].head == :call)
            throw(ReduceError("function assignment for $(expr.head) not supported"))
        else
            show_expr(io,expr.args[1])
            print(io,ListPrint()>0 ? "=" : ":=")
            print(io,match(r"[^(=$)]",string(expr.head)).match*"(")
            show_expr(io,expr.args[1])
            print(io,",")
            show_expr(io,expr.args[2])
            print(io,")")
        end
    elseif expr.head == :for
        print(io,"for ")
        show_expr(io,expr.args[1])
        show_expr(io,expr.args[2])
    elseif expr.head == :block
        lxpr = linefilter!(expr)
        if length(lxpr.args) == 1
            show_expr(io,lxpr.args[1])
        else
            print(io,"begin ")
            show_expr(io,lxpr.args[1])
            for k ∈ 2:length(lxpr.args[2:end])+1
                print(io,";")
                show_expr(io,lxpr.args[k])
            end
            print(io," end")
        end
    elseif expr.head == :function
        print(io,"procedure ")
        show_expr(io,expr.args[1])
        print(io,";")
        show_expr(io,expr.args[2])
    elseif expr.head == :return
        print(io,"return ")
        show_expr(io,expr.args[1])
    elseif expr.head == :(::)
        show_expr(io,expr.args[1])
    elseif expr.head == :macrocall
        if expr.args[1] ∈ (Symbol("@big_str"),Symbol("@int128_str"),GlobalRef(Core,Symbol("@big_str")),GlobalRef(Core,Symbol("@int128_str")))
            print(io,expr.args[end])
        else
            throw(ReduceError("Macro $(expr.args[1]) block structure not supported\n\n$expr"))
        end
    elseif expr.head == :(:)
        show_expr(io,expr.args[1])
        print(io,":")
        show_expr(io,expr.args[2])
        print(io," ")
    elseif expr.head == :vcat
        print(io,"mat((")
        for i ∈ 1:length(expr.args)-1
            show_expr(io,expr.args[i])
            print(io,"),(")
        end
        show_expr(io,expr.args[end])
        print(io,"))")
    elseif expr.head == :row
        for i ∈ 1:length(expr.args)-1
            show_expr(io,expr.args[i])
            print(io,",")
        end
        print(io,expr.args[end])
    elseif expr.head == :tuple
        ListPrint(ListPrint()+1)
        print(io,"{")
        l = length(expr.args)
        for i ∈ 1:l
            show_expr(io,expr.args[i])
            i ≠ l && print(io,",")
        end
        print(io,"}")
        ListPrint(ListPrint()-1)
    elseif expr.head == :.
        if expr.args[1] == :MathConstants
            print(io,JSymReplace(string(expr)))
        else
            throw(ReduceError("$(expr.args[1]) module scope not supported"))
        end
    elseif expr.head == :line; nothing
    else
        throw(ReduceError("Nested :$(expr.head) block structure not supported\n\n$expr"))
    end
end

const infix_ops = ["+", "-", "*", "/", "**", "^"]
isinfix(args) = replace(args,' '=>"") in infix_ops

function show_expr(io::IO, ex)
    ((ex == :nothing) | (typeof(ex) == LineNumberNode)) && (return nothing)
    if typeof(ex) <: AbstractFloat && isinf(ex)
        print(io,(r > 0 ? "" : "-")*"infinity")
        return nothing
    elseif typeof(ex) <: Irrational
        print(io,unparse_irrational(ex))
        return nothing
    elseif typeof(ex) <: Complex{Bool}
        print(io,"i")
    elseif typeof(ex) <: Complex
        show_expr(io,ex.re)
        print(io,"+")
        show_expr(io,ex.im)
        print(io,"*i")
    end
    if typeof(ex) <: Matrix
        print(io, "mat(")
        li = size(ex)[1]
        lj = size(ex)[2]
        for i ∈ 1:li
            print(io, "(")
            for j ∈ 1:lj-1
                show_expr(io,ex[i,j])
                print(io,",")
            end
            show_expr(io,ex[i,lj])
            print(io,")")
            i ≠ li && print(io,",")
        end
        print(io,")")
    elseif typeof(ex) <: Vector
        print(io,"mat(")
        l = length(ex)
        for i ∈ 1:l
            print(io,"(")
            show_expr(io,ex[i])
            print(io,")")
            i ≠ l && print(io,",")
        end
        print(io,")")
    elseif typeof(ex) <: Adjoint
        print(io,"mat((")
        l = length(ex)
        for i ∈ 1:l
            show_expr(io,ex[i])
            i ≠ l && print(io,",")
        end
        print(io,"))")
    elseif typeof(ex) <: Tuple
        ListPrint(ListPrint()+1)
        print(io,"{")
        l = length(ex)
        for i ∈ 1:l
            show_expr(io,ex[i])
            i ≠ l && print(io,",")
        end
        print(io,"}")
        ListPrint(ListPrint()-1)
    elseif typeof(ex) <: Pair
        show_expr(io,ex[1])
        print(io," => ")
        if typeof(ex[2]) <: Array
            show_expr(io,ex[2][1])
            for k ∈ 2:length(ex[2])
                print(io, k ≠ 1 ? " and " : " when ")
                show_expr(io, ex[2][k])
            end
        else
            show_expr(io,ex[2])
        end
    elseif typeof(ex) <: Array
        if length(size(ex)) > 2
            throw(ReduceError("parsing of $(typeof(ex)) not supported."))
        end
    else
        edit = IOBuffer()
        print(edit, ex)
        print(io, edit |> take! |> String |> JSymReplace)
    end
end

@inline function unparse_irrational(ex::T) where T <: Irrational
    if ex == ℯ
        return "e"
    elseif ex == π
        return "pi"
    else
        throw(ReduceError("$(typeof(ex)) not yet supported"))
    end
end

function unparse(expr::Expr)
    str = Array{String,1}(undef,0)
    io = IOBuffer()
    if expr.head == :block
        for line ∈ expr.args # block structure
            show_expr(io,line)
            push!(str,String(take!(io)))
        end
        return rtrim(str)
    else
        show_expr(io, expr)
        return push!(str,String(take!(io)))
    end
end

"""
    unfoldgen(::Symbol,::Symbol)

Parser generator that outputs code to walk and manipulate Julia expressions
"""
function unfoldgen(fun::Symbol,mode::Symbol)
    modefun = Symbol(:parse,"_",mode)
    fargs = if mode != :args
        [:(r::Union{<:Expr,<:Symbol})]
    else
        [:(r::Union{<:Expr,<:Symbol}),Expr(:...,:s)]
    end
    sargs = mode != :args ? [:(RExpr(r))] : [:(RExpr(r)),Expr(:...,:s)]
    return quote
        function $fun($(fargs...))
            if typeof(r) == Symbol
                convert(Expr,$fun($(sargs...)))
            else
                unfold(Symbol.($(string.([mode,fun])))...,$(fargs...))
            end
        end
    end
end

@noinline function unfold_expr(mode::Symbol, fun::Symbol, expr::Expr, s...; force=false)
    #expr = mode == :unary ? squash(ixpr) : ixpr
    force && return unfold_expr_force(mode,fun,expr,s...)
    if expr.head in [:call,:block,:(:)]
        if expr.args[1] == :(==) && fun == :solve
            return solve(RExpr(expr.args[2]),s...) |> parse
        else
            return unfold_expr_force(mode,fun,expr,s...)
        end
    elseif expr.head == :block
        return Expr(expr.head,unfold_expr.(mode,fun,expr.args,s...)...)
    elseif expr.head == :(:)
        return Expr(expr.head,unfold_expr_force.(mode,fun,expr.args,s...)...)
    elseif expr.head == :return
        return Expr(expr.head,unfold_expr.(mode,fun,expr.args,s...)...)
    elseif expr.head == :for
        return Expr(expr.head,expr.args[1],unfold_expr.(mode,fun,expr.args[2:end],s...)...)
    elseif occursin(r"=$",string(expr.head))
        if (typeof(expr.args[1]) == Expr) && (expr.args[1].head == :call)
            return Expr(expr.head,expr.args[1],unfold_expr_force.(mode,fun,expr.args[2:end],s...)...)
        else
            return Expr(expr.head,expr.args[1],unfold_expr(mode,fun,expr.args[2],s...;force=true))
        end
    elseif expr.head == :function
        return Expr(expr.head,expr.args[1],unfold_expr_force.(mode,fun,expr.args[2:end],s...)...)
    elseif expr.head == :(::)
        return Expr(:(::),unfold_expr(mode,fun,expr.args[1],s...),expr.args[2])
    elseif expr.head ∈ [:macrocall,:vcat,:row]
        return unfold_expr_force(mode,fun,expr,s...)
    elseif expr.head == :line
        return nothing
    else
        throw(ReduceError("Nested :$(expr.head) block structure not supported\n\n$expr"))
    end
end

function unfold_expr_force(mode::Symbol, fun::Symbol, ex, s...)
    out = nothing
    (typeof(ex) == LineNumberNode) && (return ex)
    if mode == :unary
        out = parse_unary(string(fun),RExpr(ex);be=0)
    elseif mode == :switch
        out = parse_switch(string(fun),RExpr(ex);be=0)
    elseif mode ==:args
        out = parse_args(string(fun),RExpr(ex),s...;be=0)
    else
        throw(ReduceError("Parse mode not supported."))
    end
    fun ∈ switchtex ? out : typeof(ex) <: Matrix ? mat(out) : convert(Expr, out)
end

function unfold_expr(mode::Symbol, fun::Symbol, ex, s...; force=true)
    typeof(ex) in [Nothing,LineNumberNode] ? ex : unfold_expr_force(mode,fun,ex,s...)
end

function unfold(mode::Symbol,fun::Symbol,expr::Expr,s...)
    #expr = mode == :unary ? squash(ixpr) : ixpr
    if expr.head == :block
        out = Any[]
        for line ∈ expr.args # block structure
            push!(out,unfold_expr(mode,fun,line,s...))
        end
        return Expr(:block,out...)
    else
        return unfold_expr(mode,fun,expr,s...)
    end
end

function mat(expr)
    if typeof(expr) == Expr && expr.head == :vcat
        if typeof(expr.args[1]) == Expr && expr.args[1].head == :row
            rows = length(expr.args)
            cols = length(expr.args[1].args)
            out = Array{Any,2}(undef,rows,cols)
            for k ∈ 1:rows
                for l ∈ 1:cols
                    out[k,l] = expr.args[k].args[l]
                end
            end
            return out
        else
            rows = length(expr.args)
            out = Array{Any,1}(undef,rows)
            for k ∈ 1:rows
                out[k] = expr.args[k]
            end
            return out
        end
    else
        return expr
    end
end

mat(expr,original...) = |(typeof.(original)...) <: Matrix ? mat(expr) : expr


text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

# Compiler: Julia 1.5

# Short description of this file: Implementation of Euler system with 6 target
#functions, considering η as a constant, and dealing with functions at
#denominator using an ϵ>0.

using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
using Plots, PyPlot
using DifferentialEquations, DiffEqBase, DiffEqOperators, LinearAlgebra, OrdinaryDiffEq
using Quadrature, Cubature, Cuba
using Random

Random.seed!(50)

@parameters x y η t θ
@variables u1(..), u2(..), u3(..), u4(..), u5(..), u6(..)
@derivatives Dt'~t
@derivatives Dx'~x
@derivatives Dy'~y
@derivatives Dη'~η
@derivatives Dxx''~x
@derivatives Dyy''~y

epsilon = 0.000001

#functions
u  = u1(x,y,η,t,θ)
v  = u2(x,y,η,t,θ)
w  = u3(x,y,η,t,θ)
θm = u4(x,y,η,t,θ)
qm = u5(x,y,η,t,θ)
z  = u6(x,y,η,t,θ)

#physical constants
g  = 9.81 #constant
γ  = 1.4 #air cp/cv=1.4
ρd = 1.225 #dry air density
p0 = 101325 #reference pressure 10^5 Pa
Rd = 8.31 #gas constant for dry air
ηc = 0.2 #page 8 of Advanced Research WRF V4 (http://dx.doi.org/10.5065/1dfh-6p97))

#simulation inputs
sum_of_qs = 0.1 #read from file
#dx =
#dy =
#dz =
Δx = 1 #dummy
Δy = 1 #dummy

#scale factors
dist = x^2 + y^2 #dummy for distance on earth  (MANFRIN's advice)
Dxd  = 2*x #Dx(dist)
Dyd  = 2*y #Dy(dist)
mx   = Δx/max(abs(dist),epsilon)
my   = Δy/max(abs(dist),epsilon)

#definitions
αd = 1/ρd
p  = p0*(Rd*θm/(p0*αd))^γ
c1 = 2*ηc*ηc/((1 - ηc)*(1 - ηc)*(1 - ηc))
c2 = -ηc*(4 + ηc + ηc*ηc)/((1 - ηc)*(1 - ηc)*(1 - ηc))
c3 = 2*(1 + ηc + ηc*ηc)/((1 - ηc)*(1 - ηc)*(1 - ηc))
c4 = -(1 + ηc)/((1 - ηc)*(1 - ηc)*(1 - ηc))
B  = c1 + c2*η + c3*η*η + c4*η*η*η
ps = 10 #dummy
pt = 1 #dummy
pd = B*(ps - pt) + (η - B)*(p0 - pt) + pt
μd = (ps - pt)*(c2 + 2*c3*η + 3*c4*η*η) + (p0 - pt)*(1 - c2 - 2*c3*η - 3*c4*η*η) #[μd = Dη(pd)]
D2pd = (ps - pt)*(2*c3 + 6*c4*η) + (p0 - pt)*(-2*c3 - 6*c4*η)
U  = μd*u/max(abs(my),epsilon)
V  = μd*v/max(abs(mx),epsilon)
W  = μd*w/max(abs(my),epsilon)
Θm = μd*θm
Qm = μd*qm
ϕ  = g*z
α  = αd/(1 + sum_of_qs)

#divergences
#divVu  = Dx(U*u)  + Dy(V*u)  + Dη(Ω*u)
#divVv  = Dx(U*v)  + Dy(V*v)  + Dη(Ω*v)
#divVw  = Dx(U*w)  + Dy(V*w)  + Dη(Ω*w)
#divVθm = Dx(U*θm) + Dy(V*θm) + Dη(Ω*θm)
#divV   = Dx(V)    + Dy(V)    + Dη(V)
#divVqm = Dx(U*qm) + Dy(V*qm) + Dη(Ω*qm)

divVu1 = μd*u/Δy*(u*Dxd + 2*dist*Dx(u))
divVu2 = μd/Δx*(u*v*Dyd + dist*u*Dy(v) + dist*v*Dy(u))
divVu3 = 0 #dist/Δy*(D2pd*u*ω + μd*Dη(u)*ω + μd*u*Dη(ω))
divVu  = divVu1 + divVu2 + divVu3

divVv1 = μd/Δy*(u*v*Dxd + dist*u*Dx(v) + dist*v*Dx(u))
divVv2 = μd/Δx*(v*v*Dyd + dist*2*v*Dy(v))
divVv3 = 1/max(abs(my),epsilon)*D2pd*w*v #1/my*(D2pd*w*v + μd*Dη(ω)*v + μd*ω*Dη(v))
divVv  = divVv1 + divVv2 + divVv3

divVw1 = μd/Δy*(Dxd*u*w + dist*Dx(u)*w + dist*u*Dx(w))
divVw2 = μd/Δx*(Dyd*v*w + dist*Dy(v)*w + dist*v*Dy(w))
divVw3 = 0 #1/my*(D2pd*ω*w + μd*Dη(ω)*w   + μd*ω*Dη(w))
divVw  = divVw1 + divVw2 + divVw3

divVθm1 = μd/Δy*(Dxd*u*θm + dist*Dx(u)*θm + dist*u*Dx(θm))
divVθm2 = μd/Δx*(Dyd*v*θm + dist*Dy(v)*θm + dist*v*Dy(θm))
divVθm3 = 0 #1/my*(D2pd*ω*θm + μd*Dη(ω)*θm   + μd*ω*Dη(θm))
divVθm  = divVθm1 + divVθm2 + divVθm3

divV1 = μd/Δx*(Dxd*v + dist*Dx(v))
divV2 = μd/Δx*(Dyd*v + dist*Dy(v))
divV3 = 1/max(abs(mx),epsilon)*(D2pd*v + μd*Dη(v))
divV  = divV1 + divV2 + divV3

divVqm1 = μd/Δy*(Dxd*u*qm + dist*Dx(u)*qm + dist*u*Dx(qm))
divVqm2 = μd/Δx*(Dyd*v*qm + dist*Dy(v)*qm + dist*v*Dy(qm))
divVqm3 = 0 #1/my*(D2pd*ω*qm + μd*Dη(ω)*qm   + μd*ω*Dη(qm))
divVqm  = divVqm1 +divVqm2 + divVqm3

#inner products
innerV∇ϕ = U*g*Dx(z) + V*g*Dy(z)

#geographic parametrization
Ωe = 0.0000727 #dummy #angular rotation rate of earth
ψ  = 1 #dummy #latitude on earth
re = 6371 #dummy #radius of earth
αr = 1 #dummy #local rotation angle between y-axis and the meridians

f  = 2*Ωe*sin(ψ)
e  = 2*Ωe*cos(ψ)

n = 0.0001
#external forcings
Dymx = -Δx/max(dist*dist,epsilon)*Dyd
Dxmy = -Δy/max(dist*dist,epsilon)*Dxd

FU  =  ((mx/max(abs(my),epsilon))*(f + u*(my/max(abs(mx),epsilon))*Dymx - v*Dxmy)*V - (u/re + e*cos(αr))*W)*n
FV  =  (-(my/max(abs(mx),epsilon))*(f + u*(my/max(abs(mx),epsilon))*Dymx - v*Dxmy)*U + (v/re - e*sin(αr))*W)*n
FW  =  (e*(U*cos(αr) - (mx/max(abs(my),epsilon))*V*sin(αr) + 1/re*(u*U + (mx/max(abs(my),epsilon))*v*V)))*n
FΘm = 0 #dummy
FQm = 0 #dummy

#equations (Flux-Form Euler)
eq1 = (μd/max(abs(my),epsilon)*Dt(u) + divVu + μd*α*p0*(Rd/(p0*αd))^γ*γ*θm^(γ-1)*Dx(θm) +
      (α/αd)*p0*(Rd/(p0*αd))^γ*γ*θm^(γ-1)*Dη(θm)*g*Dx(z))*n ~ FU
eq2 = (μd/max(abs(mx),epsilon)*Dt(v) + divVv + μd*α*p0*(Rd/(p0*αd))^γ*γ*θm^(γ-1)*Dy(θm) +
      (α/αd)*p0*(Rd/(p0*αd))^γ*γ*θm^(γ-1)*Dη(θm)*g*Dy(z))*n ~ FV
eq3 = (μd/max(abs(my),epsilon)*Dt(w) + divVw - g*((α/αd)*p0*(Rd/(p0*αd))^γ*γ*θm^(γ-1)*Dη(θm) -
       μd))*n     ~ FW

eq4 = (μd*Dt(θm) + divVθm)*n                   ~ FΘm #coupling with SFIRE
eq5 = divV  *n                               ~ 0.
eq6 = (g*Dt(z)  + (1.0/μd)*(innerV∇ϕ - g*W))*n ~ 0.
eq7 = (μd*Dt(qm) + divVqm  )*n                 ~ FQm  #coupling with SFIRE

#system of equations
eqs = [eq1,eq2,eq3,eq4,eq5,eq6,eq7]

# Initial and boundary conditions
x0 = 0.5
y0 = 0.5
η0 = 0.5

bcs =    [u1(x,y,η,0,θ) ~ (x - x0)^2 + (y - y0)^2 + (η - η0)^2,
          u2(x,y,η,0,θ) ~ (x - x0)^2 + (y - y0)^2 + (η - η0)^2,
          u3(x,y,η,0,θ) ~ (x - x0)^2 + (y - y0)^2 + (η - η0)^2,
          u4(x,y,η,0,θ) ~ (x - x0)^2 + (y - y0)^2 + (η - η0)^2,
          u5(x,y,η,0,θ) ~ (x - x0)^2 + (y - y0)^2 + (η - η0)^2,
          u6(x,y,η,0,θ) ~ (x - x0)^2 + (y - y0)^2 + (η - η0)^2]

#=
k = 0.4
u0 = 0.5
z0 = 0.5

bcs = [u1(x,y,η,0,θ) ~ u0/k*log(z/z0),
       u2(x,y,η,0,θ) ~ 0.,
       u3(x,y,η,0,θ) ~ 0.,
       u4(x,y,η,0,θ) ~ 0.,
       u5(x,y,η,0,θ) ~ 320,
       u6(x,y,η,0,θ) ~ 1.,
       u7(x,y,η,0,θ) ~ 10.] # physically meaningful bcs
=#

# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0),
           y ∈ IntervalDomain(0.0,1.0),
           η ∈ IntervalDomain(0.0,1.0),
           t ∈ IntervalDomain(0.0,1.0)] #Discretization

# Discretization
dx = 0.1; dy= 0.1; dη = 0.1; dt = 0.1

# Neural network
input_ = length(domains)
output = 6

n     = 16
chain = FastChain(FastDense(input_,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,output))

q_strategy = NeuralPDE.QuadratureTraining(algorithm =CubaCuhre(),reltol=1e-8,abstol=1e-8,maxiters=150)
discretization = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain,strategy = q_strategy)

pde_system = PDESystem(eqs,bcs,domains,[x,y,η,t],[u1,u2,u3,u4,u5,u6])
prob = NeuralPDE.discretize(pde_system,discretization)

listTraining = [0.0]

cb = function (p,l)
    println("Current loss is: $l")
    append!(listTraining,l)
    return false
end

res = GalacticOptim.solve(prob, ADAM(0.001), cb = cb, maxiters=1500) ##multi-step training to improve the convergence

initθ = res.minimizer
discretization2 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ; strategy = q_strategy)
initθ == discretization2.initθ
prob2 = NeuralPDE.discretize(pde_system,discretization2)
res2 = GalacticOptim.solve(prob2, ADAM(0.001), progress = true, cb = cb, maxiters=1000)

initθ2 = res2.minimizer
discretization3 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ2; strategy = q_strategy)
initθ2 == discretization3.initθ
prob3 = NeuralPDE.discretize(pde_system,discretization3)
res3 = GalacticOptim.solve(prob3, ADAM(0.001), progress = true, cb = cb, maxiters=100)

initθ3 = res3.minimizer
discretization4 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ3; strategy = q_strategy)
initθ3 == discretization4.initθ
prob4 = NeuralPDE.discretize(pde_system,discretization4)
res3 = GalacticOptim.solve(prob4, ADAM(0.0005), progress = true, cb = cb, maxiters=100)

initθ4 = res3.minimizer
discretization5 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ4; strategy = q_strategy)
initθ4 == discretization5.initθ
prob5 = NeuralPDE.discretize(pde_system,discretization5)
res5 = GalacticOptim.solve(prob5, ADAM(0.0005), progress = true, cb = cb, maxiters=300)

initθ5 = res5.minimizer
discretization6 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ5; strategy = q_strategy)
initθ5 == discretization6.initθ
prob6 = NeuralPDE.discretize(pde_system,discretization6)
res6 = GalacticOptim.solve(prob6, ADAM(0.0005), progress = true, cb = cb, maxiters=300)

initθ6 = res6.minimizer
discretization7 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ6; strategy = q_strategy)
initθ6 == discretization7.initθ
prob7 = NeuralPDE.discretize(pde_system,discretization7)
res7 = GalacticOptim.solve(prob7, ADAM(0.00025), progress = true, cb = cb, maxiters=400)

initθ7 = res7.minimizer
discretization8 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ7; strategy = q_strategy)
initθ7 == discretization8.initθ
prob8 = NeuralPDE.discretize(pde_system,discretization8)
res8 = GalacticOptim.solve(prob8, ADAM(0.00025), progress = true, cb = cb, maxiters=400)

initθ8 = res8.minimizer
discretization9 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ8; strategy = q_strategy)
initθ8 == discretization9.initθ
prob9 = NeuralPDE.discretize(pde_system,discretization9)
res9 = GalacticOptim.solve(prob9, ADAM(0.00025), progress = true, cb = cb, maxiters=500)

initθ9 = res9.minimizer
discretization10 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ9; strategy = q_strategy)
initθ9 == discretization10.initθ
prob10 = NeuralPDE.discretize(pde_system,discretization10)
res10 = GalacticOptim.solve(prob10, ADAM(0.00025), progress = true, cb = cb, maxiters=600)

initθ10 = res10.minimizer
discretization11 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ10; strategy = q_strategy)
initθ10 == discretization11.initθ
prob11 = NeuralPDE.discretize(pde_system,discretization11)
res11 = GalacticOptim.solve(prob11, ADAM(0.0005), progress = true, cb = cb, maxiters=300)

initθ11 = res11.minimizer
discretization12 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ11; strategy = q_strategy)
initθ11 == discretization12.initθ
prob12 = NeuralPDE.discretize(pde_system,discretization12)
res12 = GalacticOptim.solve(prob12, ADAM(0.0005), progress = true, cb = cb, maxiters=400)

initθ12 = res12.minimizer
discretization13 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ12; strategy = q_strategy)
initθ12 == discretization13.initθ
prob13 = NeuralPDE.discretize(pde_system,discretization13)
res13 = GalacticOptim.solve(prob13, ADAM(0.0005), progress = true, cb = cb, maxiters=250)

initθ13 = res13.minimizer
discretization14 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ13; strategy = q_strategy)
initθ13 == discretization14.initθ
prob14 = NeuralPDE.discretize(pde_system,discretization14)
res14 = GalacticOptim.solve(prob14, ADAM(0.00025), progress = true, cb = cb, maxiters=40)

initθ14 = res14.minimizer
discretization15 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ14; strategy = q_strategy)
initθ14 == discretization15.initθ
prob15 = NeuralPDE.discretize(pde_system,discretization15)
res15 = GalacticOptim.solve(prob15, ADAM(0.000001), progress = true, cb = cb, maxiters=20)

initθ15 = res15.minimizer
discretization16 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ15; strategy = q_strategy)
initθ15 == discretization16.initθ
prob16 = NeuralPDE.discretize(pde_system,discretization16)
res16 = GalacticOptim.solve(prob16, ADAM(0.000001), progress = true, cb = cb, maxiters=30)

initθ16 = res16.minimizer
discretization17 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ16; strategy = q_strategy)
initθ16 == discretization17.initθ
prob17 = NeuralPDE.discretize(pde_system,discretization17)
res17 = GalacticOptim.solve(prob17, ADAM(0.000001), progress = true, cb = cb, maxiters=20)

initθ17 = res17.minimizer
discretization18 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ17; strategy = q_strategy)
initθ17 == discretization18.initθ
prob18 = NeuralPDE.discretize(pde_system,discretization18)
res18 = GalacticOptim.solve(prob18, ADAM(0.000001), progress = true, cb = cb, maxiters=14)

q_strategy2 = NeuralPDE.QuadratureTraining(algorithm =CubaCuhre(),reltol=1e-15,abstol=1e-15,maxiters=500)

initθ18 = res18.minimizer
discretization19 = NeuralPDE.PhysicsInformedNN([dx,dy,dη,dt],chain, initθ18; strategy = q_strategy2)
initθ18 == discretization19.initθ
prob19 = NeuralPDE.discretize(pde_system,discretization19)
res19 = GalacticOptim.solve(prob19, ADAM(0.001), progress = true, cb = cb, maxiters=200)

phi = discretization18.phi

##SAVE PARAMETERS TO REPRODUCE THE RESULT WITHOUT TRAINING AGAIN
outfile = "params_euler_system_res18_2.txt"
open(outfile, "w") do f
  for i in initθ18
    println(f, i)
  end
end

##OUTPUTS
xs,ys,ηs,ts = [domain.domain.lower:dx:domain.domain.upper for (dx,domain) in zip([dx,dy,dη,dt],domains)]

u_predict  = [reshape([phi([x,y,η,t],res18.minimizer)[i] for x in xs for y in ys for η in ηs for t in ts],
              (length(xs),length(ys),length(ηs),length(ts))) for i in 1:6]

##TRAINING LOSS PLOT
Plots.plot(listTraining[3:end], yaxis=:log, ylabel="Loss function", xlabel="Epochs (ADAM)", legend=false, title="Training")
Plots.savefig("euler_training_test_18_steps.pdf")

fn = 6

##BCS PLOTS
ic(x,y) = (x - x0)^2 + (y - y0)^2
initial_cond = reshape([ic(x,y) for x in xs for y in ys], (length(xs),length(ys)))
Plots.plot(xs, ys, initial_cond, st=:surface)
Plots.plot(xs, ys, u_predict[fn][:,:,1,1], st=:surface)

##OUTPUT PLOTS
maxlim = maximum(maximum(u_predict[fn][:,:,1,t]) for t = 1:length(ts))
minlim = minimum(minimum(u_predict[fn][:,:,1,t]) for t = 1:length(ts))

result = @animate for time = 1:length(ts)
    Plots.plot(xs, ys, u_predict[fn][:,:,1,time],st=:surface,camera=(30,30), zlim=(minlim,maxlim), clim=(minlim,maxlim),
                title = string("ψ: max = ",round(maxlim, digits = 3)," min = ", round(minlim, digits = 3),"\\n"," t = ",time))
end

gif(result, "euler_6func_contour.gif", fps = 3)

ext: jl


text:
using ValueHistories, Primes, LaTeXStrings, Plots
pyplot()

_mkpath(fname) = joinpath(dirname(@__FILE__()), fname)
srand(123)

# ---------------------------------------------------------------
# README Plots

fnt = Plots.Font("ubuntu",8,:hcenter,:vcenter,0.0,
                 RGB{FixedPointNumbers.Normed{UInt8, 8}}(0.0,0.0,0.0))
default(bg_outside = RGBA(1,1,1,0),
        grid = false,
        legendfont = fnt,
        palette = :darktest,
        size = (600, 400))

# ---------------------------------------------------------------

history = QHistory(Float64)
for i = 1:100
    isprime(i) && push!(history, i, sin(.1*i))
end
p = plot(history, legend=false)

Plots.svg(_mkpath("qhistory"))

# ---------------------------------------------------------------

history = MVHistory()
for i=1:100
    x = 0.1i
    push!(history, :mysin, x, sin(x))
    push!(history, :mystring, i, "i=$i")
    isprime(i) && push!(history, :mycos, Float32(x), cos(x))
end
p = plot(history)

Plots.svg(_mkpath("mvhistory"))

max_stars_repo_name: JuliaTagBot/FileStorage

max_issues_repo_name: JuliaTagBot/FileStorage

max_forks_repo_name: JuliaTagBot/FileStorage


text:
__precompile__(true)
module DualFFT

# ---- Imported packages ---- #
using AbstractFFTs, ForwardDiff
using ForwardDiff: Dual, Partials
using AbstractFFTs: Plan, ScaledPlan
import LinearAlgebra
import FFTW

# ---- Imported base functions ---- #
import ForwardDiff: value, partials, npartials, valtype, tagtype
import AbstractFFTs: plan_fft, plan_inv, plan_bfft
import Base: *
import LinearAlgebra.mul!
import FFTW: set_timelimit, dims_howmany, unsafe_execute!, cFFTWPlan, r2rFFTWPlan, PlanPtr, FFTWPlan, ScaledPlan, destroy_plan
import AbstractFFTs: normalization, complexfloat, strides
import Random
# ---- Source files ---- #
include("complex_dual.jl")
include("abstract_dual_fft.jl")
include("dual_fft.jl")

# ---- Exported functions ---- #
export dual2array, array2dual

end


text:
"""
    decode(s)

Decode some [ETF](http://erlang.org/doc/apps/erts/erl_ext_dist.html) data.
"""
function decode(s)
    io = IOBuffer(s)
    v = first(read(io, 1))
    v == VERSION_MAGIC || error("Unknown ETF version $v")
    return decode_from_tag!(io)
end

function decode_from_tag!(io::IO)
    t = first(read(io, 1))
    haskey(DECODER_TAGS, t) || error("Unknown ETF tag $t")
    return DECODER_TAGS[t](io)
end

decode_multiple!(io::IO, n::Integer) = Any[decode_from_tag!(io) for _ in 1:n]

# Integers

decode_small_integer!(io::IO) = first(read(io, 1))
function decode_integer!(io::IO, T::Type{<:Integer}=Int32, sz::Int=4, be=true)
    n = first(reinterpret(T, read(io, sz)))
    return be ? ntoh(n) : n
end
function decode_big!(io::IO, T::Type{<:Integer}, sz::Int)
    len = decode_integer!(io, T, sz, false)
    neg = first(read(io, 1)) == 1

    # http://erlang.org/doc/apps/erts/erl_ext_dist.html#small_big_ext
    B = BigInt(256)
    n = sum(map(p -> p[2] * B^(p[1]-1), enumerate(read(io, len))))

    return neg ? -n : n
end
decode_small_big!(io::IO) = decode_big!(io, UInt8, 1)
decode_large_big!(io::IO) = decode_big!(io, UInt32, 4)

# Floats

decode_float!(io::IO) = parse(Float64, String(read(io, 31)))
decode_new_float!(io::IO) = ntoh(reinterpret(Float64, read(io, 8)))

# Tuples

decode_tuple!(io::IO, n::Integer) = Tuple(decode_multiple!(io, n))
function decode_small_tuple!(io::IO)
    len = first(read(io, 1))
    return decode_tuple!(io, len)
end
function decode_large_tuple!(io::IO)
    len = decode_integer!(io, UInt32, 4, true)
    decode_tuple!(io, len)
end

# Atoms

function decode_atom!(io::IO, T::Type{<:Integer}=UInt16, sz::Int=2)
    len = decode_integer!(io, T, sz, true)
    s = Symbol(read(io, len))

    return if s === :nil
        nothing
    elseif s === :true
        true
    elseif s === :false
        false
    else
        s
    end
end
decode_small_atom!(io::IO) = decode_atom!(io, UInt8, 1)
decode_atom_utf8!(io::IO) = decode_atom!(io, UInt16, 2)
decode_small_atom_utf8!(io::IO) = decode_atom!(io, UInt8, 1)

# Lists

decode_nil!(io::IO) = []
function decode_list!(io::IO)
    len = decode_integer!(io, UInt32, 4, true)
    xs = decode_multiple!(io, len)
    decode_from_tag!(io)
    return xs
end
function decode_string!(io::IO)
    len = decode_integer!(io, UInt16, 2, true)
    return read(io, len)
end

# Maps

function decode_map!(io::IO)
    len = decode_integer!(io, UInt32, 4, true)
    xs = decode_multiple!(io, 2 * len)
    return Dict{Any, Any}(xs[i] => xs[i+1] for i in 1:2:length(xs))
end

# Strings

function decode_binary!(io::IO)
    len = decode_integer!(io, UInt32, 4, true)
    return String(read(io, len))
end

const DECODER_TAGS = Dict(
    70 => decode_new_float!,
    # 77  => decode_bit_binary!,
    # 82  => decode_atom_cache!,
    97  => decode_small_integer!,
    98  => decode_integer!,
    99  => decode_float!,
    100 => decode_atom!,
    # 101 => decode_reference!,
    # 102 => decode_port!,
    # 103 => decode_pid!,
    104 => decode_small_tuple!,
    105 => decode_large_tuple!,
    106 => decode_nil!,
    107 => decode_string!,
    108 => decode_list!,
    109 => decode_binary!,
    110 => decode_small_big!,
    111 => decode_large_big!,
    # 112 => decode_new_fun!,
    # 113 => decode_export!,
    # 114 => decode_new_reference!,
    115 => decode_small_atom!,
    116 => decode_map!,
    # 117 => decode_fun!,
    118 => decode_atom_utf8!,
    119 => decode_small_atom_utf8!,
)


text:
using Pkg

Pkg.instantiate()

max_stars_repo_name: MathOnco/ploidyEvolution

max_issues_repo_name: MathOnco/ploidyEvolution

max_forks_repo_name: MathOnco/ploidyEvolution


text:
include("../test_fem_sa.jl")

nparts = (2,2)
prun(test_fem_sa,mpi,nparts)


text:
export TwoStageCost, two_stage_method

struct TwoStageCost{F,F2,D} <: Function
  cost_function::F
  cost_function2::F2
  estimated_solution::D
  estimated_derivative::D
end

(f::TwoStageCost)(p) = f.cost_function(p)
(f::TwoStageCost)(p,g) = f.cost_function2(p,g)

decide_kernel(kernel::CollocationKernel) = kernel
function decide_kernel(kernel::Symbol)
    if kernel == :Epanechnikov
        return EpanechnikovKernel()
    elseif kernel == :Uniform
        return UniformKernel()
    elseif kernel == :Triangular
        return TriangularKernel()
    elseif kernel == :Quartic
      return QuarticKernel()
    elseif kernel == :Triweight
      return TriweightKernel()
    elseif kernel == :Tricube
      return TricubeKernel()
    elseif kernel == :Gaussian
      return GaussianKernel()
    elseif kernel == :Cosine
      return CosineKernel()
    elseif kernel == :Logistic
      return LogisticKernel()
    elseif kernel == :Sigmoid
      return SigmoidKernel()
    else
      return SilvermanKernel()
    end
end

function construct_t1(t,tpoints)
    hcat(ones(eltype(tpoints),length(tpoints)),tpoints.-t)
end
function construct_t2(t,tpoints)
  hcat(ones(eltype(tpoints),length(tpoints)),tpoints.-t,(tpoints.-t).^2)
end
function construct_w(t,tpoints,h,kernel)
    W = @. calckernel((kernel,),(tpoints-t)/h)/h
    Diagonal(W)
end
function construct_estimated_solution_and_derivative!(data,kernel,tpoints)
  _one = oneunit(first(data))
  _zero = zero(first(data))
  e1 = [_one;_zero]
  e2 = [_zero;_one;_zero]
  n = length(tpoints)
  h = (n^(-1/5))*(n^(-3/35))*((log(n))^(-1/16))

  Wd = similar(data, n, size(data,1))
  WT1 = similar(data, n, 2)
  WT2 = similar(data, n, 3)
  x = map(tpoints) do _t
      T1 = construct_t1(_t,tpoints)
      T2 = construct_t2(_t,tpoints)
      W = construct_w(_t,tpoints,h,kernel)
      mul!(Wd,W,data')
      mul!(WT1,W,T1)
      mul!(WT2,W,T2)
      (e2'*((T2'*WT2)\T2'))*Wd,(e1'*((T1'*WT1)\T1'))*Wd
  end
  estimated_derivative = reduce(hcat,transpose.(first.(x)))
  estimated_solution = reduce(hcat,transpose.(last.(x)))
  estimated_derivative,estimated_solution
end
function construct_iip_cost_function(f,du,preview_est_sol,preview_est_deriv,tpoints)
  function (p)
      _du = PreallocationTools.get_tmp(du,p)
      vecdu = vec(_du)
      cost = zero(first(p))
      for i in 1:length(preview_est_sol)
        est_sol = preview_est_sol[i]
        f(_du,est_sol,p,tpoints[i])
        vecdu .= vec(preview_est_deriv[i]) .- vec(_du)
        cost += sum(abs2,vecdu)
      end
      sqrt(cost)
  end
end

function construct_oop_cost_function(f,du,preview_est_sol,preview_est_deriv,tpoints)
  function (p)
      cost = zero(first(p))
      for i in 1:length(preview_est_sol)
        est_sol = preview_est_sol[i]
        _du = f(est_sol,p,tpoints[i])
        cost += sum(abs2,vec(preview_est_deriv[i]) .- vec(_du))
      end
      sqrt(cost)
  end
end

get_chunksize(cs) = cs
get_chunksize(cs::Type{Val{CS}}) where CS = CS

function two_stage_method(prob::DiffEqBase.DEProblem,tpoints,data;kernel= EpanechnikovKernel(),
                          loss_func = L2Loss,mpg_autodiff = false,
                          verbose = false,verbose_steps = 100,
                          autodiff_chunk = length(prob.p))
    f = prob.f
    kernel_function = decide_kernel(kernel)
    estimated_derivative,estimated_solution = construct_estimated_solution_and_derivative!(data,kernel_function,tpoints)

    # Step - 2

    du = PreallocationTools.dualcache(similar(prob.u0), autodiff_chunk)
    preview_est_sol = [@view estimated_solution[:,i] for i in 1:size(estimated_solution,2)]
    preview_est_deriv = [@view estimated_derivative[:,i] for i in 1:size(estimated_solution,2)]
    if DiffEqBase.isinplace(prob)
      cost_function = construct_iip_cost_function(f,du,preview_est_sol,preview_est_deriv,tpoints)
    else
      cost_function = construct_oop_cost_function(f,du,preview_est_sol,preview_est_deriv,tpoints)
    end

    if mpg_autodiff
      gcfg = ForwardDiff.GradientConfig(cost_function, prob.p, ForwardDiff.Chunk{get_chunksize(autodiff_chunk)}())
      g! = (x, out) -> ForwardDiff.gradient!(out, cost_function, x, gcfg)
    else
      g! = (x, out) -> Calculus.finite_difference!(cost_function,x,out,:central)
    end
    if verbose
      count = 0 # keep track of # function evaluations
    end
    cost_function2 = function (p,grad)
      if length(grad)>0
        g!(p,grad)
      end
      loss_val = cost_function(p)
      if verbose
        count::Int += 1
        if mod(count,verbose_steps) == 0
          println("Iteration: $count")
          println("Current Cost: $loss_val")
          println("Parameters: $p")
        end
      end
      loss_val
    end

    return TwoStageCost(cost_function, cost_function2, estimated_solution, estimated_derivative)
end


text:
export zoom_generic

function zoom_generic(h :: AbstractNLPModel,
					  nlpstop :: AbstractStopping,
                  	  t₀ = 2.7, t₁ = 7.5;
					  direction = :Nwt, tol = 1e-7, verbose = false,
					  kwargs...)

    if true in (obj(h, t₀) .< obj(h, t₁))
    	tl = t₀; th = t₁
    else
    	tl = t₁; th = t₀
    end

    hl = obj(h, tl); dhl = grad(h, tl)
    hh = obj(h, th); dhh = grad(h, th)
    h₀ = obj(h, t₀); dh₀ = grad(h, t₀)

    γ = 0.8; t = t₁; tp = t₀; tqnp = t₀; hk = 0; hkm1 = 0; gkm1 = 0; hp = 0
    i = 0

    hk   = obj(h, t);    gk   = grad(h, t)
    hkm1 = obj(h, tqnp); gkm1 = grad(h, tqnp)

	OK = stop!(nlpstop)

    verbose &&
        @printf(" iter        tₗ        tₕ         t        hₗ        hₕ")
    verbose &&
        @printf("         hk         gk\n")
    verbose &&
        @printf(" %7.2e %7.2e  %7.2e  %7.2e  %7.2e %7.2e %7.2e %7.2e\n",
                i, tl[1], th[1], t[1], hl[1], hh[1], hk[1], gk[1])

	while !OK

		if true in (hk .> h₀) || true in (hl .<= hk)
			th = t; tlast = th; hlast = hh
			dhlast = dhh; hh = hk; dhh = gk
		else
			if true in (abs.(gk) .< tol)
				topt = t; iter = i
				return (topt, iter)
			elseif  true in (gk .* (th - tl) .>= 0)
				th = tl; hh = hl; dhh = dhl
			end
		tl = t; tlast = tl; hlast = hl
		dhlast = dhl; hl = hk; dhl = gk
		end

		if direction == :Nwt
			kₖ = hess(h, t)
			dN = vec(- gk ./ kₖ) #direction de Newton
		elseif direction == :Sec
			s = t .- tqnp; y = gk .- gkm1
			dN =-gk .* s ./ y
		elseif direction == :SecA
			s = t .- tqnp; y = gk .- gkm1
			Γ = 3.0 .* (gk .+ gkm1) .* s .- 6.0 .* (hk .- hkm1)
			if y .* s .+ Γ < eps(Float64) .* (s .^ 2)
				yt = y
			else
				yt = y .+ Γ ./ s
			end
			dN = -gk .* s ./ yt
		elseif direction == :Cub
			s = t .- tqnp; y = gk .- gkm1; α = -s
			z = gk .+ gkm1 .+ 3.0 .* (hk .- hkm1) ./ α
			discr = z .^ 2 .- gk .* gkm1; denom = gk .+ gkm1 .+ 2.0 .* z
			if true in (discr .> 0.0) && true in (abs.(denom) .> eps(Float64))
				w = sqrt.(discr)
				dN = -s .* (gk .+ z .+ sign.(α) .* w) ./ (denom)
			else
				dN=-gk .* s ./ y
			end
		end

		if true in ((tp .- t) .* dN .> 0.0) && true in  (dN ./ (tp .- t) .< γ)
      		tplus = t .+ dN
      		hplus = obj(h, tplus); gplus = grad(h, tplus)
      		verbose && println("N")
    	else
      		tplus = (t .+ tp) ./ 2.0
      		hplus = obj(h, tplus); gplus = grad(h, tplus)
      		verbose && println("B")
    	end

		if true in (t .> tp)
      		if true in (gplus .< 0.0)
        		tp = t; tqnp = t; t = tplus
      		else
        		tqnp = t; t = tplus
      		end
    	else
      		if true in (gplus .> 0.0)
        		tp = t; tqnp = t; t = tplus
      		else
        		tqnp = t; t = tplus
      		end
    	end

    	#mise à jour des valeurs
    	hkm1 = hk; gkm1 = gk; hk = hplus; gk = gplus
		OK = update_and_stop!(nlpstop, x = t, fx = hk, gx = gk)
    	i = i + 1
    	verbose &&
        	@printf(" %7.2e %7.2e  %7.2e  %7.2e %7.2e %7.2e %7.2e %7.2e\n",
                	i, tl[1], th[1], t[1], hl[1], hh[1], hk[1], gk[1])

    end

    return nlpstop.meta.optimal, nlpstop
end #function


text:
###############################################################################
#
# Radiative conductance
#
###############################################################################
"""
    radiative_conductance(T::FT) where {FT<:AbstractFloat}

Return the radiative conductance of leaf, given
- `T` leaf temperature
"""
function radiative_conductance(
            T::FT
) where {FT<:AbstractFloat}
    _A0::FT = 0.1579;
    _A1::FT = 0.0017;
    _A2::FT = 7.17E-6;
    _T ::FT = T - T_0(FT);

    return _A0 + _A1*T + _A2*T^2
end

###############################################################################
#
# Black body energy emittance
#
###############################################################################
"""
    black_body_emittance(T::FT) where {FT<:AbstractFloat}

Return the energy been radiated out, given
- `T` leaf temperature
"""
function black_body_emittance(
            T::FT
) where {FT<:AbstractFloat}
    return K_STEFAN(FT) * T^4
end

###############################################################################
#
# Leaf boundary layer conductance
#
###############################################################################
"""
    boundary_layer_conductance(wind::FT, width::FT) where {FT<:AbstractFloat}

Return the boundary layer conductance, given
- `wind` Wind speed
- `width` Leaf width
"""
function boundary_layer_conductance(
            wind::FT,
            width::FT
) where {FT<:AbstractFloat}
    _A::FT = FT(1.4 * 0.135)
    _B::FT = FT(0.72)

    return _A * sqrt( wind / (_B*width) )
end

###############################################################################
#
# Leaf temperature
#
###############################################################################
"""
    leaf_temperature(
                node::SPACSimple{FT},
                rad::FT,
                e_rad::FT,
                epla::FT
    ) where {FT<:AbstractFloat}
    leaf_temperature(
                node::SPACSimple{FT},
                rad::FT,
                flow::FT
    ) where {FT<:AbstractFloat}

Return leaf temperature, given
- `node` [`SPACSimple`](@ref) type struct
- `rad` Absorbed solar radiation per leaf area
- `e_rad` Emitted absorbed radiation per leaf area
- `epla` Flow rate per leaf area
- `flow` Total flow rate per basal area
"""
function leaf_temperature(
            node::SPACSimple{FT},
            rad::FT,
            e_rad::FT,
            epla::FT
) where {FT<:AbstractFloat}
    @unpack t_air,wind = node.envir;

    lambda = latent_heat_vapor(t_air) * M_H₂O(FT);
    Gr     = radiative_conductance(t_air);
    GHa    = boundary_layer_conductance(wind, node.width);
    e_lat  = lambda * epla;

    t_leaf = t_air + (rad-e_rad-e_lat) / (FT(29.3)*(Gr+GHa)) / 2;

    return t_leaf
end

function leaf_temperature(
            node::SPACSimple{FT},
            rad::FT,
            flow::FT
) where {FT<:AbstractFloat}
    lai = node.laba / node.gaba;

    e_rad = FT(0.97) * black_body_emittance(node.envir.t_air) / max(1, lai);
    epla  = flow / node.laba;

    return leaf_temperature(node, rad, e_rad, epla)
end

"""
    leaf_temperature_sunlit(
                node::SPACSimple{FT},
                rad::FT,
                f_sl::FT
    ) where {FT<:AbstractFloat}

Return leaf temperature, given
- `node` [`SPACSimple`](@ref) type struct
- `rad` Absorbed solar radiation per leaf area
- `f_sl` Total flow rate per basal area into sunlit leaves
"""
function leaf_temperature_sunlit(
            node::SPACSimple{FT},
            rad::FT,
            f_sl::FT
) where {FT<:AbstractFloat}
    e_rad  = FT(0.97) * black_body_emittance(node.envir.t_air) /
             max(1, (node.container2L).lai_sl);
    epla   = f_sl / (node.container2L).la_sl;

    return leaf_temperature(node, rad, e_rad, epla)
end

"""
    leaf_temperature_shaded(
                node::SPACSimple{FT},
                rad::FT,
                f_sh::FT
    ) where {FT<:AbstractFloat}

Return leaf temperature, given
- `node` [`SPACSimple`](@ref) type struct
- `rad` Absorbed solar radiation per leaf area
- `f_sh` Total flow rate per basal area into shaded leaves
"""
function leaf_temperature_shaded(
            node::SPACSimple{FT},
            rad::FT,
            f_sh::FT
) where {FT<:AbstractFloat}
    epla = f_sh / (node.container2L).la_sh;

    return leaf_temperature(node, rad, FT(0), epla)
end


text:
using Test
using DataFlowTasks
using DataFlowTasks: R,W,RW
using LinearAlgebra
using DataFlowTasks.TiledFactorization

include(joinpath(DataFlowTasks.PROJECT_ROOT,"test","testutils.jl"))

sch = DataFlowTasks.JuliaScheduler(200)
DataFlowTasks.setscheduler!(sch)

@testset "Julia scheduler" begin
    @testset "Fork-join" begin
        m = 10
        s = 0.1
        nw = Threads.nthreads()
        t = fetch(fork_join(m,s))
        t1 = @elapsed fetch(fork_join(m,s))
        t2 = (2+ceil(m/nw))*s
        # test that ideal vs actual time are close
        @test abs(t1-t2) < 1e-2
    end

    @testset "Stop dag worker" begin
        m = 10
        s = 0.1
        nw = Threads.nthreads()
        # stopping the dag_worker will mean finished nodes are no longer cleaned up
        DataFlowTasks.stop_dag_worker()
        fork_join(m,s)
        @test DataFlowTasks.num_nodes(sch.dag) == m+2
        # resuming dag_worker will now cleanup nodes
        DataFlowTasks.start_dag_worker()
        DataFlowTasks.sync() # wait for dag to be empty
        @test DataFlowTasks.num_nodes(sch.dag) == 0
        t1 = @elapsed fetch(fork_join(m,s))
        t2 = (2+ceil(m/nw))*s
        # test that ideal vs actual time are close
        @test abs(t1-t2) < 1e-2
    end

    @testset "Tiled cholesky factorization" begin
        m  = 100
        bsize = div(m,5)
        # create an SPD matrix
        A = rand(m,m)
        A = (A + adjoint(A))/2
        A = A + m*I
        F = TiledFactorization.cholesky(A)
        @test F.L*F.U ≈ A
    end

    @testset "Tiled lu factorization" begin
        m  = 100
        bsize = div(m,5)
        A = rand(m,m)
        F = TiledFactorization.lu(A)
        @test F.L*F.U ≈ A
    end
end


text:
@inline eps_around_one(θ::T) where T = 100sqrt(eps(one(θ)))

"""
    nlsolve!(nlsolver::AbstractNLSolver, integrator)

Solve
```math
dt⋅f(innertmp + γ⋅z, p, t + c⋅dt) + outertmp = z
```
where `dt` is the step size and `γ` and `c` are constants, and return the solution `z`.
"""
function nlsolve!(nlsolver::AbstractNLSolver, integrator, cache=nothing, repeat_step=false)
  @label REDO
  if isnewton(nlsolver)
    cache === nothing && throw(ArgumentError("cache is not passed to `nlsolve!` when using NLNewton"))
    if nlsolver.method === DIRK
      γW = nlsolver.γ * integrator.dt
    else
      γW = nlsolver.γ * integrator.dt / nlsolver.α
    end
    update_W!(nlsolver, integrator, cache, γW, repeat_step)
  end

  @unpack maxiters, κ, fast_convergence_cutoff = nlsolver

  initialize!(nlsolver, integrator)
  nlsolver.status = Divergence
  η = get_new_W!(nlsolver) ? initial_η(nlsolver, integrator) : nlsolver.ηold

  local ndz
  for iter in 1:maxiters
    nlsolver.iter = iter

    # compute next step and calculate norm of residuals
    iter > 1 && (ndzprev = ndz)
    ndz = compute_step!(nlsolver, integrator)
    if !isfinite(ndz)
      nlsolver.status = Divergence
      nlsolver.nfails += 1
      break
    end

    # check divergence (not in initial step)
    if iter > 1
      θ = ndz / ndzprev

      # When one Newton iteration basically does nothing, it's likely that we
      # are at the percision limit of floating point number. Thus, we just call
      # it convergence/divergence according to `ndz` directly.
      if abs(θ - one(θ)) <= eps_around_one(θ)
        if ndz <= one(ndz)
          nlsolver.status = Convergence
          nlsolver.nfails = 0
          break
        else
          nlsolver.status = Divergence
          nlsolver.nfails += 1
          break
        end
      end

      # divergence
      if θ > 2
        nlsolver.status = Divergence
        nlsolver.nfails += 1
        break
      end
    end

    apply_step!(nlsolver, integrator)

    # check for convergence
    iter > 1 && (η = θ / (1 - θ))
    if (iter == 1 && ndz < 1e-5) || (iter > 1 && (η >= zero(η) && η * ndz < κ))
      nlsolver.status = Convergence
      nlsolver.nfails = 0
      break
    end
  end

  if isnewton(nlsolver) && nlsolver.status == Divergence && !isJcurrent(nlsolver, integrator)
    nlsolver.status = TryAgain
    nlsolver.nfails += 1
    @goto REDO
  end

  nlsolver.ηold = η
  postamble!(nlsolver, integrator)
end

## default implementations

initialize!(::AbstractNLSolver, integrator) = nothing

initial_η(nlsolver::NLSolver, integrator) =
  max(nlsolver.ηold, eps(eltype(integrator.opts.reltol)))^(0.8)

function apply_step!(nlsolver::NLSolver{algType,iip}, integrator) where {algType,iip}
  if iip
    @.. nlsolver.z = nlsolver.ztmp
  else
    nlsolver.z = nlsolver.ztmp
  end

  nothing
end

function postamble!(nlsolver::NLSolver, integrator)
  if DiffEqBase.has_destats(integrator)
    integrator.destats.nnonliniter += nlsolver.iter

    if nlsolvefail(nlsolver)
      integrator.destats.nnonlinconvfail += 1
    end
  end
  integrator.force_stepfail = nlsolvefail(nlsolver)
  setfirststage!(nlsolver, false)
  isnewton(nlsolver) && (nlsolver.cache.firstcall = false)

  nlsolver.z
end


text: @everywhere println("yo")

max_stars_repo_name: miguelraz/CommunityServiceUNAM

max_issues_repo_name: miguelraz/CommunityServiceUNAM

max_forks_repo_name: miguelraz/CommunityServiceUNAM


text:
struct RecombinationWeights
    μ::Int
    μeff::Float64
    weights::Vector{Float64}
    positive_weights::SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}
    negative_weights::SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}
end
default_weights(λ) = [log((λ + 1)/2) - log(i) for i in 1:λ]
RecombinationWeights(λ::Int) = RecombinationWeights(default_weights(λ))
function _negative_weights_limit_sum!(w, v)
    sw = -sum(w.negative_weights)
    if sw > v
        w.negative_weights .*= v/sw
    end
    w
end
function finalize_negative_weights!(w, cov, n)
    w.negative_weights .*= 1 + cov.c1/cov.cμ
    _negative_weights_limit_sum!(w, (1 - cov.c1 - cov.cμ) / cov.cμ / n)
    _negative_weights_limit_sum!(w, 1 + 2 * _mueff(w.negative_weights)/(w.μeff + 2))
end
_mueff(w) = sum(w)^2/sum(w.^2)
function RecombinationWeights(w::Vector)
    μ = sum(w .> 0)
    positive_weights = view(w, 1:μ)
    positive_weights ./= sum(positive_weights)
    negative_weights = view(w, μ + 1:length(w))
    negative_weights ./= -sum(negative_weights)
    RecombinationWeights(μ, _mueff(positive_weights), w,
                         positive_weights, negative_weights)
end

mutable struct Sigma
    σ::Float64
    c::Float64
    d::Float64
    e::Float64
    h::Bool
    g::Int
    χn::Float64
    p::Vector{Float64}
end
function Base.show(io::IO, ::MIME"text/plain", s::Sigma)
    println(io, "Sigma")
    println(io, "  σ: $(s.σ)")
    println(io, "  c: $(s.c)")
    println(io, "  d: $(s.d)")
    println(io, "  h: $(s.h)")
end
default_cσ(n, μeff) = (μeff + 2)/(n + μeff + 3)
default_dσ(c, n, μeff) = 1 + max(0, √((μeff - 1)/(n + 1)) - 1) + c
function Sigma(σ0, n, μeff;
               c = default_cσ(n, μeff),
               d = default_dσ(c, n, μeff))
    e = √(c * (2 - c) * μeff)
    Sigma(σ0, c, d, e, true, 0, sqrt(n)*(1 - 1/(4n) + 1/(21n^2)), zeros(n))
end
function update_p!(s::Sigma, m, cov)
    s.p .= (1 - s.c) * s.p +  s.e * whiten(cov.C, m)
end
function update!(s::Sigma, m, C)
    update_p!(s, m, C)
    update!(s)
end
function update!(s::Sigma)
    norm_p = norm(s.p)
    n = length(s.p)
    s.h = norm_p^2/(1 - (1 - s.c)^(2(s.g + 1)))/n - 1 < 1 + 4 / (n + 1)
    tmp = s.c/s.d * (norm_p / s.χn - 1)
    tmp = clamp(tmp, -1, 1)
    s.σ *= exp(tmp)
    s.g += 1
    s
end

mutable struct Covariance{T}
    c::Float64
    c1::Float64
    cμ::Float64
    e::Float64
    p::Vector{Float64}
    C::T
end
function Base.show(io::IO, ::MIME"text/plain", s::Covariance)
    println(io, "Covariance")
    println(io, "  c: $(s.c)")
    println(io, "  c1: $(s.c1)")
    println(io, "  cμ: $(s.cμ)")
end
default_cc(n, μeff) = (4 + μeff/n)/(n + 4 + 2μeff/n)
default_cc1(n, μeff) = 2 / ((n + 1.3)^2 + μeff)
default_ccμ(c1, n, μeff) = min(1 - c1, 2 * (.25 + μeff - 2 + 1/μeff) / ((n + 2)^2 + μeff))
struct MEigen
    mat::Matrix{Float64}
    e::Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}
    sqrtvalues::Vector{Float64}
end
function MEigen(m)
    e = eigen(Symmetric(m))
    if e.values[1] < 0
        e.values .-= 1.1 * e.values[1]
    end
    MEigen(m, e, sqrt.(e.values))
end
function Covariance(n, μeff;
                    c = default_cc(n, μeff),
                    c1 = default_cc1(n, μeff),
                    cmu = default_ccμ(c1, n, μeff))
    Covariance(c, c1, cmu, √(c * (2 - c) * μeff),
               zeros(n),
               MEigen(diagm([exp(1e-4/n * i) for i in 0:n-1]))
              )
end
function update_p!(c::Covariance, m, h)
    c.p .*= (1 - c.c)
    if h
        c.p .+= c.e * m
    end
    c.p
end
maha_norm(c, x) = √sum(abs2, c.C.e.vectors' * x ./ c.C.sqrtvalues)
function update!(c::Covariance, m, y, perm, w, h)
    update_p!(c, m, h)
    weights = c.cμ * w.weights
    c.C.mat .*= 1 - h * c.c1 - sum(weights)
    BLAS.syr!('U', c.c1, c.p, c.C.mat)
    for (i, j) in enumerate(perm)
        w = weights[i]
        w == 0 && continue
        v = y[:, j]
        if w < 0
            w *= length(v) / (maha_norm(c, v) + 1e-9)^2
        end
        BLAS.syr!('U', w, v, c.C.mat)
    end
    c.C = MEigen(c.C.mat)
    c
end
whiten(e::MEigen, x) = e.e.vectors * ((e.e.vectors' * x) ./ e.sqrtvalues)
unwhiten(e::MEigen, x) = e.e.vectors * (e.sqrtvalues .* x)

mutable struct Parameters{V,T,N,R}
    n::Int
    λ::Int
    mean::V
    sigma::Sigma
    cov::Covariance
    weights::RecombinationWeights
    constraints::T
    noise_handling::N
    parallel_evaluation::Bool
    multi_threading::Bool
    seed::UInt
    rng::R
end
function Base.show(io::IO, ::MIME"text/plain", s::Parameters{T, N}) where {T, N}
    println(io, "CMAParameters{$(T.name.name),$(N.name.name)}")
    println(io, "  dimensions: $(s.n)")
    println(io, "  population size: $(s.λ)")
end
sigma(p::Parameters) = p.sigma.σ
default_popsize(n) = 4 + floor(Int, 3*log(n))
function Parameters(x0, σ0;
                       lower = nothing, upper = nothing,
                       constraints = _constraints(lower, upper),
                       seed = rand(UInt),
                       parallel_evaluation = false,
                       multi_threading = false,
                       noise_handling = nothing,
                       rng = MersenneTwister(UInt(seed)),
                       popsize = default_popsize(length(x0)))
    n = length(x0)
    weights = RecombinationWeights(popsize)
    sigma = Sigma(σ0, n, weights.μeff)
    cov = Covariance(n, weights.μeff)
    finalize_negative_weights!(weights, cov, n)
    Parameters(n, popsize, backtransform(constraints, x0),
                  sigma,
                  cov,
                  weights,
                  constraints,
                  noise_handling,
                  parallel_evaluation,
                  multi_threading,
                  UInt(seed),
                  rng)
end
function update!(p::Parameters, y, perm)
    sample_mean = weighted_average(y, perm, p.weights)
    p.mean .+= sigma(p) * sample_mean
    update!(p.sigma, sample_mean, p.cov)
    update!(p.cov, sample_mean, y, perm, p.weights, p.sigma.h)
    p
end
sample(p) = unwhiten(p.cov.C, randn(p.rng, p.n, p.λ))
function compute_input(p, y)
    transform!(p.constraints, sigma(p) * y .+ p.mean)
end
function evaluate(p::Parameters, f, input)
    if p.parallel_evaluation
        f(input)
    elseif p.multi_threading
        λ = size(input, 2)
        result = zeros(λ)
        Threads.@threads for i in 1:λ
            result[i] = f(@view(input[:, i]))
        end
        result
    else
        [f(@view(input[:, i])) for i in 1:size(input, 2)]
    end
end
weighted_average(y, perm, w) = sum(w.positive_weights[i] * y[:, perm[i]] for i in 1:w.μ)

struct Optimizer{P,L,S}
    p::P
    logger::L
    stop::S
end
function Base.show(io::IO, ::MIME"text/plain", s::Optimizer{P, L, S}) where {P, L, S}
    println(io, "Optimizer{$(P.name.name),$(L.name.name),$(S.name.name)}")
    print_header(s)
    print_result(s)
end
function Optimizer(x0, s0;
                   lower = nothing,
                   upper = nothing,
                   constraints = _constraints(lower, upper),
                   noise_handling = nothing,
                   popsize = default_popsize(length(x0)),
                   stop = nothing,
                   callback = (o, y, fvals, perm) -> nothing,
                   verbosity = 1,
                   parallel_evaluation = false,
                   multi_threading = false,
                   seed = rand(UInt),
                   logger = BasicLogger(x0,
                                        verbosity = verbosity,
                                        callback = callback),
                   kwargs...)
    p = Parameters(x0, s0, popsize = popsize, seed = seed,
                      lower = lower, upper = upper,
                      parallel_evaluation = parallel_evaluation,
                      multi_threading = multi_threading,
                      noise_handling = noise_handling,
                      constraints = constraints)
    Optimizer(p,
              logger,
              stop === nothing ? Stop(p.n, p.λ; kwargs...) : stop)
end

max_stars_repo_name: jbrea/CMAEvolutionStrategy

max_issues_repo_name: jbrea/CMAEvolutionStrategy

max_forks_repo_name: jbrea/CMAEvolutionStrategy


text:
using Revise
using Basins
using DifferentialEquations
using Plots

# Equations of motion: E. Ott, et al. I Physica D 76 (1994) 384-410
function forced_particle!(du, u, p, t)
γ=0.05  ; x̄ = 1.9  ; f₀=2.3  ; ω =3.5
x₀=1. ; y₀=0.;

    x=u[1]; dx = u[2]; y=u[3]; dy = u[4];

    du[1] = dx
    du[2] = -γ*dx -(-4*x*(1-x^2) + y^2) +  f₀*sin(ω*t)*x₀
    du[3] = dy
    du[4] = -γ*dy -(2*y*(x+x̄)) +  f₀*sin(ω*t)*y₀
end

df = ODEProblem(forced_particle!,rand(4),(0.0,20.0), [])
integ  = init(df, alg=Tsit5(); reltol=1e-9, save_everystep=false)

xres=300
yres=300

# range for forced pend
xg = range(0.,1.2,length=xres)
yg = range(0.,1.2,length=yres)

ω =3.5

# compute basin
#@time bsn = Basins.basin_general_ds(xg, yg, integ; dt=2*pi/ω, idxs=[1,3])

#plot(xg,yg,bsn.basin', seriestype=:heatmap)

function escape_function(u)

    if abs(u[3]) > 20 && abs(u[4])>100 && u[3]*u[4]>0
        @show integ.t
        return 2
    elseif abs(u[3]) < 1e-8 && abs(u[4])<1e-9 && u[3]*u[4]<0
        return 1
    end
    return 0
end

basin_t = zeros(length(xg),length(yg))
for (k,x) in enumerate(xg), (n,y) in enumerate(yg)
    reinit!(integ,[x,0,y,0])
    step!(integ, 2π/ω*10,true)
    while escape_function(integ.u) == 0
        step!(integ, 2π/ω,true)
        integ.t > 200 && break
    end

    basin_t[k,n]=escape_function(integ.u);
end

using JLD

@save "riddle.jld" basin_t


text:
include("plot_all.jl")
include("experiment_runtimes_table_2_and_figure_6.jl")
include("experiment_runtimes_table_1.jl")


text:
struct CronbachAlpha{T <: Real}
    alpha::T
    dropped::Vector{T}
end

function Base.show(io::IO, x::CronbachAlpha)
    @printf(io, "Cronbach's alpha for all items: %.4f\n", x.alpha)
    isempty(x.dropped) && return
    println(io, "\nCronbach's alpha if an item is dropped:")
    for (idx, val) in enumerate(x.dropped)
        @printf(io, "item %i: %.4f\n", idx, val)
    end
end

"""
    cronbachalpha(covmatrix::AbstractMatrix{<:Real})

Calculate Cronbach's alpha (1951) from a covariance matrix `covmatrix` according to
the [formula](https://en.wikipedia.org/wiki/Cronbach%27s_alpha):

```math
\\rho = \\frac{k}{k-1} (1 - \\frac{\\sum^k_{i=1} \\sigma^2_i}{\\sum_{i=1}^k \\sum_{j=1}^k \\sigma_{ij}})
```

where ``k`` is the number of items, i.e. columns, ``\\sigma_i^2`` the item variance,
and ``\\sigma_{ij}`` the inter-item covariance.

Returns a `CronbachAlpha` object that holds:

* `alpha`: the Cronbach's alpha score for all items, i.e. columns, in `covmatrix`; and
* `dropped`: a vector giving Cronbach's alpha scores if a specific item,
  i.e. column, is dropped from `covmatrix`.

# Example
```jldoctest
julia> using StatsBase

julia> cov_X = [10 6 6 6;
                6 11 6 6;
                6 6 12 6;
                6 6 6 13];

julia> cronbachalpha(cov_X)
Cronbach's alpha for all items: 0.8136

Cronbach's alpha if an item is dropped:
item 1: 0.7500
item 2: 0.7606
item 3: 0.7714
item 4: 0.7826
```
"""
function cronbachalpha(covmatrix::AbstractMatrix{<:Real})
    isposdef(covmatrix) || throw(ArgumentError("Covariance matrix must be positive definite."))
    k = size(covmatrix, 2)
    k > 1  || throw(ArgumentError("Covariance matrix must have more than one column."))
    v = vec(sum(covmatrix, dims=1))
    σ = sum(v)
    for i in axes(v, 1)
        v[i] -= covmatrix[i, i]
    end
    σ_diag = sum(i -> covmatrix[i, i], 1:k)

    alpha = k * (1 - σ_diag / σ) / (k - 1)
    if k > 2
        dropped = typeof(alpha)[(k - 1) * (1 - (σ_diag - covmatrix[i, i]) / (σ - 2*v[i] - covmatrix[i, i])) / (k - 2)
                                for i in 1:k]
    else
        # if k = 2 do not produce dropped; this has to be also
        # correctly handled in show
        dropped = Vector{typeof(alpha)}()
    end
    return CronbachAlpha(alpha, dropped)
end


text:
using TimeZones
import Dates

@testset "calendar" begin
    @testset "basic" begin
        @system SCalendar(Calendar, Controller)
        d = Dates.Date(2011, 10, 29)
        t0 = ZonedDateTime(d, tz"Asia/Seoul")
        o = :Calendar => :init => t0
        s = instance(SCalendar; config=o)
        @test s.init' == t0
        @test s.time' == t0
        @test s.date' == d
        @test s.step' == Dates.Hour(1)
        update!(s)
        @test s.time' == t0 + Dates.Hour(1)
        update!(s)
        @test s.time' == t0 + Dates.Hour(2)
    end

    @testset "stop" begin
        @system SCalendarStop(Calendar, Controller)
        t0 = ZonedDateTime(2011, 10, 29, tz"Asia/Seoul")
        t1 = t0 + Dates.Day(1)
        o = :Calendar => (init=t0, last=t1)
        s = instance(SCalendarStop; config=o)
        @test s.init' == t0
        @test s.last' == t1
        @test s.time' == t0
        @test s.step' == Dates.Hour(1)
        @test s.stop' == false
        for i in 1:24
            update!(s)
        end
        @test s.time' == t1
        @test s.stop' == true
    end

    @testset "count" begin
        @system SCalendarCount(Calendar, Controller)
        t0 = ZonedDateTime(2011, 10, 29, tz"Asia/Seoul")
        t1 = t0 + Dates.Day(1)
        n = 24
        o = :Calendar => (init=t0, last=t1)
        s = instance(SCalendarCount; config=o)
        @test s.count' == n
        r1 = simulate(SCalendarCount; config=o, stop=n)
        r2 = simulate(SCalendarCount; config=o, stop=:stop)
        r3 = simulate(SCalendarCount; config=o, stop=:count)
        @test r1 == r2 == r3
    end

    @testset "count nothing" begin
        @system SCalendarCountNothing(Calendar, Controller)
        t0 = ZonedDateTime(2011, 10, 29, tz"Asia/Seoul")
        o = :Calendar => (init=t0, last=nothing)
        s = instance(SCalendarCountNothing; config=o)
        @test s.count' === nothing
        @test s.stop' == false
    end

    @testset "count seconds" begin
        @system SCalendarCountSeconds(Calendar, Controller)
        t0 = ZonedDateTime(2011, 10, 29, tz"Asia/Seoul")
        t1 = t0 + Dates.Hour(1)
        n = 60*60
        o = (
            :Calendar => (init=t0, last=t1),
            :Clock => (step=1u"s",),
        )
        s = instance(SCalendarCountSeconds; config=o)
        @test s.count' == n
    end
end


text:
using GitWorkers

## ---------------------------------------------------------------
run_gitworker_server(;
    sys_root = "SERVER_ROOT",
    url = "REMOTE_URL",
    deb = DEB_MODE
)

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_mono_m6_opt_rho2_7(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_mono_m6_opt_rho2_7!(A_copy)
end

@inline function exp_mono_m6_opt_rho2_7!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=9
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 T2k9
    # Computing Bb2 = x*I+x*A
    coeff1=0.1980540938521463
    coeff2=0.9168156264182638
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.1980540938521463
    coeff2=0.9168156264182638
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-0.039878213296149906
    coeff2=0.3191312115175403
    coeff3=0.9305465620354068
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.30406437703941375
    coeff2=0.9357055778901153
    coeff3=0.031045948902277486
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.012054524614470144
    coeff2=-0.015473430203696856
    coeff3=0.3230094878143697
    coeff4=0.9639291036056933
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.011642655100621168
    coeff2=1.0007182791901643
    coeff3=0.05844498049251375
    coeff4=0.00182483121117615
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-0.001434162774737488
    coeff2=-0.005882107697567988
    coeff3=0.004345083585047115
    coeff4=0.013669299358486304
    coeff5=0.8211621256077539
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.8274052334396728
    coeff2=0.6467470861357121
    coeff3=0.052798886182833585
    coeff4=0.0030548445338270413
    coeff5=4.39542321609363e-5
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-0.000118740658609087
    coeff2=-0.00020274697859323616
    coeff3=0.004184759286380629
    coeff4=-0.13046443012976514
    coeff5=0.8480144820381275
    coeff6=0.8140979919427691
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-0.10696304090753009
    coeff2=0.9888245872897673
    coeff3=0.10240241498554757
    coeff4=0.00845811634632641
    coeff5=0.0002567501403605987
    coeff6=-1.308226007116608e-7
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=-4.49699371062923e-6
    coeff2=4.300156065032907e-6
    coeff3=4.354896741415178e-5
    coeff4=-0.00613552979940033
    coeff5=0.06890626823439797
    coeff6=-0.12050842503745675
    coeff7=0.9947846525992833
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=-8.044483063015104e-8
    coeff2=0.9965080883913857
    coeff3=0.11832377144773065
    coeff4=0.011665314654668413
    coeff5=0.0005336901003739011
    coeff6=1.726903525987751e-5
    coeff7=2.4795445648642135e-8
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots3)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 3
    # Computing T2k9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.9837654928593774
    coeff2=0.8167945494192191
    coeff3=0.41828764759110026
    coeff4=0.1592701070085137
    coeff5=0.03329710767309
    coeff6=0.00850052220135662
    coeff7=0.0010349274617070035
    coeff8=4.1578020415467386e-5
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    # Deallocating B7 in slot 9
    return memslots1 # Returning T2k9
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
# Resources

const resources = Dict{String, String}()

resource(f, name = basename(f)) = (@assert isfile(f); resources[name] = f)

const resroute =
  branch(req -> length(req[:path]) == 1 && haskey(resources, req[:path][1]),
         req -> d(:body => open(read, resources[req[:path][1]]),
                  :headers => Mux.fileheaders(req[:path][1])))

# Server setup

const maintp = Mustache.template_from_file(joinpath(dirname(@__FILE__), "main.html"))

app(f) = req -> render(maintp, d("id"=>Page(f).id))

function page_handler(req)
  id = try parse(req[:params][:id]) catch e @goto fail end
  haskey(pool, id) || @goto fail
  active(pool[id].value) && @goto fail

  return render(maintp, d("id"=>id))

  @label fail
  return d(:body => "Not found",
           :status => 404)
end

function ws_handler(req)
  id = try parse(req[:path][end]) catch e @goto fail end
  client = req[:socket]
  haskey(pool, id) || @goto fail
  p = pool[id].value
  active(p) && @goto fail

  p.sock = client
  @schedule @errs get(handlers(p), "init", identity)(p)
  notify(p.cb)
  while active(p)
    local data
    try
      data = read(client)
    catch e
      if isa(e, ArgumentError) && contains(e.msg, "closed")
        handle_message(p, d("type"=>"close"))
        yield() # Prevents an HttpServer task error (!?)
        return
      else
        rethrow()
      end
    end
    @errs handle_message(p, JSON.parse(String(data)))
  end
  return

  @label fail
  close(client)
end

http_default =
  mux(Mux.defaults,
      resroute,
      page(":id", page_handler),
      Mux.notfound())

ws_default =
  mux(Mux.wdefaults,
      ws_handler)

isprecompiling() = ccall(:jl_generating_output, Cint, ()) == 1

@init if !isprecompiling()
  get(ENV, "BLINK_SERVE", "true") in ("1", "true") || return
  http = Mux.http_handler(Mux.App(http_default))
  delete!(http.events, "listen")
  ws = Mux.ws_handler(Mux.App(ws_default))
  serve(Mux.Server(http, ws), port)
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using Plots

@userplot CirclePlot
@recipe function f(cp::CirclePlot)
    x, y, i = cp.args
    n = length(x)
    inds = circshift(1:n, 1 - i)
    linewidth --> range(0, 10, length = n)
    seriesalpha --> range(0, 1, length = n)
    aspect_ratio --> 1
    label --> false
    x[inds], y[inds]
end

n = 150
t = range(0, 2π, length = n)
x = sin.(t)
y = cos.(t)

anim = @animate for i ∈ 1:n
    circleplot(x, y, i)
end
gif(anim, "anim_fps15.gif", fps = 15)

# %%
dump(anim)

# %%
anim.dir

# %%
first(frame(anim), 5)

# %%
using Images
A = [load(joinpath(anim.dir, fn)) for fn in frame(anim)]
A[77]

# %%
for img in A
    IJulia.clear_output(true)
    display(img)
    sleep(1/15)
end

# %%

max_stars_repo_path: 0009/Plots.Animation object.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0009/Plots.Animation object.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0009/Plots.Animation object.jl

max_forks_repo_name: genkuroki/public


text:
import LambdaTAME: spatial_network

@testset "Alignment Drivers" begin

    A = spatial_network(20,2)
    B = spatial_network(25,2)  #different sizes is good for catching errors

    betas=[0.0,1.0]

    @suppress_out begin
        for postProcessing in [KlauAlgo(),SuccessiveKlauAlgo(),LocalSearch(verbose=false)]

            #=
            # -- TAME method tests -- #
            for method in [ΛTAME_M(),LowRankTAME_M()]
                for profile in [true,false]
                    if method === ΛTAME_M()
                    for matchingMethod in [ΛTAME_GramMatching(),ΛTAME_rankOneMatching()]
                        #@inferred align_matrices(A,B,KlauAlgo();method,matchingMethod,betas,profile)
                        #type stability if broken by profile Bool. intending to fix this by making it a type flag.
                        #@test_nothrow
                        @test_nothrow align_matrices(A,B;postProcessing,method,matchingMethod,betas,profile)
                    end
                    elseif method === LowRankTAME_M()
                        @test_nothrow align_matrices(A,B;postProcessing,method,betas,profile)
                    end
                end
            end
            =#

            # -- LowRankEigenAlign method tests -- #
            @test_nothrow align_matrices(A,B;postProcessing,method=LowRankEigenAlign_M())

        end
    end
end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME