text:
export AbstractCacheServer

"""
    AbstractCacheServer{K, ELT}
"""
abstract type AbstractCacheServer{K, ELT} end

"""
    alloc!(server, object, storage) -> server

alloc new storage on the server.
"""
function alloc!(s::AbstractCacheServer, object, storage)
    throw(MethodError(alloc!, (s, object, storage)))
end

"""
    iscacheable(server, object)

check if there is available space to storage this object's value. (if this object
was allocated on the server before.).
"""
function iscacheable(s::AbstractCacheServer, object)
    throw(MethodError(iscacheable, (s, object)))
end

"""
    iscached(server, object, [params...])

check if this object (with params) is already cached.
"""
function iscached(s::AbstractCacheServer, object, params...)
    throw(MethodError(push!, (s, object, params...)))
end

"""
    push!(server, val, object) -> server

push `val` to the storage of `object` in the `server`.
"""
function push!(s::AbstractCacheServer, val, object)
    throw(MethodError(push!, (s, val, object)))
end

"""
    update!(storage, val) -> storage
"""
function update!(storage, val)
    throw(MethodError(update!, (storage, val)))
end

"""
    pull(server, object, params...) -> value

pull `object` storage from server.
"""
function pull(server::AbstractCacheServer, object, params...)
    throw(MethodError(pull, (server, object, params...)))
end

# getindex will call pull by default
function getindex(server::AbstractCacheServer, index)
    pull(server, index)
end

"""
    delete!(server, object) -> server

delete this object from the server. (the storage will be deleted)
"""
function delete!(server::AbstractCacheServer, object)
    throw(MethodError(clear!, (server, object)))
end

"""
    clear!(server, object) -> server

clear the storage in the `server` of this `object`.
"""
function clear!(server::AbstractCacheServer, object)
    throw(MethodError(clear!, (server, object)))
end


text:
using Documenter, BibParser, BibInternal

makedocs(
    sitename = "BibParser.jl",
    authors = "Jean-François BAFFIER",
    repo="https://github.com/Humans-of-Julia/BibParser.jl/blob/{commit}{path}#L{line}",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages = [
        "BibParser" => "index.md",
        "BibTeX" => "bibtex.md",
        # "BibTeX - automa" => "bibtex_automa.md",
        "CSL-JSON" => "csl.md",
        "Internal" => "internal.md",
    ]
)

deploydocs(; repo = "github.com/Humans-of-Julia/BibParser.jl.git", devbranch = "master")

ext: jl


text:
using MPI
MPI.Init()

comm = MPI.COMM_WORLD
root = 0

r = MPI.Comm_rank(comm)

sr = MPI.Reduce(r, +, root, comm)

if MPI.Comm_rank(comm) == root
    println("sum of ranks = $sr")
end


text:
#
# SPMe with input parameters, comparing with PyBaMM
#

using PyBaMM
using SparseArrays, LinearAlgebra
using PyCall

# load model
pybamm = pyimport("pybamm")
model = pybamm.lithium_ion.SPMe(name="SPMe")
parameter_values = model.default_parameter_values
parameter_values.update(
    Dict(
        "Cation transference number" => pybamm.InputParameter("t_plus"),
        "Electrolyte conductivity [S.m-1]" => pybamm.InputParameter("kappa_e"),
    )
)
sim = pybamm.Simulation(model, parameter_values=parameter_values)

using OrderedCollections
inputs = OrderedDict{String,Float64}(
    "t_plus" => 0.4,
    "kappa_e" => 0.7,
)

# Build ODE problem
prob = get_ode_problem(sim, inputs)

# Generate data
using Sundials
sol = solve(prob, CVODE_BDF(linear_solver=:Band,jac_upper=1,jac_lower=1), reltol=1e-6, abstol=1e-6, saveat=prob.tspan[2] / 100);
V_data = get_variable(sim, sol, "Terminal voltage [V]", inputs) + 0.005*randn(size(sol.t))
t = get_variable(sim, sol, "Time [s]")

using Plots
plot(t,V_data)

# Build loss objective
using DiffEqParamEstim

loss = get_l2loss_function(sim, "Terminal voltage [V]", inputs, V_data)
cost_function = build_loss_objective(prob,CVODE_BDF(linear_solver=:Band,jac_upper=1,jac_lower=1),loss,
                                     maxiters=10000,verbose=true,saveat=sol.t)

using BenchmarkTools
@btime loss(sol)
@btime get_variable(sim, sol, "Terminal voltage [V]", inputs)
# Sanity check: plot 1-parameter cost function
vals = 0.1:0.01:0.8
cost = [cost_function([0.4,i]) for i in vals]
using Plots; plotly()
plot(vals,cost,yscale=:log10,
    xaxis = "Parameter", yaxis = "Cost", title = "1-Parameter Cost Function",
    lw = 3)

# Optimize
using Optim
result = optimize(cost_function, [0.3, 0.6], BFGS())
@show result.minimizer


text:
# Trajectory Set-Up

"""
    initializeTraj(x0::Array{Float64, 1}, xN::Array{Float64, 1}, NSteps::Int64)

Initializes the full `XU` vector (see example below) using a linear
interpolation from the initial point `x0` to the final point `xN`.

Take `x0 = [s0; v0]` and `xN = [sN; vN`] to be of size (`2n×1`). The full `XU`
vector size is `(3Nn + 2n)×1`.

For example, take a 2D system: `x0 = [sx0; sy0; vx0; vy0]` and
                               `xN = [sxN; syN; vxN, vyN]`.

The controls are `uk = [uxk; uyk]`. Thus, for `N = 4`

`XU` vector = `[x0; u0; x1; u1; x2; u2; x3; u3; x4]`

which is of size `((4 + 2) + (4 + 2) + (4 + 2) + (4 + 2) + (4)) = 28`
or equivalently `3Nn + 2n = 3(4)(2) + 2(2) = 24 + 4 = 28`, as desired.
"""
function initializeTraj(x0::Array{Float64, 1}, xN::Array{Float64, 1},
                        NSteps::Int64)
    if size(x0) != size(xN)
        error("x0 and xN must have the same size. " *
                "Currently: x0 = $x0 and xN = $xN")
    end
    if mod(size(x0, 1), 2) != 0
        error("x0 and xN must have an even number of variables. " *
              "Example, x0 = [sx0; sy0; vx0; vy0]. Currently x0 = $x0.")
    end

    # Calculate the linearly interpolated x values
    xSteps = LinRange(x0, xN, NSteps + 1)
    # Now we want to stack this as a vector in the form
    # [x0; u0; x1; ...; uN-1, xN]
    nDim = Int64(size(x0, 1) / 2)
    sizeXU = 3 * NSteps * nDim + 2 * nDim
    XUfull = zeros(sizeXU, 1)

    # println("Dimensions = $nDim")

    for k in 1:(NSteps + 1)
        iStart = 1 + (k - 1) * (3 * nDim)
        iEnd = iStart + (2 * nDim - 1)
        # println("Accessing ($iStart, $iEnd)")
        XUfull[iStart:iEnd] = xSteps[k]
    end

    return XUfull
end

function initializeTraj(x0::Array{Int64, 1}, xN::Array{Int64, 1},
                        NSteps::Int64)
    return initializeTraj(float.(x0), float.(xN), NSteps)
end

"""
    initializeTraj(x0::Array{Float64, 1}, xN::Array{Float64, 1},
               u0::Array{Float64, 1}, uN::Array{Float64, 1},
               NSteps::Int64)

Initializes the full `XU` vector (see example below) using a linear
interpolation from the initial point `x0` to the final point `xN`.

Take `x0 = [s0; v0]` and `xN = [sN; vN`] to be of size (`2n×1`). The full `XU`
vector size is `(3Nn + 2n)×1`.

For example, take a 2D system: `x0 = [sx0; sy0; vx0; vy0]` and
                               `xN = [sxN; syN; vxN, vyN]`.

The controls are `uk = [uxk; uyk]`. Thus, for `N = 4`

`XU` vector = `[x0; u0; x1; u1; x2; u2; x3; u3; x4]`

which is of size `((4 + 2) + (4 + 2) + (4 + 2) + (4 + 2) + (4)) = 28`
or equivalently `3Nn + 2n = 3(4)(2) + 2(2) = 24 + 4 = 28`, as desired.
"""
function initializeTraj(x0::Array{Float64, 1}, xN::Array{Float64, 1},
                        u0::Array{Float64, 1}, uN::Array{Float64, 1},
                        NSteps::Int64)
    if size(x0) != size(xN)
        error("x0 and xN must have the same size. " *
                "Currently: x0 = $x0 and xN = $xN")
    end
    if mod(size(x0, 1), 2) != 0
        error("x0 and xN must have an even number of variables. " *
              "Example, x0 = [sx0; sy0; vx0; vy0]. Currently x0 = $x0.")
    end

    # Calculate the linearly interpolated x values
    xSteps = LinRange(x0, xN, NSteps + 1)
    uSteps = LinRange(u0, uN, NSteps)
    # Now we want to stack this as a vector in the form
    # [x0; u0; x1; ...; uN-1, xN]
    nDim = Int64(size(x0, 1) / 2)
    sizeXU = 3 * NSteps * nDim + 2 * nDim
    XUfull = zeros(sizeXU, 1)

    # println("Dimensions = $nDim")

    for k in 1:(NSteps + 1)
        iStart = 1 + (k - 1) * (3 * nDim)
        iEnd = iStart + (2 * nDim - 1)
        # println("Accessing ($iStart, $iEnd)")
        XUfull[iStart:iEnd] = xSteps[k]

        if k ≤ NSteps
            jStart = iEnd + 1
            jEnd = jStart + (nDim - 1)
            XUfull[jStart:jEnd] = uSteps[k]
        end
    end

    return XUfull
end

"""
    initializeHoverOnlyTraj(x0::Array{Float64, 1}, xN::Array{Float64, 1},
                        uHover::Array{Float64, 1}, NSteps::Int64)

Initializes the full `XU` vector (see example below) that hovers at the initial
point and does not move. The final state is effectively appended at the end.

The control vector is set to hover only.
"""
function initializeHoverOnlyTraj(x0::Array{Float64, 1}, xN::Array{Float64, 1},
                                 uHover::Array{Float64, 1}, NSteps::Int64)
     if size(x0) != size(xN)
         error("x0 and xN must have the same size. " *
                 "Currently: x0 = $x0 and xN = $xN")
     end
     if mod(size(x0, 1), 2) != 0
         error("x0 and xN must have an even number of variables. " *
               "Example, x0 = [sx0; sy0; vx0; vy0]. Currently x0 = $x0.")
     end

     nDim = Int64(size(x0, 1) / 2)
     sizeXU = 3 * NSteps * nDim + 2 * nDim
     XUfull = zeros(sizeXU, 1)

     for k in 1:(NSteps + 1)
         iStart = 1 + (k - 1) * (3 * nDim)
         iEnd = iStart + (2 * nDim - 1)
         # println("Accessing ($iStart, $iEnd)")

         if k ≤ NSteps
             XUfull[iStart:iEnd] = x0
         else
             XUfull[iStart:iEnd] = xN
         end

         if k ≤ NSteps
             jStart = iEnd + 1
             jEnd = jStart + (nDim - 1)
             XUfull[jStart:jEnd] = uHover
         end
     end

     return XUfull
end

"""
    initializeHoverOnlyTraj(x0::Array{Float64, 1}, uHover::Array{Float64, 1},
                            NSteps::Int64)

Initializes the full `XU` vector (see example below) that hovers at the initial
point and does not move. The control vector is set to hover only.
"""
function initializeHoverOnlyTraj(x0::Array{Float64, 1},
                                 uHover::Array{Float64, 1}, NSteps::Int64)
     if mod(size(x0, 1), 2) != 0
         error("x0 and xN must have an even number of variables. " *
               "Example, x0 = [sx0; sy0; vx0; vy0]. Currently x0 = $x0.")
     end

     nDim = Int64(size(x0, 1) / 2)
     sizeXU = 3 * NSteps * nDim + 2 * nDim
     XUfull = zeros(sizeXU, 1)

     for k in 1:(NSteps + 1)
         iStart = 1 + (k - 1) * (3 * nDim)
         iEnd = iStart + (2 * nDim - 1)
         # println("Accessing ($iStart, $iEnd)")
         XUfull[iStart:iEnd] = x0

         if k ≤ NSteps
             jStart = iEnd + 1
             jEnd = jStart + (nDim - 1)
             XUfull[jStart:jEnd] = uHover
         end
     end

     return XUfull
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
# QCQP example
#    maximize x + y
#
#    s.t.  x, y >= 0
#          x^2 + y^2 <= 1
#
#    solution: (0.71, 0.71) objv = 1.414
using Xpress, Compat.Test
using Compat.SparseArrays
using Compat.LinearAlgebra

@testset "QCQP 1" begin
    model = Xpress.Model( "qcqp_01", :maximize)

    add_cvars!(model, [1., 1.], 0., Inf)

    # add_qpterms!(model, linearindices, linearcoeffs, qrowinds, qcolinds, qcoeffs, sense, rhs)
    add_qconstr!(model, [], [], [1, 2], [1, 2], [1, 1.], '<', 1.0)

    @test Xpress.get_obj(model) == [1, 1]
    @test Xpress.get_constrmatrix(model) == spzeros(1,2)
    @test Xpress.get_rhs(model) == [1]
    @test Xpress.get_lb(model) == [0, 0]
    @test Xpress.get_ub(model) == [1e20, 1e20]
    @test Xpress.get_sense(model) == [:<=]

    optimize(model)

    ans = get_optiminfo(model)
    @test ans.status_lp == :optimal

    @test isapprox(get_solution(model), [0.7071, 0.7071]; atol = 1e-3)

    @test isapprox(get_objval(model), 1.414; atol = 1e-3)
end


text:
import Composites

# Step 1: Identify Material Properties
e1 = [181e9,203e9,38.6e9,76e9]
e2 = [10.3e9,11.2e9,8.27e9,5.5e9]
g12 = [7.17e9,8.4e9,4.14e9,2.3e9]
nu12 = [0.28,0.32,0.26,0.34]
rho = zeros(Float64,4)
xt = [1500.0,3500.0,1062.0,1400.0]*1e6
xc = [1500.0,1540.0,610.0,235.0]*1e6
yt = [40.0,56.0,31.0,12.0]*1e6
yc = [246.0,150.0,118.0,53.0]*1e6
s = [68.0,98.0,72.0,34.0]*1e6
t = fill(0.0001,4)
mat = Composites.Material.(e1,e2,g12,nu12,rho,xt,xc,yt,yc,s,t)

# Step 2: Input Ply Information
matid = [4,4,4,4,4,4,4]
nply = [5,5,5,10,5,5,5]
tply = [0.0001,0.0001,0.0001,0.0001,0.0001,0.0001,0.0001]
theta = [0.0,45.0,-45.0,90.0,-45.0,45.0,0.0]
lam = Composites.Laminate(matid,nply,tply,theta)

# Step 3: Get Material Stiffness Matrix (Q)
q = Composites.getQ.(mat)

# Step 4: Determine Laminate Stiffness Matrix (A,B, and D)
A,B,D = Composites.getABD(lam,q)

# Step 5: Get Resultant Forces and Moments
force = [25000.0,0.0,0.0]
moment = [1000.0,0.0,0.0]
resultantstress = vcat(force,moment)

# Step 6: Get Resultant Strain
S = vcat(hcat(A,B),hcat(B',D))
resultantstrain = S\resultantstress

# Step 7: Solve for local ply strains and stresses in material coordinate system
lowerplystrain,upperplystrain = Composites.getplystrain(lam,resultantstrain)
lowerplystress = Composites.getplystress(lowerplystrain,q,lam)
upperplystress = Composites.getplystress(upperplystrain,q,lam)

# Step 8: Apply appropriate material failure criteria
lowermaxstress,lsfmaxstress = Composites.getmatfail(lowerplystress,mat,lam,"maxstress")
uppermaxstress,usfmaxstress = Composites.getmatfail(upperplystress,mat,lam,"maxstress")
lowertsai,lsftsai = Composites.getmatfail(lowerplystress,mat,lam,"tsaiwu")
uppertsai,usftsai = Composites.getmatfail(upperplystress,mat,lam,"tsaiwu")
lowerhashinrotem,lsfhashinrotem = Composites.getmatfail(lowerplystress,mat,lam,"hashinrotem")
upperhashinrotem,usfhashinrotem = Composites.getmatfail(upperplystress,mat,lam,"hashinrotem")


text:
# push!(LOAD_PATH,"../src/")
using Documenter, GPLinearODEMaker

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

makedocs(
    modules = [GPLinearODEMaker],
    sitename = "GPLinearODEMaker.jl",
    authors = "Christian Gilbertson",
    format = Documenter.HTML(),
    pages = [
        "Home" => "index.md",
        "gettingstarted.md",
        "User's Guide" => [
            "Kernels" => [
                "Kernel functions" => "kernel.md",
                "Adding new kernels" => "kernel_creation.md",
                ],
            "GPLinearODE struct" => "glo.md",
            "GLOM Functionality" => "nlogl.md",
            "Prior functions" => "priors.md",
        ],
        "LICENSE.md",
        hide("Indices" => "indices.md"),
        hide("Diagnostic functions" => "diagnostic.md"),
        hide("Utility functions" => "utility.md"),
    ]
)

deploydocs(
    repo = "github.com/christiangil/GPLinearODEMaker.jl.git",
    deploy_config = Documenter.GitHubActions(),
)

max_stars_repo_name: christiangil/GPLinearODEMaker

max_issues_repo_name: christiangil/GPLinearODEMaker

max_forks_repo_name: christiangil/GPLinearODEMaker


text:
using Cubature
using Statistics

function profile(ϕ, params);

    ϕmid = (ϕ[1] + ϕ[end])/2
    Δη = params.Uj * params.Lj * params.TwoΩ / params.g

    U_function(x)  =  params.Uj * sech.((x .- ϕmid) / params.Lj).^2
    LHS(x) = - (params.TwoΩ .* sin.(x) .* params.a .+ U_function(x).*tan.(x)) .* U_function(x) / params.g

    U =  U_function(ϕ)
    E = zeros((params.Nϕ+1))
    for i in 1:params.Nϕ
        result = pquadrature(LHS, ϕ[end], ϕ[i])
        E[i] =  result[1]
    end
    E = E .- mean(E)

    return U,E
end


text:
function my_f(x, y)
    return 2*x + 3*y
end;

function derivative_of_my_f(x, y)
    return ForwardDiff.derivative(t -> my_f(t, y), x)
end;

max_stars_repo_name: DomagojGalic/MyExample

max_issues_repo_name: DomagojGalic/MyExample

max_forks_repo_name: DomagojGalic/MyExample


text:
using PSSFSS
using PSSFSS.Rings: Ring
using Test

@testset "Rings" begin
    @test collect(Ring(0)) == [(0,0)]
    @test collect(Ring(1)) == [(1, -1), (1, 0), (1, 1), (-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1)]
    @test collect(Ring(2)) == [(2, -2), (2, -1), (2, 0), (2, 1), (2, 2), (-2, -2), (-2, -1), (-2, 0),
                               (-2, 1), (-2, 2), (-1, -2), (0, -2), (1, -2), (-1, 2), (0, 2), (1, 2)]
    @test length(collect(Ring(10))) == length(Ring(10)) == 80
    @test all(map(x -> max(abs(x[1]),abs(x[2])), collect(Ring(100))) .== 100)
end


text:
using StaticArrays

struct TruncatedTriangle{T}
    midpoint::SVector{3,T}
    rot::SMatrix{3,3,T,9}
    r::T
    k::T
    d_eq::T
    t::T
    vertices::SVector{6,SVector{2,T}}
end

function TruncatedTriangle(midpoint::SVector{3,T}, rot::SMatrix{3,3,T,9}, r::T, k::T, d_eq::T) where T
    L = cbrt(1/k * 8/3 * (2*r+1)^3 / ((r^2+4*r+1)*(r+1)) * 4*π/3) * (d_eq / 2)
    H = L * sqrt(T(3)) / 2
    h = (r+1)/(2*r+1) * H
    t = k * h

    a_long = 1/(2*r+1) * L
    a_short = r * a_long

    # vertices in coordinates of the triangle
    vertices = SA[
        SA[a_long/2, -H/3],
        SA[a_long/2 + a_short/2, sqrt(T(3))/2*a_short - H/3],
        SA[a_short/2, h - H/3],
        SA[-a_short/2, h - H/3],
        SA[-a_long/2 - a_short/2, sqrt(T(3))/2*a_short - H/3],
        SA[-a_long/2, -H/3]
    ]

    return TruncatedTriangle(midpoint, rot, r, k, d_eq, t, vertices)
end

struct TTGenerator{Fr,Fk,Fd}
    r::Fr
    k::Fk
    d_eq::Fd
end

function (g::TTGenerator)(midpoint, rot)
    #midpoint = L * @SVector rand(3)
    #rot = GeometryTools.random_rotation()
    TruncatedTriangle(midpoint, rot, g.r(), g.k(), g.d_eq())
end

volume(shape::TruncatedTriangle) = π/6 * shape.d_eq^3

function set_midpoint(t::TruncatedTriangle, midpoint)
    return TruncatedTriangle(midpoint, t.rot, t.r, t.k, t.d_eq, t.t, t.vertices)
end

function bounding_box(shape::TruncatedTriangle)
    lower = shape.midpoint
    upper = shape.midpoint
    for vertex in shape.vertices
        vertex_a = shape.rot * SA[vertex[1], vertex[2], -shape.t/2]
        vertex_b = shape.rot * SA[vertex[1], vertex[2], shape.t/2]
        lower = min.(lower, vertex_a, vertex_b)
        upper = max.(upper, vertex_a, vertex_b)
    end
    return lower + shape.midpoint, upper + shape.midpoint
end

function inside(shape::TruncatedTriangle, r0)
    testline(p, line) = (p[2] - line[1][2])*(line[2][1] - line[1][1]) - (line[2][2] - line[1][2])*(p[1] - line[1][1])

    # Rotate to coordinates of the triangle as used above
    r = shape.rot' * (r0 - shape.midpoint)

    # the triangle is within the x-y plane.
    return -shape.t/2 < r[3] < shape.t/2 &&
        testline(r, shape.vertices[SA[3,2]]) < 0 &&
        testline(r, shape.vertices[SA[4,3]]) < 0 &&
        testline(r, shape.vertices[SA[5,4]]) < 0 &&
        testline(r, shape.vertices[SA[5,6]]) > 0 &&
        testline(r, shape.vertices[SA[6,1]]) > 0 &&
        testline(r, shape.vertices[SA[1,2]]) > 0
#    return -shape.t/2 < r[3] < shape.t / 2 &&
#        testline(r, shape.vertices[SA[1,2]]) > 0 &&
#        testline(r, shape.vertices[SA[2,3]]) > 0 &&
#        testline(r, shape.vertices[SA[3,4]]) > 0 &&
#        testline(r, shape.vertices[SA[4,5]]) > 0 &&
#        testline(r, shape.vertices[SA[5,6]]) > 0 &&
#        testline(r, shape.vertices[SA[6,1]]) > 0 &&
end


text:
function transit_orb(t,x,dt,nsub)
# Computes a transit lightcurve, normalized to unity, as a
# function of time, t, usually in HJD or BJD.
#
# Input parameters (x) are:
# x[1] = P  (units of day)
# x[2] = inc = inclination angle
# x[3] = p = R_p/R_* = radius of planet in units of radius of star
# x[4] = t0 = mid-point of transit
# x[5] = u1 = linear limb-darkening coefficient
# x[6] = u2 = quadratic limb-darkening coefficient
# x[7] = f0 = uneclipsed flux
# x[8] = a/R_* = semi-major axis divided by R_*
# x[9] = e = eccentricity
# x[10] = omega = longitude of pericentre
# Compute time of pericentre passage:
# The true anomaly at the time of transit:
f1=1.5*pi-x[10]*pi/180.0
ecc=x[9]
tp=(x[4]+x[1]*sqrt(1.0-ecc^2)/2.0/pi*(ecc*sin(f1)/(1.0+ecc*cos(f1))
    -2.0/sqrt(1.0-ecc^2)*atan2(sqrt(1.0-ecc^2)*tan(0.5*f1),1.0+ecc)))
fluxoft = 0.0
for j=1:nsub
  tsub = t-dt/2.0+dt*(j-0.5)/nsub
  m=2.0*pi/x[1]*(tsub-tp)
  f=kepler(m,ecc)
  radius=x[8]*(1.0-ecc^2)/(1.0+ecc*cos(f))
# Now compute sky separation:
  z0=radius*sqrt(1.0-(sin(x[2]*pi/180.0)*sin(x[10]*pi/180.0+f))^2)
  if (sin(x[10]*pi/180.0+f) < 0.0) && (z0 <= (1.0+x[3]))
    if x[3] < 1.0
      fluxoft += occultquad(z0,x[5],x[6],x[3])
    else
# We'll assume that the smaller object is not limb-darkened:
      fluxoft += occultquad(z0,  0.,  0.,x[3])
    end
  else
    fluxoft += 1.0
  end
end
# Take mean of flux over sub-time steps:
flux = fluxoft*x[7]/nsub
return flux
end


text:
using Test
using DataFlowTasks
using DataFlowTasks: R,W,RW, execute_dag
using LinearAlgebra

sch = DataFlowTasks.StaticScheduler()
DataFlowTasks.setscheduler!(sch)

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

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

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

    @testset "Tiled lu factorization" begin
        m  = 1000
        bsize = div(m,5)
        A = rand(m,m)
        t = tiled_lu(A,bsize)
        execute_dag(sch)
        F = fetch(t)
        @test F.L*F.U ≈ A
    end
end


text:
export CKKSParams, FixedRational

################################################################################
#                        CKKS Scheme definition
################################################################################

struct CKKSParams <: SHEShemeParams
    # The Cypertext ring over which operations are performed
    ℛ
    relin_window::Int
    σ::Float64
end
CKKSParams(ℛ, relin_window=default_relin_window(ℛ), σ=8/sqrt(2pi)) = CKKSParams(ℛ, relin_window, σ)
scheme_name(p::Type{CKKSParams}) = "CKKS"

# From the RLWE perspective, ℛ is both the plain and ciphertext. The encoder
# takes care of the conversion to/from complex numbers.
ℛ_plain(p::CKKSParams) = p.ℛ
ℛ_cipher(p::CKKSParams) = p.ℛ

π⁻¹(params::CKKSParams, plaintext) = params.ℛ(plaintext)
π(params::CKKSParams, b) = b

𝒩(params::CKKSParams) = RingSampler(params.ℛ, DiscreteNormal(0, params.σ))
𝒢(params::CKKSParams) = RingSampler(params.ℛ, DiscreteNormal(0, params.σ))

################################################################################
#                        CKKS Scheme definition
################################################################################

"""
    FixedRational{den, T<:Integer} <: Real
Rational number type, with numerator of type T and fixed denominator `den`.
"""
struct FixedRational{denom, T}
    x::T
    Base.reinterpret(::Type{FixedRational{denom, T}}, x::T) where {T,denom} =
        new{denom, T}(x)
    Base.reinterpret(::Type{FixedRational{denom}}, x::T) where {T,denom} =
        new{denom, T}(x)
    function FixedRational{denom, T}(x::Real) where {T, denom}
        n = round(BigInt, big(x)*denom)
        if n < 0
            n = modulus(T) + n
        end
        new{denom, T}(convert(T, n))
    end
end
drop_last(::Type{FixedRational{denom, T}}) where {T, denom} = FixedRational{denom, drop_last(T)}
drop_last(::Type{FixedRational{denom}}) where {denom} = FixedRational{denom}

function Base.convert(::Type{Float64}, fr::FixedRational{denom, T}) where {T, denom}
    n = convert(Integer, fr.x)
    if n > div(modulus(T), 2)
        n = n - modulus(T)
    end
    Float64(n/denom)
end

function Base.show(io::IO, fr::FixedRational{denom}) where {denom}
    print(io, convert(Float64, fr))
end

function maybe_wide_mul(a::Integer, b::Integer)
    T = promote_type(typeof(a), typeof(b))
    c = widemul(a, b)
    isa(T, BigInt) && return Float64(c)
    c < typemax(T) && return T(c)
    isa(c, BigInt) && return Float64(c)
    return c
end
maybe_wide_mul(a, b) = a*b

function maybe_wide_sq(a::Integer, b::Integer)
    T = promote_type(typeof(a), typeof(b))
    c = big(a)^b
    isa(T, BigInt) && return Float64(c)
    c < typemax(T) && return T(c)
    isa(c, BigInt) && return Float64(c)
    return c
end
maybe_wide_sq(a, b) = a^b

Base.:^(a::Type{FixedRational{denom, T}}, b::Int64) where {T,denom} = FixedRational{maybe_wide_sq(denom, b), T}
Base.:^(a::Type{FixedRational{denom}}, b::Int64) where {T,denom} = FixedRational{maybe_wide_sq(denom, b)}
Base.:*(a::Type{FixedRational{denom, T}}, b::Number) where {T,denom} = FixedRational{b/denom, T}
Base.:*(a::Type{FixedRational{denom}}, b::Number) where {T,denom} = FixedRational{b/denom, T}
Base.:*(a::Type{FixedRational{denom1, T}}, b::Type{FixedRational{denom2, T}}) where {T,denom1,denom2} = FixedRational{maybe_wide_mul(denom1, denom2), T}
Base.:*(a::Type{FixedRational{denom1}}, b::Type{FixedRational{denom2}}) where {denom1,denom2} = FixedRational{maybe_wide_mul(denom1, denom2)}
StructArrays.createinstance(FR::Type{FixedRational{denom, T}}, x::T) where {T,denom} =
    reinterpret(FR, x)


text:
# [69](https://projecteuler.net/problem=69)

#Euler totient function looks like
#```math
#\varphi(n) = n \prod_{p|n} ( 1 - \frac{1}{p} )
#```
#
#To find the solution n* which maximizes our object function
#```math
#\text{max} \frac{n}{\varphi(n)} = \frac{1}{ \prod\limits_{p|n} (1-\frac{1}{p}) }, n \leq 1000000
#```
#
#is equivalent to
#```math
#\text{min} \prod_{p|n} (1-\frac{1}{p}), n \leq 1000000
#```
#
#Distinct prime factors $`p_i \in \{p|n\}`$ are always positive integers that are larger than 1,
#hence $`0 < 1-\frac{1}{p} < 1`$ always holds. To minimize the above object function, we need
#as many distince prime factors as possible from the number n*. Now we comprehend this problem
#as to figure out a integer n* where n* <= 1000000 and has the most distinct prime factors
#among the ingeters less or equal to itself.
#
#Let's think about this problem in the reverse direction. The most ideal integer for this problem
#should ship all possible primes, e.g. $`n^* =\prod([2,3,5,7,11,\ldots])`$. Moreover, there are infinite
#number of primes, and the constraint $`n\leq 1000000`$ is exactly telling us when we should stop
#the infinite production.
#
for i in 1:20
    @printf "%d\t%8d\t%s\n" i prod(primes(i)) "$(prod(primes(i))<1_000_000)"
end

#Output
#```
#1	       1	true
#2	       2	true
#3	       6	true
#4	       6	true
#5	      30	true
#6	      30	true
#7	     210	true
#8	     210	true
#9	     210	true
#10	     210	true
#11	    2310	true
#12	    2310	true
#13	   30030	true
#14	   30030	true
#15	   30030	true
#16	   30030	true
#17	  510510	true
#18	  510510	true    *
#19	 9699690	false
#20	 9699690	false
#```

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
using YAML, DataFrames, CSV, Plots
using Statistics
using Streamfall

here = @__DIR__
data_path = joinpath(here, "../test/data/campaspe/")

# Load and generate stream network
network = YAML.load_file(joinpath(data_path, "campaspe_network.yml"))
sn = create_network("Example Network", network)

# Load climate data
date_format = "YYYY-mm-dd"
climate_data = CSV.File(joinpath(data_path, "climate/climate_historic.csv"),
                          comment="#",
                          dateformat=date_format) |> DataFrame

dam_level_fn = joinpath(data_path, "dam/historic_levels_for_fit.csv")
dam_releases_fn = joinpath(data_path, "dam/historic_releases.csv")
hist_dam_levels = CSV.File(dam_level_fn, dateformat=date_format) |> DataFrame
hist_dam_releases = CSV.File(dam_releases_fn, dateformat=date_format) |> DataFrame

# Combine recorded extractions and releases together
# hist_dam_releases[:, "406000_releases_[ML]"] .+ hist_dam_releases[:, "406000_extractions_[ML]"]
hist_dam_releases[:, "406000_extractions_[ML]"] = hist_dam_releases[:, "406000_releases_[ML]"]

# Subset to same range
climate_data, hist_dam_levels, hist_dam_releases = Streamfall.align_time_frame(climate_data,
                                                                               hist_dam_levels,
                                                                               hist_dam_releases)

climate = Climate(climate_data, "_rain", "_evap")

@info "Running example stream..."

reset!(sn)

dam_id, dam_node = sn["406000"]
Streamfall.run_node!(sn, dam_id, climate; extraction=hist_dam_releases)

h_data = hist_dam_levels[:, "Dam Level [mAHD]"]
n_data = dam_node.level

nnse_score = Streamfall.NNSE(h_data, n_data)
nse_score = Streamfall.NSE(h_data, n_data)
rmse_score = Streamfall.RMSE(h_data, n_data)

@info "Obj Func Scores:" rmse_score nnse_score nse_score

nse = round(nse_score, digits=4)
rmse = round(rmse_score, digits=4)

# Displaying results and saving figure
# plot(h_data,
#      legend=:bottomleft,
#      title="Calibrated IHACRES\n(RMSE: $(rmse); NSE: $(nse))",
#      label="Historic", xlabel="Day", ylabel="Dam Level [mAHD]")

# display(plot!(n_data, label="IHACRES"))

# savefig("calibrated_example.png")


text:
using MetaGen
using Gen

num_videos = 2
num_frames = 2

objects_observed = Matrix{Array{Detection2D}}(undef, num_videos, num_frames)
camera_trajectories = Matrix{Camera_Params}(undef, num_videos, num_frames)
for v=1:num_videos
	f = 1

	labels = [2, 5]
	xs = [50.1, 225.5]
	ys = [60.2, 173.1]
	temp = Array{Detection2D}(undef, length(labels))
	for i = 1:length(labels)
		label = labels[i]
		x = xs[i]
		y = ys[i]
		temp[i] = (x, y, label) #so now I have an array of detections
	end
	#turn that array of detections into an array of an array of detections sorted by receptive_field
	#temp_sorted_into_rfs = map(rf -> filter(p -> within(p, rf), temp), receptive_fields)
	objects_observed[v, f] = temp

	c = Camera_Params(camera_location = Coordinate(0.01,0.02,0.001), camera_focus = Coordinate(1.0,1.0,1.0))
	camera_trajectories[v, f] = c

	#################################
	#frame 2
	for f = 2:2

		labels = [2]
		xs = [50.1]
		ys = [60.2]
		temp = Array{Detection2D}(undef, length(labels))
		for i = 1:length(labels)
			label = labels[i]
			x = xs[i]
			y = ys[i]
			temp[i] = (x, y, label) #so now I have an array of detections
		end
		#turn that array of detections into an array of an array of detections sorted by receptive_field
		#temp_sorted_into_rfs = map(rf -> filter(p -> within(p, rf), temp), receptive_fields)
		#objects_observed[v, f] = temp_sorted_into_rfs
		objects_observed[v, f] = temp

		c = Camera_Params(camera_location = Coordinate(0.01,0.02,0.001), camera_focus = Coordinate(1.0,1.0,1.0))
		camera_trajectories[v, f] = c
	end
end

cm = choicemap()
for v = 1:num_videos
	# cm[:videos => v => :v_matrix => :lambda_fa => 1 => :fa] = 0.000001
	# cm[:videos => v => :v_matrix => :lambda_fa => 2 => :fa] = 0.000001
	# cm[:videos => v => :v_matrix => :lambda_fa => 3 => :fa] = 0.000001
	# cm[:videos => v => :v_matrix => :lambda_fa => 4 => :fa] = 0.000001
	# cm[:videos => v => :v_matrix => :lambda_fa => 5 => :fa] = 0.1 #super high false alarm rate for category 5
	# cm[:videos => v => :v_matrix => :miss_rate => 1 => :miss] = 0.000001 #tiny miss rate for category 1
	# cm[:videos => v => :v_matrix => :miss_rate => 2 => :miss] = 0.01
	# cm[:videos => v => :v_matrix => :miss_rate => 3 => :miss] = 0.01
	# cm[:videos => v => :v_matrix => :miss_rate => 4 => :miss] = 0.99 #huge miss rate for category 4
	# cm[:videos => v => :v_matrix => :miss_rate => 5 => :miss] = 0.01
	for f = 1:num_frames
		camera_params = camera_trajectories[v, f]
		cm[:videos => v => :frame_chain => f => :camera => :camera_location_x] = camera_params.camera_location.x
		cm[:videos => v => :frame_chain => f => :camera => :camera_location_y] = camera_params.camera_location.y
		cm[:videos => v => :frame_chain => f => :camera => :camera_location_z] = camera_params.camera_location.z
		cm[:videos => v => :frame_chain => f => :camera => :camera_focus_x] = camera_params.camera_focus.x
		cm[:videos => v => :frame_chain => f => :camera => :camera_focus_y] = camera_params.camera_focus.y
		cm[:videos => v => :frame_chain => f => :camera => :camera_focus_z] = camera_params.camera_focus.z

		#for rf = 1:num_receptive_fields
			#println("objects_observed[v, f][rf] ", objects_observed[v, f][rf])
			#println("type ", typeof(objects_observed[v, f][rf]))
		cm[:videos => v => :frame_chain => f => :observations_2D] = convert(Array{Any, 1}, objects_observed[v, f])
		#end
	end
	#def
end

gt_trace,_ = Gen.generate(main, (num_videos, num_frames), cm);
#println(gt_trace)
gt_choices = get_choices(gt_trace)

#they run fine when main isn't (static)
new_trace,_ = regenerate(gt_trace, select(:videos => 1 => :init_scene));
gt_choices = get_choices(new_trace)
regenerate(gt_trace, select(:videos => 2 => :v_matrix => :lambda_fa => 1 => :fa));

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
using Test

using DesignAutomation

@testset "A* example A" begin
    # node = DesignAutomation.airport_dfs(3, "../src/data/A.csv", 10)
    # print(node.matrix)

    node = DesignAutomation.airport_a_star(3, "../src/data/A.csv", 8)

    expected = [
        -1  -1  -1  -1  -1
        -1   0   0   0  -1
        -1   1   2   0  -1
        -1   1   0   0  -1
        -1   1   2   0  -1
        -1   1   0   0  -1
        -1   1   2   0  -1
        -1   0   0   0  -1
        -1  -1  -1  -1  -1
    ]

    @test node.matrix == expected
    @test node.g == 4
    @test node.h == 0
    @test node.depth == 8
end

@testset "A* example B" begin
    node = DesignAutomation.airport_a_star(4, "../src/data/B.csv", 8)

    expected = [
        -1  -1  -1  -1  -1  -1  -1
        -1   0   0   1   0   0  -1
        -1   0   2   1   2   0  -1
        -1   0   0   1   0   0  -1
        -1   0   2   1   2   0  -1
        -1   0   0   0   0   0  -1
        -1  -1  -1  -1  -1  -1  -1
    ]

    @test node.matrix == expected
    @test node.g == 3
    @test node.h == 0
    @test node.depth == 8
end

@testset "A* example C" begin
    node = DesignAutomation.airport_a_star(6, "../src/data/C.csv", 9)

    expected = [
        -1  -1  -1  -1  -1  -1  -1  -1  -1
        -1   0   0   0   0   0   0   0  -1
        -1   0   0   2   0   2   0   0  -1
        -1   0   2   1   1   1   2   0  -1
        -1   0   0   2   0   2   0   0  -1
        -1   0   0   0   0   0   0   0  -1
        -1  -1  -1  -1  -1  -1  -1  -1  -1
    ]

    @test node.matrix == expected
    @test node.g == 2
    @test node.h == 0
    @test node.depth == 9
end

@testset "A* example D" begin
    node = DesignAutomation.airport_a_star(9, "../src/data/D.csv", 18)

    expected = [
        -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
        -1   0   0   0   0   0   0   0   0  -1
        -1   1   2   0   0   2   0   0   0  -1
        -1   1   0   2   0   1   2   0   0  -1
        -1   0   1   1   1   1   0   0   0  -1
        -1   0   2   0   2   0   1   2   0  -1
        -1   0   0   0   0   2   1   0   0  -1
        -1   0   0   0   0   0   2   0   0  -1
        -1   0   0   0   0   0   0   0   0  -1
        -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
    ]

    @test node.matrix == expected
    @test node.g == 8
    @test node.h == 0
    @test node.depth == 18
end

@testset "A* example E" begin
    node = DesignAutomation.airport_a_star(9, "../src/data/E.csv", 19)

    expected = [
        -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
        -1   1   0   0  -1  -1  -1  -1  -1  -1
        -1   1   2   0   0  -1  -1  -1  -1  -1
        -1   1   0   2   0   0  -1  -1  -1  -1
        -1   0   1   1   1   2   0  -1  -1  -1
        -1   0   2   0   1   0   2   0  -1  -1
        -1   0   0   2   1   1   1   2   0  -1
        -1   0   0   0   2   0   2   0   0  -1
        -1   0   0   0   0   0   0   0   0  -1
        -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
    ]

    @test node.matrix == expected
    @test node.g == 9
    @test node.h == 0
    @test node.depth == 19
end

@testset "A* example F" begin
    node = DesignAutomation.airport_a_star(11, "../src/data/F.csv", 23)

    expected = [
        -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
        -1   0   1   0   0   0  -1  -1  -1   0   0   0  -1
        -1   0   0   1   2   0  -1  -1   0   0   0   0  -1
        -1   0   2   1   0   0  -1  -1   0   2   0   0  -1
        -1   0   0   1   2   0  -1  -1   0   1   2   0  -1
        -1   0   2   1   0   2   0  -1   1   0   0   0  -1
        -1   0   0   0   1   1   1   1   1   2   0   0  -1
        -1   0   0   0   2   0   2   0   2   0   0   0  -1
        -1   0   0   0   0   0   0   0   0   0   0   0  -1
        -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
    ]

    @test node.matrix == expected
    @test node.g == 11
    @test node.h == 0
    @test node.depth == 23
end

max_stars_repo_name: camirmas/DesignAutomation

max_issues_repo_name: camirmas/DesignAutomation

max_forks_repo_name: camirmas/DesignAutomation


text:
using Documenter

include("../src/NewsLookout.jl")

makedocs(
    sitename = "NewsLookout",
    format = Documenter.HTML(),
    modules = [NewsLookout]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(repo = "github.com/sandeep-sandhu/NewsLookout.jl.git")


text:
using SpecialSets

set_context!(Context([
    Associative.([
        +,
        *,
        &,
        |,
        ⊻,
    ]);
    Commutative.([
        +,
        &,
        |,
        ⊻,
    ]);
    Closure(+, TypeSet(Number))
    Closure(+, TypeSet(Int))
    Closure(-, TypeSet(Number))
    Closure(-, TypeSet(Int))
    Closure(*, TypeSet(Number))
    Closure(*, TypeSet(Int))
    Signature(inv, [TypeSet(Number)], TypeSet(Number))
    Signature(/, [TypeSet(Number), TypeSet(Number)], TypeSet(Float64))
    Signature(^, [Positive, TypeSet(Number)], Positive)
    Signature(^, [Zero, TypeSet(Number)], Set([0, 1]))
    Signature(^, [TypeSet(Number), Even], Nonnegative)
    Signature(^, [TypeSet(Real), TypeSet(Real)], TypeSet(Real))
    Signature(^, [TypeSet(Number), TypeSet(Number)], TypeSet(Number))
    Signature(abs, [TypeSet(Number)], Nonnegative)
    Signature(sqrt, [TypeSet(Real)], Nonnegative)
    Signature(sin, [TypeSet(Real)], GreaterThan{Number}(-1, true) ∩ LessThan{Number}(1, true))
    Signature(cos, [TypeSet(Real)], GreaterThan{Number}(-1, true) ∩ LessThan{Number}(1, true))
    Signature(tan, [TypeSet(Number)], TypeSet(Number))
    Signature(log, [TypeSet(Number)], TypeSet(Float64))
    Signature(diff, [TypeSet(Number), TypeSet(Number)], TypeSet(Number))  # FIXME
    Closure(&, TypeSet(Bool))
    Closure(|, TypeSet(Bool))
    Signature(!, [TypeSet(Bool)], TypeSet(Bool))
]))


text:
module VarBCFiles

using Dates, Printf, DataStructures, DataFrames
import Base: read, write, show, size, getindex, ==, merge!, length, filter!, keys
import DataFrames: DataFrame
export VarBC, VarBCRecord, increment, label,keys, readParamIncrement

include("VarBCRecord.jl")

struct VarBC  # <: AbstractDict{String, VarBCRecord}
    datetime::DateTime
    version::String
    header1::String
    header2::String
    records::OrderedDict{String,VarBCRecord}
end

length(a::VarBC) = length(a.records)
show(io::IO, ::MIME"text/plain", a::VarBC) = println(io, "VarBC with $(length(a)) records for $(a.datetime)")

"""
    filter!(f,a)

Update VarBCRecords in `a`, removing elements for which `f` is `false`. The function `f` is passed `key=>value` pairs

# Example
```jldoctest

julia> filter!(p-> p.second.label == "sfcobs", a)

```
"""
filter!(f,a::VarBC) = filter!(f,a.records)

"""
    merge!(a,b)

Merges the VarBC Records from `b` into `a`. If a key is present in both `a` and `b` the value from `b` is used.
Header information in `a` is not taken from `b`
"""
merge!(a::VarBC,b::VarBC) = merge!(a.records,b.records)

"""
    a == b

Returns `true` if `a` and `b` are equal
"""
function ==(a::VarBC,b::VarBC)
    a.datetime == b.datetime &&
    a.header1 == b.header1 &&
    b.header2 == b.header2 &&
    a.version == a.version &&
    a.records == a.records
end

"""
    keys(a)

Returns keys
"""
keys(a::VarBC) = keys(a.records)

"""
    read(filename,VarBC)

returns a object of type `VarBC`
"""
function read(fname::String,::Type{VarBC})
    io = open(fname)
    version = readline(io); @assert version == "VARBC_cycle.version006"
    header1 = readline(io)
    header2 = readline(io)
    str, date, time = split(header1)
    numrecord, dummy1 = split(header2)
    dt = DateTime("$date$(lpad(time,6,"0"))","yyyymmddHHMMSS")
    out = VarBC(dt,version,header1, header2,OrderedDict{String,VarBCRecord}())
    while !eof(io)
       push!(out.records, read(io,VarBCRecord))
    end
    close(io)
    return out
end

"""
    readParamIncrement(filename,key)

returns namedTuple of paramIncrement
"""
function readParamIncrement(fname::String,key)
    io = open(fname)
    version = readline(io);
    @assert version == "VARBC_cycle.version006"
    header1 = readline(io)
    header2 = readline(io)
    str, date, time = split(header1)
    dt = DateTime("$date$(lpad(time,6,"0"))","yyyymmddHHMMSS")
    while !eof(io)
       rec = read(io,VarBCRecord)
       val = last(rec)
       if key== first(rec)
         return ( ;datetime=dt, zip(Symbol.("p",val.predcs),val.params-val.param0)...)
       end
    end
    close(io)
end

"""
    write(fname,a)

Write VarBC struct to file
"""
function write(fname::String,a::VarBC)
    io = open(fname,"w")
    write(io,a.version,"\n")
    write(io,a.header1,"\n")
    write(io,a.header2,"\n")
    i=1
    for (key,rec) in a.records
        write(io,rec,key,i)
        i=i+1
    end
    close(io)
end

end # module


text:
using Documenter, EconUtils

makedocs(
    format = :html,
    sitename = "EconUtils.jl",
    pages = [
        "index.md",
        "GettingStarted.md",
        "API.md",
        "Examples.md",
        "References.md"
    ]
)

deploydocs(
    deps = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
    repo = "github.com/JuliaEconometrics/EconUtils.jl.git",
    julia  = "0.7.0-DEV.2279"
)


text:
function findpivot(ss::StandardSimplexe{T}; verbose::Bool = false)::Tuple{Int, Int} where T
    if isdegenerated(ss)
        return BlandFindpivot(ss, verbose = verbose)
    end
    M = ss.M
    m = size(M, 1) - 1
    c = @view M[end, 1:end-1]
    b = @view M[1:end-1, end]
    entering = argmin(c)
    verbose && @show entering
    @assert c[entering] < 0 "simplexe is in an optimal state" #ss was optimal
    val = T(Inf)
    leaving = -1
    for k in 1:m
        if (ss.freelines[k]) && (M[k, entering] > 0) && ((tobeat = b[k]/M[k, entering]) < val)
            leaving = k
            val = tobeat
        end
    end
    verbose && @show leaving
    return leaving, entering
end
function BlandFindpivot(ss::StandardSimplexe{T}; verbose::Bool = verbose)::Tuple{Int, Int} where T
    M = ss.M
    m = size(M, 1) - 1
    c = @view M[end, 1:end-1]
    b = @view M[1:end-1, end]
    entering = findfirst(t->t<zero(T), c)
    verbose && @show entering
    @assert c[entering] < 0 "simplexe is in an optimal state" #ss was optimal
    val = T(Inf)
    leaving = -1
    for k in 1:m
        if (ss.freelines[k]) && (M[k, entering] > 0) && ((tobeat = b[k]/M[k, entering]) < val)
            leaving = k
            val = tobeat
        end
    end
    verbose && @show leaving
    return leaving, entering
end
function isdegenerated(ss::StandardSimplexe)::Bool
    M = ss.M
    vTD = @view M[1:end-1, end]
    return any(iszero, vTD[ss.freelines])
end

function phase2!(ss::StandardSimplexe{T}; verbose::Bool = false, kmax = 10000)::Status where T
    k = 0
    m2, n2 = size(ss.M)
    m, n = m2 - 1,  n2 -1
    while !isOptimal(ss) && (k+=1) <= kmax
        verbose && @show k
        verbose && println(ss)
        i, j = findpivot(ss, verbose = verbose)
        if i == -1
            ss.status = Unbounded
            break
        end
        pivot!(ss.M, i, j, ss.Binv, freelines = ss.freelines)
        ss.b_idx[i] = j
    end
    if isOptimal(ss)
        ss.status = Optimal
        ss.xstar = zeros(T, n)
        freebasis = ss.b_idx[ss.freelines]
        xB = ss.M[[ss.freelines; false], end]
        ss.xstar[freebasis] = xB
        ss.vstar = -ss.M[end, end]
    end

    status = ss.status
    verbose && (println("$status Simplexe");  println(ss))
    verbose && @show ss.xstar
    verbose && @show ss.vstar
    ss.status
end

max_stars_repo_name: fbastin/Jasmin

max_issues_repo_name: fbastin/Jasmin

max_forks_repo_name: fbastin/Jasmin


text:
#!/usr/bin/julia -f

import FunctionWrappers
import FunctionWrappers: FunctionWrapper
using Base.Test

immutable CallbackF64
    f::FunctionWrapper{Float64,Tuple{Int}}
end
(cb::CallbackF64)(v) = cb.f(v)
gen_closure(x) = y->x + y

@testset "As field" begin
    f1 = @inferred CallbackF64(identity)
    f2 = @inferred CallbackF64(sin)
    f3 = @inferred CallbackF64(gen_closure(2))
    @test typeof(f1) == typeof(f2) == typeof(f3)
    @test @inferred(f1(1)) === 1.0
    @test @inferred(f2(1)) === sin(1)
    @test @inferred(f3(1)) === 3.0
end

typealias F64AnyFunc FunctionWrapper{Float64,Tuple{Any}}

@testset "Any input" begin
    f1 = @inferred F64AnyFunc(identity)
    f2 = @inferred F64AnyFunc(gen_closure(2))
    @test typeof(f1) === typeof(f2)
    @test @inferred(f1(1)) === 1.0
    @test @inferred(f1(1.0)) === 1.0
    @test @inferred(f1(1f0)) === 1.0
    @test @inferred(f2(1)) === 3.0
end

typealias F64F64Func FunctionWrapper{Float64,Tuple{Float64}}

@testset "Convert" begin
    f1 = @inferred F64F64Func(sin)
    f2 = @inferred F64F64Func(f1)
    f3 = @inferred convert(F64F64Func, f1)
    @test f1 === f2
    @test f1 === f3
end

typealias NumberAnyFunc FunctionWrapper{Number,Tuple{Any}}

@testset "Abstract Return" begin
    @test NumberAnyFunc(sin)(1) === sin(1)
    @test NumberAnyFunc(identity)(1) === 1
end

@testset "Precompile" begin
    @test FunctionWrappers.identityAnyAny(1) === 1
end


text:
"""
    Field(nx::Int64, ny::Int64, dx::Float64, dy::Float64, data::Matrix{Float64})

Temperature field type. nx and ny are the dimensions of the field.
The array data contains also ghost layers, so it will have dimensions
[nx+2, ny+2]
"""
mutable struct Field{T<:AbstractArray}
    nx::Int64
    ny::Int64
    # Size of the grid cells
    dx::Float64
    dy::Float64
    # The temperature values in the 2D grid
    data::T
end

# outer constructor with default cell sizes and initialized data
Field(nx::Int64, ny::Int64, data) = Field{typeof(data)}(nx, ny, 0.01, 0.01, data)

# extend deepcopy to new type
Base.deepcopy(f::Field) = Field(f.nx, f.ny, f.dx, f.dy, deepcopy(f.data))

"""
    initialize(rows::Int, cols::Int)

Initialize two temperature field with (nrows, ncols) number of
rows and columns.
"""
function initialize(nrows = 1000, ncols = 1000, arraytype = Matrix)

    data = arraytype(zeros(nrows+2, ncols+2))
    previous = Field(nrows, ncols, data)

    # generate a specific field with boundary conditions
    generate_field!(previous)
    current = Base.deepcopy(previous)

    return previous, current
end

"""
    generate_field!(field0::Field)

Generate a temperature field.  Pattern is disc with a radius
of nx / 6 in the center of the grid. Boundary conditions are
(different) constant temperatures outside the grid.
"""
function generate_field!(field::Field)
    # Square of the disk radius
    radius2 = (field.nx / 6.0)^2

    for j = 1:field.ny+2
        for i = 1:field.nx+2
            ds2 = (i - field.nx / 2)^2 + (j - field.ny / 2)^2
            if ds2 < radius2
                field.data[i,j] = 5.0
            else
                field.data[i,j] = 65.0
            end
        end
    end

    # Boundary conditions
    field.data[:,1] .= 20.0
    field.data[:,field.ny+2] .= 70.0
    field.data[1,:] .= 85.0
    field.data[field.nx+2,:] .= 5.0
end


text:
# One action, random +1/-1 observation, one timestep long, obs-dependent +1/-1 reward every time:
#If my agent can learn the value in (1.) but not this one - meaning it can learn a constant
#reward but not a predictable one! - it must be that backpropagation through my network is broken.
import AlphaZero.GI

struct GameSpec <: GI.AbstractGameSpec end

mutable struct GameEnv <: GI.AbstractGameEnv
  terminated :: Bool
end

const INITSTATE = (; terminated=false)

GI.init(::GameSpec, state=INITSTATE) = GameEnv(state.terminated)

GI.spec(::GameEnv) = GameSpec()

GI.two_players(::GameSpec) = false

const ACTIONS = collect(1:1)

GI.actions(::GameSpec) = ACTIONS

GI.actions_mask(g::GameEnv) = [true]

GI.current_state(g::GameEnv) = (; terminated=g.terminated)

GI.white_playing(g::GameEnv) = true

GI.game_terminated(g::GameEnv) = g.terminated

GI.white_reward(g::GameEnv) = rand([-1.,1.])

function GI.play!(g::GameEnv, action)
  g.terminated = true
end

GI.vectorize_state(::GameSpec, state) = Float32[!state.terminated]


text:
# used for constructing appropraite dict from YAML object.
function construct_type_map(t::Symbol, constructor::YAML.Constructor,
                            node::YAML.MappingNode)
    mapping = _symbol_dict(YAML.construct_mapping(constructor, node))
    mapping[:tag] = t
    mapping
end

function construct_type_map(t::Symbol, constructor::YAML.Constructor,
                            node::YAML.SequenceNode)
    YAML.construct_sequence(constructor, node)
end

const yaml_types = let
    pairs = [("!Cartesian", :Cartesian),
            ("!Smolyak", :Smolyak),
            ("!Random", :Random),
             ("!Normal", :Normal),
             ("!MarkovChain", :MarkovChain),
             ("!Product", :Product),
             ("!PoissonProcess", :PoissonProcess),
             ("!DeathProcess", :DeathProcess),
             ("!AgingProcess", :AgingProcess),
             ("!VAR1", :VAR1)]
    Dict{AbstractString,Function}([(t, (c, n) -> construct_type_map(s, c, n))
                           for (t, s) in pairs])
end

mutable struct SModel{ID} <: ASModel{ID}
    data::Dict{Symbol,Any}
end

function SModel(url::AbstractString)
    if match(r"(http|https):.*", url) != nothing
        res = HTTP.request("GET", url)
        txt = (String(res.body))
        data = _symbol_dict(YAML.load(txt, yaml_types))
    else
        data = _symbol_dict(YAML.load_file(url, yaml_types))
    end
    id = gensym(basename(url))
    return SModel{id}(_symbol_dict(data))
end

function get_symbols(model::ASModel)
    syms = model.data[:symbols]
    symbols = OrderedDict{Symbol,Vector{Symbol}}()
    for k in keys(syms)
        symbols[Symbol(k)] = [Symbol(e) for e in syms[k]]
    end
    return symbols
end

function get_variables(model::ASModel)
    symbols = get_symbols(model)
    vars = cat(values(symbols)...; dims=1)
    dynvars = setdiff(vars, symbols[:parameters] )
end

function get_name(model::ASModel)
    get(model.data, :name, "modeldoesnotwork")
end

### TODO: look for duplicates
eq_to_expr(eq::Symbol) = eq
function eq_to_expr(eq::Expr)
   if !(eq.head==:call && eq.args[1]==:(==))
      eq
   else
      rhs = eq.args[3]
      lhs = eq.args[2]
      :($rhs-$lhs)
   end
end

function get_equations(model::ASModel)

    eqs = model.data[:equations]
    recipe = RECIPES[:dtcc]
    _symbols = get_symbols(model)

    # prep equations: parse to Expr
    _eqs = OrderedDict{Symbol,Vector{Expr}}()
    # for k in Dolo.keys(recipe[:specs])
    for k in Dolo.keys(eqs)
      if !(k in Dolo.keys(recipe[:specs]))
           ii=0
           dist = zeros(length(Dolo.keys(recipe[:specs])))
           for kk in Dolo.keys(recipe[:specs])
             ii+=1
             dist[ii]=compare(Hamming(), string(kk), string(k))
           end
           kk = collect(keys(Dolo.keys(recipe[:specs]).dict))[findmax(dist)[2]]
           msg = string("You are using a wrong name for a sub-section title.`$(k)` is not accepted. ",
                        "Maybe you meant `$(kk)`?")
           error(msg)
      end
      (k in [:arbitrage,]) && continue
      these_eq = get(eqs, (k), [])
      # verify that we have at least 1 equation if section is required
      if !get(recipe[:specs][k], :optional, false)
          length(these_eq) == 0 && error("equation section $k required")
      end
      # finally pass in the expressions
      _eqs[k] = Expr[_to_expr(eq) for eq in these_eq]
    end
    # end
    # handle the arbitrage, arbitrage_exp, controls_lb, and controls_ub
    if haskey(recipe[:specs], :arbitrage) && (:arbitrage in keys(eqs))
        c_lb, c_ub, arb = _handle_arbitrage(eqs[:arbitrage],
                                            _symbols[:controls])
        _eqs[:arbitrage] = arb
        _eqs[:controls_lb] = c_lb
        _eqs[:controls_ub] = c_ub
    end

    defs = get_definitions(model)
    dynvars = cat(get_variables(model), [k[1] for k in keys(defs)]...; dims=1)

    for eqtype in keys(_eqs)
        _eqs[eqtype] = [sanitize(eq,dynvars) for eq in _eqs[eqtype]]
    end

    return _eqs
end

function get_infos(model::ASModel)
    get(model.data, :infos, Dict())
end

function get_options(model::ASModel; options=Dict())
    opts = get(model.data, :options, Dict())
    return rmerge(opts, options)
end

function get_definitions(model::ASModel)::OrderedDict{Tuple{Symbol,Int},SymExpr}
    # parse defs so values are Expr
    defs = get(model.data,:definitions, Dict())
    dynvars = cat(get_variables(model), keys(defs)...; dims=1)
    _defs = OrderedDict{Tuple{Symbol,Int},SymExpr}(
        [ (k,0) => sanitize(_to_expr(v), dynvars) for (k, v) in defs]
    )
    return _defs
end

function get_calibration(model::ASModel)
    calib = get(model.data, :calibration, Dict())
    # prep calib: parse to Expr, Symbol, or Number
    _calib  = OrderedDict{Symbol,Union{Expr,Symbol,Number}}()
    # add calibration for definitions
    for k in keys(calib)
        _calib[Symbol(k)] = _expr_or_number(calib[k])
    end
    # so far _calib is a symbolic calibration
    calibration = solve_triangular_system(_calib)
    symbols = get_symbols(model)
    return ModelCalibration(calibration, symbols)
end

function get_domain(model::ASModel)
    domain = deepcopy(get(model.data, :domain, Dict()))
    calib = get_calibration(model)
    # TODO deal with empty dict and make robust construction
    states = get_symbols(model)[:states]
    min = [eval_with(calib, domain[(k)][1]) for k in states]
    max = [eval_with(calib, domain[(k)][2]) for k in states]
    return Domain(states, min, max)
end

function get_grid(model::ASModel; options=Dict())
    domain = get_domain(model)
    d = length(domain.states)
    options = get_options(model; options=Dict(:grid=>options))
    grid_dict = options[:grid]
    if :type in keys(grid_dict)
        grid_dict[:tag] = grid_dict[:type]
    end
    if grid_dict[:tag] == :Cartesian
        orders = get(grid_dict, :orders, [20 for i=1:d])
        grid = CartesianGrid{d}(domain.min, domain.max, orders)
        if length(orders)!=length(model.calibration[:states])
            msg = string("Check the dimension of the matrix given in the yaml file, section: options-grid-orders. ",
                         "Expected to be of dimension $([length(model.calibration[:states])])")
            error(msg)
        end
    elseif grid_dict[:tag] == :Smolyak
        mu = get(grid_dict, :mu, 3)
        grid = SmolyakGrid{d}(domain.min, domain.max, mu)
    elseif grid_dict[:tag] == :Random
        n = get(grid_dict, :N, 200)
        grid = RandomGrid{d}(domain.min, domain.max, n)
    else
        error("Unknown grid type.")
    end
    return grid
end

function get_exogenous_old(model::ASModel)
    exo_dict = get(model.data,:exogenous,Dict{Symbol,Any}())
    if length(exo_dict)==0
        exo_dict = get(model.data[:options], :exogenous, Dict{Symbol,Any}())
    end
    exo = _symbol_dict(exo_dict)
    calib = get_calibration(model)
    exogenous = _build_exogenous_entry(exo, calib)
    return exogenous
end

#

function get_exogenous(model::AModel)
    rdata = model.rdata
    exo_dict = model.rdata[:exogenous]
    cond = !(exo_dict.tag=="tag:yaml.org,2002:map")
    if cond
        # old style exogenous block
        return get_exogenous_old(model)
    else
        syms = cat( [[Symbol(strip(e)) for e in split(k, ",")] for k in keys(exo_dict)]..., dims=1)
        expected = model.symbols[:exogenous]
        if (syms != expected)
            msg = string("In 'exogenous' section, shocks must be declared in the same order as shock symbol. Found: ", syms, ". Expected: ", expected, ".")
            throw(ErrorException(msg))
        end
        calibration = model.calibration.flat
        processes = []
        for k in keys(exo_dict)
            v = exo_dict[k]
            p = Dolang.eval_node(v, calibration, minilang, FromGreek())
            push!(processes, p)
        end
        return ProductProcess(processes...)
    end
end

function set_calibration!(model::ASModel, key::Symbol, value::Union{Real,Expr, Symbol})
    model.data[:calibration][key] = value
end

mutable struct Model{ID} <: AModel{ID}

    data::Dict{Symbol,Any}
    rdata::YAML.Node
    name::String           # weakly immutable
    filename::String
    symbols::OrderedDict{Symbol,Array{Symbol,1}}        # weakly immutable
    equations      # weakly immutable
    definitions    # weakly immutable
    # functions      # weakly immutable
    factories::Dict
    calibration::ModelCalibration
    exogenous
    domain
    grid
    options

    function Model{ID}(data::Dict{Symbol,Any}, rdata::YAML.Node; print_code::Bool=false, filename="<string>") where ID

        model = new{ID}(data, rdata)
        model.name = get_name(model)
        model.filename = filename
        model.symbols = get_symbols(model)
        model.equations = get_equations(model)
        model.definitions = get_definitions(model)
        model.factories = Dict()
        model.calibration = get_calibration(model)
        model.exogenous = get_exogenous(model)
        model.domain = get_domain(model)
        model.options = get_options(model)
        model.grid = get_grid(model)

        factories = Dict{Symbol,Dolang.FlatFunctionFactory}()
        for (eqtype,equations) in model.equations
            spec = Dolo.RECIPES[:dtcc][:specs][eqtype]
            spec_eqs = spec[:eqs]
            symbols = model.symbols
            arguments = OrderedDict{Symbol,Union{Vector{Symbol},Vector{Tuple{Symbol,Int}}}}([
                Symbol(s)=>[(e,t) for e in symbols[Symbol(sg)]]
                for (sg,t,s) in spec_eqs
            ])
            arguments[:p] = [e[1] for e in arguments[:p]]
            tdefs = Dolo.get_definitions(model)
            if :target in keys(spec)
                sg,t,s = spec[:target]
                targets = [(e,t) for e in symbols[Symbol(sg)]]
                lhs = [eq.args[2] for eq in equations]
                rhs = [eq.args[3] for eq in equations]
                expressions = OrderedDict(zip(targets,rhs))
                @assert Dolang.stringify.(lhs) == Dolang.stringify.(targets)
            else
                expressions = OrderedDict( [Symbol("out_",i)=>eq_to_expr(eq) for (i,eq) in enumerate(equations)])
            end
            ff = Dolang.FlatFunctionFactory(expressions, arguments, tdefs; funname=eqtype )
            factories[eqtype] = ff
            code = Dolang.gen_generated_gufun(ff; dispatch=typeof(model))
            print_code && println(code)
            Core.eval(Dolo, code)
        end

        model.factories = factories

        # TEMP: until we have a better method in Dolang
        # Create definitions
        vars = cat(model.symbols[:exogenous], model.symbols[:states], model.symbols[:controls]; dims=1)
        args = OrderedDict(
            :past => [(v,-1) for v in vars],
            :present => [(v,0) for v in vars],
            :future => [(v,1) for v in vars],
            :params => model.symbols[:parameters]
        )
        fff = Dolang.FlatFunctionFactory(model.definitions, args, typeof(model.definitions)())
        code = Dolang.gen_generated_gufun(fff;dispatch=typeof(model), funname=:evaluate_definitions)
        Core.eval(Dolo,code)

        return model
    end

end

_numeric_mod_type(::Model{ID}) where {ID} = Model{ID}

function Base.show(io::IO, model::Model)
    print(io, "Model")
end

#### import functions

function Model(url::AbstractString; print_code=false)
    # it looks like it would be cool to use the super constructor ;-)
    if match(r"(http|https):.*", url) != nothing
        res = HTTP.request("GET", url)
        txt = String(res.body)
    else
        txt = read(open(url), String)
    end
    data = _symbol_dict(load(txt, yaml_types))
    rdata = Dolang.yaml_node_from_string(txt)
    id = gensym()
    fname = basename(url)
    return Model{id}(data, rdata; print_code=print_code, filename=fname)
end

function set_calibration!(model::Model, key::Symbol, value::Union{Real,Expr, Symbol})
    model.data[:calibration][key] = value
    model.calibration = get_calibration(model)
    model.exogenous = get_exogenous(model)
    model.domain = get_domain(model)
    # model.options = get_options(model) # we don't substitute calib here
    model.grid = get_grid(model)
    model.calibration
end

function set_calibration!(model::Model, values::AbstractDict{Symbol,T}) where T
    for (key,value) in values
        model.data[:calibration][key] = value
    end
    model.calibration = get_calibration(model)
    model.exogenous = get_exogenous(model)
    model.domain = get_domain(model)
    # model.options = get_options(model) # we don't substitute calib here
    model.grid = get_grid(model)
    model.calibration
end

function set_calibration!(model::Model; kwargs...)
    set_calibration!(model, Dict(kwargs))
end

"""
Imports the model from a yaml file specified by the `url` input
parameter, and returns the corresponding `Dolo` model object.
"""
function yaml_import(url; print_code::Bool=false)
    Model(url; print_code=print_code)
end

function features(model::Dolo.ASModel)
    features = Dict{Symbol, Bool}()
    features[:one_state] = (length(model.symbols[:states])==1)
    features[:one_control] = (length(model.symbols[:controls])==1)
    features[:one_dimensional] = features[:one_state] && features[:one_state]
    rhs_transitions = [eq.args[3] for eq in model.equations[:transition]]
    it = Dolang.IncidenceTable( model.equations[:transition] )
    features[:nonstochastic_transitions] = (length( intersect(it.by_date[0], model.symbols[:exogenous] ) )==0)
    it = Dolang.IncidenceTable( [model.equations[:controls_lb]; model.equations[:controls_lb]] )
    features[:bounds_are_constant] = (length(it.by_var) == 0)
    return features
end


text:
@testset "ok" begin

@test parse_raw("a 1") == Any[Any[(:a,1)]]

@test parse_raw("a 1.0") == Any[Any[(:a,1.0)]]

@test parse_raw("a \"b\"") == Any[Any[(:a,"b")]]

@test parse_raw(" a 1") == Any[Any[(:a,1)]]

@test parse_raw(" a 1.0") == Any[Any[(:a,1.0)]]

@test parse_raw(" a \"b\"") == Any[Any[(:a,"b")]]

@test parse_raw("a 1 ") == Any[Any[(:a,1)]]

@test parse_raw("a 1.0 ") == Any[Any[(:a,1.0)]]

@test parse_raw("a \"b\" ") == Any[Any[(:a,"b")]]

@test parse_raw("a [b 1]") == Any[Any[(:a,Any[(:b,1)])]]

@test parse_raw(" a [b 1]") == Any[Any[(:a,Any[(:b,1)])]]

@test parse_raw("a [ b 1]") == Any[Any[(:a,Any[(:b,1)])]]

@test parse_raw("a [b 1 ]") == Any[Any[(:a,Any[(:b,1)])]]

@test parse_raw("a [b 1] ") == Any[Any[(:a,Any[(:b,1)])]]

@test parse_raw("a [b [c 1]]") == Any[Any[(:a,Any[(:b,Any[(:c,1)])])]]

@test parse_raw("a [b [c 1 d 2.0 e \"3\"]]") == Any[Any[(:a,Any[(:b,Any[(:c,1),(:d,2.0),(:e,"3")])])]]

end


text:
@testset "trace" begin

global_tracer_provider(
    TracerProvider(
        span_processor=MultiSpanProcessor(
            SimpleSpanProcessor(
                InMemorySpanExporter()
            )
        )
    )
)

tracer = get_tracer("test")

with_span("test_async", tracer) do
    @sync for i in 1:5
        @async with_span("asyn_span_$i", tracer) do
            current_span()["my_id"] = i
        end
    end
end

end


text:
using Markdown

struct SupportedMethods
    names::Vector{Symbol}
end

SupportedMethods() = SupportedMethods([
    n for n in names(PyBase; all=true) if occursin(r"^__[a-z]+__$", string(n))
])

function Base.show(io::IO, ::MIME"text/markdown", sm::SupportedMethods)
    println(io, "| **Supported method** |")
    println(io, "|:-- |")
    for n in sm.names
        println(io, "| ", _pymodelref(n), " |")
    end
end

function Markdown.MD(sm::SupportedMethods)
    io = IOBuffer()
    show(io, "text/markdown", sm)
    seek(io, 0)
    return Markdown.parse(io)
end

# Base.show(io::IO, sm::SupportedMethods) = show(io, "text/plain", sm)

Base.show(io::IO, mime::MIME"text/plain", sm::SupportedMethods) =
    _show(io, mime, sm)
Base.show(io::IO, mime::MIME"text/html", sm::SupportedMethods) =
    _show(io, mime, sm)

_show(io, mime, sm) = show(io, mime, Markdown.MD(sm))


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script to compile data for some descriptive statistics
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, CSV, XLSX, LinearAlgebra, Statistics, StatsBase

# other scripts
dir_home = "x:/VIVES/1-Personal/Florian/git/2022_01_comext/src/"
include(dir_home * "functions.jl")

# location of data input/output (io)
dir_io = "C:/Users/u0148308/data/comext/"
dir_dropbox = "C:/Users/u0148308/Dropbox/BREXIT/"

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

# Notes:
#   - import raw data, initial cleaning
#   - subset and compute data for figure
#   - loop over months and add to DataFrame
#   - export data
#   - import data and plot figures (different script)

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

# timespan
years = string.(2001:2021)
months = lpad.(1:12, 2, '0')

# EU27
ctrys_EU27 = ["AT", "BE", "BG", "CY", "CZ", "DE", "DK", "EE", "ES", "FI", "FR", "GR", "HR", "HU",
             "IE", "IT", "LT", "LU", "LV", "MT", "NL", "PL", "PT", "RO", "SE", "SI", "SK"]

# initialize dataframes
df_BE = DataFrame(DECLARANT_ISO=String[], PARTNER_ISO=String[], TRADE_TYPE=String[], PRODUCT_NC=String[], FLOW=String[], PERIOD=Int64[],
                VALUE_IN_EUROS=Union{Missing, Float64}[], QUANTITY_IN_KG=Union{Missing, Float64}[],
                PRODUCT_SITC=String[], PRODUCT_CPA2008=String[], PRODUCT_CPA2002=String[], PRODUCT_CPA2_1=String[],
                PRODUCT_BEC=String[], PRODUCT_BEC5=String[], PRODUCT_SECTION=String[])

df_fig1 = DataFrame(DECLARANT_ISO=String[], PARTNER_ISO=String[], FLOW=String[], PERIOD=Int64[],
                    VALUE_IN_EUROS=Float64[], QUANTITY_IN_KG=Float64[], PRICE_INDEX=Float64[])

df_fig2 = DataFrame(DECLARANT_ISO=String[], PARTNER_ISO=String[], FLOW=String[], PERIOD=Int64[], PRODUCT_NC_digits=String[],
                    VALUE_IN_EUROS=Float64[], QUANTITY_IN_KG=Float64[], PRICE_INDEX=Float64[])

df_fig3 = DataFrame(DECLARANT_ISO=String[], PARTNER_ISO=String[], FLOW=String[], PERIOD=Int64[],
                    VALUE_SHARE=Float64[], QUANTITY_SHARE=Float64[])

df_tab1 = DataFrame(DECLARANT_ISO=String[], PARTNER_ISO=String[], FLOW=String[], PRODUCT_NC_digits=String[], PERIOD=Int64[],
                    VALUE_IN_EUROS=Float64[], QUANTITY_IN_KG=Float64[])

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

only_BE = true
if only_BE == true
    println(" - Using individual country data for computations. \n")
else
    println(" - Using raw data from Eurostat for computations. \n")
end

for i in years
    for j in months

        # import and clean data
        if only_BE == true
            df = import_ctry(i, j, "BE")
        else
            df = initial_cleaning(i, j) # if computed from raw data with all DECLARANT
            df = append_WORLD(df)
            df = append_EU(df, ctrys_EU27)

            # possibly include in initial cleaning
            subset!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> !ismissing(x)))
            subset!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> x != 0))
            transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(Float64), renamecols=false)

        end

        # append data for figures/tables
        append!(df_BE, df) # try to produce one file for BE
        append!(df_fig1, data_fig1(df, ["BE"], ["GB", "EU", "DE", "FR", "NL"]))
        append!(df_fig2, data_fig2(df, ["BE"], ["GB"], 2))
        append!(df_fig3, data_fig3(df, ["BE"], ["GB", "EU", "DE", "FR", "NL"]))
        append!(df_tab1, tab1(df, ["BE"], ["GB", "EU", "DE", "FR", "NL"], 2))

        println(" ✓ Data for $j/$i has been successfully added. \n")

    end
end

# export Dropbox
CSV.write(dir_dropbox * "rawdata/comext/BE/" * "df_BE" * ".csv", df_BE)
CSV.write(dir_dropbox * "results/" * "df_fig1" * ".csv", df_fig1)
CSV.write(dir_dropbox * "results/" * "df_fig2" * ".csv", df_fig2)
CSV.write(dir_dropbox * "results/" * "df_fig3" * ".csv", df_fig3)
CSV.write(dir_dropbox * "results/" * "df_tab1" * ".csv", df_tab1)

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


text:
#===============================================================================

    AmpOutput.jl

    amplification output formats per step or ramp:
    long (full) and short (for conversion to JSON)

    Author: Tom Price
    Date:   June 2019

===============================================================================#

import DataStructures.OrderedDict
import StaticArrays: SArray, SMatrix, SVector
import Ipopt: IpoptSolver #, NLoptSolver

#===============================================================================
    Field definitions >>
===============================================================================#

const AMPLONGOUTPUT_FIELD_DEFS = [
    Field(:raw_data,                     Array{<: Real,3}), ## fr_ary3
    ## computed in amp_analysis()
    Field(:background_subtracted_data,   Array{<: Real,3}), ## mw_ary3
    Field(:norm_data,                    SArray{S,<: Real} where {S}), ## wva_data
    Field(:k_deconv,                     DeconvolutionMatrices), ## k4dcv
    Field(:deconvoluted_data,            Array{<: Real,3}), ## dcvd_ary3
    Field(:rbbs_ary3,                    Array{Float_T,3}), ## calibrated_data
    # cq_method,                           Symbol
    ## for ct method
    Field(:ct_fluos,                     SVector{C,Float_T} where {C}),
    ## computed by fit_model!()
    Field(:bl_fit,                       SMatrix{W,C,Union{AmpModelFit,Symbol}} where {W,C}),
    Field(:bl_notes,                     SMatrix{W,C,Vector{String}} where {W,C}),
    Field(:blsub_fluos,                  Array{Float_T,3}),
    Field(:blsub_fluos_flb,                  Array{Float_T,3}),
    Field(:quant_fit,                    SMatrix{W,C,Union{AmpModelFit,Symbol}} where {W,C}),
    Field(:quant_status,                 SMatrix{W,C,Symbol} where {W,C}),
    Field(:coefs,                        Array{Float_T,3}),
    Field(:d0,                           SMatrix{W,C,Float_T} where {W,C}),
    Field(:quant_fluos,                  Array{Float_T,3}),
    Field(:dr1_pred,                     Array{Float_T,3}),
    Field(:dr2_pred,                     Array{Float_T,3}),
    Field(:dr1_pred1,                     Array{Float_T,3}),
    Field(:dr2_pred1,                     Array{Float_T,3}),
    Field(:max_dr1,                      SMatrix{W,C,Float_T} where {W,C}),
    Field(:max_dr2,                      SMatrix{W,C,Float_T} where {W,C}),
    Field(:cyc_vals_4cq,                 SMatrix{W,C,OrderedDict{Symbol,Float_T}} where {W,C}),
    Field(:eff_vals_4cq,                 SMatrix{W,C,OrderedDict{Symbol,Float_T}} where {W,C}),
    Field(:cq_raw,                       SMatrix{W,C,Float_T} where {W,C}),
    Field(:cq,                           SMatrix{W,C,Float_T} where {W,C}),
    Field(:eff,                          SMatrix{W,C,Float_T} where {W,C}),
    Field(:cq_fluo,                      SMatrix{W,C,Float_T} where {W,C}),
    ## computed in amp_analysis()
    Field(:qt_fluos,                     SMatrix{W,C,Float_T} where {W,C}),
    Field(:max_qt_fluo,                  Float_T),
    ## computed by report_cq!()
    Field(:max_bsf,                      SMatrix{W,C,Float_T} where {W,C}),
    Field(:scaled_max_bsf,               SMatrix{W,C,Float_T} where {W,C}),
    Field(:scaled_max_dr1,               SMatrix{W,C,Float_T} where {W,C}),
    Field(:scaled_max_dr2,               SMatrix{W,C,Float_T} where {W,C}),
    Field(:why_NaN,                      SMatrix{W,C,String} where {W,C}),
    ## allelic discrimination output
    Field(:assignments_adj_labels_dict,  OrderedDict{Symbol,Vector{String}}),
    Field(:agr_dict,                     OrderedDict{Symbol,AssignGenosResult})]

const AMPSHORTOUTPUT_FIELDNAMES = [
    :rbbs_ary3,
    :blsub_fluos,
    :blsub_fluos_flb,
    :dr1_pred,
    :dr2_pred,
    :dr1_pred1,
    :dr2_pred1,
    :quant_status,
    :coefs,
    :cq,
    :d0,
    :ct_fluos,
    :assignments_adj_labels_dict]

#===============================================================================
    struct generation >>
===============================================================================#

abstract type AmpOutput end

SCHEMA = AMPLONGOUTPUT_FIELD_DEFS
@make_struct_from_SCHEMA AmpLongOutput AmpOutput true
@make_constructor_from_SCHEMA AmpLongOutput

SCHEMA = subset_schema(
    AMPLONGOUTPUT_FIELD_DEFS,
    AMPSHORTOUTPUT_FIELDNAMES)
@make_struct_from_SCHEMA AmpShortOutput AmpOutput true
@make_constructor_from_SCHEMA AmpShortOutput

#===============================================================================
    constructors >>
===============================================================================#

function AmpOutput(
    ::Type{Val{long}},
    i                           ::AmpInput,
    background_subtracted_data  ::Array{<: Real,3},
    k_deconv                    ::DeconvolutionMatrices,
    deconvoluted_data           ::Array{<: Real,3},
    norm_data                   ::SArray{S,R} where {S <: Tuple, R <: Real},
    norm_wells                  ::AbstractVector{Symbol},
    calibrated_data             ::Array{<: Real,3},
    ct_fluos                    ::AbstractVector;
)
    NaN_array2        = amp_init(i, NaN_T)
    # fitted_init       = amp_init(i, FIT[i.amp_model]())
    empty_vals_4cq    = amp_init(i, OrderedDict{Symbol, Float_T}())
    AmpLongOutput(
        i.raw_data, ## formerly fr_ary3
        background_subtracted_data, ## formerly mw_ary3
        norm_data, ## formerly wva_data
        k_deconv, ## formerly k4dcv
        deconvoluted_data, ## formerly dcvd_ary3
        calibrated_data, ## formerly rbbs_ary3
        # cq_method,
        ## ct_fluos
        SVector{i.num_channels, Float_T}(ct_fluos), ## ct_fluos
        ## model fit
        fitted_init, ## bl_fit,
        amp_init(i, Vector{String}()), ## bl_notes
        calibrated_data, ## blsub_fluos
        calibrated_data,
        fitted_init, ## quant_fit,
        amp_init(i, :not_fitted), ## quant_status
        fill(NaN_T, 1, i.num_wells, i.num_channels), ## coefs ## NB size = 1 for 1st dimension may not be correct for the chosen model
        NaN_array2, ## d0
        calibrated_data, ## quant_fluos,
        calibrated_data, ## dr1_pred
        calibrated_data, ## dr2_pred
        calibrated_data, ## dr1_pred
        calibrated_data, ## dr2_pred
        NaN_array2, ## max_dr1
        NaN_array2, ## max_dr2
        empty_vals_4cq, ## cyc_vals_4cq
        empty_vals_4cq, ## eff_vals_4cq
        NaN_array2, ## cq_raw
        NaN_array2, ## cq
        NaN_array2, ## eff
        NaN_array2, ## cq_fluo
        ## set_qt_fluos!()
        NaN_array2, ## qt_fluos
        NaN_array2, ## max_qt_fluo
        ## report_cq!()
        NaN_array2, ## max_bsf
        NaN_array2, ## scaled_max_bsf
        NaN_array2, ## scaled_max_dr1
        NaN_array2, ## scaled_max_dr2
        amp_init(i, ""), ## why_NaN
        ## allelic discrimination
        OrderedDict{Symbol, Vector{String}}(), ## assignments_adj_labels_dict
        OrderedDict{Symbol, AssignGenosResult}()) ## agr_dict
end ## constructor

function AmpOutput(
    ::Type{Val{short}},
    i                           ::AmpInput,
    background_subtracted_data  ::Array{<: Real,3},
    k_deconv                    ::DeconvolutionMatrices,
    deconvoluted_data           ::Array{<: Real,3},
    norm_data                   ::SArray{S,R} where {S <: Tuple, R <: Real},
    norm_well_nums              ::AbstractVector,
    calibrated_data             ::Array{<: Real,3},
    ct_fluos                    ::AbstractVector;
    reporting                   ::Function = roundoff(JSON_DIGITS) ## reporting function
)
    NaN_array2 = amp_init(i, NaN_T)
    cd = reporting.(calibrated_data)
    AmpShortOutput(
        cd, ## formerly rbbs_ary3
        cd, ## blsub_fluos
        cd,
        cd, ## dr1_pred
        cd, ## dr2_pred
        cd, ## dr1_pred
        cd, ## dr2_pred
        amp_init(i, :not_fitted), ## quant_status
        fill(NaN_T, 1, i.num_wells, i.num_channels),
        NaN_array2, ## d0
        NaN_array2, ## cq
        ct_fluos, ## ct_fluos
        OrderedDict{Symbol, AssignGenosResult}()) ## agr_dict
end ## constructor

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
(0 != incx && 0 != incy && 0 != incz) || throw(ArgumentError("zero increment"))
(0 < ix && 0 < iy && 0 < iz) || throw(BoundsError())
ix+(n-1)*abs(incx) <= length(x) || throw(BoundsError())
iy+(n-1)*abs(incy) <= length(y) || throw(BoundsError())
iz+(n-1)*abs(incz) <= length(z) || throw(BoundsError())


text:
include("./staterequest.jl")
include("./static export.jl")
include("./connections.jl")


text:
#!/usr/bin/env julia

using Weber

version = v"0.0.3"
sid,trial_skip,adapt = @read_args("Frequency Discrimination ($version).",
                                  adapt=[:levitt,:bayes])

const atten_dB = 30
const n_trials = 60
const feedback_delay = 750ms

isresponse(e) = iskeydown(e,key"p") || iskeydown(e,key"q")

const standard_freq = 1kHz
const standard = @> tone(standard_freq,100ms) ramp attenuate(atten_dB)

function one_trial(adapter)
  first_lower = rand(Bool)
  resp = response(adapter,key"q" => "first_lower",key"p" => "second_lower",
                  correct=(first_lower ? "first_lower" : "second_lower"))

  signal() = @> tone((1-delta(adapter))*standard_freq,100ms) begin
    ramp
    attenuate(atten_dB)
  end
  stimuli = first_lower? [signal,standard] : [standard,signal]

  [moment(feedback_delay,play,stimuli[1]),
   show_cross(),
   moment(900ms,play,stimuli[2]),
   moment(100ms + 300ms,display,
          "Was the first [Q] or second sound [P] lower in pitch?"),
   resp,await_response(isresponse)]
end

experiment = Experiment(
  skip=trial_skip,
  columns = [
    :sid => sid,
    :condition => "example",
    :version => version,
    :standard => standard_freq
  ]
)

setup(experiment) do
  addbreak(moment(record,"start"))

  addbreak(instruct("""

    On each trial, you will hear two beeps. Indicate which of the two beeps you
heard was lower in pitch. Hit 'Q' if the first beep was lower, and 'P' if the
second beep was lower.
"""))

  if adapt == :levitt
    adapter = levitt_adapter(down=3,up=1,min_delta=0,max_delta=1,
                             big=2,little=sqrt(2),mult=true)
  else
    adapter = bayesian_adapter(min_delta = 0,max_delta = 0.95)
  end

  @addtrials let a = adapter
    for trial in 1:n_trials
      addtrial(one_trial(a))
    end

    function threshold_report()
      mean,sd = estimate(adapter)
      thresh = round(mean,3)*standard_freq
      thresh_sd = round(sd,3)*standard_freq

      # define this string during run time when we know
      # what the threshold estimate is.
      "Threshold $(thresh)Hz (SD: $thresh_sd)\n"*
      "Hit spacebar to continue..."
    end

    addbreak(moment(display,threshold_report,clean_whitespace=false),
             await_response(iskeydown(key":space:")))
  end

end

run(experiment)


text:
using PowerDynamics
using NetworkDynamics
using LaTeXStrings
include("plotting.jl")
include("ThirdOrderEq.jl")
include("ThirdOrderEqGovernorIEEEG1.jl")
include("ThirdOrderEqExciterIEEEDC1A.jl")

line_list = []
append!(line_list, [StaticLine(from = 1, to = 2, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 2, to = 3, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 1, to = 3, Y = -1im / 0.02)])

node_list = []
append!(node_list, [SlackAlgebraic(U = 1, Y_n = 1)])
append!(node_list, [ThirdOrderEq(H = 5.0, P = 1.0, D = 0.1, Ω = 50, E_f = 0.2, X_d_dash = 0.111, T_d_dash = 0.1, X_q_dash = 0.103, X_d = 0.1)])
append!(node_list, [ThirdOrderEq(H = 5.0, P = 1.0, D = 0.1, Ω = 50, E_f = 0.2,X_d_dash = 0.111, T_d_dash = 0.1, X_q_dash = 0.103, X_d = 0.1)])

powergrid = PowerGrid(node_list, line_list)

operationpoint = find_operationpoint(powergrid)

result = simulate(
    Perturbation(2, :ω, Inc(0.1)),
    powergrid,
    operationpoint,
    timespan = (0.0, 10),
)

plot_res(result, powergrid, 2)

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
struct Normal <: Distribution{Float64} end
struct BroadcastedNormal <: Distribution{Array{Float64}} end

function broadcast_shapes_or_crash(args...)
    # This expression produces a DimensionMismatch exception if the shapes are
    # incompatible.  There seems to currently be no ready-made exception-free
    # way to check this, and I think it makes more sense to wait until that
    # capability exists in Julia proper than to re-implement it here.
    s = Base.Broadcast.broadcast_shape(map(size, args)...)
end

function assert_has_shape(x, expected_shape; msg="Shape assertion failed")
    if size(x) != expected_shape
        throw(DimensionMismatch(string(msg,
                                       " Expected shape: $expected_shape",
                                       " Actual shape: $(size(x))")))
    end
    nothing
end

"""
    normal(mu::Real, std::Real)

Samples a `Float64` value from a normal distribution.
"""
const normal = Normal()

"""
    broadcasted_normal(mu::AbstractArray{<:Real, N1},
                       std::AbstractArray{<:Real, N2}) where {N1, N2}

Samples an `Array{Float64, max(N1, N2)}` of shape
`Broadcast.broadcast_shapes(size(mu), size(std))` where each element is
independently normally distributed.  This is equivalent to (a reshape of) a
multivariate normal with diagonal covariance matrix, but its implementation is
more efficient than that of the more general `mvnormal` for this case.

The shapes of `mu` and `std` must be broadcast-compatible.

If all args are 0-dimensional arrays, then sampling via
`broadcasted_normal(...)` returns a `Float64` rather than properly returning an
`Array{Float64, 0}`.  This is consistent with Julia's own inconsistency on the
matter:

```jldoctest
julia> typeof(ones())
Array{Float64,0}

julia> typeof(ones() .* ones())
Float64
```
"""
const broadcasted_normal = BroadcastedNormal()

function logpdf(::Normal, x::Real, mu::Real, std::Real)
    var = std * std
    diff = x - mu
    -(diff * diff)/ (2.0 * var) - 0.5 * log(2.0 * pi * var)
end

function logpdf(::BroadcastedNormal,
                x::Union{AbstractArray{<:Real}, Real},
                mu::Union{AbstractArray{<:Real}, Real},
                std::Union{AbstractArray{<:Real}, Real})
    assert_has_shape(x, broadcast_shapes_or_crash(mu, std);
                     msg="Shape of `x` does not agree with the sample space")
    var = std .* std
    diff = x .- mu
    sum(-(diff .* diff) ./ (2.0 * var) .- 0.5 * log.(2.0 * pi * var))
end

function logpdf_grad(::Normal, x::Real, mu::Real, std::Real)
    precision = 1. / (std * std)
    diff = mu - x
    deriv_x = diff * precision
    deriv_mu = -deriv_x
    deriv_std = -1. / std + (diff * diff) / (std * std * std)
    (deriv_x, deriv_mu, deriv_std)
end

function logpdf_grad(::BroadcastedNormal,
                     x::Union{AbstractArray{<:Real}, Real},
                     mu::Union{AbstractArray{<:Real}, Real},
                     std::Union{AbstractArray{<:Real}, Real})
    assert_has_shape(x, broadcast_shapes_or_crash(mu, std);
                     msg="Shape of `x` does not agree with the sample space")
    precision = 1.0 ./ (std .* std)
    diff = mu .- x
    deriv_x = sum(diff .* precision)
    deriv_mu = sum(-deriv_x)
    deriv_std = sum(-1.0 ./ std .+ (diff .* diff) ./ (std .* std .* std))
    (deriv_x, deriv_mu, deriv_std)
end

random(::Normal, mu::Real, std::Real) = mu + std * randn()
is_discrete(::Normal) = false

function random(::BroadcastedNormal,
                mu::Union{AbstractArray{<:Real}, Real},
                std::Union{AbstractArray{<:Real}, Real})
    broadcast_shape = broadcast_shapes_or_crash(mu, std)
    mu .+ std .* randn(broadcast_shape)
end

(::Normal)(mu, std) = random(Normal(), mu, std)
(::BroadcastedNormal)(mu, std) = random(BroadcastedNormal(), mu, std)

has_output_grad(::Normal) = true
has_argument_grads(::Normal) = (true, true)

has_output_grad(::BroadcastedNormal) = true
has_argument_grads(::BroadcastedNormal) = (true, true)

export normal
export broadcasted_normal


text:
using Sockets

function matgen(n, seed)
    tmp = seed / n / n
    [tmp * (i - j) * (i + j - 2) for i = 1:n, j = 1:n]
end

function mul(a, b)
    c = zeros(size(a, 1), size(b, 2))
    M, N = size(c)
    K = size(b, 1)
    @inbounds Threads.@threads for n ∈ 1:N
        for k ∈ 1:K
            @simd ivdep for m ∈ 1:M
                @fastmath c[m, n] += a[m, k] * b[k, n]
            end
        end
    end
    return c
end

function calc(n)
    n = round(Int, n / 2) * 2
    a = matgen(n, 1.0)
    b = matgen(n, 2.0)
    c = mul(a, b)
    c[Int(n / 2)+1, Int(n / 2)+1]
end

function notify(msg)
    try
        socket = connect("localhost", 9001)
        write(socket, msg)
        close(socket)
    catch
        # standalone usage
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    n = length(ARGS) > 0 ? parse(Int, ARGS[1]) : 100

    left = calc(101)
    right = -18.67
    if abs(left - right) > 0.1
        println(stderr, "$(left) != $(right)")
        exit(1)
    end

    notify("Julia (no BLAS)\t$(getpid())")
    t = time()
    results = calc(n)
    elapsed = time() - t
    notify("stop")

    println(results)
    println("time: $(elapsed) s")
end

max_stars_repo_name: clemenswasser/benchmarks

max_issues_repo_name: clemenswasser/benchmarks

max_forks_repo_name: clemenswasser/benchmarks


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

include("model_ConnectivityDestination.jl")
include("model_ConnectivityHop.jl")
include("model_ConnectivityInformation.jl")
include("model_ConnectivityIssue.jl")
include("model_ConnectivityParameters.jl")
include("model_ConnectivitySource.jl")
include("model_EffectiveNetworkSecurityRule.jl")
include("model_FlowLogInformation.jl")
include("model_FlowLogProperties.jl")
include("model_FlowLogStatusParameters.jl")
include("model_IssueContext.jl")
include("model_NetworkInterfaceAssociation.jl")
include("model_NetworkWatcherListResult.jl")
include("model_NetworkWatcherPropertiesFormat.jl")
include("model_NextHopParameters.jl")
include("model_NextHopResult.jl")
include("model_PacketCapture.jl")
include("model_PacketCaptureFilter.jl")
include("model_PacketCaptureListResult.jl")
include("model_PacketCaptureParameters.jl")
include("model_PacketCaptureQueryStatusResult.jl")
include("model_PacketCaptureResult.jl")
include("model_PacketCaptureStorageLocation.jl")
include("model_QueryTroubleshootingParameters.jl")
include("model_Resource.jl")
include("model_RetentionPolicyParameters.jl")
include("model_SecurityGroupNetworkInterface.jl")
include("model_SecurityGroupViewParameters.jl")
include("model_SecurityGroupViewResult.jl")
include("model_SecurityRuleAssociations.jl")
include("model_SubResource.jl")
include("model_SubnetAssociation.jl")
include("model_Topology.jl")
include("model_TopologyAssociation.jl")
include("model_TopologyParameters.jl")
include("model_TopologyResource.jl")
include("model_TroubleshootingDetails.jl")
include("model_TroubleshootingParameters.jl")
include("model_TroubleshootingProperties.jl")
include("model_TroubleshootingRecommendedActions.jl")
include("model_TroubleshootingResult.jl")
include("model_VerificationIPFlowParameters.jl")
include("model_VerificationIPFlowResult.jl")
include("model_NetworkWatcher.jl")
include("model_PacketCaptureResultProperties.jl")
include("model_SecurityRule.jl")


text:
using LinearAlgebra

A = [
   1 2
   2 1
]

e = eigen(A)
V = e.vectors
@show Diagonal(e.values)
V⁻¹ = inv(V)
@show V⁻¹ * A * V

F = rand(3,5)
s = svd(F)

F |> display
s.U * Diagonal(s.S) * s.Vt |> display

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
struct CouetteFlow{T <: Real, Int <: Integer} <: FluidFlowProblem
    rho_0::T
    u_max::T
    ν::T
    NX::Int
    NY::Int
    domain_size::Tuple{T, T}
end

function CouetteFlow(
    ν = 1.0 / 6.0,
    scale = 2,
    NX = 5 * scale + 0,
    NY = NX,
    domain_size = (1.0, 1.0),
)
    # # Since our velocity profile is constant over x it is sufficient to only take
    # # 3 nodes in the x direction
    NX = 1
    u_max = 0.015 / scale
    u_max = 0.01 / scale
    # u_max = sqrt(0.001) / scale

    return CouetteFlow(1.0, u_max, ν, NX, NY, domain_size)
end

function density(
    q::Quadrature,
    problem::CouetteFlow,
    x::T,
    y::T,
    timestep::Real = 0.0,
) where {T <: Real}
    return 1.0
end

function pressure(
    q::Quadrature,
    problem::CouetteFlow,
    x::T,
    y::T,
    timestep::Real = 0.0,
) where {T <: Real}
    return 1.0
end

function velocity(problem::CouetteFlow, x::T, y::T, timestep::Real = 0.0) where {T <: Real}
    return [
        y
        0.0
    ]
end
function velocity_gradient(
    problem::CouetteFlow,
    x::T,
    y::T,
    timestep::Real = 0.0,
) where {T <: Real}
    u_x = 0.0
    v_y = 0.0
    u_y = 1.0
    v_x = 0.0

    return [u_x v_x; u_y v_y]
end
function force(
    problem::CouetteFlow,
    x::Int,
    y::Int,
    time::Real = 0.0,
) where {Int <: Integer}
    return [
        0.0
        0.0
    ]
end

"""
Apply a bounce back to the bottom of the domain and a moving wall to the
top of the domain
"""
boundary_conditions(problem::CouetteFlow) = [
    BounceBack(South(), 1:(problem.NX), 1:(problem.NY)),
    MovingWall(North(), 1:(problem.NX), 1:(problem.NY), [problem.u_max, 0]),
]


text:
export ess, actime

# Effective sample size (ESS)
actypes = (:bm, :imse, :ipse)

function ess(c::MCMCChain, pars::Ranges=1:ncol(c.samples); vtype::Symbol=:imse, args...)
  @assert in(vtype, actypes) "Unknown ESS type $vtype"

  return nrow(c.samples)*var(c, pars; vtype=:iid)./var(c, pars; vtype=vtype, args...)
end

ess(c::MCMCChain, par::Real; vtype::Symbol=:imse, args...) = ess(c, par:par; vtype=vtype, args...)

# Integrated autocorrelation time
function actime(c::MCMCChain, pars::Ranges=1:ncol(c.samples); vtype::Symbol=:imse, args...)
  @assert in(vtype, actypes) "Unknown integrated autocorrelation time type $vtype"

  return var(c, pars; vtype=vtype, args...)./var(c, pars; vtype=:iid)
end

actime(c::MCMCChain, par::Real; vtype::Symbol=:imse, args...) = actime(c, par:par; vtype=vtype, args...)


text:
module Samantha

using Parameters
using OnlineStats
using Random
import Base: ==, hash
import UUIDs: UUID, uuid4

# Core files
include("imports-exports.jl")
include("abstracts.jl")
include("util.jl")
include("container.jl")
include("defaults.jl")
include("neurons/neurons.jl")
include("synapses/synapses.jl")
include("layer.jl")
include("patchclamp.jl")
include("agent.jl")
include("interface.jl")

# Standard Library
# TODO: include("Stdlib.jl")

# External Optional Dependencies
# TODO: include("external/External.jl")

end