text:
# Autogenerated wrapper script for Ccluster_jll for armv7l-linux-gnueabihf
export libccluster

using Arb_jll
using FLINT_jll
using GMP_jll
using MPFR_jll
JLLWrappers.@generate_wrapper_header("Ccluster")
JLLWrappers.@declare_library_product(libccluster, "libccluster.so.1.1.7")
function __init__()
    JLLWrappers.@generate_init_header(Arb_jll, FLINT_jll, GMP_jll, MPFR_jll)
    JLLWrappers.@init_library_product(
        libccluster,
        "lib/libccluster.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Patches

using AWS
using Dates
using HTTP
using JSON
using GitHub
using Mocking
using OrderedCollections: LittleDict

version = v"1.1.0"
status = 200
headers = Pair[
    "x-amz-id-2" => "x-amz-id-2",
    "x-amz-request-id" => "x-amz-request-id",
    "Date" => "Tue, 16 Jun 2020 21:29:18 GMT",
    "x-amz-bucket-region" => "us-east-1",
    "Content-Type" => "application/xml",
    "Transfer-Encoding" => "chunked",
    "Server" => "AmazonS3"
]

body = """
    <?xml version=\"1.0\" encoding=\"UTF-8\"?>\n
    <ListBucketResult xmlns=\"http://s3.amazonaws.com/doc/2006-03-01/\">
        <Name>sample-bucket</Name>
        <Prefix></Prefix>
        <Marker></Marker>
        <MaxKeys>1000</MaxKeys>
        <IsTruncated>false</IsTruncated>
        <Contents>
            <Key>test.txt</Key>
            <LastModified>2020-06-16T21:37:34.000Z</LastModified>
            <ETag>"d41d8cd98f00b204e9800998ecf8427e"</ETag>
            <Size>0</Size>
            <Owner>
                <ID>id</ID>
                <DisplayName>matt.brzezinski</DisplayName>
            </Owner>
            <StorageClass>STANDARD</StorageClass>
        </Contents>
    </ListBucketResult>
    """

response = HTTP.Messages.Response()

web_access_key = "web_identity_access_key"
web_secret_key = "web_identity_secret_key"
web_sesh_token = "web_session_token"

function _response!(; version::VersionNumber=version, status::Int64=status, headers::Array=headers, body::String=body)
    response.version = version
    response.status = status
    response.headers = headers
    response.body = Vector{UInt8}(body)

    return response
end

_aws_http_request_patch = @patch function AWS._http_request(request::Request)
    return response
end

_cred_file_patch = @patch function dot_aws_credentials_file()
    return ""
end

_config_file_patch = @patch function dot_aws_config_file()
    return ""
end

_web_identity_patch = @patch function AWS._http_request(request)
    creds = Dict(
        "AccessKeyId" => web_access_key,
        "SecretAccessKey" => web_secret_key,
        "SessionToken" => web_sesh_token,
        "Expiration" => string(now(UTC))
    )

    result = Dict("AssumeRoleWithWebIdentityResult" => Dict("Credentials" => creds))

    return HTTP.Response(200, ["Content-Type" => "text/json", "charset" => "utf-8"], body=json(result))
end

_github_tree_patch = @patch function tree(repo, tree_obj; kwargs...)
    if tree_obj == "master"
        return Tree("test-sha", HTTP.URI(), [Dict("path"=>"apis", "sha"=>"apis-sha")], false)
    else
        return Tree("test-sha", HTTP.URI(), [Dict("path"=>"test-2020-01-01.normal.json")], false)
    end
end

end

ext: jl


text:
using Luxor
using Dates

function render(pic::Pic, scale::Int64, max_width::Int64=500, max_height::Int64=500, generation::Int64=0)
    fname = "evolisa.png"

    Drawing(max_width, max_height, fname)
    background("black")
    sethue("black")

    begin
        for i in 1:length(pic.polygons)
            @layer render(pic.polygons[i], scale)
        end

    end

    finish()

    img = load(fname)
    img
end

function render(pol::Polygon, scale::Int64)
    pts = pol.pts
    scaled_points = map((p::Pt) -> Point(round(Int, p.x*scale), round(Int, p.y*scale)), pts)

    brush = pol.brush
    (red, green, blue, alpha) = (brush.r/255, brush.g/255, brush.b/255, brush.alpha/100)
    setcolor(red, green, blue, alpha)

    poly(scaled_points, :fill, close=true)
end

max_stars_repo_name: llrt/evolisa-jl

max_issues_repo_name: llrt/evolisa-jl

max_forks_repo_name: llrt/evolisa-jl


text:
# Monte Carlo expectations

function loglikelihood_mc(a::Distribution, b::Distribution, N::Integer)
    mean(logpdf.(b, rand(a, N)))
end

function loglikelihood_mc(a::HMM, b::HMM, τ::Integer, N::Integer)
    mean([
        HMMBase.loglikelihood(b, rand(a, τ), logl = true, robust = true)
        for _ in 1:N
    ])
end


text:
"""
    create_info(args...; kwargs...)

A type-preserving wrapper around `mne.create_info`.

See [`PyMNE.mne.create_info`](@ref) for the associated MNE docstring.
"""
function create_info(args...; kwargs...)
    # `mne.create_info(...)` on its own gives us back a Julia `Dict`, but we actually
    # want the Python object, so we have to use `pycall`
    return pycall(mne.create_info, PyObject, args...; kwargs...)
end

"""
    get_info(::PyObject)

Extract an `Info` property from an MNE object while preserving Python type.
"""
get_info(o::PyObject) = o."info"

# TODO: add a converting for Julia Dict to Python Info <: Dict


text:
# function constraint_frequency_stab_OPF_MP_PL(pm::_PM.AbstractPowerModel)
#     base_nws = pm.setting["base_list"]
#     cont_nws = pm.setting["Cont_list"]
#     ev_syncarea = pm.setting["syncarea"]
#
#     for (base, nw, br) in cont_nws
#         load = _PM.ref(pm, nw, :load, 1) # get index later on
#         syncarea  = _PM.ids(pm, nw, :arcs_reserves_syn)
#         Phvdcoaux = _PM.var(pm, nw)[:Phvdcoaux]
#         Phvdccaux = _PM.var(pm, nw)[:Phvdccaux]
#         Pconv1 = _PM.var(pm, base, :pconv_tf_fr)
#         Pconv2 = _PM.var(pm, nw, :pconv_tf_fr)
#         Zb1 = _PM.var(pm, nw)[:zb1]
#
#         # for i in syncarea # when there is more than one synchornous areas
#         i = ev_syncarea
#         JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) / load["pd"])
#         JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) / load["pd"])
#         JuMP.@constraint(pm.model,  Phvdccaux[i]  == 0)
#         JuMP.@constraint(pm.model,  Pconv1[3] - Pconv2[3] == Pconv1[4] - Pconv2[4] )
#         JuMP.@constraint(pm.model,  Pconv1[3]== Pconv1[4])
#
#         JuMP.set_upper_bound(Phvdcoaux[i], 72/load["pd"])
#         JuMP.set_upper_bound(Phvdccaux[i], 72/load["pd"])
#         Pg = _PM.var(pm, nw, :Pgg, i)/load["pd"]
#         Pf = _PM.var(pm, nw, :Pff, i)/load["pd"]
#         reserves = _PM.ref(pm, nw, :reserves)
#         bi_bp = Dict([((reserves["syncarea"]),i ) for (i,reserves) in _PM.ref(pm, nw, :reserves)])
#         Td = reserves[bi_bp[i]]["Td"]
#         Tg= reserves[bi_bp[i]]["Tg"]
#         Tf = reserves[bi_bp[i]]["Tf"]
#         H = reserves[bi_bp[i]]["H"]
#         Tcl = reserves[bi_bp[i]]["Tcl"]
#         M = 1 # largest possible Pl=0.15 for now
#
#         z1 = _PM.var(pm, nw, :z1, i)
#         z4 = _PM.var(pm, nw, :z4, i)
#         z11 = _PM.var(pm, nw, :z11, i)
#         z41 = _PM.var(pm, nw, :z41, i)
#         z12 = _PM.var(pm, nw, :z12, i)
#         z42 = _PM.var(pm, nw, :z42, i)
#         e = 1e-06
#         p2r = (50/(2*H))
#
#         ############## frequency minimum in interval 1 ##########################################
#          JuMP.@constraint(pm.model,  Phvdcoaux[i]  +e <=   Pf + (Pg/Tg)*Tf + M*(1-z11))
#          JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + (Pg/Tg)*Tf - M*(z11))
#          JuMP.@constraint(pm.model, 0 <=   Phvdcoaux[i]  + M*(1-z12))
#          JuMP.@constraint(pm.model,  0  >= e +  Phvdcoaux[i]  - M*(z12))
#          JuMP.@constraint(pm.model,  0 <= z11 + z12 - 2*z1)
#          JuMP.@constraint(pm.model,  z11 + z12 - 2*z1 <= 1)
#
#          k11 = _PM.var(pm, nw, :k11, i)
#          k12 = _PM.var(pm, nw, :k12, i)
#
#          JuMP.@constraint(pm.model,k11 ==  (H/50)*(Pg/Tg + Pf/Tf) )
#          JuMP.@constraint(pm.model,k12 ==  (Phvdcoaux[i])/2)
#          JuMP.@constraint(pm.model,k12^2 <=  k11)
#
#          ############## frequency minimum in interval 2 ##########################################
#          JuMP.@constraint(pm.model,  Phvdcoaux[i]  + e <=   Pf + Pg +  M*(1-z41))
#          JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + Pg   - M*(z41))
#          JuMP.@constraint(pm.model, Pf + (Pg/Tg)*(Tf)   <=   Phvdcoaux[i]  + M*(1-z42))
#          JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*(Tf)   >= e +  Phvdcoaux[i]  - M*(z42))
#          JuMP.@constraint(pm.model,  0 <= z41 + z42 - 2*z4)
#          JuMP.@constraint(pm.model,  z41 + z42 - 2*z4 <= 1)
#
#          k41 = _PM.var(pm, nw, :k41, i)
#          k42 = _PM.var(pm, nw, :k42, i)
#          k43 = _PM.var(pm, nw, :k43, i)
#
#          JuMP.@constraint(pm.model, k41 ==  (H/50 - Pf*Tf/4) )
#          JuMP.@constraint(pm.model, k42 ==  (Pg/Tg) )
#          JuMP.@constraint(pm.model, k43 ==  (Phvdcoaux[i] - Pf )/2 )
#          JuMP.@constraint(pm.model,[k41/sqrt(2), k42/sqrt(2), k43, 0] in JuMP.RotatedSecondOrderCone() )
#
#          JuMP.@constraint(pm.model, z1 + z4 == 1)
#      end
# end
#
#
# function constraint_frequency_stab_OPF_MP_FSNS(pm::_PM.AbstractPowerModel)
#     base_nws = pm.setting["base_list"]
#     cont_nws = pm.setting["Cont_list"]
#     ev_syncarea = pm.setting["syncarea"]
#
#     for (base, nw, br) in cont_nws
#         load = _PM.ref(pm, nw, :load, 1)
#         syncarea  = _PM.ids(pm, nw, :arcs_reserves_syn)
#         Phvdcoaux = _PM.var(pm, nw)[:Phvdcoaux]
#         Phvdccaux = _PM.var(pm, nw)[:Phvdccaux]
#         Pconv1 = _PM.var(pm, base, :pconv_tf_fr)
#         Pconv2 = _PM.var(pm, nw, :pconv_tf_fr)
#         Zb1 = _PM.var(pm, nw)[:zb1]
#
#         # for i in syncarea
#         i = ev_syncarea
#         conv_conn = _PM.ref(pm, nw, :bus_arcs_conv, i)
#
#         if pm.setting["FSprotection"] == true
#             # display("FSprotection")
#             Mmax = 72/ load["pd"]
#             # for k in conv_conn
#             #     display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= - (Pconv1[k]) /  load["pd"]) )
#             #     display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  <= - (Pconv1[k]) / load["pd"] + Mmax*(1-Zb1[k]) ) )
#             #     display(JuMP.@constraint(pm.model,  Phvdccaux[i]  >= - Pconv2[k] /  load["pd"] -  Mmax*(1-Zb1[k]) ) )
#             #     display(JuMP.@constraint(pm.model,  Phvdccaux[i]  <= - Pconv2[k] / load["pd"] +   Mmax*(1-Zb1[k]) ) )
#             # end
#             # JuMP.@constraint(pm.model,  sum(Zb1[k] for k in conv_conn)==1 )
#
#             display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_convs_dirct, i)) / load["pd"]) )
#             display(JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_convs_dirct, i)) /load["pd"]) )
#             # display(JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) /load["pd"] + sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) /load["pd"] ) )
#             JuMP.@constraint(pm.model,  Pconv1[3] - Pconv2[3] == Pconv1[4] - Pconv2[4] )
#             JuMP.@constraint(pm.model,  Pconv1[3]== Pconv1[4])
#
#          elseif pm.setting["NSprotection"] == true
#             # display("NSprotection")
#             # display(load["pd"])
#             display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) /load["pd"]))
#             display(JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) / load["pd"]))
#             JuMP.@constraint(pm.model,  Pconv1[3] - Pconv2[3] == Pconv1[4] - Pconv2[4])
#             JuMP.@constraint(pm.model,  Pconv1[3]== Pconv1[4])
#
#         end
#
#         JuMP.set_upper_bound(Phvdcoaux[i], 72/load["pd"])
#         JuMP.set_upper_bound(Phvdccaux[i], 72/load["pd"])
#
#         Pg = _PM.var(pm, nw, :Pgg, i)/load["pd"]
#         Pf = _PM.var(pm, nw, :Pff, i)/load["pd"]
#
#         reserves = _PM.ref(pm, nw, :reserves)
#         bi_bp = Dict([((reserves["syncarea"]),i ) for (i,reserves) in _PM.ref(pm, nw, :reserves)])
#         Td = reserves[bi_bp[i]]["Td"]
#         Tg= reserves[bi_bp[i]]["Tg"]
#         Tf = reserves[bi_bp[i]]["Tf"]
#         Tcl = reserves[bi_bp[i]]["Tcl"]
#         H = reserves[bi_bp[i]]["H"]
#
#         z1 = _PM.var(pm, nw, :z1, i)
#         z2 = _PM.var(pm, nw, :z2, i)
#         z3 = _PM.var(pm, nw, :z3, i)
#         z4 = _PM.var(pm, nw, :z4, i)
#         z11 = _PM.var(pm, nw, :z11, i)
#         z21 = _PM.var(pm, nw, :z21, i)
#         z31 = _PM.var(pm, nw, :z31, i)
#         z41 = _PM.var(pm, nw, :z41, i)
#         z12 = _PM.var(pm, nw, :z12, i)
#         z22 = _PM.var(pm, nw, :z22, i)
#         z32 = _PM.var(pm, nw, :z32, i)
#         z42 = _PM.var(pm, nw, :z42, i)
#         e = 1e-06
#         p2r = (50/(2*H))
#         M = 1
#         # M = ub(Pf + Pg +  Phvdccaux[i])
#        ############## frequency minimum in interval 1 ##########################################
#          JuMP.@constraint(pm.model,  Phvdcoaux[i]  +e <=   Pf + (Pg/Tg)*Tf + M*(1-z11))
#          JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + (Pg/Tg)*Tf - M*(z11))
#          JuMP.@constraint(pm.model, 0 <=   Phvdcoaux[i]  + M*(1-z12))
#          JuMP.@constraint(pm.model,  0  >= e +  Phvdcoaux[i]  - M*(z12))
#          JuMP.@constraint(pm.model,  0 <= z11 + z12 - 2*z1)
#          JuMP.@constraint(pm.model,  z11 + z12 - 2*z1 <= 1)
#
#          k11 = _PM.var(pm, nw, :k11, i)
#          k12 = _PM.var(pm, nw, :k12, i)
#
#          JuMP.@constraint(pm.model,k11 ==  (H/50)*(Pg/Tg + Pf/Tf) )
#          JuMP.@constraint(pm.model,k12 ==  (Phvdcoaux[i])/2)
#          JuMP.@constraint(pm.model,k12^2 <=  k11)
#
#          ############## frequency minimum in interval 2 ##########################################
#          JuMP.@constraint(pm.model,  Phvdcoaux[i]  +e  <=  Pf + (Pg/Tg)*Tcl + M*(1-z21))
#          JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + (Pg/Tg)*Tcl - M*(z21))
#          JuMP.@constraint(pm.model, Pf + (Pg/Tg)*Tf <=   Phvdcoaux[i]  + M*(1-z22))
#          JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*Tf  >= e +  Phvdcoaux[i]  - M*(z22))
#          JuMP.@constraint(pm.model,  0 <= z21 + z22 - 2*z2)
#          JuMP.@constraint(pm.model,  z21 + z22 - 2*z2 <= 1)
#
#          k21 = _PM.var(pm, nw, :k21, i)
#          k22 = _PM.var(pm, nw, :k22, i)
#          k23 = _PM.var(pm, nw, :k23, i)
#
#          JuMP.@constraint(pm.model,k21 ==  (H/50 - Pf*Tf/4))
#          JuMP.@constraint(pm.model,k22 ==  Pg/Tg)
#          JuMP.@constraint(pm.model,k23 ==  (Phvdcoaux[i] - Pf)/2)
#
#          JuMP.@constraint(pm.model,[k21/sqrt(2), k22/sqrt(2), k23, 0] in JuMP.RotatedSecondOrderCone() )
#
#          ############## frequency minimum in interval 3 ##########################################
#          JuMP.@constraint(pm.model,  Phvdcoaux[i]  +e  <=   Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux[i]  + M*(1-z31))
#          JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux[i]  - M*(z31))
#          JuMP.@constraint(pm.model, Pf + (Pg/Tg)*Tcl <=   Phvdcoaux[i]  + M*(1-z32))
#          JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*Tcl  >= e+  Phvdcoaux[i]  - M*(z32))
#          JuMP.@constraint(pm.model,  0 <= z31 + z32 - 2*z3)
#          JuMP.@constraint(pm.model,  z31 + z32 - 2*z3 <= 1)
#
#          k31 = _PM.var(pm, nw, :k31, i)
#          k32 = _PM.var(pm, nw, :k32, i)
#          k33 = _PM.var(pm, nw, :k33, i)
#
#          JuMP.@constraint(pm.model, k31 ==  (H/50 - Pf*Tf/4 + Phvdccaux[i]*Tcl^2/(Td*4) ) )
#          JuMP.@constraint(pm.model, k32 ==  (Pg/Tg + Phvdccaux[i]/Td) )
#          JuMP.@constraint(pm.model, k33 ==  (Phvdcoaux[i] - Pf + Phvdccaux[i]*Tcl/Td)/2 )
#
#          JuMP.@constraint(pm.model,[k31/sqrt(2), k32/sqrt(2), k33, 0] in JuMP.RotatedSecondOrderCone() )
#
#         ############## frequency minimum in interval 4 ##########################################
#         JuMP.@constraint(pm.model,  Phvdcoaux[i]  + e <=   Pf + Pg +  Phvdccaux[i]  + M*(1-z41))
#         JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + Pg +  Phvdccaux[i]  - M*(z41))
#         JuMP.@constraint(pm.model, Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux[i]  <=   Phvdcoaux[i]  + M*(1-z42))
#         JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux[i]   >= e +  Phvdcoaux[i]  - M*(z42))
#         JuMP.@constraint(pm.model,  0 <= z41 + z42 - 2*z4)
#         JuMP.@constraint(pm.model,  z41 + z42 - 2*z4 <= 1)
#
#         k41 = _PM.var(pm, nw, :k41, i)
#         k42 = _PM.var(pm, nw, :k42, i)
#         k43 = _PM.var(pm, nw, :k43, i)
#
#         JuMP.@constraint(pm.model, k41 ==  (H/50 - Pf*Tf/4 - Phvdccaux[i]*(Td + 2*Tcl)/4) )
#         JuMP.@constraint(pm.model, k42 ==  (Pg/Tg) )
#         JuMP.@constraint(pm.model, k43 ==  (Phvdcoaux[i] - Pf - Phvdccaux[i])/2 )
#         JuMP.@constraint(pm.model,[k41/sqrt(2), k42/sqrt(2), k43, 0] in JuMP.RotatedSecondOrderCone() )
#
#          JuMP.@constraint(pm.model, z1 +z2 + z3+ z4 == 1)
#      end
# end
#
# function constraint_frequency_stab_TNEP(pm::_PM.AbstractPowerModel)
#     if pm.setting["Permanentloss"] == true
#         constraint_frequency_stab_TNEP_MP_PL(pm)
#     elseif pm.setting["FSprotection"] == true || pm.setting["NSprotection"] == true
#         constraint_frequency_stab_TNEP_MP_FSNS(pm)
#     end
# end
#
#
#
# function constraint_frequency_stab_OPF_MP_PLdim(pm::_PM.AbstractPowerModel)
#     base_nws = pm.setting["base_list"]
#     cont_nws = pm.setting["Cont_list"]
#     ev_syncarea = pm.setting["syncarea"]
#
#     for (base, nw, br) in cont_nws
#         load = _PM.ref(pm, nw, :load, 1) # get index later on
#         syncarea  = _PM.ids(pm, nw, :arcs_reserves_syn)
#         Phvdcoaux = _PM.var(pm, nw)[:Phvdcoaux]
#         Phvdccaux = _PM.var(pm, nw)[:Phvdccaux]
#         Pconv1 = _PM.var(pm, base, :pconv_tf_fr)
#         Pconv2 = _PM.var(pm, nw, :pconv_tf_fr)
#         Zb1 = _PM.var(pm, nw)[:zb1]
#
#         # for i in syncarea # when there is more than one synchornous areas
#         i = ev_syncarea
#         display("connstraint: post_acdcscopf_nocl_PLdim")
#         JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) / load["pd"])
#         JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) / load["pd"])
#         JuMP.@constraint(pm.model,  Phvdccaux[i]  == 0)
#
#         JuMP.set_upper_bound(Phvdcoaux[i], 72/load["pd"])
#         JuMP.set_upper_bound(Phvdccaux[i], 72/load["pd"])
#         Pg = _PM.var(pm, nw, :Pgg, i)/load["pd"]
#         Pf = _PM.var(pm, nw, :Pff, i)/load["pd"]
#         reserves = _PM.ref(pm, nw, :reserves)
#         bi_bp = Dict([((reserves["syncarea"]),i ) for (i,reserves) in _PM.ref(pm, nw, :reserves)])
#         Td = reserves[bi_bp[i]]["Td"]
#         Tg= reserves[bi_bp[i]]["Tg"]
#         Tf = reserves[bi_bp[i]]["Tf"]
#         H = reserves[bi_bp[i]]["H"]
#
#         M = 1 # largest possible Pl=0.15 for now
#
#         z1 = _PM.var(pm, nw, :z1, i)
#         z4 = _PM.var(pm, nw, :z4, i)
#         z11 = _PM.var(pm, nw, :z11, i)
#         z41 = _PM.var(pm, nw, :z41, i)
#         z12 = _PM.var(pm, nw, :z12, i)
#         z42 = _PM.var(pm, nw, :z42, i)
#         e = 1e-06
#         p2r = (50/(2*H))
#
#         ############## frequency minimum in interval 1 ##########################################
#          JuMP.@constraint(pm.model,  Phvdcoaux[i]  +e <=   Pf + (Pg/Tg)*Tf + M*(1-z11))
#          JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + (Pg/Tg)*Tf - M*(z11))
#          JuMP.@constraint(pm.model, 0 <=   Phvdcoaux[i]  + M*(1-z12))
#          JuMP.@constraint(pm.model,  0  >= e +  Phvdcoaux[i]  - M*(z12))
#          JuMP.@constraint(pm.model,  0 <= z11 + z12 - 2*z1)
#          JuMP.@constraint(pm.model,  z11 + z12 - 2*z1 <= 1)
#
#          k11 = _PM.var(pm, nw, :k11, i)
#          k12 = _PM.var(pm, nw, :k12, i)
#
#          JuMP.@constraint(pm.model,k11 ==  (H/50)*(Pg/Tg + Pf/Tf) )
#          JuMP.@constraint(pm.model,k12 ==  (Phvdcoaux[i])/2)
#          JuMP.@constraint(pm.model,k12^2 <=  k11)
#
#          ############## frequency minimum in interval 2 ##########################################
#          JuMP.@constraint(pm.model,  Phvdcoaux[i]  + e <=   Pf + Pg +  M*(1-z41))
#          JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + Pg   - M*(z41))
#          JuMP.@constraint(pm.model, Pf + (Pg/Tg)*(Tf)   <=   Phvdcoaux[i]  + M*(1-z42))
#          JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*(Tf)   >= e +  Phvdcoaux[i]  - M*(z42))
#          JuMP.@constraint(pm.model,  0 <= z41 + z42 - 2*z4)
#          JuMP.@constraint(pm.model,  z41 + z42 - 2*z4 <= 1)
#
#          k41 = _PM.var(pm, nw, :k41, i)
#          k42 = _PM.var(pm, nw, :k42, i)
#          k43 = _PM.var(pm, nw, :k43, i)
#
#          JuMP.@constraint(pm.model, k41 ==  (H/50 - Pf*Tf/4) )
#          JuMP.@constraint(pm.model, k42 ==  (Pg/Tg) )
#          JuMP.@constraint(pm.model, k43 ==  (Phvdcoaux[i] - Pf )/2 )
#          JuMP.@constraint(pm.model,[k41/sqrt(2), k42/sqrt(2), k43, 0] in JuMP.RotatedSecondOrderCone() )
#
#          JuMP.@constraint(pm.model, z1 + z4 == 1)
#      end
# end


text:
export LowStorageRungeKutta2N
export LSRK54CarpenterKennedy, LSRK144NiegemannDiehlBusch, LSRKEulerMethod

"""
    LowStorageRungeKutta2N(f, RKA, RKB, RKC, Q; dt, t0 = 0)

This is a time stepping object for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

The constructor builds a low-storage Runge-Kutta scheme using 2N
storage based on the provided `RKA`, `RKB` and `RKC` coefficient arrays.

The available concrete implementations are:

  - [`LSRK54CarpenterKennedy`](@ref)
  - [`LSRK144NiegemannDiehlBusch`](@ref)
"""
mutable struct LowStorageRungeKutta2N{T, RT, AT, Nstages} <: AbstractODESolver
    "time step"
    dt::RT
    "time"
    t::RT
    "rhs function"
    rhs!
    "Storage for RHS during the LowStorageRungeKutta update"
    dQ::AT
    "low storage RK coefficient vector A (rhs scaling)"
    RKA::NTuple{Nstages, RT}
    "low storage RK coefficient vector B (rhs add in scaling)"
    RKB::NTuple{Nstages, RT}
    "low storage RK coefficient vector C (time scaling)"
    RKC::NTuple{Nstages, RT}

    function LowStorageRungeKutta2N(
        rhs!,
        RKA,
        RKB,
        RKC,
        Q::AT;
        dt = 0,
        t0 = 0,
    ) where {AT <: AbstractArray}

        T = eltype(Q)
        RT = real(T)

        dQ = similar(Q)
        fill!(dQ, 0)

        new{T, RT, AT, length(RKA)}(RT(dt), RT(t0), rhs!, dQ, RKA, RKB, RKC)
    end
end

function LowStorageRungeKutta2N(
    spacedisc::AbstractSpaceMethod,
    RKA,
    RKB,
    RKC,
    Q::AT;
    dt = 0,
    t0 = 0,
) where {AT <: AbstractArray}
    rhs! =
        (x...; increment) ->
            SpaceMethods.odefun!(spacedisc, x..., increment = increment)
    LowStorageRungeKutta2N(rhs!, RKA, RKB, RKC, Q; dt = dt, t0 = t0)
end

"""
    dostep!(Q, lsrk::LowStorageRungeKutta2N, p, time::Real,
            [slow_δ, slow_rv_dQ, slow_scaling])

Use the 2N low storage Runge--Kutta method `lsrk` to step `Q` forward in time
from the current time `time` to final time `time + getdt(lsrk)`.

If the optional parameter `slow_δ !== nothing` then `slow_rv_dQ * slow_δ` is
added as an additionall ODE right-hand side source. If the optional parameter
`slow_scaling !== nothing` then after the final stage update the scaling
`slow_rv_dQ *= slow_scaling` is performed.
"""
function dostep!(
    Q,
    lsrk::LowStorageRungeKutta2N,
    p,
    time,
    slow_δ = nothing,
    slow_rv_dQ = nothing,
    in_slow_scaling = nothing,
)
    dt = lsrk.dt

    RKA, RKB, RKC = lsrk.RKA, lsrk.RKB, lsrk.RKC
    rhs!, dQ = lsrk.rhs!, lsrk.dQ

    rv_Q = realview(Q)
    rv_dQ = realview(dQ)

    groupsize = 256

    for s in 1:length(RKA)
        rhs!(dQ, Q, p, time + RKC[s] * dt, increment = true)

        slow_scaling = nothing
        if s == length(RKA)
            slow_scaling = in_slow_scaling
        end
        # update solution and scale RHS
        event = Event(device(Q))
        event = update!(device(Q), groupsize)(
            rv_dQ,
            rv_Q,
            RKA[s % length(RKA) + 1],
            RKB[s],
            dt,
            slow_δ,
            slow_rv_dQ,
            slow_scaling;
            ndrange = length(rv_Q),
            dependencies = (event,),
        )
        wait(device(Q), event)
    end
end

@kernel function update!(dQ, Q, rka, rkb, dt, slow_δ, slow_dQ, slow_scaling)
    i = @index(Global, Linear)
    @inbounds begin
        if slow_δ !== nothing
            dQ[i] += slow_δ * slow_dQ[i]
        end
        Q[i] += rkb * dt * dQ[i]
        dQ[i] *= rka
        if slow_scaling !== nothing
            slow_dQ[i] *= slow_scaling
        end
    end
end

"""
    LSRKEulerMethod(f, Q; dt, t0 = 0)

This function returns a [`LowStorageRungeKutta2N`](@ref) time stepping object
for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This method uses the LSRK2N framework to implement a simple Eulerian forward time stepping scheme for the use of debugging.

### References

"""
function LSRKEulerMethod(
    F,
    Q::AT;
    dt = nothing,
    t0 = 0,
) where {AT <: AbstractArray}
    T = eltype(Q)
    RT = real(T)

    RKA = (RT(0),)

    RKB = (RT(1),)

    RKC = (RT(0),)

    LowStorageRungeKutta2N(F, RKA, RKB, RKC, Q; dt = dt, t0 = t0)
end

"""
    LSRK54CarpenterKennedy(f, Q; dt, t0 = 0)

This function returns a [`LowStorageRungeKutta2N`](@ref) time stepping object
for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This uses the fourth-order, low-storage, Runge--Kutta scheme of Carpenter
and Kennedy (1994) (in their notation (5,4) 2N-Storage RK scheme).

### References

    @TECHREPORT{CarpenterKennedy1994,
      author = {M.~H. Carpenter and C.~A. Kennedy},
      title = {Fourth-order {2N-storage} {Runge-Kutta} schemes},
      institution = {National Aeronautics and Space Administration},
      year = {1994},
      number = {NASA TM-109112},
      address = {Langley Research Center, Hampton, VA},
    }
"""
function LSRK54CarpenterKennedy(
    F,
    Q::AT;
    dt = 0,
    t0 = 0,
) where {AT <: AbstractArray}
    T = eltype(Q)
    RT = real(T)

    RKA = (
        RT(0),
        RT(-567301805773 // 1357537059087),
        RT(-2404267990393 // 2016746695238),
        RT(-3550918686646 // 2091501179385),
        RT(-1275806237668 // 842570457699),
    )

    RKB = (
        RT(1432997174477 // 9575080441755),
        RT(5161836677717 // 13612068292357),
        RT(1720146321549 // 2090206949498),
        RT(3134564353537 // 4481467310338),
        RT(2277821191437 // 14882151754819),
    )

    RKC = (
        RT(0),
        RT(1432997174477 // 9575080441755),
        RT(2526269341429 // 6820363962896),
        RT(2006345519317 // 3224310063776),
        RT(2802321613138 // 2924317926251),
    )

    LowStorageRungeKutta2N(F, RKA, RKB, RKC, Q; dt = dt, t0 = t0)
end

"""
    LSRK144NiegemannDiehlBusch((f, Q; dt, t0 = 0)

This function returns a [`LowStorageRungeKutta2N`](@ref) time stepping object
for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This uses the fourth-order, 14-stage, low-storage, Runge--Kutta scheme of
Niegemann, Diehl, and Busch (2012) with optimized stability region

### References

    @article{niegemann2012efficient,
      title={Efficient low-storage Runge--Kutta schemes with optimized stability regions},
      author={Niegemann, Jens and Diehl, Richard and Busch, Kurt},
      journal={Journal of Computational Physics},
      volume={231},
      number={2},
      pages={364--372},
      year={2012},
      publisher={Elsevier}
    }
"""
function LSRK144NiegemannDiehlBusch(
    F,
    Q::AT;
    dt = 0,
    t0 = 0,
) where {AT <: AbstractArray}
    T = eltype(Q)
    RT = real(T)

    RKA = (
        RT(0),
        RT(-0.7188012108672410),
        RT(-0.7785331173421570),
        RT(-0.0053282796654044),
        RT(-0.8552979934029281),
        RT(-3.9564138245774565),
        RT(-1.5780575380587385),
        RT(-2.0837094552574054),
        RT(-0.7483334182761610),
        RT(-0.7032861106563359),
        RT(0.0013917096117681),
        RT(-0.0932075369637460),
        RT(-0.9514200470875948),
        RT(-7.1151571693922548),
    )

    RKB = (
        RT(0.0367762454319673),
        RT(0.3136296607553959),
        RT(0.1531848691869027),
        RT(0.0030097086818182),
        RT(0.3326293790646110),
        RT(0.2440251405350864),
        RT(0.3718879239592277),
        RT(0.6204126221582444),
        RT(0.1524043173028741),
        RT(0.0760894927419266),
        RT(0.0077604214040978),
        RT(0.0024647284755382),
        RT(0.0780348340049386),
        RT(5.5059777270269628),
    )

    RKC = (
        RT(0),
        RT(0.0367762454319673),
        RT(0.1249685262725025),
        RT(0.2446177702277698),
        RT(0.2476149531070420),
        RT(0.2969311120382472),
        RT(0.3978149645802642),
        RT(0.5270854589440328),
        RT(0.6981269994175695),
        RT(0.8190890835352128),
        RT(0.8527059887098624),
        RT(0.8604711817462826),
        RT(0.8627060376969976),
        RT(0.8734213127600976),
    )

    LowStorageRungeKutta2N(F, RKA, RKB, RKC, Q; dt = dt, t0 = t0)
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
using Pkg, DrWatson

begin
	using MonteCarloMeasurements
    using BSplines
	using StanQuap
	using StatisticalRethinking
end

df = CSV.read(sr_datadir("cherry_blossoms.csv"), DataFrame; missingstring = "NA")
df = dropmissing(df, :doy)

num_knots = 15
knot_list = quantile(df.year, range(0, 1, length = num_knots))
basis = BSplineBasis(4, knot_list)
B = basismatrix(basis, df.year)

stan4_7 = "
data {
    int n;
    int k;
    int doy[n];
    matrix[n, k] B;
}
parameters {
    real a;
    vector[k] w;
    real<lower=0> sigma;
}
transformed parameters {
    vector[n] mu;
    mu = a + B * w;
}
model {
    for (i in 1:n) {
        doy[i] ~ normal(mu[i], sigma);
    }
    a ~ normal(100, 10);
    w ~ normal(0, 10);
    sigma ~ exponential(1);
}
";

data = (n = size(B, 1), k = size(B, 2), doy = df.doy, B = B)
init = (a = 100.0, sigma = 20.0)
q4_7s, m4_7s, o4_7s = stan_quap("m4.7s", stan4_7; data, init)

if !isnothing(m4_7s)
      part4_7s = read_samples(m4_7s, :particles)
      part4_7s |> display
      println()
      chns = read_samples(m4_7s)
      chns |> display
end


text:
module C

using Base.Test

# Error free transformations, as described in
# [LaLo05] "Solving Triangular Systems More Accurately and Efficiently", Ph. Langlois and N. Louvet
#
# Notation:
#
#    ⊕ ≡ floating point add
#    ⊖ ≡ floating point subtract
#    ⊗ ≡ floating point multiply
#    ⊘ ≡ floating point divide

function twosum{T<:AbstractFloat}(a::T, b::T)
    x = a+b
    z = x-a
    y = (a-(x-z)) + (b-z)
    (x,y)
end

# Split a floating point number into high and low bits
function split(a::Float64)
    z = a*134217729 #(2^27 + 1)
    x = z - (z-a)
    y = a - x
    (x,y)
end

function fast_split(a::Float64)
    ah = Float64(Float32(a))
    ah, a - ah
end

twoproduct(a,b) = twoproduct(promote(a,b)...)
function twoproduct{T}(a::T,b::T)
    x = a*b
    ah,al = split(a)
    bh,bl = split(b)
    y = al*bl - (((x - ah*bh) - al*bh) - ah*bl)
    (x,y)
end

approx_twoproduct(a,b) = approx_twoproduct(promote(a,b)...)
function approx_twoproduct{T}(a::T,b::T)
    x = a*b
    ah,al = fast_split(a)
    bh,bl = fast_split(b)
    y = al*bl - (((x - ah*bh) - al*bh) - ah*bl)
    (x,y)
end

function approxtwodiv(a,b)
    x = a/b
    v,w = twoproduct(x,b)
    y = (a-v-w)/b
    (x,y)
end

@testset "twosum" begin
    for i=1:1000
        a,b = rand(2)
        x,y = twosum(a,b)

        @test abs(y) < eps(x)
        @test big(a)+big(b) == big(x)+big(y)
    end
end

@testset "twoproduct" begin
    for i=1:1000
        a,b = rand(2)
        x,y = twoproduct(a,b)

        @test abs(y) < eps(x)
        @test big(a)*big(b) == big(x)+big(y)
    end
end

@testset "approxtwodiv" begin
    for i=1:1000
        a,b = rand(2)
        x,y = approxtwodiv(a,b)

        @test abs(y) < eps(x)
        @test abs(big(a)/big(b) - (big(x)+big(y))) < max(eps(a),eps(b))
    end
end

# Want correctly rounded version of
#
#   a * (N-i)/(N-1) + b * (i-1)/(N-1)
#
# Rearrange:
#
# (Ah + Al) * (N-i) + (Bh + Bl) * (i-1)
#
#   R ≡ c + (Mh + Ml) * i
#
#   We have:
#
#       abs(Ml) < eps(Mh) ⇒ abs(Ml⊗i) < 2*eps(Mh⊗i)
#       ... FIXME
#
#   xh,xl = twoproduct(Mh, i)
#   y = xl + Ml*i
#   zh,zl = twosum(c, xh)
#
#

immutable LinSpace{T} <: Range{T}
    a::T
    b::T
    c::T
    Mh::T
    Ml::T
    N::Int
end

macro myshow(ex)
    s = string(ex)
    quote
        println($s * " = " * string(map(x->@sprintf("%a",x), $(esc(ex)))))
    end
end

function LinSpace{T<:AbstractFloat}(a::T,b::T,N::Integer)
    M = (big(b) - big(a)) / (N-1)
    Mh = Float64(M)
    Ml = Float64(M - Mh)
    #@myshow Mh,Ml
    #=
    if abs(a) < abs(b)
        a,b = b,a
    end
    =#
    LinSpace(a,b, a, Mh, Ml, N)
end

LinSpace{T}(a::T,b::T,N) = LinSpace(promote(float(a),float(b))..., N)

Base.size(L::LinSpace) = (L.N,)
Base.linearindexing(L::LinSpace) = Base.LinearFast()

function Base.getindex(L::LinSpace, i::Int)
    # xh+xl = Mh*i
    # yh+yl = Ml*i
    # zh+zl ≈ (Mh+Ml)*(i-1)
    # rh+rl = c + zh
    xh,xl = twoproduct(L.Mh, i-1)
    yh,yl = twoproduct(L.Ml, i-1)
    zh,zl = xh, xl + yh
    rh,rl = twosum(L.c, zh)

    rh + (rl + zl)
    #=
    xl += L.Ml*i
    zh,zl = twosum(L.c, xh)
    zl += xl
    zh + zl
    =#
end

Base.step(L::LinSpace) = L.Mh

function bigindex(L::LinSpace, i::Int)
    big(L.a)*(big(L.N-i)/(L.N-1)) + big(L.b)*(big(i-1)/(L.N-1))
end

function verify(a, b, N, L)
    for i= 1:N
        x = bigindex(L, i)
        d = L[i] - Float64(x)
        if d != 0
            # @show (x - Float64(x))/eps(Float64(x))
            return L, i
        end
    end
    return nothing
end

badly_rounded_cases = []
for i=1:1000
    a = rand()
    b = rand()
    N = rand(1:10000)
    ret = verify(a,b,N, LinSpace(a,b,N))
    if ret !== nothing
        println("(LinSpace($a, $b, $N), $(ret[2]))")
        push!(badly_rounded_cases, ret)
    end
end

L,i = badly_rounded_cases[1]
@myshow bigindex(L,i)
@myshow L[i] - bigindex(L,i)

#=
pf(x) = @printf("%a\n", x)

# Float64 - 53 bits significand
# Float32 - 24 bits significand

a = 0x1.1111111111111p0
b = 0x1.1111111111111p0

af = Float64(Float32(a))
bf = Float64(Float32(b))
ar = a - af
br = b - bf

m  = af*bf # exact
mr = af*ar + bf*br + ar*br # approximate, but with extra bits

#pf(m)
#pf(mr)

@show BigFloat(a) * BigFloat(b)
@show BigFloat(m) + BigFloat(mr)
=#

end

L = C.L


text:
#=

This code is based on routines found in the scipy python library.

The license for scipy is included below:

Copyright (c) 2001, 2002 Enthought, Inc.
All rights reserved.

Copyright (c) 2003-2016 SciPy Developers.
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

  a. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
  b. Redistributions in binary form must reproduce the above copyright
     notice, this list of conditions and the following disclaimer in the
     documentation and/or other materials provided with the distribution.
  c. Neither the name of Enthought nor the names of the SciPy Developers
     may be used to endorse or promote products derived from this software
     without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS
BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.

=#
mutable struct ConvergenceError <: Exception
    msg::AbstractString
end

"""
Given a function `f` and an initial guessed range `x1` to `x2`, the routine
expands the range geometrically until a root is bracketed by the returned values
`x1` and `x2` (in which case zbrac returns true) or until the range becomes
unacceptably large (in which case a `ConvergenceError` is thrown).

##### Arguments

- `f::Function`: The function you want to bracket
- `x1::T`: Initial guess for lower border of bracket
- `x2::T`: Initial guess ofr upper border of bracket
- `;ntry::Int(50)`: The maximum number of expansion iterations
- `;fac::Float64(1.6)`: Expansion factor (higher ⟶ larger interval size jumps)

##### Returns

- `x1::T`: The lower end of an actual bracketing interval
- `x2::T`: The upper end of an actual bracketing interval

##### References

This method is `zbrac` from numerical recipies in C++

##### Exceptions

- Throws a `ConvergenceError` if the maximum number of iterations is exceeded
"""
function expand_bracket(f::Function, x1::T, x2::T;
                        ntry::Int=50, fac::Float64=1.6) where T<:Number

    # x1 <= x2 || throw(ArgumentError("x1 must be less than x2"))

    f1 = f(x1)
    f2 = f(x2)

    for j=1:ntry
        if f1*f2 < 0.0
            return x1, x2
        end

        if abs(f1) < abs(f2)
            x1 += fac*(x1 - x2)
            f1 = f(x1)
        else
            x2 += fac*(x2 - x1)
            f2 = f(x2)
        end
    end

    throw(ConvergenceError("failed to find bracket in $ntry iterations"))
end

expand_bracket(f::Function, x1::T; ntry::Int=50, fac::Float64=1.6) where {T<:Number} =
    expand_bracket(f, 0.9x1, 1.1x1; ntry=ntry, fac=fac)

"""
Given a function `f` defined on the interval `[x1, x2]`, subdivide the
interval into `n` equally spaced segments, and search for zero crossings of the
function. `nroot` will be set to the number of bracketing pairs found. If it is
positive, the arrays `xb1[1..nroot]` and `xb2[1..nroot]` will be filled
sequentially with any bracketing pairs that are found.

##### Arguments

- `f::Function`: The function you want to bracket
- `x1::T`: Lower border for search interval
- `x2::T`: Upper border for search interval
- `n::Int(50)`: The number of sub-intervals to divide `[x1, x2]` into

##### Returns

- `x1b::Vector{T}`: `Vector` of lower borders of bracketing intervals
- `x2b::Vector{T}`: `Vector` of upper borders of bracketing intervals

##### References

This is `zbrack` from Numerical Recepies Recepies in C++
"""
function divide_bracket(f::Function, x1::T, x2::T, n::Int=50) where T<:Number
    x1 <= x2 || throw(ArgumentError("x1 must be less than x2"))

    xs = linspace(x1, x2, n)
    dx = xs[2] - xs[1]

    x1b = T[]
    x2b = T[]

    f1 = f(x1)

    for x in xs[2:end]
        f2 = f(x)

        if f1*f2 <= 0.0
            push!(x1b, x-dx)
            push!(x2b, x)
        end
        f1 = f2
    end

    return x1b, x2b
end

__zero_docstr_arg_ret = """
##### Arguments

- `f::Function`: The function you want to bracket
- `x1::T`: Lower border for search interval
- `x2::T`: Upper border for search interval
- `;maxiter::Int(500)`: Maximum number of bisection iterations
- `;xtol::Float64(1e-12)`: The routine converges when a root is known to lie
  within `xtol` of the value return. Should be >= 0. The routine modifies this to
  take into account the relative precision of doubles.
- `;rtol::Float64(2*eps())`:The routine converges when a root is known to lie
  within `rtol` times the value returned of the value returned. Should be ≥ 0

##### Returns

- `x::T`: The found root

##### Exceptions

- Throws an `ArgumentError` if `[x1, x2]` does not form a bracketing interval
- Throws a `ConvergenceError` if the maximum number of iterations is exceeded

"""

## Bisection

"""
Find the root of the `f` on the bracketing inverval `[x1, x2]` via bisection.

$__zero_docstr_arg_ret

##### References

Matches `bisect` function from scipy/scipy/optimize/Zeros/bisect.c
"""
function bisect(f::Function, x1::T, x2::T; maxiter::Int=500,
                xtol::Float64=1e-12, rtol::Float64=2*eps()) where T<:AbstractFloat

    tol = xtol + rtol*(abs(x1) + abs(x2))

    f1, f2 = f(x1), f(x2)

    if f1 * f2 > 0
        throw(ArgumentError("Root must be bracketed by [x1, x2]"))
    end

    # maybe we got lucky and either x1 or x2 is a root
    if f1 == 0.0
        return x1
    end

    if f2 == 0.0
        return x2
    end

    dm = x2 - x1

    for i=1:maxiter
        dm *= 0.5
        xm = x1 + dm
        fm = f(xm)

        # move bracketing interval up if sign(f(xm)) == sign(f(x1))
        if fm*f1 >= 0.0
            x1 = xm
        end

        if fm == 0.0 || abs(dm) < tol
            return xm
        end
    end

    throw(ConvergenceError("Failed to converge in $maxiter iterations"))
end

## Brent's algorithm
abstract type BrentExtrapolation end

struct BrentQuadratic <: BrentExtrapolation
end

struct BrentHyperbolic <: BrentExtrapolation
end

@inline evaluate(::BrentQuadratic, fcur, fblk, fpre, dpre, dblk) =
    -fcur*(fblk*dblk - fpre*dpre) / (dblk*dpre*(fblk - fpre))

@inline evaluate(::BrentHyperbolic, fcur, fblk, fpre, dpre, dblk) =
    -fcur*(fblk - fpre)/(fblk*dpre - fpre*dblk)

function _brent_body(BE::BrentExtrapolation, f::Function,
                     xa::T, xb::T, maxiter::Int=500,
                     xtol::Float64=1e-12,
                     rtol::Float64=2*eps()) where T<:AbstractFloat
    xpre, xcur = xa, xb
    xblk = fblk = spre = scur = 0.0

    fpre = f(xpre)
    fcur = f(xcur)

    if fpre*fcur > 0
        throw(ArgumentError("Root must be bracketed by [xa, xb]"))
    end

    # maybe we got lucky and x1 or x2 is a root of f
    if fpre == 0.0
        return xpre
    end

    if fcur == 0.0
        return xcur
    end

    for i=1:maxiter
        if fpre*fcur < 0
            xblk = xpre
            fblk = fpre
            spre = scur = xcur - xpre
        end

        if abs(fblk) < abs(fcur)
            xpre = xcur
            xcur = xblk
            xblk = xpre

            fpre = fcur
            fcur = fblk
            fblk = fpre
        end

        tol = xtol + rtol*abs(xcur)
        sbis = (xblk - xcur)/2;

        if fcur == 0.0 || abs(sbis) < tol
            return xcur
        end

        if abs(spre) > tol && abs(fcur) < abs(fpre)
            if xpre == xblk
                # interpolate
                stry = -fcur*(xcur - xpre)/(fcur - fpre)
            else
                # extrapolate
                dpre = (fpre-fcur) / (xpre-xcur)
                dblk = (fblk-fcur) / (xblk-xcur)
                stry = evaluate(BE, fcur, fblk, fpre, dpre, dblk)
            end

            if 2*abs(stry) < min(abs(spre), 3*abs(sbis) - tol)
                # good short step
                spre = scur
                scur = stry
            else
                # bisect
                spre = sbis
                scur = sbis
            end
        else
            # bisect
            spre = sbis
            scur = sbis
        end

        xpre = xcur
        fpre = fcur

        if abs(scur) > tol
            xcur += scur
        else
            xcur += copysign(tol, sbis)
        end

        fcur = f(xcur)
    end

    throw(ConvergenceError("Failed to converge in $maxiter iterations"))
end

"""
Find the root of the `f` on the bracketing inverval `[x1, x2]` via brent's algo.

$__zero_docstr_arg_ret

##### References

Matches `brentq` function from scipy/scipy/optimize/Zeros/bisectq.c
"""
function brent(f::Function, xa::T, xb::T; maxiter::Int=500,
               xtol::Float64=1e-12, rtol::Float64=2*eps()) where T<:AbstractFloat
    _brent_body(BrentQuadratic(), f, xa, xb, maxiter, xtol, rtol)
end

"""
Find a root of the `f` on the bracketing inverval `[x1, x2]` via modified brent

This routine uses a hyperbolic extrapolation formula instead of the standard
inverse quadratic formula. Otherwise it is the original Brent's algorithm, as
implemented in the `brent` function.

$__zero_docstr_arg_ret

##### References

Matches `brenth` function from scipy/scipy/optimize/Zeros/bisecth.c
"""
function brenth(f::Function, xa::T, xb::T; maxiter::Int=500,
                xtol::Float64=1e-12, rtol::Float64=2*eps()) where T<:AbstractFloat
    _brent_body(BrentHyperbolic(), f, xa, xb, maxiter, xtol, rtol)
end

## Ridder's method

"""
Find a root of the `f` on the bracketing inverval `[x1, x2]` via ridder algo

$__zero_docstr_arg_ret

##### References

Matches `ridder` function from scipy/scipy/optimize/Zeros/ridder.c
"""
function ridder(f::Function, xa::T, xb::T; maxiter::Int=500,
                xtol::Float64=1e-12, rtol::Float64=2*eps()) where T<:AbstractFloat
    tol = xtol + rtol*(abs(xa) + abs(xb))

    fa, fb = f(xa), f(xb)

    if fa * fb > 0
        throw(ArgumentError("Root must be bracketed by [xa, xb]"))
    end

    # maybe we got lucky and either xa or xb is a root
    if fa == 0.0
        return xa
    end

    if fb == 0.0
        return xb
    end

    for i=1:maxiter
        dm = 0.5*(xb - xa)
        xm = xa + dm
        fm = f(xm)
        dn = sign(fb - fa)*dm*fm/sqrt(fm*fm - fa*fb)
        xn = xm - sign(dn) * min(abs(dn), abs(dm) - .5*tol)
        fn = f(xn)

        if fn*fm < 0.0
            xa, fa, xb, fb = xn, fn, xm, fm
        elseif fn*fa < 0.0
            xb, fb = xn, fn
        else
            xa, fa = xn, fn
        end

        if fn == 0.0 || abs(xb - xa) < tol
            return xn
        end
    end

    throw(ConvergenceError("Failed to converge in $maxiter iterations"))
end


text:
# https://projecteuler.net/problem=15
# answer: 137846528820

f20 = factorial(big(20))
f40 = factorial(big(40))
println(Int(f40 / f20 / f20))

max_stars_repo_name: the-gigi/project-euler

max_issues_repo_name: the-gigi/project-euler

max_forks_repo_name: the-gigi/project-euler


text:
using FreqTables
using Base.Test

x = repeat(["a", "b", "c", "d"], outer=[100]);
# Values not in order to test discrepancy between index and levels with CategoricalArray
y = repeat(["D", "C", "A", "B"], inner=[10], outer=[10]);

tab = @inferred freqtable(x)
@test tab == [100, 100, 100, 100]
@test names(tab) == [["a", "b", "c", "d"]]
tab = @inferred freqtable(y)
@test tab == [100, 100, 100, 100]
@test names(tab) == [["A", "B", "C", "D"]]
tab = @inferred freqtable(x, y)
@test tab == [30 20 20 30;
              30 20 20 30;
              20 30 30 20;
              20 30 30 20]
@test names(tab) == [["a", "b", "c", "d"], ["A", "B", "C", "D"]]

tab =freqtable(x, y,
               subset=1:20,
               weights=repeat([1, .5], outer=[10]))
@test tab == [2.0 3.0
              1.0 1.5
              3.0 2.0
              1.5 1.0]
@test names(tab) == [["a", "b", "c", "d"], ["C", "D"]]

using CategoricalArrays
cx = CategoricalArray(x)
cy = CategoricalArray(y)

tab = @inferred freqtable(cx)
@test tab == [100, 100, 100, 100]
@test names(tab) == [["a", "b", "c", "d"]]
tab = @inferred freqtable(cy)
@test tab == [100, 100, 100, 100]
@test names(tab) == [["A", "B", "C", "D"]]
tab = @inferred freqtable(cx, cy)
@test tab == [30 20 20 30;
              30 20 20 30;
              20 30 30 20;
              20 30 30 20]
@test names(tab) == [["a", "b", "c", "d"], ["A", "B", "C", "D"]]

using Missings
const ≅ = isequal
mx = Array{Union{String, Missing}}(x)
my = Array{Union{String, Missing}}(y)
mx[1] = missing
my[[1, 10, 20, 400]] = missing

mcx = categorical(mx)
mcy = categorical(my)

tab = freqtable(mx)
tabc = freqtable(mcx)
@test tab == tabc == [99, 100, 100, 100, 1]
@test names(tab) ≅ names(tabc) ≅ [["a", "b", "c", "d", missing]]
tab = freqtable(my)
tabc = freqtable(mcy)
@test tab == tabc == [100, 99, 99, 98, 4]
@test names(tab) ≅ names(tabc) ≅ [["A", "B", "C", "D", missing]]
tab = freqtable(mx, my)
tabc = freqtable(mcx, mcy)
@test tab == tabc == [30 20 20 29 0;
                      30 20 20 29 1;
                      20 30 30 20 0;
                      20 29 29 20 2;
                      0   0  0  0 1]
@test names(tab) ≅ names(tabc) ≅ [["a", "b", "c", "d", missing],
                                  ["A", "B", "C", "D", missing]]

tab = freqtable(mx, skipmissing=true)
tabc = freqtable(mcx, skipmissing=true)
@test tab == tabc == [99, 100, 100, 100]
@test names(tab) ≅ names(tabc) ≅ [["a", "b", "c", "d"]]
tab = freqtable(my, skipmissing=true)
tabc = freqtable(mcy, skipmissing=true)
@test names(tab) ≅ names(tabc) ≅ [["A", "B", "C", "D"]]
@test tab == tabc == [100, 99, 99, 98]
tab = freqtable(mx, my, skipmissing=true)
tabc = freqtable(mcx, mcy, skipmissing=true)
@test tab == tabc == [30 20 20 29;
                      30 20 20 29;
                      20 30 30 20;
                      20 29 29 20]

using DataFrames, CSV
iris = CSV.read(joinpath(Pkg.dir("DataFrames"), "test/data/iris.csv"), categorical=false);
iris[:LongSepal] = iris[:SepalLength] .> 5.0
tab = freqtable(iris, :Species, :LongSepal)
@test tab == [28 22
               3 47
               1 49]
@test names(tab) == [["setosa", "versicolor", "virginica"], [false, true]]
tab = freqtable(iris, :Species, :LongSepal, subset=iris[:PetalLength] .< 4.0)
@test tab == [28 22
               3  8]
@test names(tab) == [["setosa", "versicolor"], [false, true]]

# Issue #5
@test freqtable([Set(1), Set(2)]) == [1, 1]
@test freqtable([Set(1), Set(2)], [Set(1), Set(2)]) == eye(2)


text:
# Tests analytic differentiation on transit_poly.jl:
#include("../src/transit_poly.jl")
# include("../src/transit_poly_struct.jl")
#include("../src/dJv_seriesdk.jl")
#=using PyPlot
if VERSION >= v"0.7"
  using Test
  if !(eval(:(@isdefined skip_plots)))
    skip_plots = false
  end
else
  using Base.Test
  if !(isdefined(:skip_plots))
    skip_plots = false
  end
end
=#
if !(eval(:(@isdefined skip_plots)))
  skip_plots = true
end

# Now, carry out finite-difference derivative computation:
function transit_poly_grad_num(r::T,b::T,u_n::Array{T,1},dq0::T) where {T <: Real}
  dq = big(dq0)
# Make BigFloat versions of r, b & u_n:
  r_big = big(r)
  b_big = big(b); u_big = big.(u_n)
# Compute flux to BigFloat precision:
  tp=transit_poly(r_big,b_big,u_big)
# Now, compute finite differences:
#  tp_grad_big= zeros(BigFloat,2+length(u_n))
  tp_grad_big= zeros(BigFloat,2)
  tp_plus = transit_poly(r_big+dq,b_big,u_big)
  tp_minus = transit_poly(r_big-dq,b_big,u_big)
  tp_grad_big[1] = (tp_plus-tp_minus)*.5/dq
  if r == 1.0 && length(u_n) > 2
    tp_grad_big[1] = (tp_plus-tp)/dq
  end
  tp_plus = transit_poly(r_big,b_big+dq,u_big)
  if b_big > dq
    tp_minus = transit_poly(r_big,b_big-dq,u_big)
    width = 2dq
  else
    tp_minus = tp
    width = dq
  end
  tp_grad_big[2] = (tp_plus-tp_minus)/width
#  for i=1:length(u_n)
#    u_tmp = copy(u_big); u_tmp[i] += dq
#    tp_plus = transit_poly(r_big,b_big,u_tmp)
#    u_tmp[i] -= 2dq
#    tp_minus = transit_poly(r_big,b_big,u_tmp)
#    tp_grad_big[2+i] = (tp_plus-tp_minus)*.5/dq
#  end
return convert(Float64,tp),convert(Array{Float64,1},tp_grad_big)
end

include("test_transit_poly_gradient_case.jl")

@testset "transit_poly_gradient" begin
# No limb-darkening:
#global u_n = [0.0]
#global test_name = "uniform"
#include("test_transit_poly_gradient_case.jl")
#fgrid0 = fgrid
fgrid0 = test_transit_poly_gradient_case([0.0],"uniform")
#read(STDIN,Char)
# Linear limb-darkening:
#global u_n = [1.0]
#global test_name = "linear"
#include("test_transit_poly_gradient_case.jl")
#fgrid1 = fgrid
fgrid1 = test_transit_poly_gradient_case([1.0],"linear")
#read(STDIN,Char)
# Quadratic limb-darkening:
#global u_n = [2.0,-1.0]
#global test_name = "quadratic"
#include("test_transit_poly_gradient_case.jl")
#fgrid2 = fgrid
fgrid2 = test_transit_poly_gradient_case([2.0,-1.0],"quadratic")
#read(STDIN,Char)
# Cubic limb-darkening:
#global u_n = [3.0,-3.0,1.0]
#global test_name = "cubic"
#include("test_transit_poly_gradient_case.jl")
#fgrid3 = fgrid
fgrid3 = test_transit_poly_gradient_case([3.0,-3.0,1.0],"cubic")
fgrid4 = test_transit_poly_gradient_case(ones(4)*0.1,"4th order")
fgrid5 = test_transit_poly_gradient_case(ones(5)*0.1,"5th order")
# Now, for arbitrary limb-darkening:
#global u_n = ones(10)*0.1
#global test_name = "10th order"
#include("test_transit_poly_gradient_case.jl")
#fgridn = fgrid
fgridn = test_transit_poly_gradient_case(ones(10)*0.1,"10th order")

if ~skip_plots
  r0 = [0.01,0.1,0.5,0.999999,1.0,1.000001,2.0,10.,100.0]
  r0_name =["0.01","0.1","0.5","0.99","1","1.01","2","10","100"]
  nb = 50

  epsilon = 1e-9; delta = 1e-3
  fig,axes = subplots(3,3)
  for i=1:length(r0)
    global r=r0[i]
    ax = axes[i]
    if r < 1.0
      bgrid = abs.([linearspace(0.0,epsilon,nb); linearspace(epsilon,delta,nb); linearspace(delta,r-delta,nb);
       -logarithmspace(log10(delta),log10(epsilon),nb) .+ r; linearspace(r-epsilon,r+epsilon,nb); logarithmspace(log10(epsilon),log10(delta),nb) .+ r;
       linearspace(r+delta,1-r-delta,nb); -logarithmspace(log10(delta),log10(epsilon),nb) .+ (1-r); linearspace(1-r-epsilon,1-r+epsilon,nb);
       logarithmspace(log10(epsilon),log10(delta),nb) .+ (1-r); linearspace(1-r+delta,1+r-delta,nb); -logarithmspace(log10(delta),log10(epsilon),nb) .+ (1+r);linearspace(1+r-epsilon,1+r-1e-13,nb)])
       nticks = 14
       xticknames=[L"$10^{-13}$",L"$10^{-12}$",L"$10^{-3}$",L"$r-10^{-3}$",L"$r-10^{-12}$",L"$r+10^{-12}$",L"$r+10^{-3}$",
       L"$1-r-10^{-3}$",L"$1-r-10^{-12}$",L"$1-r+10^{-12}$",L"$1-r+10^{-3}$",L"$1+r-10^{-3}$",L"$1+r-10^{-12}$",L"$1+r-10^{-13}$"]
    else
      bgrid = [r-1+1e-13; logarithmspace(log10(epsilon),log10(delta),nb) .+ (r-1); linearspace(r-1+delta,r-delta,nb);
       -logarithmspace(log10(delta),log10(epsilon),nb) .+ r; linearspace(r-epsilon,r+epsilon,nb); logarithmspace(log10(epsilon),log10(delta),nb) .+ r;
       linearspace(r+delta,r+1-delta,nb); -logarithmspace(log10(delta),log10(epsilon),nb) .+ (r+1); r+1-1e-13]
       nticks = 8
       xticknames=[L"$r-1+10^{-13}$",L"$r-1+10^{-3}$",L"$r-10^{-3}$",L"$r-10^{-12}$",L"$r+10^{-12}$",L"$r+10^{-3}$",
       L"$r+1-10^{-3}$",L"$r+1-10^{-13}$"]
    end
    nbgrid = length(bgrid)
    ax[:plot](bgrid,fgrid0[i,1:nbgrid],label="Uniform")
    ax[:plot](bgrid,fgrid1[i,1:nbgrid],label="Linear")
    ax[:plot](bgrid,fgrid2[i,1:nbgrid],label="Quadratic")
    ax[:plot](bgrid,fgrid3[i,1:nbgrid],label="Cubic")
    ax[:plot](bgrid,fgrid4[i,1:nbgrid],label="N=4")
    ax[:plot](bgrid,fgrid5[i,1:nbgrid],label="N=5")
    ax[:plot](bgrid,fgridn[i,1:nbgrid],label="N=10")
    if mod(i,3) == 0
      ax[:set_xlabel]("b")
    end
    if i <= 3
      ax[:set_ylabel]("F")
    end
    ax[:legend](loc="upper left",fontsize=6)
    ax[:set_title](string("r = ",r0[i]),fontsize=6)
  end
end
end


text:
using SMC, Base.Test
srand(125)

## Linear Gaussian

dx, dy = 10, 6
A = randn(dx,dx); A /= 0.9norm(A)
B = randn(dy,dx); B /= 1.1norm(B)
Q = randn(dx,dx); Q *= Q'; Q += 0.5*eye(dx); Q += Q'; Q /= 5
R = randn(dy,dy); R *= R'; R += 0.5*eye(dy); R += R'; R /= 20

lg  = LinearGaussian(A,B,Q,R)
hmm = HMM(lg)
x0  = randn(dx)

K = 50
N = 500
(states, observations) = generate(lg, x0, K)

srand(155)
x00 = x0+randn(dx)/5
@time kf = kalmanfilter(lg, observations, x00, eye(dx))
println("KF    : $(norm(kf.means-states)/norm(states))")
srand(521)
@time ks = kalmansmoother(lg, observations, kf)
println("KS    : $(norm(ks.means-states)/norm(states))")

srand(155)
@time (psf, ess) = particlefilter(hmm, observations, N, bootstrapprop(lg))

@test length(psf)==K

pfm  = mean(psf)
pfmm = zeros(dx,K)
for k in 1:K
    pfmm[:,k] = pfm[k]
end

@test norm(pfmm-states)/norm(states) < 1.2
println("PF    : $(norm(pfmm-states)/norm(states))")

srand(521)
@time pslffbs = particlesmoother_lffbs(hmm, psf, 15)

pfm  = mean(pslffbs)
pfmm = zeros(dx,K)
for k in 1:K
    pfmm[:,k] = pfm[k]
end
@test norm(pfmm-states)/norm(states) < 0.9
println("PSLFBS: $(norm(pfmm-states)/norm(states))")

srand(521)
@time (pslbbis, ess) = particlesmoother_lbbis(hmm, observations,
                                              psf, bootstrapprop(lg))

psm4  = mean(pslbbis)
psmm4 = zeros(dx,K)
for k in 1:K
    psmm4[:,k] = psm4[k]
end

@test norm(psmm4-states)/norm(states) < 1.
println("PSBISL: $(norm(psmm4-states)/norm(states))")

srand(521)
@time (psllbbis, ess) = particlesmoother_llbbis(hmm, observations, psf,
                                                5round(Int,log(N)),
                                                bootstrapprop(lg))

psm5  = mean(psllbbis)
psmm5 = zeros(dx,K)
for k in 1:K
    psmm5[:,k] = psm5[k]
end

@test norm(psmm5-states)/norm(states) < 1.0
println("PSBISLL: $(norm(psmm5-states)/norm(states))")

srand(521)
@time pss_fh = particlesmoother_fearnhead_lg(lg, observations, psf)

psm  = mean(pss_fh)
psmm = zeros(dx,K)
for k in 1:K
    psmm[:,k] = psm[k]
end

@test norm(psmm - states)/norm(states) < 1.5 # quite poor...
println("PS_FH  : $(norm(psmm-states)/norm(states))")


text:
#MODULE: SERVICES
import HttpCommon
import DataFrames
import CSV

using DataFrames: DataFrame

export Response, read, dataframe!, stats!, warnings, drilldown, pivot

#GLOBALS
const STATS        = "X-Clarus-Stats"
const WARNINGS     = "X-Clarus-Messages"
const MESSAGES     = "messages"

#REQUEST STATS
const GRID_TOTAL   = "GridTotal"
const REQUEST_ID   = "RequestId"
const GRID_ID      = "GridId"
const CALC_TIME    = "CalcTime"
const TASK_TIME    = "TaskTime"
const NUM_MTM      = "NumMTM"
const NUM_TRADES   = "NumTrades"
const NUM_WARNINGS = "NumWarnings"

#END-GLOBALS

mutable struct Response
  httpresponse::HttpCommon.Response
  stats::Nullable{Dict{String,String}}
  dataframe::Nullable{DataFrame}
  warnings::Nullable{Any}
  Response(r::HttpCommon.Response) = new(r,Nullable{Dict{String,String}}(),Nullable{DataFrames.DataFrame}())
end

function text(r::Response)
  return String(r.httpresponse.data)
end

function data(r::Response)
  return r.httpresponse.data
end

function headers(r::Response)
  return r.httpresponse.headers
end

function status(r::Response)
  return r.httpresponse.status
end

function read(r::Response,  sink::Type=DataFrame, ; dateformat::Dates.DateFormat=Dates.ISODateFormat)
  return CSV.read(IOBuffer(r.httpresponse.data),sink;delim = ',', dateformat=dateformat)
end

function Base.print(io::IO, r::Response)
  return Base.showall(dataframe!(r))
end

function dataframe!(r::Response)
  if isnull(r.dataframe)
    r.dataframe = Nullable(read(r,DataFrame))
  end
  return get(r.dataframe)
end

function Base.show(io::IO, r::Response)
  return Base.show(dataframe!(r))
end

function stats!(r::Response)
  if isnull(r.stats)
    head = r.httpresponse.headers
    if STATS in keys(head)
      d = Dict()
      for (i,s) in enumerate(split(head[STATS],';'))
          pair = split(s,'=')
          push!(d, strip(pair[1])=> strip(pair[2]))
      end
      r.stats = Nullable(d)
    else
      return Dict()
    end
  end
  return get(r.stats)
end

#getters

function stat!(r::Response,stat::String)
  return get(stats!(r),stat,nothing)
end

function total!(r::Response)
  return stat!(r, GRID_TOTAL)
end

function requestid!(r::Response)
  return stat!(r,REQUEST_ID)
end

function gridid!(r::Response)
  return stat!(r,GRID_ID)
end
function warnings!(r::Response)
  if isnull(r.warnings)
    r.warnings = Nullable(get(r.httpresponse.headers,WARNINGS,nothing))
  end
  return r.warnings
end

function pivot(grid::Response,rowAxis="Currency",colAxis="SubType",ccy="USD",view="Latest",output=credentials.default_outputtype)
  return Response(gridrequest(("GridId"=>gridid!(grid),"Row"=>rowAxis,"Col"=>colAxis,"reportCcy"=>ccy,"View"=>view),output))
end

function drilldown(grid::Response,row="Total",col="Total",view="Default",output=credentials.default_outputtype)
  return Response(gridrequest(("GridId"=>gridid!(grid),"DrilldownRow"=>row,"DrilldownCol"=>col,"DrilldownView"=>view),output))
end

function gridrequest(params,output=credentials.default_outputtype)
  headers = ("Content-Type"=>"application/json","User-Agent"=>credentials.user_agent)
  return Requests.post(url(UTIL_SERVICE,"Grid",output),json = Dict(params),headers = Dict(headers))
end


text:
# Copyright (c) 2019 Guido Kraemer
# Copyright (c) 2019 Simon Castano
# Distributed under the MIT software license, see the accompanying
# file COPYING or http://www.opensource.org/licenses/mit-license.php.

include("../script/Script.jl")
include("../script/Witness.jl")
include("Outpoint.jl")
include("TxIn.jl")
include("TxOut.jl")

abstract type Transaction end

"""
Bitcoin transactions are broadcast between peers in a serialized byte format,
called raw format. It is this form of a transaction which is SHA256(SHA256())
hashed to create the TXID and, ultimately, the merkle root of a block
containing the transaction—making the transaction format part of the consensus
rules.

A raw transaction has the following top-level format:

- Transaction version number; currently version 1 or 2. Programs creating
transactions using newer consensus rules may use higher version numbers.
Version 2 means that BIP 68 applies.
- A marker which MUST be a 1-byte zero value: `0x00` (BIP 141)
- A flag which MUST be a 1-byte non-zero value: `0x01` (BIP 141)
- Transaction inputs
- Transaction outputs
- A time (Unix epoch time) or block number (BIP 68)

A transaction may have multiple inputs and outputs, so the `TxIn` and `TxOut`
structures may recur within a transaction.
"""
struct Tx <: Transaction
    version     :: UInt32
    marker      :: UInt8
    flag        :: UInt8
    inputs      :: Vector{TxIn}
    outputs     :: Vector{TxOut}
    witnesses   :: Vector{Witness}
    locktime    :: UInt32
end

"""
    Tx(io::IOBuffer) -> Tx

Parse an `IOBuffer` to a `Tx`
"""
function Tx(io::IOBuffer)

    version         =   ltoh(read(io, UInt32))

    x               =   CompactSizeUInt(io).value
    x == zero(x)    ?   segwit = true : segwit = false

    if segwit
        marker, flag    =   x, read(io, UInt8)
        @assert flag    ==  0x01
        txin_count      =   CompactSizeUInt(io).value
    else
        marker, flag    =   0xff, 0xff
        txin_count      =   x
    end

    inputs  = TxIn[TxIn(io)     for i ∈ 1:txin_count]
    outputs = TxOut[TxOut(io)   for i ∈ 1:CompactSizeUInt(io).value]

    if segwit
        witness_count = txin_count
        @assert witness_count > 0
        witnesses = Witness[Witness(io) for i ∈ 1:witness_count]
    else
        witnesses = [Witness()]
    end

    locktime = ltoh(read(io, UInt32))

    return Tx(version, marker, flag, inputs, outputs, witnesses, locktime)
end

function Base.show(io::IO, tx::Tx)
    println(io, "Transaction:      " * bytes2hex(reverse(hash256(tx))))
    println(io, "  Version:        " * string(tx.version,           base = 10))
    println(io, "  Marker:         " * string(tx.marker,            base = 10))
    println(io, "  Flag:           " * string(tx.flag,              base = 10))
    println(io, "  Input counter:  " * string(length(tx.inputs),    base = 10))
    println(io, "  Output counter: " * string(length(tx.outputs),   base = 10))
    println(io, "  Lock time:      " * string(tx.locktime,          base = 10))
end

"""
    serialize(tx::Tx) -> Vector{UInt8}

Returns the byte serialization of the transaction
"""
function serialize(tx::Tx)
    result = bytes(tx.version, len=4, little_endian=true)

    (tx.marker, tx.flag) == (0xff, 0xff) ? bip141 = false : bip141 = true
    bip141 ? append!(result, [tx.marker, tx.flag]) : nothing

    l = CompactSizeUInt(length(tx.inputs))
    append!(result, serialize(l))
    for input in tx.inputs
        append!(result, serialize(input))
    end

    l = CompactSizeUInt(length(tx.outputs))
    append!(result, serialize(l))
    for output in tx.outputs
        append!(result, serialize(output))
    end

    if bip141
        for i ∈ 1:length(tx.inputs)
            append!(result, serialize(tx.witnesses[i]))
        end
    end
    append!(result, bytes(tx.locktime, len=4, little_endian=true))
    return result
end

total_output(tx::Tx) = sum(x -> x.amount, tx.outputs)

"""
    iscoinbase(tx::Tx) -> Bool

Returns whether this transaction is a coinbase transaction or not
"""
function iscoinbase(tx::Tx)
    outpoint = tx.inputs[1].prevout
    length(tx.inputs) == 1 && outpoint.txid == fill(0x00, 32) && outpoint.index == 0xffffffff
end

"""
    coinbase_height(tx::Tx) ->

Returns the height of the block this coinbase transaction is in
Returns an `AssertionError` if `tx` isn't a coinbase transaction
"""
function coinbase_height(tx::Tx)
    @assert iscoinbase(tx) "This is not a coinbase transaction"
    height_bytes = tx.inputs[1].scriptsig.data[1]
    return to_int(height_bytes, little_endian=true)
end


text:
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.
#
# This file depends upon the Sudoku module being defined first.
using Test
using Random

function solve_random_puzzles()
    rng = MersenneTwister(123456)
    for i = 1:50
        solution, puzzle = Sudoku.random_puzzle(3,rng,1000,25)
        @test Sudoku.valid_puzzle(Sudoku.as_values(solution))
        # Exclusively rules-based solution
        puzzle1 = deepcopy(puzzle)
        iteration, uncertainty = Sudoku.iterative_solve!(puzzle1,1)
        @test iteration > 0
        @test uncertainty == 0
        @test Sudoku.valid_puzzle(Sudoku.as_values(puzzle1))
        for j = 1:length(puzzle1.grid)
            @test Sudoku.get_value(puzzle1.grid[j]) == Sudoku.get_value(solution.grid[j])
        end
        # Repeat using exclusively rules-evaluation in backtrack
        puzzle2 = deepcopy(puzzle)
        result_puzzle2 = Sudoku.backtrack_solve(puzzle2,1,1,2)
        @test length(result_puzzle2) == 1
        @test Sudoku.valid_puzzle(Sudoku.as_values(result_puzzle2[1]))
        for j = 1:length(result_puzzle2[1].grid)
            @test Sudoku.get_value(result_puzzle2[1].grid[j]) == Sudoku.get_value(solution.grid[j])
        end
        # Repeat using logic and not actually backtracking
        result_puzzle3 = Sudoku.backtrack_solve(puzzle2,1,2,2)
        @test length(result_puzzle3) == 1
        @test Sudoku.valid_puzzle(Sudoku.as_values(result_puzzle3[1]))
        for j = 1:length(result_puzzle3[1].grid)
            @test Sudoku.get_value(result_puzzle3[1].grid[j]) == Sudoku.get_value(solution.grid[j])
        end
    end
end

function solve_example_puzzles()
    # Puzzle example from
    # https://en.wikipedia.org/wiki/Sudoku accessed 1/1/2022
    # Puzzle image license is CC0 (Public Domain Dedication) per
    # https://en.wikipedia.org/wiki/Sudoku#/media/File:Sudoku_Puzzle_by_L2G-20050714_standardized_layout.svg
    puzzle_values = [
        5 3 0 0 7 0 0 0 0;
        6 0 0 1 9 5 0 0 0;
        0 9 8 0 0 0 0 6 0;
        8 0 0 0 6 0 0 0 3;
        4 0 0 8 0 3 0 0 1;
        7 0 0 0 2 0 0 0 6;
        0 6 0 0 0 0 2 8 0;
        0 0 0 4 1 9 0 0 5;
        0 0 0 0 8 0 0 7 9]
    puzzle1 = Sudoku.SolvablePuzzle(3)
    Sudoku.assign_values!(puzzle1,puzzle_values)
    iterations, uncertainty = Sudoku.iterative_solve!(puzzle1,2)
    @test uncertainty == 0
    @test Sudoku.satisfies(puzzle_values,Sudoku.as_values(puzzle1))
    @test Sudoku.valid_puzzle(Sudoku.as_values(puzzle1))

    # Repeat w/ simple backtracking
    puzzle2 = Sudoku.SolvablePuzzle(3)
    Sudoku.assign_values!(puzzle2,puzzle_values)
    results = Sudoku.backtrack_solve(puzzle2,8,81,2)
    @test length(results) == 1
    @test Sudoku.satisfies(puzzle_values,Sudoku.as_values(results[1]))
    @test Sudoku.valid_puzzle(Sudoku.as_values(results[1]))
    @test Sudoku.as_values(results[1]) == Sudoku.as_values(puzzle1)

    # Puzzle example from
    # https://sourceforge.net/projects/winsudoku/ accessed Jan 30, 2022
    # Example taken from BSD-licensed project
    puzzle_values = [
        0 2 1 0 6 0 0 0 8;
        8 0 0 0 0 0 7 5 0;
        0 0 0 0 0 1 0 0 0;
        0 0 4 9 0 3 0 0 0;
        9 0 0 1 0 0 0 3 0;
        1 0 0 0 0 8 0 0 0;
        0 0 0 0 0 0 5 0 0;
        0 7 0 0 0 9 0 0 6;
        0 8 3 7 0 0 1 0 4]
    puzzle = Sudoku.SolvablePuzzle(3)
    Sudoku.assign_values!(puzzle,puzzle_values)
    # Backtrack solve using
    #   * up to 8-fold compound rules to simplify backtracking
    #   * up to 57 recursive steps (24 knowns)
    #   * solution cap of 2 to detect a unique solution
    results = Sudoku.backtrack_solve(puzzle,8,81-24,2)
    @test length(results) == 1
    if length(results) == 1
        @test Sudoku.valid_puzzle(Sudoku.as_values(results[1]))
        @test Sudoku.satisfies(Sudoku.as_values(puzzle),Sudoku.as_values(results[1]))
    end
end

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia


text:
"""
Analyze filtered fitnesses for emid data. Code created as part of the revision of
   SDERGM paper
"""

## Load dataj
using Utilities,AReg,StaticNets,JLD,MLBase,StatsBase#,DynNets
using PyCall; pygui(:qt); using PyPlot

#estimate and save for half the datase (after LTRO) or whole data?
halfPeriod = false
fold_Path =  "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/"
loadFilePartialName = "Weekly_eMid_Data_from_"
halfPeriod ? periodEndStr =  "2012_03_12_to_2015_02_27.jld" : periodEndStr =
                                                "2009_06_22_to_2015_02_27.jld"
@load(fold_Path*loadFilePartialName*periodEndStr, AeMidWeekly_T,banksIDs,
      inactiveBanks,
      YeMidWeekly_T,weekInd,datesONeMid ,degsIO_T,strIO_T)

N_test = 50

# single snapshot estimates
matY_T = YeMidWeekly_T[1:N_test,1:N_test,3:end]
 matA_T = matY_T.>0
    N = length(matA_T[1,:,1])
    T = size(matA_T)[3]
    Ttrain = 100#round(Int, T/2) #70 106 #
    threshVar = 0.00#5
    degsIO_T = [sumSq(matA_T,2);sumSq(matA_T,1)]
 allFitSS =  StaticNets.estimate( StaticNets.SnapSeqNetDirBin1(degsIO_T);
            identPost = false,identIter= true )

#plot(allFitSS[:,1])

#Estimate gas model on  windows of length Ttrain
model = SdErgmDirBin1(degsIO_T[:,1:Ttrain],"FISHER-DIAG", "ONE_PAR_ALL")
gasParEstOnTrain,~ = estimateTarg(model;SSest = allFitSS[:,1:Ttrain] )

#f(x) = DynNets.score_driven_filter_or_dgp(model, x;groupsInds = model.groupsInds)[2]
GBA = 1
allpar0 = array_2_vec_all_par(model, gasParEstOnTrain)
BA_re_hat = allpar0[end-1:end]
function UnrestrAB(vecRePar::Array{<:Real,1}) #restrict a vector of only A and B
      diag_B_Re = vecRePar[1:GBA]
      diag_A_Re = vecRePar[GBA+1:end]
      diag_B_Un = log.(diag_B_Re ./ (1 .- diag_B_Re ))
      diag_A_Un = log.(diag_A_Re)
      vecUnPar =  [diag_B_Un; diag_A_Un]
      return vecUnPar
end

BA_un_hat =  UnrestrAB(BA_re_hat)
w_hat = allpar0[1:2*N]
#%%
function loglike(Model::DynNets.SdErgmDirBin1, degs_t::Array{<:Real,1},
                        f_t::Array{<:Real,1})

    thetas_mat_t_exp, exp_mat_t = StaticNets.expMatrix2(StaticNets.fooErgmBin1,f_t)
    exp_deg_t = sum(exp_mat_t,dims = 2)

    loglike_t = sum(f_t.*degs_t) -  sum(UpperTriangular(log.(1 .+ thetas_mat_t_exp))) #  sum(log.(1 + thetas_mat_t_exp))
    return loglike_t
 end

using ForwardDiff
obsT = model.obsT
T_train = size(obsT)[2]
f_T(x) = logLike_T(model::SdErgmDirBin1, obsT, x)
# @time ForwardDiff.hessian(f_t, allpar0)
#@time FiniteDiff.finite_difference_hessian(f_T, allpar0)

covA0 =  ForwardDiff.hessian(f_T, allpar0)./T_train
# estimate outer product of scores
vec_of_f_t =[x -> logLike_t(model::SdErgmDirBin1, obsT[:, 1:t], x) for t in 1:Ttrain]
# run once to precompile the functions
[f(allpar0) for f in vec_of_f_t]
vec_g_t = [ForwardDiff.gradient(f, allpar0) for f in vec_of_f_t[2:end]]
global covB0 = zeros(size(covA0))
 for g_t in vec_g_t
      global covB0 += g_t * g_t'
 end
 global covB0 = covB0./ T_train
#%%
diagCorrect = 1e-3
covMatHat = pinv(covA0) * covB0 * pinv(covA0)
      covMatHat = Symmetric(covMatHat)

covMatHat_corr = (covMatHat .+ Diagonal(diagCorrect.*ones(size(covMatHat)[1])))

parDistr = MvNormal(allpar0, diag(covMatHat_corr))

samp = rand(parDistr, 2000)
indsToKeep = (0 .< samp[end-1, :] .< 1) .& (samp[end, :] .>0)
sum(indsToKeep)
#hist(samp[end,:])
N_samp_par = 1000
allpar_sample = samp[:, indsToKeep][:, 1:N_samp_par]

#%%
ftotIO_0 = meanSq(estimateSnapSeq(model; degsIO_T = obsT[:, 1:5]), 2)
ftot_T_est, _ = score_driven_filter_or_dgp(model, allpar0,ftotIO_0 = ftotIO_0)
ftot_T, _ = score_driven_filter_or_dgp(model, allpar_sample[:,3],  ftotIO_0 = ftotIO_0)
plot(ftot_T_est')
#plot(ftot_T')

Nterms = 2*N

global gasFiltPar_conf = zeros(N_samp_par, Nterms, T_train)
      for i= 1:N_samp_par
       # @show(i)
       global gasFiltPar_conf[i,:,:], _ = score_driven_filter_or_dgp(model, allpar_sample[:,i],  ftotIO_0 = ftotIO_0)
      end

indsToRem = [any(isnan.(gasFiltPar_conf[i, : , :])) for i in 1:N_samp_par ]
sum(indsToRem)
gasFiltPar_conf = gasFiltPar_conf[ .!indsToRem, : , : ]

using Statistics
# Compute confidence bands
quant_vals = [0.95, 0.05]
 conf_band = zeros(length(quant_vals),T_train, Nterms)
 for t=1:T_train
     for k=1:Nterms
      filt_t = gasFiltPar_conf[:,k,t]
      conf_band[:,t,k] = Statistics.quantile(filt_t[.!isnan.(filt_t)],quant_vals)
    end
 end

figure()
      legTex = ["DGP";"SDERGM"; "95%"]
      node = 11
      plot(1:T_train,ftot_T_est[node, :],"r")
      plt.fill_between(1:T_train, conf_band[1, :, node], y2 =conf_band[2,:,node],color =(0.9, 0.2 , 0.2, 0.1)  )#, color='b', alpha=.1)

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
module CoreModule

include("ProgramConstants.jl")
include("Dataset.jl")
include("OptionsStruct.jl")
include("Equation.jl")
include("Operators.jl")
include("Options.jl")

import .ProgramConstantsModule:
    CONST_TYPE,
    MAX_DEGREE,
    BATCH_DIM,
    FEATURE_DIM,
    RecordType,
    SRConcurrency,
    SRSerial,
    SRThreaded,
    SRDistributed
import .DatasetModule: Dataset
import .OptionsStructModule: Options
import .EquationModule: Node, copy_node, string_tree, print_tree
import .OptionsModule: Options
import .OperatorsModule:
    plus,
    sub,
    mult,
    square,
    cube,
    pow,
    div,
    log_abs,
    log2_abs,
    log10_abs,
    log1p_abs,
    sqrt_abs,
    acosh_abs,
    neg,
    greater,
    greater,
    relu,
    logical_or,
    logical_and,
    gamma,
    erf,
    erfc,
    atanh_clip

end


text:
"""
    fm = mergespec_spaceonly(ra,f,kmax)

Spectral merging of coarse field `ra` and fine field `f` at wavenumber `kmax`.
"""
function mergespec_spaceonly(ra,f,kmax)

(nx,ny)=size(f);
(nax,nay)=size(ra);

DFTr=zeros(Complex{Float64},nax,nax);
DFTf=zeros(Complex{Float64},nx,nx);

#fft di tutto il campo f
DFTf=fft(f)

#fft di tutto il campo r
DFTr=fft(ra);

DFTr=fftshift(DFTr); #centro la fft
DFTf=fftshift(DFTf); #centro la fft
DFTfm=zeros(Complex{Float64},nx,nx);
DFTr2=zeros(Complex{Float64},nax+1,nax+1);

#DFTr[nax+1,:,:]=NaN;
#DFTr[:,nax+1,:]=NaN;
DFTr2[1:nax,1:nax]=DFTr[:,:]
DFTr2[nax+1,1:nax]=conj(DFTr[1,:]);
DFTr2[1:nax,nax+1]=conj(DFTr[:,1]);

kmax2=kmax^2;

# We need to fix the phases of the large scale field!
# The first pixel of the large field is centered in dxl/2 with dxl=2pi/nax
# The first pixel of the fine field is centered in dxs/2 with dxf=2pi/nax
ddx=2*pi/nax/2-2*pi/nx/2;

 for j=1:nx
     for i=1:nx
            kx=-div(nx,2)+i-1;
            ky=-div(nx,2)+j-1;
            k2=(kx^2+ky^2);
            ir=div(nax,2)+1+kx;
            jr=div(nax,2)+1+ky;
            if(k2<=kmax2)
                DFTfm[i,j]=DFTr2[ir,jr]*exp(-1im*ddx*kx-1im*ddx*ky);;
            else
                DFTfm[i,j]=DFTf[i,j];
            end
     end
end

DFTfm=ifftshift(DFTfm);
fm=ifft(DFTfm);
fm=real(fm);

end


text:
# Autogenerated wrapper script for Zellij_jll for i686-linux-gnu
export zellij

JLLWrappers.@generate_wrapper_header("Zellij")
JLLWrappers.@declare_executable_product(zellij)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        zellij,
        "bin/zellij",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export FxRates

#TODO. Dividir la descarga y la generaración y grabación del dataframe

function ECB_FxRates(ruta::String,from::Date)

    file = HTTP.get("https://www.ecb.europa.eu/stats/eurofxref/eurofxref-hist.zip").body
    r = ZipFile.Reader(IOBuffer(file))

    FxHist = CSV.File(r.files[1], dateformat="yyyy-mm-dd",delim=',',quotechar='"',missingstring="N/A") |> DataFrame

    FxHist = FxHist[FxHist[:,1].>=from,:]

    CSV.write(ruta*"/tDivisa.csv",FxHist,delim=";",decimal=',')

    FxHist
end


text:
"""
forwardsub(L,b)

Solve the lower-triangular linear system with matrix `L` and
right-hand side vector `b`.
"""
function forwardsub(L,b)

n = size(L,1)
x = zeros(n)
x[1] = b[1]/L[1,1]
for i = 2:n
    s = sum( L[i,j]*x[j] for j=1:i-1 )
    x[i] = ( b[i] - s ) / L[i,i]
end

return x
end

"""
backsub(U,b)

Solve the upper-triangular linear system with matrix `U` and
right-hand side vector `b`.
"""
function backsub(U,b)

n = size(U,1)
x = zeros(n)
x[n] = b[n]/U[n,n]
for i = n-1:-1:1
    s = sum( U[i,j]*x[j] for j=i+1:n )
    x[i] = ( b[i] - s ) / U[i,i]
end

return x
end

"""
lufact(A)

Compute the LU factorization of square matrix `A`, returning the
factors.
"""
function lufact(A)

n = size(A,1)
L = diagm(0=>ones(n))  # ones on main diagonal, zeros elsewhere
U = float(copy(A))

# Gaussian elimination
for j = 1:n-1
  for i = j+1:n
    L[i,j] = U[i,j] / U[j,j]   # row multiplier
    U[i,j:n] -= L[i,j]*U[j,j:n]
  end
end

return L,triu(U)
end

max_stars_repo_name: tobydriscoll/fnc

max_issues_repo_name: tobydriscoll/fnc

max_forks_repo_name: tobydriscoll/fnc


text:
export segy_write_append

function segy_write_append(file::String, block::SeisBlock)

    # Open buffer for writing
    s = open(file, "r+")

    # Write FileHeader
    check_fileheader(s, block.fileheader)

    # Write Data
    ns,ntraces = size(block.data)
    for t in 1:ntraces
	write_trace(seekend(s), block, t)
    end

    close(s)
end


text:
using BenchmarkTools, ModelingToolkit

include(joinpath(pwd(), "src/models/car.jl"))

x = zeros(model.n)
u = zeros(model.m)
w = zeros(model.d)
h = 0.0
t = 0
q = zeros(model.nq)

@benchmark M_func($model, $q)
@benchmark C_func($model, $q, $q)
@benchmark fd($model, $x, $x, $u, $w, $h, $t)

n = 15
A = Diagonal(@SVector ones(n))

L(q) = transpose(q) * A * q

f(q) = ForwardDiff.gradient(L, q)
g(q) = 2.0 * A * q
x = rand(n)
f(x) == g(x)
@benchmark f($x)
@benchmark g($x)

# generate fast functions
@variables x_sym[1:n]

F = g(x_sym);
fast = eval(ModelingToolkit.build_function(F, x_sym,
            parallel=ModelingToolkit.MultithreadedForm())[1])
fast_sparsity = ModelingToolkit.sparsejacobian(F, x_sym)
∇fast! = eval(ModelingToolkit.build_function(fast_sparsity, x_sym,
            parallel=ModelingToolkit.MultithreadedForm())[2])
∇fast = similar(fast_sparsity, Float64)

y = ones(n)
@benchmark fast!($y, $x)
@benchmark fast($x)

l = L(x_sym);
fast = eval(ModelingToolkit.build_function([l], x_sym,
            parallel=ModelingToolkit.MultithreadedForm())[1])
fast_sparsity = ModelingToolkit.sparsejacobian([l], x_sym)
dfast = eval(ModelingToolkit.build_function(fast_sparsity, x_sym,
            parallel=ModelingToolkit.MultithreadedForm())[1])

l2 = dfast(x_sym)

fast2 = eval(ModelingToolkit.build_function(simplify.(l2), x_sym,
            parallel=ModelingToolkit.MultithreadedForm())[1])
fast2_sparsity = ModelingToolkit.sparsejacobian(simplify.(l2), x_sym)
∇fast2! = eval(ModelingToolkit.build_function(fast2_sparsity, x_sym,
            parallel=ModelingToolkit.MultithreadedForm())[2])
∇fast2 = similar(fast2_sparsity, Float64)

# pendulum

L(q, q̇) = 0.5 * 1.0 * q̇^2.0 - 1.0 * 9.81 * 1.0 * cos(q)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
using ArgParse
using JSON
using DataStructures
using NCDatasets

function runOneCmd(cmd)
    println(">> ", string(cmd))
    run(cmd)
end

function pleaseRun(cmd)
    if isa(cmd, Array)
        for i = 1:length(cmd)
            runOneCmd(cmd[i])
        end
    else
        runOneCmd(cmd)
    end
end

println("""
This file produce an netcdf file with only `time` and `time_bound` variables
for users to append variables to other files.

Notice: This program uses no-leap calendar (365 days fixed)

""")

s = ArgParseSettings()
@add_arg_table s begin

    "--output"
        help = "The file you want to produce."
        arg_type = String
        required = true

    "--years"
        help = "How many years?"
        arg_type = Int64
        default = 1

end

parsed = DataStructures.OrderedDict(parse_args(ARGS, s))

JSON.print(parsed, 4)

doy = 365

t    = zeros(Float64,    doy)
bnds = zeros(Float64, 2, doy)

for d=1:doy*parsed["years"]
    bnds[1, d] = d-1
    bnds[2, d] = d
    t[d] = (bnds[1, d] + bnds[2, d]) / 2.0
end

Dataset(parsed["output"], "c") do ds

    defDim(ds, "time", Inf)
    defDim(ds, "d2", 2)

    defVar(ds, "time", t, ("time", ), ; attrib = Dict(
        "long_name" => "time",
        "bounds"    => "time_bound",
        "calendar"  =>  "noleap",
        "units"     => "days since 0001-01-01 00:00:00",
    ))

    defVar(ds, "time_bound", bnds, ("d2", "time"), ; attrib = Dict(
        "long_name" => "boundaries for time-averaging interval",
        "units"     => "days since 0001-01-01 00:00:00",
    ))

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# Write to text / REPL
testshow(A::TimeAxisArray) = show(IOBuffer(),MIME("text/plain"),A)
map(testshow, TimeAxisArray[A1d1, A1d2, A2d1, A3d, A4d])

# Read to / write from file

testfile = "test.csv"

writedlm(testfile, A1d1)
@test A1d1 == readtimeaxisarray(testfile, datetime(), header=false)

writedlm(testfile, A2d1)
@test A2d1 == readtimeaxisarray(testfile, datetime(), namedaxes=true)

writedlm(testfile, A4d)
@test A4d == readtimeaxisarray(testfile, datetime(), namedaxes=true, headlines=3)

rm(testfile)

max_stars_repo_name: GordStephen/TimeAxisArrays

max_issues_repo_name: GordStephen/TimeAxisArrays

max_forks_repo_name: GordStephen/TimeAxisArrays


text:
module DCDdouble

using Parameters

export load_dcd, natoms, nframes, timestep, hascell, cell, positions, timepassed

"""
    load_dcd(f::String)

Load a dcd trajectory from file `f`.
"""
function load_dcd(f::String)
    io = open(f)
    load_dcd(io)
end

"""
    load_dcd(io::IO)

Load a dcd trajectory from stream `io`.
"""
function load_dcd(io::IO)
    checkmagic(io)
    File(io)
end

"""
    checkmagic(io::IO)

Check whether the stream `io` is a valid dcd file by looking at its magic bytes.
"""
function checkmagic(io::IO)
    magic = read(io, 8)
    magic == UInt8[0x54, 0x00, 0x00, 0x00, 0x43, 0x4f, 0x52, 0x44] || throw("invalid magic bytes")
    seekstart(io)
end

"""
    @with_kw struct File

Stores information about a dcd trajectory file.
"""
@with_kw struct File
    io::IO
    natoms::Int64
    nframes::Int64
    δt::Float64
    hascell::Bool
    nbytes_header::Int64
    nbytes_frame::Int64
end

"""
    File(io::IO)

Construct a [`DCD.File`](@ref) object from an IO stream `io`.
"""
function File(io::IO)
    header = Header(io)
    File(
        io=io,
        natoms=natoms(header),
        nframes=nframes(header),
        δt=timestep(header),
        hascell=hascell(header),
        nbytes_header=nbytes_header(header),
        nbytes_frame=nbytes_frame(header)
        )
end

"""
    iostream(f::File)

Get the [`IOStream`](@ref) associated with the [`DCD.File`](@ref) object.
"""
iostream(f::File) = f.io

"""
    natoms(f::File)

Get the number of atoms.
"""
natoms(f::File) = f.natoms

"""
    nframes(f::File)

Get the number of frames.
"""
nframes(f::File) = f.nframes

"""
    timestep(f::File)

Get the time passed between frames.
"""
timestep(f::File) = f.δt

"""
    hascell(f::File)

Check if the dcd trajectory file has unit cell information.
"""
hascell(f::File) = f.hascell

"""
    nbytes_header(f::File)

Get the number of bytes of the header.
"""
nbytes_header(f::File) = f.nbytes_header

"""
    nbytes_frame(f::File)

Get the number of bytes of a single frame.
"""
nbytes_frame(f::File) = f.nbytes_frame

"""
    Base.eltype(::Type{File})

[`DCD.File`](@ref) objects are iterators over [`DCD.Frame`](@ref) objects.
"""
Base.eltype(::Type{File}) = Frame

function Base.getindex(f::File, i::Int64)
    1 <= i <= f.nframes || throw(BoundsError(f, i))
    return Frame(f, i)
end
Base.getindex(f::File, i::Number) = f[convert(Int64, i)]
Base.getindex(f::File, I) = [f[i] for i in I]
Base.firstindex(f::File) = 1
Base.lastindex(f::File) = f.nframes
Base.length(f::File) = f.nframes

function Base.iterate(f::File, frame::Int64=1)
    if frame > nframes(f)
        return nothing
    else
        return Frame(f, frame), frame + 1
    end
end

"""
    @with_kw struct Header

Stores information from the dcd trajectory file header.
"""
@with_kw struct Header
    nframes::Int32
    firststep::Int32
    stepincrement::Int32
    laststep::Int32
    δt::Float32
    hascell::Int32
    natoms::Int32
    nbytes::Int64
end

"""
    header(io::IO)

Construct a [`DCD.Header`](@ref) object from an [`IOStream`](@ref).
"""
function Header(io::IO)
    skip(io, 8)
    nframes = read(io, Int32)
    firststep = read(io, Int32)
    stepincrement = read(io, Int32)
    laststep = read(io, Int32)
    skip(io, 20)
    δt = read(io, Float32)
    hascell = read(io, Int32)
    skip(io, 216)
    natoms = read(io, Int32)
    skip(io, 4)
    nbytes = position(io)

    Header(
        nframes=nframes,
        firststep=firststep,
        stepincrement=stepincrement,
        laststep=laststep,
        δt=δt,
        hascell=hascell,
        natoms=natoms,
        nbytes=nbytes
        )
end

"""
    natoms(h::Header)

Get the number of atoms.
"""
natoms(h::Header)::Int64 = h.natoms

"""
    nframes(h::Header)

Get the number of frames.
"""
nframes(h::Header)::Int64 = h.nframes

"""
    timestep(h::Header)

Get the time passed between two subsequent frames.
"""
timestep(h::Header)::Float64 = h.δt * h.stepincrement

"""
    getcell(h::Header)

Check whether the dcd trajectory file has unit cell information.
"""
hascell(h::Header)::Bool = h.hascell

"""
    nbytes_header(h::Header)

Get the number of bytes of the header.
"""
nbytes_header(h::Header)::Int64 = h.nbytes

"""
    nbytes_frame(h::Header)

Get the number of bytes of a frame.
"""
nbytes_frame(h::Header)::Int64 = hascell(h) ? natoms(h) * 12 + 24 + 56 : natoms(h) * 12 + 24
# Syntax note: [Bool expression ? Action1 if True : Action2 if False]

"""
    seekframe(f::File, index::Int64)

Move the file's [`IOStream`](@ref) to the position of the indicated frame.
"""
function seekframe(f::File, index::Int64)
    1 <= index <= f.nframes || throw(BoundsError(f, index))
    pos = nbytes_header(f) + (index - 1) * nbytes_frame(f)
    seek(f.io, pos)
end

"""
    @with_kw struct Frame

Stores simulation time, unit cell information and positions.
"""
@with_kw struct Frame
    time::Float64
    cell::Union{Vector{Float64},Nothing}
    positions::Array{Float64,2}
end

"""
    Frame(f::File, index::Int64)

Construct a [`DCD.Frame`](@ref) object from a [`DCD.File`](@ref) `f` with the specified `index`.
"""
function Frame(f::File, index::Int64)
    seekframe(f, index)
    time = index * timestep(f)
    io = iostream(f)
    if hascell(f)
        skip(io, 4)
        cell = Vector{Float64}(undef, 6)
        read!(io, cell)
        skip(io, 4)
    else
        cell = nothing
    end
    positions = Array{Float64,2}(undef, 3, natoms(f))
    for dim in 1:3
        skip(io, 4)
        for iatom in 1:natoms(f)
            positions[dim, iatom] = read(io, Float32)
        end
        skip(io, 4)
    end

    Frame(time=time, cell=cell, positions=positions)
end

"""
    timepassed(f::Frame)

Get the simulation time that has passed since the beginning of the simulation.
"""
timepassed(f::Frame) = f.time

"""
    positions(f::Frame)

Get the current positions.
"""
positions(f::Frame) = f.positions

function cell(f::Frame)
    cell = zeros(6)
    cell[1] = f.cell[1]
    cell[2] = f.cell[3]
    cell[3] = f.cell[6]
    cell[4] = 90.0 - asind(f.cell[5])
    cell[5]  = 90.0 - asind(f.cell[4])
    cell[6] = 90.0 - asind(f.cell[2])
    cell
end

end # module


text:
using Test

using GRIB

@testset "from_grib_date_time" begin
    result = CfGRIB.from_grib_date_time(20160706, 1944)

    @test result == 1467834240
end

@testset "test_build_valid_time" begin
    forecast_reference_time = [0]
    forecast_period = [0]

    vt_dims, vt_data = CfGRIB.build_valid_time(
        forecast_reference_time,
        forecast_period
    )

    @test vt_dims == ()
    @test size(vt_data) == ()

    forecast_reference_time = [0, 31536000]
    forecast_period = 0

    vt_dims, vt_data = CfGRIB.build_valid_time(
        forecast_reference_time,
        forecast_period
    )

    @test vt_dims == ("time", )
    @test size(vt_data) == (2, )

    forecast_reference_time = 0
    forecast_period = [0, 12, 24, 36]

    vt_dims, vt_data = CfGRIB.build_valid_time(
        forecast_reference_time,
        forecast_period
    )

    @test vt_dims == ("step", )
    @test size(vt_data) == (4,)
    # @test np.allclose((data - data[..., :1]) / 3600, forecast_period)

    forecast_reference_time = [0, 31536000]
    forecast_period = [0, 12, 24, 36]

    vt_dims, vt_data = CfGRIB.build_valid_time(
        forecast_reference_time,
        forecast_period
    )

    @test vt_dims == ("time", "step")

    @test size(vt_data) == (4, 2)
    # @test np.allclose((data - data[..., :1]) / 3600, forecast_period)
end


text:
"""
    abstract type CPLayerVertex end

A vertex that contains the CP "true" object.
"""
abstract type CPLayerVertex end

"""
    abstract type FixedEdgesVertex <: CPLayerVertex end

A vertex whose edges won't be modified during the solving.
"""
abstract type FixedEdgesVertex <: CPLayerVertex end

"""
    struct ValueVertex <: CPLayerVertex

A vertex corresponding to the value of a domain, that is connected to
a variable only if the value is in the domain of that variable.
"""
struct ValueVertex <: CPLayerVertex
    value       ::Int
end

"""
    struct ConstraintVertex <: CPLayerVertex

A vertex corresponding to a constraint, that is connected to
a variable only if the constraint affects that variable.
"""
struct ConstraintVertex <: FixedEdgesVertex
    constraint  ::Constraint
end

"""
    struct VariableVertex <: CPLayerVertex

A vertex corresponding to a variable, that is connected to
a constraint only if that constraint affects the variable,
and connected to a value only if that value is in the domain
of the variable.
"""
struct VariableVertex <: FixedEdgesVertex
    variable    ::AbstractVar
end

"""
    struct CPLayerGraph <: AbstractGraph{Int}

Graph representing the current status of the CPModel.
It is a tripartite graph, linking 3 types of nodes: constraints, variables and values.
A constraint is connected to a variable if the constraint affects the variable.
A variable is connected to a value if the value is in the variable's domain.

Since the relations between constraints and variables are not modified during the solving,
they are stored in an inner graph, `fixedEdgesGraph`, that won't be edited after creation.
On the contrary, domains get pruned during the solving. To always keep an up-to-date representation,
none of the value-variables edges are stored, and everytime they're needed, they are gotten from
the `CPModel` directly, hence the "layer" in the name of the struct.

Each vertex is indexed by an integer, in a specific order. You first have the constraints,
then the variables and finally the values.
"""
struct CPLayerGraph <: LightGraphs.AbstractGraph{Int}
    cpmodel                     ::Union{CPModel, Nothing}
    idToNode                    ::Array{CPLayerVertex}
    nodeToId                    ::Dict{CPLayerVertex, Int}
    fixedEdgesGraph             ::LightGraphs.Graph
    numberOfConstraints         ::Int
    numberOfVariables           ::Int
    numberOfValues              ::Int
    totalLength                 ::Int

    """
        function CPLayerGraph(cpmodel::CPModel)

    Create the graph corresponding to the CPModel.
    The graph gets linked to `cpmodel` and does not need to get updated by anyone when domains are pruned.
    """
    function CPLayerGraph(cpmodel::CPModel)
        variables = Set{AbstractVar}(values(cpmodel.variables))
        valuesOfVariables = branchable_values(cpmodel)
        numberOfConstraints = length(cpmodel.constraints)
        numberOfValues = length(valuesOfVariables)

        variableConnections = Tuple{AbstractVar, AbstractVar}[]

        # Take into account IntVarViews that are only declared in constraints
        for constraint in cpmodel.constraints
            for constraintVar in variablesArray(constraint)
                while typeof(constraintVar) <: Union{IntVarView, BoolVarView}
                    push!(variables, constraintVar)

                    # Storing variable connections
                    push!(variableConnections, (constraintVar, constraintVar.x))

                    constraintVar = constraintVar.x
                end
            end
        end

        variables = collect((variables))

        # We sort the variables by their id to get a consistent order
        sort!(variables; by=(x) -> x.id)

        numberOfVariables = length(variables)

        totalLength = numberOfConstraints + numberOfVariables + numberOfValues
        nodeToId = Dict{CPLayerVertex, Int}()
        idToNode = Array{CPLayerVertex}(undef, totalLength)

        # Filling constraints
        for i in 1:numberOfConstraints
            idToNode[i] = ConstraintVertex(cpmodel.constraints[i])
            nodeToId[ConstraintVertex(cpmodel.constraints[i])] = i
        end

        # Filling variables
        for i in 1:numberOfVariables
            idToNode[numberOfConstraints + i] = VariableVertex(variables[i])
            nodeToId[VariableVertex(variables[i])] = numberOfConstraints + i
        end

        # Filling values
        for i in 1:numberOfValues
            idToNode[numberOfConstraints + numberOfVariables + i] = ValueVertex(valuesOfVariables[i])
            nodeToId[ValueVertex(valuesOfVariables[i])] = numberOfConstraints + numberOfVariables + i
        end

        fixedEdgesGraph = Graph(numberOfConstraints + numberOfVariables)
        for id in 1:numberOfConstraints
            constraint = idToNode[id].constraint
            varArray = variablesArray(constraint)
            for x in varArray
                if x.id != "SEAPEARL_one"
                    add_edge!(fixedEdgesGraph, id, nodeToId[VariableVertex(x)])
                end
            end
        end

        for (x1, x2) in variableConnections
            v1, v2 = VariableVertex(x1), VariableVertex(x2)
            add_edge!(fixedEdgesGraph, nodeToId[v1], nodeToId[v2])
        end

        return new(cpmodel, idToNode, nodeToId, fixedEdgesGraph, numberOfConstraints, numberOfVariables, numberOfValues, totalLength)
    end

    """
        function CPLayerGraph()

    Create an empty graph, needed to implement the `zero` function for the LightGraphs.jl interface.
    """
    function CPLayerGraph()
        return new(nothing, CPLayerVertex[], Dict{CPLayerVertex, Int}(), Graph(0), 0, 0, 0, 0)
    end
end


text:
#=
Created on Thurday 20 Feb 2010
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Plot the results of the honeybee spillover.
=#

using StatsPlots, HypothesisTests, BSON

BSON.@load "results/spillover.bson" Ps_during Ps_after M M_a M_b M_ab kl_ot_ab  kl_ot_a kl_ot_b kl_ot_softmax kl_ot_neutral

# tests
# -----

# test if kl_ot_ab < kl_ot_neutral
@show wsr = SignedRankTest(kl_ot_neutral, kl_ot_ab)
@show pvalue(wsr, tail=:right)

# plotting
# --------

p_mut = boxplot(kl_ot_neutral, label="neutral", legend=:none)
boxplot!(kl_ot_ab, label="OT (both)")
boxplot!(kl_ot_a, label="OT (plants fixed)")
boxplot!(kl_ot_b, label="OT (poll. fixed)")
boxplot!(kl_ot_softmax, label="OT (none)")
ylabel!("Kullback-Leibler (nats)")
xticks!(1:5, ["neutral", "OT (both)", "OT (plants fixed)", "OT (poll. fixed)", "OT (none)"])
title!("Predicted versus observed pollination\nduring honeybee spillover")

savefig(p_mut, "figures/honeybee_spillover.svg")

# Same figure, but on training
using OTSIN

kl_ot_ab = Float64[]
kl_ot_a = Float64[]
kl_ot_b = Float64[]
kl_ot_softmax = Float64[]
kl_ot_neutral = Float64[]

for P in Ps_after
    a, b = marginals(P)
    P = P[a.>0, b.>0]
    M = M_ab[a.>0, b.>0]
    a, b = marginals(P)
    # optimal transport both
    Q = optimaltransport(M, a, b)
    push!(kl_ot_ab, KL(P, Q))
    Q .= optimaltransport(M, a, nothing)
    push!(kl_ot_a, KL(P, Q))
    Q .= optimaltransport(M, nothing, b)
    push!(kl_ot_b, KL(P, Q))
    Q .= optimaltransport(M, nothing, nothing)
    push!(kl_ot_softmax, KL(P, Q))
    Q .= optimaltransport(a, b)
    push!(kl_ot_neutral, KL(P, Q))
end

p_mut = boxplot(kl_ot_neutral, label="neutral", legend=:none)
boxplot!(kl_ot_ab, label="OT (both)")
boxplot!(kl_ot_a, label="OT (plants fixed)")
boxplot!(kl_ot_b, label="OT (poll. fixed)")
boxplot!(kl_ot_softmax, label="OT (none)")
ylabel!("Kullback-Leibler (nats)")
xticks!(1:5, ["neutral", "OT (both)", "OT (plants fixed)", "OT (poll. fixed)", "OT (none)"])
title!("Predicted versus observed pollination\nafter honeybee spillover")

savefig(p_mut, "figures/honeybee_spillover_train.svg")

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN