text:
""" Data structure for undirected wiring diagrams.
"""
module UndirectedWiringDiagrams
export UndirectedWiringDiagram, UntypedUWD, TypedUWD,
  outer_box, box, junction, nboxes, njunctions, boxes, junctions,
  ports, ports_with_junction, junction_type, port_type,
  add_box!, add_junction!, add_junctions!, set_junction!, add_wire!,
  add_wires!, singleton_diagram, cospan_diagram, junction_diagram,
  ocompose, substitute

using ...Present, ...CategoricalAlgebra.CSets, ...CategoricalAlgebra.Limits
using ...CategoricalAlgebra.FinSets: FinSet, FinFunction
using ...Theories: dom, codom, compose, ⋅, id, oplus
import ..DirectedWiringDiagrams: box, boxes, nboxes, add_box!, add_wire!,
  add_wires!, singleton_diagram, ocompose, substitute

# Data types
############

@present TheoryUWD(FreeSchema) begin
  Box::Ob
  Port::Ob
  OuterPort::Ob
  Junction::Ob

  box::Hom(Port, Box)
  junction::Hom(Port, Junction)
  outer_junction::Hom(OuterPort, Junction)
end

const UndirectedWiringDiagram = const AbstractUWD = AbstractACSetType(TheoryUWD)
const UntypedUWD = CSetType(TheoryUWD, index=[:box, :junction, :outer_junction])

@present TheoryTypedUWD <: TheoryUWD begin
  Type::Data

  port_type::Attr(Port, Type)
  outer_port_type::Attr(OuterPort, Type)
  junction_type::Attr(Junction, Type)

  compose(junction, junction_type) == port_type
  compose(outer_junction, junction_type) == outer_port_type
end

const AbstractTypedUWD = AbstractACSetType(TheoryTypedUWD)
const TypedUWD = ACSetType(TheoryTypedUWD,
                           index=[:box, :junction, :outer_junction])

function (::Type{UWD})(nports::Int) where UWD <: AbstractUWD
  d = UWD()
  add_parts!(d, :OuterPort, nports)
  return d
end

function (::Type{UWD})(port_types::AbstractVector) where UWD <: AbstractUWD
  d = UWD()
  nports = length(port_types)
  add_parts!(d, :OuterPort, nports, outer_port_type=port_types)
  return d
end

UndirectedWiringDiagram(nports::Int) = UntypedUWD(nports)
UndirectedWiringDiagram(port_types::AbstractVector{T}) where T =
  TypedUWD{T}(port_types)

# Imperative interface
######################

# Accessors
#----------

outer_box(::AbstractACSet) = 0
box(d::AbstractACSet, args...) = subpart(d, args..., :box)
junction(d::AbstractACSet, args...; outer::Bool=false) =
  subpart(d, args..., outer ? :outer_junction : :junction)

function junction(d::AbstractACSet, port::Tuple{Int,Int})
  box, nport = port
  box == outer_box(d) ?
    junction(d, nport, outer=true) : junction(d, ports(d, box)[nport])
end

nboxes(d::AbstractACSet) = nparts(d, :Box)
njunctions(d::AbstractACSet) = nparts(d, :Junction)
boxes(d::AbstractACSet) = parts(d, :Box)
junctions(d::AbstractACSet) = parts(d, :Junction)

ports(d::AbstractACSet; outer::Bool=false) =
  parts(d, outer ? :OuterPort : :Port)
ports(d::AbstractACSet, box) =
  box == outer_box(d) ? parts(d, :OuterPort) : incident(d, box, :box)
ports_with_junction(d::AbstractACSet, junction; outer::Bool=false) =
  incident(d, junction, outer ? :outer_junction : :junction)

junction_type(d::AbstractACSet, args...) = subpart(d, args..., :junction_type)
port_type(d::AbstractACSet, args...; outer::Bool=false) =
  subpart(d, args..., outer ? :outer_port_type : :port_type)

function port_type(d::AbstractACSet, port::Tuple{Int,Int})
  box, nport = port
  box == outer_box(d) ?
    port_type(d, nport, outer=true) : port_type(d, ports(d, box)[nport])
end

# Mutators
#---------

add_box!(d::AbstractUWD; data...) = add_part!(d, :Box; data...)

function add_box!(d::AbstractUWD, nports::Int; data...)
  box = add_box!(d; data...)
  ports = add_parts!(d, :Port, nports, box=fill(box, nports))
  box
end

function add_box!(d::AbstractUWD, port_types::AbstractVector; data...)
  box = add_box!(d; data...)
  nports = length(port_types)
  ports = add_parts!(d, :Port, nports, box=fill(box, nports),
                     port_type=port_types)
  box
end

add_junction!(d::AbstractUWD; data...) = add_part!(d, :Junction; data...)
add_junction!(d::AbstractUWD, type; data...) =
  add_part!(d, :Junction; junction_type=type, data...)

add_junctions!(d::AbstractUWD, njunctions::Int; data...) =
  add_parts!(d, :Junction, njunctions; data...)
add_junctions!(d::AbstractUWD, types::AbstractVector; data...) =
  add_parts!(d, :Junction, length(types); junction_type=types, data...)

function set_junction!(d::AbstractUWD, port, junction; outer::Bool=false)
  if has_subpart(d, :junction_type)
    ptype, jtype = port_type(d, port, outer=outer), junction_type(d, junction)
    all(ptype .== jtype) || error(
      "Domain error: port type $ptype and junction type $jtype do not match")
  end
  set_subpart!(d, port, outer ? :outer_junction : :junction, junction)
end
set_junction!(d::AbstractUWD, junction; kw...) =
  set_junction!(d, :, junction; kw...)

function set_junction!(d::AbstractUWD, port::Tuple{Int,Int}, junction)
  box, nport = port
  if box == outer_box(d)
    set_junction!(d, nport, junction, outer=true)
  else
    set_junction!(d, ports(d, box)[nport], junction)
  end
end

""" Wire together two ports in an undirected wiring diagram.

A convenience method that creates and sets junctions as needed. Ports are only
allowed to have one junction, so if both ports already have junctions, then the
second port is assigned the junction of the first. The handling of the two
arguments is otherwise symmetric.

FIXME: When both ports already have junctions, the two junctions should be
*merged*. To do this, we must implement `merge_junctions!` and thus also
`rem_part!`.
"""
function add_wire!(d::AbstractUWD, port1::Tuple{Int,Int}, port2::Tuple{Int,Int})
  j1, j2 = junction(d, port1), junction(d, port2)
  if j1 > 0
    set_junction!(d, port2, j1)
  elseif j2 > 0
    set_junction!(d, port1, j2)
  else
    j = has_subpart(d, :junction_type) ?
      add_junction!(d, port_type(d, port1)) : add_junction!(d)
    set_junction!(d, port1, j)
    set_junction!(d, port2, j)
  end
end
add_wire!(d, wire::Pair) = add_wire!(d, first(wire), last(wire))

function add_wires!(d::AbstractUWD, wires)
  for wire in wires
    add_wire!(d, wire)
  end
end

# Other constructors
#-------------------

function singleton_diagram(::Type{UWD}, port_types;
                           data...) where UWD <: AbstractUWD
  d = UWD(port_types)
  add_box!(d, port_types; data...)
  js = add_junctions!(d, port_types)
  set_junction!(d, js)
  set_junction!(d, js, outer=true)
  return d
end

""" Undirected wiring diagram defined by a cospan.

The wiring diagram has a single box. The ports of this box, the outer ports, the
junctions, and the connections between them are defined by the cospan. Thus,
this function generalizes [`singleton_diagram`](@ref).

See also: [`junction_diagram`](@ref).
"""
function cospan_diagram(::Type{UWD}, inner::FinFunction{Int,Int},
                        outer::FinFunction{Int,Int}, junction_types=nothing;
                        data...) where UWD <: AbstractUWD
  @assert codom(inner) == codom(outer)
  if isnothing(junction_types)
    junction_types = length(codom(inner))
  end
  cospan_diagram(UWD, collect(inner), collect(outer), junction_types; data...)
end
function cospan_diagram(::Type{UWD}, inner::AbstractVector{Int},
                        outer::AbstractVector{Int}, junction_types;
                        data...) where UWD <: AbstractUWD
  d = UWD(map_port_types(junction_types, outer))
  add_box!(d, map_port_types(junction_types, inner); data...)
  add_junctions!(d, junction_types)
  set_junction!(d, inner)
  set_junction!(d, outer, outer=true)
  return d
end

""" Undirected wiring diagram with no boxes, only junctions.

See also: [`singleton_diagram`](@ref), [`cospan_diagram`](@ref).
"""
function junction_diagram(::Type{UWD}, outer::FinFunction{Int,Int},
                          junction_types=nothing) where UWD <: AbstractUWD
  if isnothing(junction_types)
    junction_types = length(codom(outer))
  end
  junction_diagram(UWD, collect(outer), junction_types)
end
function junction_diagram(::Type{UWD}, outer::AbstractVector{Int},
                          junction_types) where UWD <: AbstractUWD
  d = UWD(map_port_types(junction_types, outer))
  add_junctions!(d, junction_types)
  set_junction!(d, outer, outer=true)
  return d
end

map_port_types(::Int, g) = length(g)
map_port_types(types::AbstractVector, g) = types[g]

# Operadic interface
####################

function ocompose(f::AbstractUWD, gs::AbstractVector{<:AbstractUWD})
  substitute(f, boxes(f), gs)
end
function ocompose(f::AbstractUWD, i::Int, g::AbstractUWD)
  substitute(f, [i], [g])
end

function substitute(f::UWD, indices::AbstractVector{Int},
                    gs::AbstractVector{UWD}) where {UWD <: AbstractUWD}
  @assert length(indices) == length(gs)
  h = UWD()
  copy_parts!(h, f, OuterPort=:)

  # Copy boxes from original diagram and diagrams to be substituted.
  gs_index = Vector{Union{Int,Nothing}}(nothing, nboxes(f))
  gs_index[indices] = 1:length(indices)
  for (i,j) in enumerate(gs_index)
    if isnothing(j)
      copy_parts!(h, f, Box=[i], Port=ports(f, i))
    else
      g = gs[j]
      copy_parts!(h, g, Box=boxes(g), Port=flat(ports(g, boxes(g))))
    end
  end

  # Add junctions obtained from pushout of attributed junction sets.
  Junctions = ACSetTableType(UWD, :Junction)
  f_junctions, g_junctions, apex = Junctions(), Junctions(), Junctions()
  copy_parts_only!(f_junctions, f)
  for g in gs
    copy_parts_only!(g_junctions, g)
    copy_parts_only!(apex, g, Junction=junction(g, outer=true))
  end
  f_inner = ACSetTransformation(apex, f_junctions, Junction=flat(
    junction(f, ports(f, i)) for (i,j) in enumerate(gs_index) if !isnothing(j)))
  g_outer = ACSetTransformation(apex, g_junctions, Junction=oplus([
    FinFunction(junction(g, outer=true), njunctions(g)) for g in gs ]))
  colim = pushout(f_inner, g_outer)
  copy_parts_only!(h, ob(colim))

  # Assign junctions to all ports.
  f_leg, g_leg = (leg[:Junction] for leg in colim)
  gs_incl = legs(coproduct([ FinSet(njunctions(g)) for g in gs ]))
  set_junction!(h, map(f_leg, junction(f, outer=true)), outer=true)
  set_junction!(h, flatmap(enumerate(gs_index)) do (i,j)
    if isnothing(j)
      map(f_leg, junction(f, ports(f, i)))
    else
      map(gs_incl[j] ⋅ g_leg, junction(gs[j]))
    end
  end)
  return h
end

flat(x) = reduce(vcat, x, init=Int[])
flatmap(f, xs...) = mapreduce(f, vcat, xs..., init=Int[])

end


text:
# JOIN clause.

mutable struct JoinClause <: AbstractSQLClause
    over::Union{SQLClause, Nothing}
    joinee::SQLClause
    on::SQLClause
    left::Bool
    right::Bool
    lateral::Bool

    JoinClause(;
               over = nothing,
               joinee,
               on,
               left = false,
               right = false,
               lateral = false) =
        new(over, joinee, on, left, right, lateral)
end

JoinClause(joinee; over = nothing, on, left = false, right = false, lateral = false) =
    JoinClause(over = over, joinee = joinee, on = on, left = left, right = right, lateral = lateral)

JoinClause(joinee, on; over = nothing, left = false, right = false, lateral = false) =
    JoinClause(over = over, joinee = joinee, on = on, left = left, right = right, lateral = lateral)

"""
    JOIN(; over = nothing, joinee, on, left = false, right = false, lateral = false)
    JOIN(joinee; over = nothing, on, left = false, right = false, lateral = false)
    JOIN(joinee, on; over = nothing, left = false, right = false, lateral = false)

A `JOIN` clause.

# Examples

```jldoctest
julia> c = FROM(:p => :person) |>
           JOIN(:l => :location,
                on = OP("=", (:p, :location_id), (:l, :location_id)),
                left = true) |>
           SELECT((:p, :person_id), (:l, :state));

julia> print(render(c))
SELECT
  "p"."person_id",
  "l"."state"
FROM "person" AS "p"
LEFT JOIN "location" AS "l" ON ("p"."location_id" = "l"."location_id")
```
"""
JOIN(args...; kws...) =
    JoinClause(args...; kws...) |> SQLClause

dissect(scr::Symbol, ::typeof(JOIN), pats::Vector{Any}) =
    dissect(scr, JoinClause, pats)

function PrettyPrinting.quoteof(c::JoinClause, ctx::QuoteContext)
    ex = Expr(:call, nameof(JOIN), quoteof([c.joinee, c.on], ctx)...)
    if c.left
        push!(ex.args, Expr(:kw, :left, c.left))
    end
    if c.right
        push!(ex.args, Expr(:kw, :right, c.right))
    end
    if c.lateral
        push!(ex.args, Expr(:kw, :lateral, c.lateral))
    end
    if c.over !== nothing
        ex = Expr(:call, :|>, quoteof(c.over, ctx), ex)
    end
    ex
end

rebase(c::JoinClause, c′) =
    JoinClause(over = rebase(c.over, c′),
               joinee = c.joinee, on = c.on, left = c.left, right = c.right, lateral = c.lateral)

ext: jl


text:
function getIndices()
    dates = Date(1920,1,1):Date(2015,9,25)

    ## measure time
    t0 = time()

    ## procIds = addprocs(3)

    ## @everywhere using DataFrames
    ## @everywhere using TimeData
    ## @everywhere using EconDatasets

indexSymb = ["^DJI",
"^DJT",
"^DJU",
"^DJA",
"^TV.N",
"^NYA",
"^NUS",
"^NIN",
"^NWL",
"^NTM",
"^IXIC",
"^TV.O",
"^NDX",
"^IXBK",
"^IXFN",
"^IXF",
"^IXID",
"^IXIS",
"^IXK",
"^IXTR",
"^IXUT",
"^NBI",
"^GSPC",
"^OEX",
"^MID",
"^SML",
"^SPSUPX",
"^XAX",
"^IIX",
"^NWX",
"^XMI",
"^PSE",
"^SOXX",
"^RUI",
"^RUT",
"^RUA",
"^DOT",
"^DWC",
"^BATSK",
"^TYX",
"^TNX",
"^FVX",
"^IRX",
"^DJC",
"^XAU",
"^NSEBANK",
"^CRSLDX",
"^CNXIT",
"^CRSMID",
"^NSEI",
"^NSMIDCP",
"^AORD",
"^SSEC",
"^HSI",
"^BSESN",
"^JKSE",
"^KLSE",
"^N225",
"^NZ50",
"^STI",
"^KS11",
"^TWII",
"^ATX",
"^BFX",
"^FCHI",
"^GDAXI",
"^OSEAX",
"^MIBTEL",
"^SSMI",
"^FTSE",
"FPXAA.PR",
"MICEXINDEXCF.ME",
"GD.AT",
"^CCSI",
"^TA100",
"^MERV",
"^BVSP",
"^GSPTSE",
"^MXX",
"^GSPC",
"BSE-100.BO",
"BSE-200.BO",
"BSE-500.BO",
"BSE-AUTO.BO",
"BSE-BANK.BO",
"BSE-CD.BO",
"BSE-FMCG.BO",
"BSE-HC.BO",
"BSE-IT.BO",
"BSE-METAL.BO",
"BSE-MIDCAP.BO",
"BSE-OILGAS.BO",
"^BSESN",
"BSE-SMLCAP.BO",
    "BSE-TECK.BO"]

## download indices
indexPrices = readYahooAdjClose(dates, indexSymb)

valsTn = convert(Timenum, indexPrices)

## delete columns without observations
##------------------------------------

onlyNAs = x -> all(isna(x))

## find columns without observations
onlyNANames = Symbol[]
for col in eachcol(valsTn.vals)
    if onlyNAs(col[2])
        push!(onlyNANames, col[1])
    end
end

## delete columns
for nam in onlyNANames
    delete!(valsTn.vals, nam)
end

## save result
##------------

pathToStore = joinpath(Pkg.dir("EconDatasets"), "data", "Indices.csv")
writeTimedata(pathToStore, valsTn)

## rmprocs(procIds)

t1 = time()
elapsedTime = t1-t0
mins, secs = divrem(elapsedTime, 60)

println("elapsed time: ", int(mins), " minutes, ", ceil(secs), "
seconds")

end


text:
# This file was generated, do not modify it. # hide
unwanted = [:peripheral_pulse, :nasogastric_tube, :nasogastric_reflux,
        :nasogastric_reflux_ph, :feces, :abdomen, :abdomcentesis_appearance, :abdomcentesis_total_protein]
select!(datac, Not(unwanted));


text:
module solve_dist
function solve()
    return 1;
end
end


text:
export TrivialMeasure

struct TrivialMeasure <: PrimitiveMeasure end

@primitive TrivialMeasure

sampletype(::TrivialMeasure) = Nothing


text:
# Control script for calculating CPCP.
#
# Hongyang Zhou, hyzhou@umich.edu 11/18/2019

include("fieldline.jl")
################################################################################
# Main

# Step 1: create seeding points
createSeeds("seeds.txt", "../data/")

# Step 2: field line tracing (in Paraview)
pvpython = `pvpython ../scripts/streamtrace_save.py`;
run(pvpython)

# Step 3: find magnetopause boundary
file_stream = "../data/streamline.txt";
x, y, z = findBoundary(file_stream, true);
plot(x,y);
saveBoundary(x, y, z, "../data/boundary.txt");

# Step 4: interpolate onto the boundary points (in Paraview)
pvpython = `pvpython ../scripts/Interpolate_On_Boundary_batch.py`;
run(pvpython)

# Step 5: 2D line integral
file_integration = "../data/boundary_value.txt";
ϕ = integrate_along_boundary(file_integration, true)

file_integration = "../data/boundary_value_hall.txt";
ϕ_hall= integrate_along_boundary_hall(file_integration, true)


text:
using DiffEqBase: set_t!, set_u!, set_ut!
using OrdinaryDiffEq

println("First")
# set_X!(integrator, integrator.X) should not change the result.
@testset "Trivial $setter ($alg, inplace=$iip)" for alg in [RK4, Trapezoid],
    setter in [set_t!, set_u!, set_ut!],
    iip in [false, true]
  if iip
    f = (du,u,p,t) -> (du .= 2u)
    prob = ODEProblem{iip}(f,[1/2],(0.0,2.0))
  else
    f = (u,p,t) -> (2u)
    prob = ODEProblem{iip}(f,1/2,(0.0,2.0))
  end
  t_half = 1.0

  integrator1 = init(prob, alg(); save_everystep=false)
  integrator2 = init(prob, alg(); save_everystep=false)

  step!(integrator1, t_half)  # "inexact" stepping w/o tstops
  if setter === set_t!
    set_t!(integrator1, integrator1.t)
  elseif setter === set_u!
    set_u!(integrator1, integrator1.u)
  else
    set_ut!(integrator1, integrator1.u, integrator1.t)
  end
  solve!(integrator1)

  solve!(integrator2)

  @test integrator1.t == integrator2.t
  if alg === Trapezoid
    rtol = integrator1.opts.reltol * 100
    atol = integrator1.opts.abstol
    @test integrator1.u ≈ integrator2.u  rtol=rtol atol=atol
  else
    @test integrator1.u == integrator2.u
  end
end

println("Second")

@testset "Resolve with $setter ($alg, inplace=$iip)" for alg in [RK4, Trapezoid],
    setter in [set_t!, set_ut!],
    iip in [false, true]
  if iip
    f = (du,u,p,t) -> (du .= 2u * cos(2π * t))
    prob1 = ODEProblem{iip}(f,[1/2],(0.0,1.0))
  else
    f = (u,p,t) -> (2u * cos(2π * t))
    prob1 = ODEProblem{iip}(f,1/2,(0.0,1.0))
  end
  prob2 = remake(prob1; tspan=(0.0,2.0))
  integrator1 = init(prob1, alg(); save_everystep=false)
  integrator2 = init(prob2, alg(); save_everystep=false)

  solve!(integrator1)
  if setter === set_t!
    set_t!(integrator1, 0)
  else
    set_ut!(integrator1, integrator1.u, 0)
  end
  solve!(integrator1)

  solve!(integrator2)

  rtol = integrator1.opts.reltol
  atol = integrator1.opts.abstol
  if alg === Trapezoid
    rtol *= 100
  end
  @test integrator1.u ≈ integrator2.u  rtol=rtol atol=atol
end

@testset "set_proposed_dt!" begin
  # 1) Setup problem
  f = (du,u,p,t) -> (du[1] = p)
  u0 = [1.0]
  tspan = (0.0,5.0)
  p = -1
  prob = ODEProblem(f,u0,tspan,p)

  # 2) Setup callback
  condition(u,t,integrator) = u[1] <= 0.0
  # function that sets dt that works with adaptive solvers:
  function affect!(integrator)
    integrator.p = 1
    set_proposed_dt!(integrator,0.5)
  end
  cb = DiscreteCallback(condition,affect!;save_positions=(true,true))

  sol = solve(prob,RK4(), dt = 0.1,adaptive = false, callback=cb)
  @test diff(sol.t)[1] == 0.1
  @test diff(sol.t)[13] == 0.5

  sol    = solve(prob,RK4(), dt = 0.1,adaptive = true, callback=cb)
  @test diff(sol.t)[1] == 0.1
  @test diff(sol.t)[4] == 0.5
end


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

include("common_io.jl")
using LinearAlgebra

struct Wishart
  gamma::Float64
  m::Int
end

function logsumexp(x)
  mx = maximum(x)
  log.(sum(exp.(x .- mx))) .+ mx
end

sumsq(v) = sum(abs2, v)

function ltri_unpack(D, LT)
  d=length(D)
  make_row(r::Int, L) = hcat(reshape([ L[i] for i=1:r-1 ],1,r-1), D[r], zeros(1,d-r))
  row_start(r::Int) = div((r-1)*(r-2),2)
  inds(r) = row_start(r) .+ (1:r-1)
  vcat([ make_row(r, LT[inds(r)]) for r=1:d ]...)
end

function get_Q(d,icf)
  ltri_unpack(exp.(icf[1:d]),icf[d+1:end])
end

# Gradient helpers
function pack(alphas,means,icf)
  [alphas[:];means[:];icf[:]]
end

function unpack(d,k,packed)
  alphas = reshape(packed[1:k],1,k)
  off = k
  means = reshape(packed[(1:d*k) .+ off],d,k)
  icf_sz = div(d*(d + 1),2)
  off += d*k
  icf = reshape(packed[off+1:end],icf_sz,k)
  (alphas,means,icf)
end

function log_gamma_distrib(a, p)
  out = 0.25 * p * (p - 1) * 1.1447298858494002 #convert(Float64, log(pi))
	for j in 1:p
    out += (logabsgamma(a + 0.5*(1 - j)))[1]
  end
	out
end

function log_wishart_prior(wishart::Wishart, sum_qs, Qs, icf)
  p = size(Qs[1],1)
  n = p + wishart.m + 1
  C = n*p*(log(wishart.gamma) - 0.5*log(2)) - log_gamma_distrib(0.5*n, p)

  frobenius = 0.
  for Q in Qs
    frobenius += sum(abs2,diag(Q))
  end
  frobenius += sum(abs2,icf[d+1:end,:])
	0.5*wishart.gamma^2 * frobenius - wishart.m*sum(sum_qs) - k*C
end

# IO
function read_gmm_instance(fn,replicate_point)
  fid = open(fn)
  lines = readlines(fid)
  close(fid)
  line=split(lines[1]," ")
  d = parse(Int,line[1])
  k = parse(Int,line[2])
  n = parse(Int,line[3])
  icf_sz = div(d*(d + 1),2)
  off = 1

  alphas = zeros(Float64,1,k)
  for i in 1:k
    alphas[i] = parse(Float64,lines[i+off])
  end
  off += k

  means = zeros(Float64,d,k)
  for ik in 1:k
    line=split(lines[ik+off]," ")
    for id in 1:d
      means[id,ik] = parse(Float64,line[id])
    end
  end
  off += k

  icf = zeros(Float64,icf_sz,k)
  for ik in 1:k
    line=split(lines[ik+off]," ")
    for i in 1:icf_sz
      icf[i,ik] = parse(Float64,line[i])
    end
  end
  off += k

  if replicate_point
    x_ = zeros(Float64,d,1)
    line=split(lines[1+off]," ")
    for id in 1:d
      x_[id] = parse(Float64,line[id])
    end
    x = repeat(x_,1,n)
    off += 1
  else
    x = zeros(Float64,d,n)
    for ix in 1:n
      line=split(lines[ix+off]," ")
      for id in 1:d
        x[id,ix] = parse(Float64,line[id])
      end
    end
    off += n
  end
  line = split(lines[1+off]," ")
  wishart = Wishart(parse(Float64,line[1]),parse(Int,line[2]))
  (alphas,means,icf,x,wishart)
end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using ContactImplicitMPC
using LinearAlgebra
using Quaternions

@show Threads.nthreads()

include("continuous_policy_v2.jl")

# ## Simulation
s = get_simulation("point_foot_quadruped", "flat_3D_lc", "flat")
model = s.model
env = s.env

# ## Reference Trajectory
ref_traj = deepcopy(get_trajectory(s.model, s.env,
	joinpath(module_dir(), "src/dynamics/point_foot_quadruped/gaits/inplace_trot_v4.jld2"),
    # joinpath(module_dir(), "src/dynamics/point_foot_quadruped/gaits/stand_euler_v0.jld2"),
    load_type = :split_traj_alt));

H = ref_traj.H
h = ref_traj.h

# ## MPC setup
N_sample = 5
H_mpc = 10
h_sim = h / N_sample
H_sim = 3000
κ_mpc = 2.0e-4

v0 = 0.0
obj = TrackingVelocityObjective(model, env, H_mpc,
    v = [Diagonal(1e-3 * [[1,1,1]; 1e+3*[1,1,1]; fill([1,1,1], 4)...]) for t = 1:H_mpc],
	q = [relative_state_cost(1e-0*[1e-2,1e-2,1], 3e-1*[1,1,1], 1e-0*[0.2,0.2,1]) for t = 1:H_mpc],
	u = [Diagonal(3e-3 * vcat(fill([1,1,1], 4)...)) for t = 1:H_mpc],
	v_target = [1/ref_traj.h * [v0;0;0; 0;0;0; v0;0;0; v0;0;0; v0;0;0; v0;0;0] for t = 1:H_mpc],)

p = ci_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
	mode = :configuration,
	ip_opts = InteriorPointOptions(
					undercut = 5.0,
					κ_tol = κ_mpc,
					r_tol = 1.0e-4, # TODO maybe relax this parameter
					diff_sol = true,
					solver = :empty_solver,
					max_time = 1e5),
    n_opts = NewtonOptions(
        r_tol = 3e-5,
        max_time=1.0e-2,
		solver=:ldl_solver,
        threads=false,
        verbose=false,
        max_iter = 5),
    mpc_opts = CIMPCOptions(
		# live_plotting=true
		));

# ## Disturbances
idx_d1 = 20
idx_d2 = idx_d1 + 350
idx_d3 = idx_d2 + 350
idx_d4 = idx_d3 + 350
idx_d5 = idx_d4 + 350
idx_d6 = idx_d5 + 350
idx = [idx_d1, idx_d2, idx_d3, idx_d4, idx_d5, idx_d6]
# impulses = [[0,0,10.0], [0,0,-10.0], [0,5,0.0], [0,-5,0.0], [5,0,0.0], [-5,0,0.0]]
impulses = [1.0 * [15.0,25.0,30.0], [0,0,0.0], [0,0,0.0], [0,0,0.0], [0,0,0.0], [0,0,0.0]]
d = impulse_disturbances(impulses, idx)

w = [[0.0,0.0,0.0] for i=1:H_sim/N_sample]
# w = [[0,0.07,0.0] for i=1:H_sim/N_sample]
# w = [[0,0.0,0.0] for i=1:H_sim/N_sample]
# d = open_loop_disturbances(w, N_sample)

# ## Initial conditions
q1_sim, v1_sim = initial_conditions(ref_traj);

# ## Simulator
# u_ref = deepcopy(ref_traj.u[1])
# p = open_loop_policy([u_ref for i=1:1000], N_sample=N_sample)
# p = open_loop_policy([(fill(ref_traj.u, 100)...)...], N_sample=N_sample)
# sim = simulator(s, H_sim, h=h_sim, dist=d);
sim = simulator(s, H_sim, h=h_sim, policy=p, dist=d);

using BenchmarkTools
# ## Simulate
q1_sim0 = deepcopy(q1_sim)
# q1_sim0[4:6] += [0.02, 0.04, 0.02]
Δx = 0.0*[0.1, 0.2, 0.1]
q1_sim0[1:3] += Δx
q1_sim0[7:9] += Δx
q1_sim0[10:12] += Δx
q1_sim0[13:15] += Δx
q1_sim0[16:18] += Δx

RoboDojo.simulate!(sim, q1_sim0, v1_sim)

# # ## Visualizer
vis = ContactImplicitMPC.Visualizer()
ContactImplicitMPC.open(vis)

# ## Visualize
set_light!(vis)
set_floor!(vis, grid=true)
set_background!(vis)
anim = visualize!(vis, model, sim.traj.q; Δt=h_sim)

# u_force = impulses[1]
# u_force ./= norm(u_force)
# u_force .*= 0.35
# force_vis = ArrowVisualizer(vis[:force])
# setobject!(force_vis, MeshPhongMaterial(color=RGBA(1.0, 0.0, 0.0, 1.0)))
# settransform!(force_vis,
#             Point(sim.traj.q[1][1], sim.traj.q[1][2], sim.traj.q[1][3]),
#             Vec(u_force[1], u_force[2], u_force[3]),
#             shaft_radius=0.035,
#             max_head_radius=0.075)
# setvisible!(vis[:force], true)
# for t = 1:length(sim.traj.q)
#     MeshCat.atframe(anim, t) do
#         t < 75 ? setvisible!(vis[:force], true) : setvisible!(vis[:force], false)
#     end
# end
# MeshCat.setanimation!(vis, anim)

# 12 / 6000 == 0.01 / 5

# thr = 1.0e-3
# N = 18
# H_ref = length(ref_traj.q)
# foot1_ref = [([[q[6 .+ 3] > thr ? 0 : 1 for q in ref_traj.q] for t = 1:N]...)...]
# foot2_ref = [([[q[9 .+ 3] > thr ? 0 : 1 for q in ref_traj.q] for t = 1:N]...)...]
# foot3_ref = [([[q[12 .+ 3] > thr ? 0 : 1 for q in ref_traj.q] for t = 1:N]...)...]
# foot4_ref = [([[q[15 .+ 3] > thr ? 0 : 1 for q in ref_traj.q] for t = 1:N]...)...]
# time_ref = range(0, stop=h * (N * H_ref - 1), length=N * H_ref)

plt = plot(layout=(4, 1));

plt = plot!(plt, time_ref, foot1_ref,
    linetype=:steppost, color=:orange, width=3.0, label="", yaxis=nothing, subplot=1)
plt = plot!(plt, time_ref, foot2_ref,
    linetype=:steppost, color=:lightgreen, width=3.0, label="", yaxis=nothing, subplot=2)
plt = plot!(plt, time_ref, foot3_ref,
    linetype=:steppost, color=:cyan, width=3.0, label="", yaxis=nothing, subplot=3)
plt = plot!(plt, time_ref, foot4_ref,
    linetype=:steppost, color=:magenta, width=3.0, label="", yaxis=nothing, subplot=4)

# foot1_sim = [q[6 .+ 3] > thr ? 0.01 : 0.99 for q in sim.traj.q]
# foot2_sim = [q[9 .+ 3] > thr ? 0.01 : 0.99 for q in sim.traj.q]
# foot3_sim = [q[12 .+ 3] > thr ? 0.01 : 0.99 for q in sim.traj.q]
# foot4_sim = [q[15 .+ 3] > thr ? 0.01 : 0.99 for q in sim.traj.q]
# time_sim = range(0, stop=h_sim * (H_sim - 1), length=H_sim)

# plt = plot!(plt, time_sim, foot1_sim[2:end-1],
#     linetype=:steppost, color=:black, width=2.0, label="", yaxis=nothing, subplot=1)
# plt = plot!(time_sim, foot2_sim[2:end-1],
#     linetype=:steppost, color=:black, width=2.0, label="", yaxis=nothing, subplot=2)
# plt = plot(plt, time_sim, foot3_sim[2:end-1],
#     linetype=:steppost, color=:black, width=2.0, label="", yaxis=nothing, subplot=3)
# plt = plot!(plt, time_sim, foot4_sim[2:end-1],
#     linetype=:steppost, color=:black, width=2.0, label="", yaxis=nothing, subplot=4, xlabel="time (s)")

# display(plt)
# savefig(plt, "/home/taylor/Downloads/point_foot_quadruped_push_recovery.png")

# # # ## Timing result
# # # Julia is [JIT-ed](https://en.wikipedia.org/wiki/Just-in-time_compilation) so re-run the MPC setup through Simulate for correct timing results.
process!(sim.stats, N_sample) # Time budget
H_sim * h_sim / sum(sim.stats.policy_time) # Speed ratio
plot(sim.stats.policy_time, xlabel="timestep", ylabel="mpc time (s)", label="", linetype=:steppost)

# # convert_frames_to_video_and_gif("point_foot_push_recovery")

# sim.policy.ref_traj
# for i = 1:100
# 	@show 1 + Int(floor((i-1)/N_sample)) % 10
# end

# sim.traj
# sim.policy.ref_traj


text:
function available_queue_families(device)
    qf_count = Ref{UInt32}(0)
    vkGetPhysicalDeviceQueueFamilyProperties(device, qf_count, C_NULL)
    # vkGetPhysicalDeviceQueueFamilyProperties(C_NULL, qf_count, C_NULL)
    qf_props = Array{VkQueueFamilyProperties}(undef, qf_count[])
    vkGetPhysicalDeviceQueueFamilyProperties(device, qf_count, qf_props)
    return qf_props
end

function select_queue_family(physical_device; predicate = qf->qf.queueFlags & VK_QUEUE_GRAPHICS_BIT)
    avail_qf = available_queue_families(physical_device)
    for (index, qf) in enumerate(avail_qf)
        # @warn Bool(predicate(qf))
        if Bool(predicate(qf))
            return index
        end
    end
    throw(ErrorException("No suitable queue family found"))
end

function create_queue_info(queue_family_index; queue_count = 1, queue_priorities = [1.0])
    cqueue_info = Container(Ref(VkDeviceQueueCreateInfo(VK_STRUCTURE_TYPE_DEVICE_QUEUE_CREATE_INFO, C_NULL, 0, queue_family_index, queue_count, Base.unsafe_convert(Ptr{Float32}, Array{Float32,1}(queue_priorities)))), [queue_priorities])

    cqueue_info
end


text:
"""
**`Helium.readhe`** -*Function*.

    `Helium.readhe(source)` => `Matrix`

Parse a binary file into a matrix. The binary files .he format includes a
header that describes the size and the data type of the matrix, as well as the
endianness of the file.
The first 56 bytes contains:
[number of rows, number of columns, data type, endianness,
 if column names exists, if row names exists, number of column of supplement]

"""
function readhe(heFile::String)
        # Get data properties
        matInfo = readheader(heFile)

        # Pre-allocate matrix
        mat = Array{matInfo.datatype}(undef, matInfo.numrows, matInfo.numcols);

        open(heFile) do file
            seek(file, 56);
            read!(file, mat)
        end

        # Check endianness
        dictEndian = Dict("Little-endian"=>ntoh,
                          "Big-endian"=>ltoh,
                          "4030201"=>"Little-endian",
                          "1020304"=>"Big-endian");
        hostendian =  dictEndian[string(ENDIAN_BOM, base = 16)]

        if hostendian != matInfo.endianness
                mat .= dictEndian[hostendian].(mat)
        end

        return mat
end

"""
**`Helium.getcolnames`** -*Function*.

    `Helium.getcolnames(source)` => `Array{String, 1}`

Return an array of string containing the column names associated to the data.

"""
function getcolnames(heFile::String)
        # Get data properties
        matInfo = readheader(heFile)

        if matInfo.hascolnames
                # Read COL NAMES
                global cNames =  Array{String,1}[]
                open(heFile, "r") do file
                        seek(file, (7 + matInfo.numrows*matInfo.numcols)*matInfo.bytesize)
                        cNames  = readline(file)
                end
                cNames = string.(split(cNames, ","));

                return cNames
        else
                error("No column names in the file!!!")
        end
end

"""
**`Helium.getrownames`** -*Function*.

    `Helium.getrownames(source)` => `Array{String, 1}`

Return an array of string containing the row names associated with the data.

"""
function getrownames(heFile::String)
        # Get data properties
        matInfo = readheader(heFile)

        if matInfo.hasrownames
                # Read ROW NAMES
                global rNames =  Array{String,1}[]
                open(heFile, "r") do file
                        seek(file, (7 + matInfo.numrows*matInfo.numcols)*matInfo.bytesize)
                        if matInfo.hascolnames
                                rNames= readline(file)
                        end
                        rNames  = readline(file)
                end
                rNames = string.(split(rNames, ","));

                return rNames
        else
                error("No row names in the file!!!")
        end
end

"""
**`Helium.getsupp`** -*Function*.

    `Helium.getsupp(source)` => `Array{String, 2}`

Return an array of string containing the supplement data associated with the
matrix data.

"""
function getsupp(heFile::String)
        # Get data properties
        matInfo = readheader(heFile)

        if matInfo.numcolsupp != 0
                # Read ROW NAMES
                global supp =  Array{String}(undef, matInfo.numrows+matInfo.hascolnames, matInfo.numcolsupp)
                global i = 1
                open(heFile, "r") do file
                        seek(file, (7 + matInfo.numrows*matInfo.numcols)*matInfo.bytesize)

                        if matInfo.hascolnames
                                tmp = readline(file)
                        end

                        if matInfo.hasrownames
                                tmp = readline(file)
                        end

                        while !eof(file)
                                tmp = readline(file)
                                supp[:, i] = string.(split(tmp, ","));
                                i += 1
                        end
                end

                supp #return supp
        else
                error("No supplement in the file!!!")
        end
end

"""
**`Helium.readheader`** -*Function*.

    `Helium.readheader(source)` => `NamedTuple{(:numrows, :numcols,
                                                :datatype, :endianness,
                                                 :hascolnames, :hasrownames,
                                                 :bytesize),
                                                 Tuple{Int64,Int64,
                                                       DataType,String,
                                                       Bool,Bool,
                                                       Int64}}`

Return a tuple containing respectively number of rows, number of columns,
data type, endianness, if there is column names , and if if there is row
names.

"""
function readheader(heFile::String)

        try
           realpath(abspath(heFile))
        catch e
           error("File name or path incorrect!!!")
        end

        heFile =  realpath(abspath(heFile))

        # Generate dictionary for reading type and endianness
        dictCtrl =  Dict("f64"=>Float64,
                    "e64"=>Int64,
                    "f32"=>Float32,
                    "e32"=>Int32,
                    "b00"=>Bool,
                    "4030201"=>"Little-endian",
                    "1020304"=>"Big-endian");

        dictDataSize =  Dict("f64"=>8,
                        "e64"=>8,
                        "f32"=>4,
                        "e32"=>4,
                        "b00"=>1);

        # Read control information
        matCtrl = Array{Int64}(undef, 1 ,7);
        open(heFile) do file
                read!(file, matCtrl)
        end

        # Check endianness
        dictEndian = Dict("Little-endian"=>ntoh,
                          "Big-endian"=>ltoh);

        if ENDIAN_BOM != matCtrl[4]
                hostendian =  dictCtrl[string(ENDIAN_BOM, base = 16)]
                matCtrl .= dictEndian[hostendian].(matCtrl)
        end

        return (numrows = matCtrl[1], numcols = matCtrl[2],
                datatype = dictCtrl[string(matCtrl[3], base = 16)],
                endianness = dictCtrl[string(matCtrl[4], base = 16)],
                hascolnames = Bool(matCtrl[5]), hasrownames = Bool(matCtrl[6]),
                numcolsupp = matCtrl[7],
                bytesize = dictDataSize[string(matCtrl[3], base = 16)])

end


text:
"""
Custom number types for extended-precision computations.
"""
module CustomTypes

export myI, myF, myC, tomyI, tomyF, tomyC

#------------------------------------
const myI = Int64
const myF = Float64
const myC = ComplexF64
#------------------------------------
#setprecision(1024)
#const myI = Int64
#const myF = BigFloat
#const myC = Complex{BigFloat}
#------------------------------------

function tomyI(v)::myI
   return convert(myI,parse(myF,"$v"))
end

function tomyF(v)::myF
   return parse(myF,"$v")
end

function tomyC(v)::myC
   return parse(myC,"$v")
end

end


text:
### SELECTED DATA FILE INDEXES REQUIRED FOR TGFb MODEL ONLY

# INDEX VALUES CHANGE WITH UPDATES TO MODEL.NET

initial_cond_update_array = [];
rate_constant_update_array = [];
control_parameter_update_array = [];

# initial conditions
push!(initial_cond_update_array,[vegf_initial_cond_index,0.0]) # VEGF

# rates
# VEGF dosing
push!(rate_constant_update_array,[vegf_forward_rate_index,4.0]) # addition_VEGF: [] =([])=> VEGF
push!(rate_constant_update_array,[vegf_reverse_rate_index,0.1]) # addition_VEGF: VEGF =([])=> [] (reverse)

# push modifications to RAM
#UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
#UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)


text:
using GLVisualize
include(GLVisualize.dir("examples", "ExampleRunner.jl"))
using ExampleRunner
importall ExampleRunner
using GLAbstraction, GLWindow, Colors
using FileIO, GeometryTypes, Reactive, Images

function image_url(path)
    path, img = splitdir(path)
    path, imgfolder = splitdir(path)
    path, pr = splitdir(path)
    path, package = splitdir(path)
    path, ci = splitdir(path)
    "https://github.com/SimiDCI/GLVisualizeCI.jl/blob/master/reports/$ci/$package/$pr/$imgfolder/$img?raw=true"
end

function create_mosaic(io, folder, width = 150)
    images = filter(x-> endswith(x, ".jpg"), readdir(folder))
    for im in images
        println(io, """<img src="$(image_url(joinpath(folder, im)))"
            alt="$(im)" width=$(width)px"/>
        """)
    end
end

if haskey(ENV, "CI_REPORT_DIR")
    full_folder = ENV["CI_REPORT_DIR"]
    recording = true
else
    recording = false
    full_folder = ""
end

files = [
    "introduction/rotate_robj.jl",
    "introduction/screens.jl",
    "plots/3dplots.jl",
    "plots/lines_scatter.jl",
    "plots/hybrid.jl",
    "camera/camera.jl",
    "gui/color_chooser.jl",
    "gui/image_processing.jl",
    "gui/buttons.jl",
    "gui/fractal_lines.jl",
    "gui/mandalas.jl",
    "plots/drawing.jl",
    "interactive/graph_editing.jl",
    "interactive/mario_game.jl",
    "text/text_particle.jl",
]

map!(x-> GLVisualize.dir("examples", x), files, files)
files = union(files, ExampleRunner.flatten_paths(GLVisualize.dir("examples")))
#push!(files, GLVisualize.dir("test", "summary.jl"))
files = unique(normpath.(files))
# Create an examplerunner, that displays all examples in the example folder, plus
# a runtest specific summary.
resolution = (200, 200)
rootscreen = glscreen(
    resolution = resolution, visible = false,
    focus = false, debugging = true
)
config = ExampleRunner.RunnerConfig(
    record = false,
    record_image = true,
    files = files,
    thumbnail = false,
    screencast_folder = full_folder,
    resolution = resolution,
    rootscreen = rootscreen
)
window = config.window
isdir(full_folder)

if recording
    imgpath = joinpath(full_folder, "images")
    isdir(imgpath) || mkdir(imgpath)
end

for path in config.files
    isopen(config.rootscreen) || break
    try
        println(basename(path))
        test_module = ExampleRunner._test_include(path, config)
        ExampleRunner.display_msg(test_module, config)
        GLWindow.poll_glfw()
        GLWindow.reactive_run_till_now()
        yield()
        render_frame(config.rootscreen)
        swapbuffers(config.rootscreen)
        if recording
            name = basename(path)[1:end-3]
            name = joinpath(imgpath, name * ".jpg")
            GLWindow.screenshot(config.rootscreen, path = name)
        end
    catch e
        #failed[i] = true
        bt = catch_backtrace()
        ex = CapturedException(e, bt)
        showerror(STDERR, ex)
        config[:success] = false
        config[:exception] = ex
    finally
        yield()
        empty!(window)
        window.color = RGBA{Float32}(1,1,1,1)
        window.clear = true
        GLVisualize.empty_screens!()
        empty!(window.cameras)
        GLVisualize.add_screen(window) # make window default again!
        gc()
    end
    yield()
end

failures = filter(config.attributes) do k, dict
    !dict[:success]
end
if recording
    open(joinpath(full_folder, "report.md"), "w") do io
        println(io, "### Test Images:")
        create_mosaic(io, imgpath)
        if !isempty(failures)
            println(io, "### Failures:")
            for (k, dict) in failures
                println(io, "file: $k")
                Base.showerror(io, "$(dict[:exception])")
                println("\n")
            end
        else
            println(io, "No failures! :)")
        end
    end
end

isempty(failures) || error("Tests did not pass!")


text:
#
# Data for the simulation
#

using Parameters

@with_kw struct Data

  eps :: Float64 = 5.0
  sig :: Float64 = 0.5

  cutoff :: Float64 = 2.0
  side :: Float64 = 100.

  N :: Int64 = 10_000

  # Data for computation
  eps4 :: Float64 = 4*eps
  sig6 :: Float64 = sig^6
  sig12 :: Float64 = sig^12

end

#function Data(eps,sig,cutoff,side,N::Int)
#  eps4 = 4*eps
#  sig6 = sig^6
#  sig12 = sig6^2
#  return Data(eps,sig,cutoff,side,N,eps4,sig6,sig12)
#end
#
#function Data(N::Int)
#  eps = 5.0
#  sig = 0.5
#  cutoff = 2.
#  side = 100.
#  return Data(eps,sig,cutoff,side,N)
#end


text:
module MLStyle

# Flags
export @use, use, @used
# Match Implementation
export @match, gen_match
# DataTypes
export @data

# Pervasive Patterns
export Many, Do
# Active Patterns
export @active
# Extensibilities
export def_pattern, def_app_pattern, def_gapp_pattern, mk_pattern, mk_app_pattern, mk_gapp_pattern, def_record, def_active_pattern
# Exceptions
export PatternUnsolvedException, InternalException, SyntaxError, UnknownExtension, @syntax_err
# Syntax Sugars
export @as_record
export @λ, gen_lambda
export @when, @otherwise, gen_when

# convenient modules
export Modules

include("Err.jl")
using MLStyle.Err

include("Extension.jl")
using MLStyle.Extension

include("Internal/Toolz.jl")

include("Render.jl")

include("MatchCore.jl")
using MLStyle.MatchCore

include("Infras.jl")
using MLStyle.Infras

include("Pervasives.jl")
using MLStyle.Pervasives

include("Qualification.jl")

include("TypeVarExtraction.jl")

include("StandardPatterns/TypeVarDecons.jl")
include("StandardPatterns/Active.jl")
using MLStyle.Active

include("Record.jl")
using MLStyle.Record

include("DataType.jl")
using MLStyle.DataType

include("StandardPatterns/Uncomprehensions.jl")

include("StandardPatterns/LambdaCases.jl")
using MLStyle.LambdaCases

include("StandardPatterns/WhenCases.jl")
using MLStyle.WhenCases

include("Modules/Modules.jl")

end # module


text:
# This file contains macros and functions
# used in the tests of this package

"""
Create a temporary file
"""
macro file()
    return esc(:((file, io) = mktemp()))
end

"""
Create an instance of the API
"""
macro instantiate()
    return esc(:(s = API()))
end

"""
Construct a function for constructing an
exception passing only a type and an ending
"""
macro construct_exceptions(prefix::AbstractString)
    return esc(quote
        """
        Construct an exception passing only a type and an ending
        """
        function construct_exceptions(
            type::AbstractString,
            ending::AbstractString
        )
            return eval(
                Meta.parse(
                    string($prefix, type, "_", ending)
                )
            )
        end
    end)
end

"""
Construct an array of exceptions passing only a type and a tuple of endings
"""
function construct_exceptions(
    type::AbstractString,
    endings::AbstractString...
)
    return construct_exceptions.(type, endings)
end

"""
Retrieve error messages from exceptions:
"""
function retrieve_messages(exceptions::Union{DataType, Tuple})
    return sprint.(showerror, exceptions)
end

"""
Append a call to a file and retrieve an error message from the exception:
"""
function retrieve_messages(exception::DataType, file::AbstractString)
    return sprint(showerror, exception(file))
end

"""
Append calls to files and retrieve error messages from exceptions:
"""
function retrieve_messages(exceptions::Tuple, file::AbstractString)
    return [ sprint(showerror, exception(file)) for exception in exceptions ]
end


text:
# enable type inference for constructor
Linear(::Type{T}, ::Type{I}, weight::W, bias::B) where {T, I, W, B} =
    Linear{T, W, B, I}(weight, bias, zeros(weight), bias === nothing ? nothing : zeros(bias))

Linear(weight::WeightType{T}, bias::BiasType{T}; nbatch=1) where {T} =
    Linear(T, linear_input_type(weight, nbatch), weight, bias)
Linear(weight::WeightType{T}; nbatch=1) where {T} =
    Linear(weight, nothing; nbatch=nbatch)

import Base: size

size(op::Linear) = size(op.weight)
size(op::Linear, d::Int) = reverse(size(op.weight))[d]

function size(op::Linear, d::Symbol)
    if d == :input
        return size(op, 1)
    elseif d == :output
        return size(op, 2)
    else
        # TODO: custom exception
        throw(ErrorException("Invalid tag, Linear layer's
        size tag should be :input or :output"))
    end
end

function zerograd!(op::Linear)
    fill!(op.grad_weight, 0)
    op.grad_bias === nothing ? fill!(op.bias, 0) : nothing
    return op
end

function initialize!(op::Linear)
    stdv = 1 / sqrt(size(op.weight, 1))
    rand!(op.weight)
    op.weight = stdv * (2 * op.weight - 1) # stdv * uniform(-1, 1)

    if op.bias !== nothing
        rand!(op.bias)
        op.bias = stdv * (2 * op.bias - 1)
    end
    return op
end

# Initialization
function Linear(::Type{T}, in_features::Int, out_features::Int;bias=true, nbatch=1) where T
    op = Linear(zeros(T, out_features, in_features), bias ? zeros(T, out_features) : nothing; nbatch=nbatch)
    return initialize!(op)
end

Linear(::Type{T}, shape::Tuple{Int, Int};bias=true, nbatch=1) where T =
    Linear(T, shape...; bias=bias, nbatch=nbatch)

Linear(shape::Tuple{Int, Int};bias=true, nbatch=1) =
    Linear(Float64, shape; bias=bias, nbatch=nbatch)
Linear(shape::Int...;bias=true, nbatch=1) =
    Linear(shape;bias=bias, nbatch=nbatch)

# TODO: use static array when the layer is small, benchmark needed for this


text:
export expectation, computehermite, unravel_index, rieszmap, extractmoments

# expectation of a function with respect to an atomic measure
function expectation(μ::AbstractMeasureLike, f::Function)
    sum(δ -> δ.weight * f(δ.center), μ.atoms)
end

# computation of Hermite polynomials' coefficients
function computehermite(maxdegree::Int64)
    a = zeros(Int64, maxdegree+1, maxdegree+1)
    a[1, 1] = 1
    a[2, 2] = 2

    for n=2:maxdegree
        a[n+1, 1] = - a[n, 2]
        a[n+1, n] = 2 * a[n, n-1]
        a[n+1, n+1] = 2 * a[n, n]
        for k=2:maxdegree
            a[n+1, k] = 2 * a[n, k-1] - k * a[n, k+1]
        end
    end

    return a
end

function unravel_index(x, n)
    return (div(x - 1, n) + 1, mod(x - 1, n) + 1)
end

# it is implicitly assumed that a given monomial cannot be found twice in y
function rieszmap(y::Array{MT, 1}, φ::AbstractPolynomialLike) where {MT<:AbstractMoment}
    sum(moment_value(y[i]) * coefficient(φ, monomial(y[i])) for i=1:length(y))
end

function extractmoments(M::MomentMatrix)
    n = M.Q.n
    σ, monos = sortmonovec(vcat([[M.x[i] * M.x[j] for j=1:n] for i=1:n]...))
    [moment(M.Q[unravel_index(σ[i], n)...], monos[i]) for i=1:length(σ)]
end


text:
using StatsBase
using Base.Test

##### weighted var & std

x = rand(10)
wv = weights(rand(10))
m = mean(x, wv)

@test var(x, wv)           ≈ sum(abs2.(x .- m), wv) ./ sum(wv)
@test var(x, wv; mean=0)   ≈ sum(abs2.(x), wv) ./ sum(wv)
@test var(x, wv; mean=1.0) ≈ sum(abs2.(x .- 1.0), wv) ./ sum(wv)

@test std(x, wv)           ≈ sqrt(var(x, wv))
@test std(x, wv; mean=0)   ≈ sqrt(var(x, wv; mean=0))
@test std(x, wv; mean=1.0) ≈ sqrt(var(x, wv; mean=1.0))

(m, v) = mean_and_var(x)
@test m == mean(x)
@test v == var(x)

(m, s) = mean_and_std(x)
@test m == mean(x)
@test s == std(x)

(m, v) = mean_and_var(x, wv)
@test m == mean(x, wv)
@test v == var(x, wv)

(m, s) = mean_and_std(x, wv)
@test m == mean(x, wv)
@test s == std(x, wv)

x = rand(5, 6)
w1 = rand(5)
w2 = rand(6)
wv1 = weights(w1)
wv2 = weights(w2)
m1 = mean(x, wv1, 1)
m2 = mean(x, wv2, 2)

@test var(x, wv1, 1; mean=0) ≈ sum(abs2.(x) .* w1, 1) ./ sum(wv1)
@test var(x, wv2, 2; mean=0) ≈ sum(abs2.(x) .* w2', 2) ./ sum(wv2)

@test var(x, wv1, 1; mean=m1) ≈ sum(abs2.(x .- m1) .* w1, 1) ./ sum(wv1)
@test var(x, wv2, 2; mean=m2) ≈ sum(abs2.(x .- m2) .* w2', 2) ./ sum(wv2)

@test var(x, wv1, 1) ≈ sum(abs2.(x .- m1) .* w1, 1) ./ sum(wv1)
@test var(x, wv2, 2) ≈ sum(abs2.(x .- m2) .* w2', 2) ./ sum(wv2)

@test std(x, wv1, 1)          ≈ sqrt.(var(x, wv1, 1))
@test std(x, wv2, 2)          ≈ sqrt.(var(x, wv2, 2))
@test std(x, wv1, 1; mean=0)  ≈ sqrt.(var(x, wv1, 1; mean=0))
@test std(x, wv2, 2; mean=0)  ≈ sqrt.(var(x, wv2, 2; mean=0))
@test std(x, wv1, 1; mean=m1) ≈ sqrt.(var(x, wv1, 1; mean=m1))
@test std(x, wv2, 2; mean=m2) ≈ sqrt.(var(x, wv2, 2; mean=m2))

for d in 1:2
    (m, v) = mean_and_var(x, d)
    @test m == mean(x, d)
    @test v == var(x, d)

    (m, s) = mean_and_std(x, d)
    @test m == mean(x, d)
    @test s == std(x, d)
end

(m, v) = mean_and_var(x, wv1, 1)
@test m == mean(x, wv1, 1)
@test v == var(x, wv1, 1)

(m, v) = mean_and_var(x, wv2, 2)
@test m == mean(x, wv2, 2)
@test v == var(x, wv2, 2)

(m, s) = mean_and_std(x, wv1, 1)
@test m == mean(x, wv1, 1)
@test s == std(x, wv1, 1)

(m, s) = mean_and_std(x, wv2, 2)
@test m == mean(x, wv2, 2)
@test s == std(x, wv2, 2)

##### skewness & kurtosis

wv = weights(ones(5) * 2.0)

@test skewness(1:5)             ≈  0.0
@test skewness([1, 2, 3, 4, 5]) ≈  0.0
@test skewness([1, 2, 2, 2, 5]) ≈  1.1731251294063556
@test skewness([1, 4, 4, 4, 5]) ≈ -1.1731251294063556

@test skewness([1, 2, 2, 2, 5], wv) ≈ 1.1731251294063556

@test kurtosis(1:5)             ≈ -1.3
@test kurtosis([1, 2, 3, 4, 5]) ≈ -1.3
@test kurtosis([1, 2, 3, 3, 2]) ≈ -1.1530612244897953

@test kurtosis([1, 2, 3, 4, 5], wv) ≈ -1.3

##### general moments

x = collect(2.0:8.0)
@test moment(x, 2) ≈ sum((x .- 5).^2) / length(x)
@test moment(x, 3) ≈ sum((x .- 5).^3) / length(x)
@test moment(x, 4) ≈ sum((x .- 5).^4) / length(x)
@test moment(x, 5) ≈ sum((x .- 5).^5) / length(x)

@test moment(x, 2, 4.0) ≈ sum((x .- 4).^2) / length(x)
@test moment(x, 3, 4.0) ≈ sum((x .- 4).^3) / length(x)
@test moment(x, 4, 4.0) ≈ sum((x .- 4).^4) / length(x)
@test moment(x, 5, 4.0) ≈ sum((x .- 4).^5) / length(x)

w = weights([1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0])
x2 = collect(2.0:6.0)
@test moment(x, 2, w) ≈ sum((x2 .- 4).^2) / 5
@test moment(x, 3, w) ≈ sum((x2 .- 4).^3) / 5
@test moment(x, 4, w) ≈ sum((x2 .- 4).^4) / 5
@test moment(x, 5, w) ≈ sum((x2 .- 4).^5) / 5


text:
#module Core
  export is_cons, assp, ext_s, call_fresh, length, equals, list, nil, Nil, bind, mplus, mzero,
    Var, var, conj, disj, car, cdr, take, take_all, vars_equal, cons

immutable Nil
end
nil() = Nil()

immutable Var
  index :: Int
end

is_var(x) = isa(x, Var)
vars_equal(x1 :: Var, x2 :: Var) = x1.index == x2.index

car(x :: Pair) = first(x)
cdr(x :: Pair) = last(x)
car(x :: Array) = x[1]
cdr(x :: Array) = x[2 : end]

is_cons(d) = isa(d, Pair) && isa(cdr(d), Pair)
is_pair(x) = isa(x, Pair) # && !isa(cdr(x), Pair)

cons(x, y) = Pair(x, y)

list() = nil()
function list(values...)
  Pair(values[1], list(values[2:end]...))
end

is_procedure(elt) = isa(elt, Function)

assp(func :: Function, alist :: Nil) = nil()
function assp(func :: Function, alist :: Pair)
  first_pair = car(alist)
  first_pair == nil() && return nil()
  first_value = car(first_pair)
  func(first_value) && return first_pair
  assp(func, cdr(alist))
end

function walk(u, s)
  if is_var(u)
    pr = assp(v -> vars_equal(u, v), s)
    pr != nil() && return walk(cdr(pr), s)
  end
  u
end

ext_s(x, v, s) = occurs(x,v,s) ? s : cons(Pair(x, v), s)

mzero = nil()
unit = s_c :: Pair -> cons(s_c, mzero)

function unify(u, v, s)
  u = walk(u, s)
  v = walk(v, s)
  is_var(u) && is_var(v) && vars_equal(u, v) && return s
  is_var(u) && return ext_s(u, v, s)
  is_var(v) && return ext_s(v, u, s)
  if is_pair(u) && is_pair(v)
    s = unify(car(u), car(v), s)
    s != nil() && begin
      return unify(cdr(u), cdr(v), s)
    end
  else
    u == v && return s
  end
  nil()
end

# Call function f with a fresh variable.
function call_fresh(f)
  s_c -> begin
    c = cdr(s_c)
    if isa(f, Expr)
      f = eval(f)
    end
    f(Var(c))(Pair(car(s_c), c + 1))
  end
end

disj(g1, g2) = s_c :: Pair -> mplus(g1(s_c), g2(s_c))
conj(g1, g2) = s_c :: Pair -> bind(g1(s_c), g2)

function equals(u, v)
  s_c :: Pair -> begin
    s = unify(u, v, car(s_c))
    s != mzero ? unit(Pair(s, cdr(s_c))) : mzero
  end
end

mplus(goals1 :: Nil, goals2) = goals2
mplus(goals1 :: Function, goals2) = () -> mplus(goals2, goals1())
mplus(goals1 :: Pair, goals2) = cons(car(goals1), mplus(goals2, cdr(goals1)))

bind(stream :: Nil, goal :: Any) = mzero
bind(stream :: Function, goal) = () -> bind(stream(), goal)
bind(stream :: Pair, goal :: Any) = mplus(goal(car(stream)), bind(cdr(stream), goal))

import Base.string

string(p :: Pair) = last(p) != nil() ? Base.string("(", p[1], " . ", p[2],")") : Base.string("(", p[1],")")
string(v :: Var) = Base.string("#(", v.index, ")")
string(v :: Nil) = Base.string("#(nil)")

import Base.print
print(io :: IO, p :: Pair) = print(io, string(p))
print(io :: IO, v :: Var)  = print(io, string(v))
print(io :: IO, v :: Nil)  = print(io, string(v))
import Base.show
show(io :: IO, p :: Pair) = print(io, p)
show(io :: IO, v :: Var) = print(io, v)
show(io :: IO, v :: Nil) = print(io, v)

#call_fresh and so on
################################
############### macro definitions
################################
export @Zzz, @fresh, @conj_, @disj_, @conde
export @fresh_helper #this gets exported for testing purpose

macro Zzz(g)
  :(s_c -> () -> $(esc(g))(s_c))
end

macro conj_(g0, g...)
  if (isempty(g))
    return :(@Zzz($(esc(g0))))
  else
    local t = [:($(esc(gg))) for gg in g]
    quote
      conj(@Zzz($(esc(g0))), @conj_($(t...)))
    end
  end
end

macro disj_(g0, g...)
  if (isempty(g))
    return :(@Zzz($(esc(g0))))
  else
    local t = [:($(esc(gg))) for gg in g]
    quote
      disj(@Zzz($(esc(g0))), @disj_($(t...)))
    end
  end
end

macro conde(g...)
  local cc = [ begin
      if isa(gg, Expr)
        if gg.head == :tuple
          [ :($(esc(a))) for a in gg.args]
        else
        [ :($(esc(gg))) ]
        end
      else
        [ :($(esc(gg))) ]
        #todo: add other types check, this is in case when gg is Symbol
      end
    end
    for gg in g ]
  local values = [ :(@conj_($(c...))) for c in cc]

  :( @disj_($(values...)) )
end

macro fresh_helper(g0, vars...)
  if (isempty(vars))
    return :($(esc(g0)))
  else
    local vars0 = vars[1]
    if length(vars) > 1
      local vars1 = vars[2:end]
      local exp = :($(esc(Expr(:->, :($vars0), :(@fresh_helper($g0, $(vars1...)))))))
      quote
        call_fresh($exp)
      end
    else
      local exp2 = :($(esc(Expr(:->, vars0, g0))))
      return :( call_fresh($exp2) )
    end
  end
end

macro fresh(vars, g0, g...)
  if isempty(g)
    local c = begin
      if isa(vars, Symbol)
        [:($(esc(eval(QuoteNode(vars)))))]
      elseif isa(vars, Expr)
        [:($(esc(eval(isa(v, Symbol) ? QuoteNode(v) : v)))) for v in vars.args]
      else
        []
      end
    end
    #println(c)
    :(@fresh_helper($(esc(g0)), $(c...) ))
  else
    local glist = [:($(esc(gg))) for gg in g]
    :(@fresh($(esc(vars)), @conj_($(esc(g0)), $(glist...)) ) )
  end
end

############# 5.2
function pull(s)
  #println(string("pull ", s))
  is_procedure(s) ? pull(s()) : s
end

function take_all(s)
  s = pull(s)
  #println(s)
  s == nil() && return []
  [car(s), take_all(cdr(s))...]
end

function take(n, s)
  n == 0 &&  return []
  s = pull(s)
  s == nil() && return []
  [car(s), take(n - 1, cdr(s))...]
end

################
### reification
###############

reify_name(n :: Int) = string("_.", n)

function map(f :: Function, p :: Pair)
  p == nil() && return nil()
  Pair(f(car(p)), map(f, cdr(p)))
end

function reify_s(v, s)
  length(n :: Nil) = 0
  length(p:: Pair) = cdr(p) == nil() && car(p) == nil() ? 0 : (car(p)!= nil() ? 1: 1 + length(cdr(p)))
  v = walk(v, s)
  #println(s)
  is_var(v) && return cons(Pair(v, reify_name(length(s))), s)
  is_pair(v) && return reify_s(cdr(v), reify_s(car(v), s))
  s
end

function walk_star(v, s)
  v = walk(v, s)
  is_var(v) && return v
  is_pair(v) && return cons(walk_star(car(v), s), walk_star(cdr(v), s))
  v
end

mk_reify(s_cs) = Base.map(reify_state_1st_var, s_cs)

export reify_state_1st_var, mk_reify

function reify_state_1st_var(s_c :: Pair)
  v = walk_star(Var(0), car(s_c))
  walk_star(v, reify_s(v, nil()))
end

empty_state = Pair(nil(), 0)

call_empty_state(g) = g(empty_state)

#############3 run macros
export @run, @run_star

#todo: fix im lazy now to do it
macro run(n, vars, g0, g...)
  #println("in run ", vars, " ", g0)
  "this should eval expression"
  local t = [:($(esc(gg))) for gg in g]
  :(mk_reify(take($n, call_empty_state(@fresh($(esc(vars)), $(esc(g0)), $(t...))))))
end

macro run_star(vars, g0, g...)
  "this should eval expression too"
  local t = [:($(esc(gg))) for gg in g]
  :(mk_reify(take_all(call_empty_state(@fresh($(esc(vars)), $(esc(g0)), $(t...))))))
end

function occurs(x, v, s)
  v = walk(v, s)
  is_var(v) && return vars_equal(v, x)
  is_pair(v) && (occurs(x, car(v), s) || occurs(x, cdr(v), s))
end
#next is module end
#end


text:
module BayesianTransGaussian

using Reexport
using Distributions
using Plots
using Roots
using Distances

include("./BTGModel.jl")

@reexport using .BTGModel

# TODO Finish Configuring Documenter.jl to generate documentation
# TODO When publishing ready, use Registrator.jl
# TODO Optimize code using Profile.jl and BenchmarkTools.jl
# TODO Give Cameron feedback on organization/interface
# TODO Pick which transforms/kernels we want to implement
# TODO Add any covariate transforms we want (polynomial, orthopolynomial for 1d, induced point)
# TODO Add tests to test/runtests.jl, if we want to do unit testing

# See BTGTest.jl for a basic example
# I recommend using Revise.jl (see workflow tips section in Julia docs)

"""
    predictquantile(mdl, q, x)

The inverse cumulative distribution function. The value at location `x` at the `q`th quantile.
"""
function predictquantile(mdl::Model, q, x)
    d = exp.(-pairwise(Euclidean(), mdl.X, reshape(x, 1, length(x)), dims=1))
    startpoint = sum(d .* mdl.Y) / sum(d)
    weights = computeweights(mdl)
    return find_zero(y -> btgdistribution(mdl, x, y, weights) - q, startpoint)
end

"""
    predictmedian(mdl, x)

The value at location `x` at the 0.5th quantile.
"""
predictmedian(mdl::Model, x) = predictquantile(mdl, 0.5, x)

"""
    equalinterval(mdl, density, x)

The equal tailed `p`-credible interval of the value at location `x`.
"""
function equalinterval(mdl::Model, p, x)
    hw = p / 2
    lq = predictquantile(mdl, 0.5 - hw, x)
    m = predictmedian(mdl, x)
    uq = predictquantile(mdl, 0.5 + hw, x)
    return lq, m, uq
end

"""
    predictmode(mdl, x)

The value with the highest probability density at the location `x`.
"""
function predictmode(mdl::Model, x)
    # TODO
    return 1.0
end

"""
    narrowinterval(mdl, p, x)

The narrowest `p`-credible interval of the value at location `x`.
"""
function narrowinterval(mdl::Model, p, x)
    # TODO Need Derivatives for this (?)
    return 0.0, 1.0, 2.0
end

"""
    paramterestimate(mdl)

Gives MAP point estimates of the transformation parameters λ and kernel parameters θ.
"""
function parameterestimate(mdl::Model)
    # TODO Need Derivatives for this
    return (1.0,), (1.0,)
end

# TODO For all plots, consider writing Plots.jl recipes rather than plot functions

"""
    plotbtg(mdl, range, resolution)

For 1D inputs, plots the mode and narrowest intervals at the locations in `range` with
the specified `resolution`.
"""
function plotbtg(mdl::Model, p, range, resolution)
    # TODO
end

"""
    plotdensity(mdl, x, resolution)

Plots the probability density of the values at location `x` with the specified `resolution`.
"""
function plotdensity(mdl::Model, p, x, resolution)
    lq, m, uq = equalinterval(mdl, p, x)
    w = (uq - lq) * 0.1
    r = range(lq - w, stop=uq + w, length=resolution)
    weights = computeweights(mdl)
    plot(y -> btgdensity(mdl, x, y, weights), r) # TODO more output options, proper formatting
    vline!([lq, m, uq])
end

function plotdistribution(mdl::Model, p, x, resolution)
    lq, m, uq = equalinterval(mdl, p, x)
    w = (uq - lq) * 0.1
    r = range(lq - w, stop=uq + w, length=resolution)
    weights = computeweights(mdl)
    plot(y -> btgdistribution(mdl, x, y, weights), r) # TODO more output options, proper formatting
end

"""
    crossvalidate(mdl)

The LOOCV error of given BTGModel.
"""
function crossvalidate(mdl::Model)
    # TODO
    return 1.0
end

# TODO Hyperparameter grid search using crossvalidate?

export
    predictmedian,
    predictmode,
    narrowinterval,
    equalinterval,
    parameterestimate,
    plotbtg,
    plotdensity,
    plotdistribution,
    crossvalidate

end # module

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
"""
Set of functions for generating signals.
"""

function random_walk(process_noise,
                    measurement_noise,
                    mean_state0,
                    precision_state0;
                    time_horizon=100)
    "Generate data according to a random walk"

    # Dimensionality
    d = size(process_noise, 1)

    # Preallocate
    y = zeros(d, time_horizon)
    x = zeros(d, time_horizon+1)

    # Initialize state
    x[:, 1] = sqrt(inv(precision_state0))*randn(d) + mean_state0

    for t = 1:time_horizon

        # Evolve state
        x[:, t+1] = sqrt(process_noise)*randn(d) + x[:, t]

        # Observe
        y[:, t] = sqrt(measurement_noise)*randn(d) + x[:, t+1]

    end
    return y, x
end

function LGDS(A, C, Q, R, x0; T=100)
    """
    Generate data according to a linear Gaussian dynamical system.

    Note: assumes univariate observations
    """

    # Dimensionality
    Dx = size(A,1)
    Dy = size(C,1)

    # Check for Dy = 1
    if Dy != 1; error("Assumes univariate observations"); end

    # Preallocate
    x = zeros(Dx,T)
    y = zeros(T,)

    # Initialize state
    x_kmin1 = x0

    for k = 1:T

        # Evolve state
        x[:,k] = A*x_kmin1 .+ cholesky(Q).U*randn(Dx,)

        # Observe
        y[k,:] = C*x[:,k] .+ sqrt(R)*randn(Dy,)[1]

        # Update previous state
        x_kmin1 = x[:,k]

    end
    return y, x
end

function NLGDS(transition_function,
               emission_function,
               process_noise,
               measurement_noise,
               state0_params;
               time_horizon=100)
    "Generate data according to a nonlinear Gaussian dynamical system"

    # Dimensionality
    d = size(process_noise, 1)

    # Preallocate
    x = zeros(d, time_horizon+1)
    y = zeros(d, time_horizon)

    # Prior parameters
    mean0, var0 = state0_params

    # Initialize state
    x[:, 1] = sqrt(var0)*randn(d) .+ mean0

    for t = 1:time_horizon

        # Evolve state
        x[:, t+1] = sqrt(process_noise)*randn(d) .+ transition_function(x[:,t])

        # Observe
        y[:, t] = sqrt(measurement_noise)*randn(d) .+ emission_function(x[:,t+1])

    end
    return y, x
end

max_stars_repo_name: wmkouw/Bayesian-smoothers

max_issues_repo_name: wmkouw/Bayesian-smoothers

max_forks_repo_name: wmkouw/Bayesian-smoothers


text:
using InfiniteOpt, Ipopt, Plots
using KNITRO

exponential_relaxation(z, δ) = exp(one(typeof(z)) - z / δ) - one(typeof(z)) - log(δ)
function relaxed_log_barrier(z; δ = 0.3f0)
    return max(z > δ ? -log(z) : exponential_relaxation(z, δ), zero(typeof(z)))
end
# Set the SEIR parameters
γ = 0.303
β = 0.727
N = 100
ξ = 0.2

# Set the domain information
i_max = 0.06
t0 = 0
tf = 50

# Set the intial condition values
e0 = 1 / N
i0 = 0
r0 = 0
s0 = 1 - 1 / N

model = InfiniteModel(KNITRO.Optimizer)

@infinite_parameter(model, t ∈ [t0, tf], num_supports = 201,
    derivative_method = OrthogonalCollocation(2))

@variable(model, s ≥ 0, Infinite(t))
@variable(model, e ≥ 0, Infinite(t))
@variable(model, i ≥ 0, Infinite(t))
@variable(model, r ≥ 0, Infinite(t))
@variable(model, 0.2 ≤ u ≤ 0.8, Infinite(t), start = 0.2)

@objective(model, Min, ∫(u, t))

# Define the initial conditions
@constraint(model, s(0) == s0)
@constraint(model, e(0) == e0)
@constraint(model, i(0) == i0)
@constraint(model, r(0) == r0)

# Define the SEIR equations
@constraint(model, s_constr, ∂(s, t) == -(1 - u) * β * s * i)
@constraint(model, e_constr, ∂(e, t) == (1 - u) * β * s * i - ξ * e)
@constraint(model, i_constr, ∂(i, t) == ξ * e - γ * i)
@constraint(model, r_constr, ∂(r, t) == γ * i)

# Define the infection rate limit
@constraint(model, imax_constr, i ≤ i_max)

print(model)

optimize!(model)

r_opt = value(r, ndarray = true) * 100 # make the population fractions into percentages
s_opt = value(s, ndarray = true) * 100
i_opt = value(i, ndarray = true) * 100
e_opt = value(e, ndarray = true) * 100
u_opt = value(u)
@show obj_opt = objective_value(model)
ts = value(t)

p = plot(ts, i_opt, label = "i(t)")
plot!(ts, e_opt, label = "e(t)")
ylabel!("Pop. (%)")
xlabel!("Time (Days)")
display(p)

p1 = plot(ts, u_opt, label = "u(t)")
xlabel!("Time (Days)")
ylabel!("Distancing Ratio")
display(p1)


text:
@testset "plots_std.jl" begin
    n = 5000

    x = rand(n)

    ybm = rand.(GeneralizedExtremeValue.(x, 1.0, 0.1))
    fmbm = MaximumLikelihoodEVA(BlockMaxima(Variable("y", ybm), locationcov = [Variable("x", x)]), [0.0, 1.0, 1.0, 0.1])

    yte = rand.(GeneralizedPareto.(x, 0.1))
    fmte = MaximumLikelihoodEVA(ThresholdExceedance(Variable("y", yte), logscalecov = [Variable("x", x)]), [0.0, 1.0, 0.1])

    ys = rand(GeneralizedExtremeValue(0.0, 1.0, 0.1), n)
    fms = MaximumLikelihoodEVA(BlockMaxima(Variable("y", ys)), [0.0, 1.0, 0.1])

    @testset "standardize(y, μ, σ, ξ)" begin
        # Simple standardized values
        std̂ = Extremes.standardize(1, 0.0, 1.0, 0.1)

        @test std̂ ≈ 0.9531 atol = 0.0001

    end

    @testset "standardize(fm)" begin
        # Function returns vector with n values
        std̂ = Extremes.standardize(fmbm)

        @test length(std̂) == n

    end

    @testset "standardize(fm)" begin
        # Function returns vector with n values
        std̂ = Extremes.standardize(fmte)

        @test length(std̂) == n

    end

    @testset "standarddist(fm)" begin
        # BlockMaxima standard distribution is Gumbel
        dist = Extremes.standarddist(fmbm.model)

        @test dist == Gumbel()

    end

    @testset "standarddist(fm)" begin
        # ThresholdExceedance standard distribution is Exponential
        dist = Extremes.standarddist(fmte.model)

        @test dist == Exponential()

    end

    @testset "probplot_std_data(fm)" begin
        dfbm = probplot_std_data(fmbm)
        dfte = probplot_std_data(fmte)

        # Returns a dataframe with n values in column Model with BlockMaxima
        @test length(dfbm[:, :Model]) == n

        # Returns a dataframe with n values in column Empirical with BlockMaxima
        @test length(dfbm[:, :Empirical]) == n

        # Returns a dataframe with n values in column Model with ThresholdExceedance
        @test length(dfte[:, :Model]) == n

        # Returns a dataframe with n values in column Empirical with ThresholdExceedance
        @test length(dfte[:, :Empirical]) == n

    end

    @testset "qqplot_std_data(fm)" begin
        dfbm = qqplot_std_data(fmbm)
        dfte = qqplot_std_data(fmte)

        # Returns a dataframe with n values in column Model with BlockMaxima
        @test length(dfbm[:, :Model]) == n

        # Returns a dataframe with n values in column Empirical with BlockMaxima
        @test length(dfbm[:, :Empirical]) == n

        # Returns a dataframe with n values in column Model with ThresholdExceedance
        @test length(dfte[:, :Model]) == n

        # Returns a dataframe with n values in column Empirical with ThresholdExceedance
        @test length(dfte[:, :Empirical]) == n

    end

end


text:
#=
Su Doku (Japanese meaning number place) is the name given to a popular puzzle concept.
Its origin is unclear, but credit must be attributed to Leonhard Euler who invented a similar, and much more difficult,
puzzle idea called Latin Squares.

The objective of Su Doku puzzles, however, is to replace the blanks (or zeros)
in a 9 by 9 grid in such that each row, column, and 3 by 3 box contains each of the digits 1 to 9.
Below is an example of a typical starting puzzle grid and its solution grid.

[ Two Puzzles Here]

A well constructed Su Doku puzzle has a unique solution and can be solved by logic,
although it may be necessary to employ "guess and test" methods in order to eliminate options (there is much contested opinion over this).
The complexity of the search determines the difficulty of the puzzle;
the example above is considered easy because it can be solved by straight forward direct deduction.

The 6K text file, sudoku.txt (right click and 'Save Link/Target As...'),
contains fifty different Su Doku puzzles ranging in difficulty, but all with unique solutions (the first puzzle in the file is the example above).

By solving all fifty puzzles find the sum of the 3-digit numbers found in the top left corner of each solution grid;
for example, 483 is the 3-digit number found in the top left corner of the solution grid above.
=#
include("projecteulerutils.jl")
using DelimitedFiles

function valuesinrow(grid, row)
  a = Array{Int64, 1}()
  for x in 1:9
    if grid[row, x] != 0 push!(a, grid[row, x]) end
  end
  a
end

function valuesincolumn(grid, column)
  a = Array{Int64, 1}()
  for y in 1:9
    if grid[y, column] != 0 push!(a, grid[y, column]) end
  end
  a
end

function valuesincurrentcell(grid, x, y)
  result = Array{Int64, 1}()
  startx = round(Int64, (3 * ceil(x / 3)) - 2)
  starty = round(Int64, (3 * ceil(y / 3)) - 2)
  for x in startx:startx+2, y in starty:starty+2
    if grid[y, x] != 0 push!(result, grid[y, x]) end
  end
  result
end

function emptycell(grid)
  for x in 1:9, y in 1:9
    if grid[y, x] == 0 return x, y end
  end
end

function solve(grid)
  if all(x -> x != 0, grid) # grid is full
    return true
  else
    x, y = emptycell(grid)
    for value in setdiff(collect(1:9), append!([], [valuesincolumn(grid, x); valuesinrow(grid, y); valuesincurrentcell(grid, x, y)])) # valid values for cell
      grid[y, x] = value
      if solve(grid)
        return true
      end
      grid[y, x] = 0
    end
  end
  false
end

function calc()
  x, total = readdlm(dirname(@__FILE__()) * "/../Resources/p96.txt", ',', Int64), 0
  for i in 1:9:450
    grid = x[i:i+8, 1:9]; solve(grid)
    total += nd(grid[1, 1:3])
  end
  total
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
cd(@__DIR__)
using Pkg
Pkg.activate(".")
Pkg.instantiate()
using JSON, Statistics, GitHub, Base64, SHA, Downloads, Dates, CairoMakie

include("benchmark-library.jl")

ctx = try
     github_context()
catch e
    @warn "Not authorized" exception=e
    # bad credentials because PR isn't from a contributor
    exit()
end

project = "current-pr-project/"
# It seems, that between julia versions, the manifest must be deleted to not get problems
isdir(project) && rm(project; force=true, recursive=true)
mkdir(project)

Pkg.activate(project)
Pkg.develop([(;path="../../MakieCore"), (;path="../../"), (;path="../../CairoMakie")])
this_pr = BenchInfo(
    project=project,
    branch="current-pr",
    commit=current_commit()
)

Pkg.activate(".")

fig = run_benchmarks(ctx, [GitHub.branch(ctx.repo, "master"), this_pr]; force=true, n=5)


text:
using Polynomials, Roots, JLD, PyPlot
using JFVM
include("peteng-master/functions/rel_perms_real.jl")

pce = 1e6
swc = 0.12
sor = 0.08

figure()
sw_plot = collect(linspace(0,1,1000))
plot(sw_plot, pc_drain.(sw_plot, pce, swc), sw_plot, pc_imb.(sw_plot, pce, swc, sor))
title("Capillary Pressure, "L"$p_c$"" [Pa]")
xlabel(L"$S_w$"" [-]")
ylabel(L"$p_c$"" [Pa]")
legend(["Drainage", "Imbibition"])
savefig("Capillary Pressure/method1.png")

figure()
sw_plot = collect(linspace(0,1,1000))
plot(sw_plot, pc_drain2.(sw_plot, pce, swc), sw_plot, pc_imb2.(sw_plot, pce, swc, sor))
title("Capillary Pressure, "L"$p_c$"" [Pa]")
xlabel(L"$S_w$"" [-]")
ylabel(L"$p_c$"" [Pa]")
legend(["Drainage", "Imbibition"])
savefig("Capillary Pressure/method2.png")

figure()
sw_plot = collect(linspace(0,1,1000))
plot(sw_plot, pc_drain3.(sw_plot, pce, swc), sw_plot, pc_imb3.(sw_plot, pce, swc, sor))
title("Capillary Pressure, "L"$p_c$"" [Pa]")
xlabel(L"$S_w$"" [-]")
ylabel(L"$p_c$"" [Pa]")
legend(["Drainage", "Imbibition"])
savefig("Capillary Pressure/method3.png")

#Van-Genuchten
a_drain = 0.066
a_imb = 0.6
m_drain = 0.6
m_imb = 0.5
M_drain = 0.26
M_imb = 0.2

k = 1e-12
ϕ = 0.4
sigma = 25.0*1000 #N/m
theta = 0.785

J = (S, S_wc, S_or, a, m, M) -> a.*(((S-S_wc)./(1.0-S_or-S_wc)).^(-1.0./m)-1.0).^(1.0-m)-M
p_c = (S, S_wc, S_or, a, m, M) -> J.(S, S_wc, S_or, a, m, M).*sigma.*cos(theta).*sqrt(ϕ/k)

sw = collect(linspace(swc+0.02,1.0-sor,1000))

figure()
plot(sw, p_c.(sw, swc, sor, a_drain, m_drain, M_drain), sw, p_c.(sw, swc, sor, a_imb, m_imb, M_imb))
title("Capillary Pressure, "L"$p_c$"" [Pa]")
xlabel(L"$S_w$"" [-]")
ylabel(L"$p_c$"" [Pa]")
legend(["Drainage", "Imbibition"])
savefig("Capillary Pressure/Van_Genuchten.png")

max_stars_repo_name: FrederikLehn/coupled-solvers

max_issues_repo_name: FrederikLehn/coupled-solvers

max_forks_repo_name: FrederikLehn/coupled-solvers


text:
function expr_arr_to_block(exprs)
  block = :(begin end)
  foreach(expr -> push!(block.args, expr), exprs)
  block
end

# used in parsing
isblock(x) = length(x) == 1 && x[1] isa Expr && x[1].head == :block
function flatten_expr!(x)
    isb = isblock(x)
    if isb
        x = x[1]
        filter!(z->z isa Symbol || z.head != :line, x.args)
        x = (x.args...)
    end
    x
end


text:
# ---
# title: 91. Decode Ways
# id: problem91
# author: Indig
# date: 2021-06-02
# difficulty: Medium
# categories: String, Dynamic Programming
# link: <https://leetcode.com/problems/decode-ways/description/>
# hidden: true
# ---
#
# A message containing letters from `A-Z` is being encoded to numbers using the
# following mapping:
#
#
#
#     'A' -> 1
#     'B' -> 2
#     ...
#     'Z' -> 26
#
#
# Given a **non-empty** string containing only digits, determine the total
# number of ways to decode it.
#
# The answer is guaranteed to fit in a **32-bit** integer.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "12"
#     Output: 2
#     Explanation: It could be decoded as "AB" (1 2) or "L" (12).
#
#
# **Example 2:**
#
#
#
#     Input: s = "226"
#     Output: 3
#     Explanation: It could be decoded as "BZ" (2 26), "VF" (22 6), or "BBF" (2 2 6).
#
#
# **Example 3:**
#
#
#
#     Input: s = "0"
#     Output: 0
#     Explanation: There is no character that is mapped to a number starting with '0'. We cannot ignore a zero when we face it while decoding. So, each '0' should be part of "10" --> 'J' or "20" --> 'T'.
#
#
# **Example 4:**
#
#
#
#     Input: s = "1"
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 100`
#   * `s` contains only digits and may contain leading zero(s).
#
#
## @lc code=start
using LeetCode

function num_decodings(s)
    len = length(s)
    dp = fill(0, len + 1)
    dp[1] = 1
    for i in 1:len
        s[i] != '0' && (dp[i + 1] += dp[i])
        (i > 1 && s[i - 1] != '0' && parse(Int, s[i-1:i]) <= 26) && (dp[i + 1] += dp[i - 1])
    end
    dp[end]
end
## @lc code=end


text:
type CuModule
    ptr::Ptr{Void}

    function CuModule(ptr)
        m = new(ptr)
        finalizer(m, cuModuleUnload)
        m
    end
end

function CuModule(image::Vector{UInt8})
    p = Ptr{Void}[0]
    cuModuleLoadData(p, image)
    #cuModuleLoadDataEx(p, image, 0, CUjit_option[], Ptr{Void}[])
    CuModule(p[1])
end

Base.unsafe_convert(::Type{Ptr{Void}}, m::CuModule) = m.ptr


text:
# This script provides an example of simulating changepoint data,
# using PELT to find these changepoints, and plotting the results
#
# Requires Gadfly plotting package

using Changepoints, Distributions
using Gadfly

n = 1000           # Sample size
λ = 100            # Frequencey of changes

# Generate a sample with changing
μ, σ = 1.0, Uniform(2.0, 15.0)
y, cps = @changepoint_sampler n λ Normal(μ, σ)
pelt_cps, cost = @PELT y Normal(μ, ?)

println(pelt_cps)
println(cps)

p = plot(y, pelt_cps)


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

struct __enum_Syntax <: ProtoEnum
    SYNTAX_PROTO2::Int32
    SYNTAX_PROTO3::Int32
    __enum_Syntax() = new(0,1)
end #struct __enum_Syntax
const Syntax = __enum_Syntax()

mutable struct Option <: ProtoType
    name::AbstractString
    value::_Any
    Option(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Option

mutable struct EnumValue <: ProtoType
    name::AbstractString
    number::Int32
    options::Base.Vector{Option}
    EnumValue(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct EnumValue

struct __enum_Field_Kind <: ProtoEnum
    TYPE_UNKNOWN::Int32
    TYPE_DOUBLE::Int32
    TYPE_FLOAT::Int32
    TYPE_INT64::Int32
    TYPE_UINT64::Int32
    TYPE_INT32::Int32
    TYPE_FIXED64::Int32
    TYPE_FIXED32::Int32
    TYPE_BOOL::Int32
    TYPE_STRING::Int32
    TYPE_GROUP::Int32
    TYPE_MESSAGE::Int32
    TYPE_BYTES::Int32
    TYPE_UINT32::Int32
    TYPE_ENUM::Int32
    TYPE_SFIXED32::Int32
    TYPE_SFIXED64::Int32
    TYPE_SINT32::Int32
    TYPE_SINT64::Int32
    __enum_Field_Kind() = new(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18)
end #struct __enum_Field_Kind
const Field_Kind = __enum_Field_Kind()

struct __enum_Field_Cardinality <: ProtoEnum
    CARDINALITY_UNKNOWN::Int32
    CARDINALITY_OPTIONAL::Int32
    CARDINALITY_REQUIRED::Int32
    CARDINALITY_REPEATED::Int32
    __enum_Field_Cardinality() = new(0,1,2,3)
end #struct __enum_Field_Cardinality
const Field_Cardinality = __enum_Field_Cardinality()

mutable struct Field <: ProtoType
    kind::Int32
    cardinality::Int32
    number::Int32
    name::AbstractString
    type_url::AbstractString
    oneof_index::Int32
    packed::Bool
    options::Base.Vector{Option}
    json_name::AbstractString
    default_value::AbstractString
    Field(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Field
const __fnum_Field = Int[1,2,3,4,6,7,8,9,10,11]
meta(t::Type{Field}) = meta(t, ProtoBuf.DEF_REQ, __fnum_Field, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct _Enum <: ProtoType
    name::AbstractString
    enumvalue::Base.Vector{EnumValue}
    options::Base.Vector{Option}
    source_context::SourceContext
    syntax::Int32
    _Enum(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct _Enum

mutable struct _Type <: ProtoType
    name::AbstractString
    fields::Base.Vector{Field}
    oneofs::Base.Vector{AbstractString}
    options::Base.Vector{Option}
    source_context::SourceContext
    syntax::Int32
    _Type(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct _Type

export Syntax, _Type, Field_Kind, Field_Cardinality, Field, _Enum, EnumValue, Option


text:
using Pkg

if VERSION >= v"1.5"
    Pkg.develop(path="../SwapLiterals")
end
Pkg.test("SwapLiterals")


text:
using Documenter, MultilevelEstimators

makedocs(modules = [MultilevelEstimators],
		 sitename = "MultilevelEstimators.jl",
         format = Documenter.HTML(
                                  prettyurls = false,#!("local" in ARGS),
                                 ),
         authors = "Pieterjan Robbe",
         pages = Any[
                  "Home" => "index.md",
                  "Example" => "example.md",
                  "Manual" => "manual.md",#["Index" => "index.md",
                              # "IndexSet" => "index_set.md",
                              # "SampleMethod" => "sample_method.md",
                              # "Distribution" => "distribution.md",
                              # "Estimator" => "estimator.md",
                              # "History" => "history.md"
                              #]
                 ]
         )

deploydocs(
    repo = "github.com/PieterjanRobbe/MultilevelEstimators.jl.git",
)


text:
@enum OpCode::UInt8 OpConstant OpAdd OpSub OpMul OpDiv OpEqual OpNotEqual OpLessThan OpGreaterThan OpMinus OpBang OpTrue OpFalse OpNull OpPop OpJump OpJumpNotTruthy OpGetGlobal OpSetGlobal OpGetLocal OpSetLocal OpGetBuiltin OpGetFree OpSetFree OpGetOuter OpSetOuter OpCurrentClosure OpArray OpHash OpClosure OpIndex OpCall OpBreak OpContinue OpReturnValue OpReturn OpIllegal

struct Instructions
    codes::Vector{UInt8}
end

Base.length(ins::Instructions) = length(ins.codes)
Base.getindex(ins::Instructions, i::Int) = ins.codes[i]
Base.getindex(ins::Instructions, r::UnitRange{Int}) = Instructions(ins.codes[r])
Base.setindex!(ins::Instructions, val::UInt8, i::Int) = ins.codes[i] = val
Base.lastindex(ins::Instructions) = lastindex(ins.codes)
Base.vcat(is::Vararg{Instructions}) = Instructions(vcat(map(x -> x.codes, is)...))
Base.append!(ins::Instructions, other::Instructions) = append!(ins.codes, other.codes)
Base.splice!(ins::Instructions, r::UnitRange{Int}) = splice!(ins.codes, r)
Base.string(ins::Instructions) = begin
    out = IOBuffer()
    i = 1
    while i <= length(ins)
        def = lookup(ins[i])
        if isnothing(def)
            write(out, "ERROR: unknown opcode: $(ins[i])\n")
            i += 1
            continue
        end

        operands, n = read_operands(def, ins[i+1:end])
        write(out, @sprintf("%04d %s\n", i - 1, format_instruction(def, operands)))
        i += n + 1
    end

    return String(take!(out))
end

struct Definition
    name::String
    operand_widths::Vector{Int}
end

const OPCODE_STRINGS = Dict{OpCode,String}(
    OpAdd => "+",
    OpSub => "-",
    OpMul => "*",
    OpDiv => "/",
    OpEqual => "==",
    OpNotEqual => "!=",
    OpLessThan => "<",
    OpGreaterThan => ">",
    OpMinus => "-",
    OpBang => "!",
)

const DEFINITIONS = Dict{OpCode,Definition}(
    OpConstant => Definition("OpConstant", [2]),
    OpAdd => Definition("OpAdd", []),
    OpSub => Definition("OpSub", []),
    OpMul => Definition("OpMul", []),
    OpDiv => Definition("OpDiv", []),
    OpEqual => Definition("OpEqual", []),
    OpNotEqual => Definition("OpNotEqual", []),
    OpLessThan => Definition("OpLessThan", []),
    OpGreaterThan => Definition("OpGreaterThan", []),
    OpMinus => Definition("OpMinus", []),
    OpBang => Definition("OpBang", []),
    OpTrue => Definition("OpTrue", []),
    OpFalse => Definition("OpFalse", []),
    OpNull => Definition("OpNull", []),
    OpPop => Definition("OpPop", []),
    OpJump => Definition("OpJump", [2]),
    OpJumpNotTruthy => Definition("OpJumpNotTruthy", [2]),
    OpGetGlobal => Definition("OpGetGlobal", [2]),
    OpSetGlobal => Definition("OpSetGlobal", [2]),
    OpGetLocal => Definition("OpGetLocal", [1]),
    OpSetLocal => Definition("OpSetLocal", [1]),
    OpGetBuiltin => Definition("OpGetBuiltin", [1]),
    OpGetFree => Definition("OpGetFree", [1]),
    OpSetFree => Definition("OpSetFree", [1]),
    OpGetOuter => Definition("OpGetOuter", [1, 1, 1]),
    OpSetOuter => Definition("OpSetOuter", [1, 1, 1]),
    OpCurrentClosure => Definition("OpCurrentClosure", []),
    OpArray => Definition("OpArray", [2]),
    OpHash => Definition("OpHash", [2]),
    OpClosure => Definition("OpClosure", [2, 1]),
    OpIndex => Definition("OpIndex", []),
    OpCall => Definition("OpCall", [1]),
    OpBreak => Definition("OpBreak", []),
    OpContinue => Definition("OpContinue", []),
    OpReturnValue => Definition("OpReturnValue", []),
    OpReturn => Definition("OpReturn", []),
)

lookup(op::UInt8) = Base.get(DEFINITIONS, OpCode(op), nothing)

function format_instruction(def::Definition, operands::Vector{Int})
    operand_count = length(operands)

    if operand_count != length(def.operand_widths)
        return "ERROR: operand length mismatch. Expected $(length(def.operand_widths)), got $operand_count instead.\n"
    end

    return def.name * join(map(x -> " " * string(x), operands))
end

function make(op::OpCode, operands::Vararg{Int})::Instructions
    if op ∉ keys(DEFINITIONS)
        return Instructions(UInt8[])
    end

    def = DEFINITIONS[op]
    codes = UInt8[Integer(op)]

    for (operand, width) in zip(operands, def.operand_widths)
        if width == 2
            push!(codes, UInt8(operand >> 8))
            push!(codes, UInt8(operand & 0xFF))
        elseif width == 1
            push!(codes, UInt8(operand))
        end
    end

    return Instructions(codes)
end

function read_operands(def::Definition, ins::Instructions)::Tuple{Vector{Int},Int}
    operands = Int[]
    offset = 0
    for width in def.operand_widths
        if width == 2
            push!(operands, read_uint16(ins, offset + 1))
        elseif width == 1
            push!(operands, ins[offset+1])
        end
        offset += width
    end
    return operands, offset
end

function read_uint16(ins::Instructions, pos::Int)::UInt16
    return UInt16(ins.codes[pos]) << 8 | ins.codes[pos+1]
end


text:
using LibGEOS
using Plots
gr() # 図が早くて綺麗に書けるおまじない

mutable struct MYPolygon{T}
    vertexes::Matrix{T}

    n::Int64
    center::Vector{T}
    area::Float64
end

# 多角形の重心
function calc_center(X, Y)
    A = zero(X[1])
    for i in 1:length(X)-1
        A += X[i] * Y[i+1] - X[i+1] * Y[i]
    end
    A += X[end] * Y[1] - X[1] * Y[end]
    A /= 2

    Cx = Cy = 0
    for i in 1:length(X)-1
        Cx += (X[i] + X[i+1]) * (X[i] * Y[i+1] - X[i+1] * Y[i])
        Cy += (Y[i] + Y[i+1]) * (X[i] * Y[i+1] - X[i+1] * Y[i])
    end
    Cx += (X[end] + X[1]) * (X[end] * Y[1] - X[1] * Y[end])
    Cy += (Y[end] + Y[1]) * (X[end] * Y[1] - X[1] * Y[end])

    return [Cx/(6A), Cy/(6A)]
end

# Figure の宣言用関数
# X x座標の配列, Y y座標の配列
function MYPolygon(X::Array, Y::Array, T=Float64)
    MP = MYPolygon{T}(
        hcat(X, Y)',
        length(X),
        calc_center(X, Y),
        LibGEOS.area(vertex2LibGEOS(X, Y))
    )
    return move(MP, -MP.center...)
end

# M (N*2)行列
function MYPolygon(M::Matrix, T=Float64)
    return MYPolygon(M[1, :], M[2, :], T)
end

function Base.copy(P::MYPolygon{T}) where T
    return MYPolygon{T}(P.vertexes, P.n, P.center)
end

# Polygon を (x, y) ベクトルだけ動かす
function move!(P::MYPolygon{T}, x::T, y::T) where T
    P.vertexes .+= [x, y]
    P.center += [x, y]
    return nothing
end

function move(P::MYPolygon{T}, x::T, y::T) where T
    newP = deepcopy(P)
    move!(newP, x, y)
    return newP
end

# Polygon を θ だけ回転させる (重心中心)
function rotate!(P::MYPolygon{T}, θ::T) where T
    R = [cos(θ) -sin(θ); sin(θ) cos(θ)]

    P.vertexes.-= P.center
    P.vertexes = R * P.vertexes
    P.vertexes .+= P.center
    return nothing
end
# Polygon を θ だけ回転させる (重心中心)
function rotate(P::MYPolygon{T}, θ::T) where T
    newP = deepcopy(P)
    rotate!(newP, θ)
    return newP
end

# Polygon を θ だけ回転させる (頂点中心)
function rotate!(P::MYPolygon{T}, θ::T, n::Int) where T
    if n > P.n
        throw(DomainError(n, "id of vertex must be less or equal than the number of vertexes of P"))
    end

    R = [cos(θ) -sin(θ); sin(θ) cos(θ)]

    P.vertexes.-= P.vertexes[:, n]
    P.vertexes = R * P.vertexes
    P.vertexes .+= P.vertexes[:, n]
    return nothing
end
# Polygon を θ だけ回転させる (頂点中心)
function rotate(P::MYPolygon{T}, θ::T, n::Int) where T
    newP = deepcopy(P)
    rotate!(newP, θ, n)
    return newP
end

# 2次元ベクトルの外積
function cross(p::Array{T, 1}, q::Array{T, 1}) where T
    return p[1] * q[2] - q[1] * p[2]
end

# 線分ABと点Pの距離
function distance_seg2point(A::Array{Float64, 1}, B::Array{Float64, 1}, P::Array{Float64, 1})
    a = B[1] - A[1]
    b = B[2] - A[2]
    if a == 0 && b == 0
        return norm(A - P)
    else
        Q = zeros(Float64, 2)
        Q[1] = (a^2 * P[1] + b^2 * A[1] + a * b * (P[2] - A[2])) / (a^2 + b^2)
        Q[2] = (a^2 * A[2] + b^2 * P[2] + a * b * (P[1] - A[1])) / (a^2 + b^2)

        t = 0
        if a != 0
            t = (B[1] - Q[1]) / a
        else
            t = (B[2] - Q[2]) / b
        end

        if 0 ≤ t ≤ 1
            return norm(Q - P)
        else
            return min(norm(A - P), norm(B - P))
        end
    end
end

# 図形を倍率dだけ拡大する
function expand!(P::MYPolygon{T}, d::T) where T
    for i in 1:P.n
        # 重心から頂点へのベクトル
        v = P.vertexes[:, i] - P.center
        P.vertexes[:, i] = v * d + P.center
    end
end

function expand(P::MYPolygon{T}, d::T) where T
    newP = deepcopy(P)
    expand!(newP, d)
    return newP
end

# 図形を描画する
function display(Polygons::MYPolygon...; center=false, vertex=false)
    xmax = ymax = -floatmax(Float64)
    xmin = ymin = floatmax(Float64)

    pl = plot(size=(400, 400))

    for P in Polygons
        # 表示範囲
        xmax = maximum(P.vertexes[1, :]) > xmax ? maximum(P.vertexes[1, :]) : xmax
        xmin = minimum(P.vertexes[1, :]) < xmin ? minimum(P.vertexes[1, :]) : xmin
        ymax = maximum(P.vertexes[2, :]) > ymax ? maximum(P.vertexes[2, :]) : ymax
        ymin = minimum(P.vertexes[2, :]) < ymin ? minimum(P.vertexes[2, :]) : ymin

        # 図形が閉じるように終点を追加
        V = hcat(P.vertexes, P.vertexes[:, 1])
        plot!(pl, V[1, :], V[2, :], fill=0)
        if center != false
            scatter!(pl, [P.center[1]], [P.center[2]], marker=center)
        end
        if vertex != false
            scatter!(pl, V[1, :], V[2, :], marker=vertex)
        end
    end
    if xmax - xmin > ymax - ymin
        xmax += (xmax - xmin)/2
        xmin -= (xmax - xmin)/2

        ycent = (ymax + ymin)/2
        ymax = ycent + (xmax - xmin)/2
        ymin = ycent - (xmax - xmin)/2
    else
        ymax += (ymax - ymin)/2
        ymin -= (ymax - ymin)/2

        xcent = (xmax + xmin)/2
        xmax = xcent + (ymax - ymin)/2
        xmin = xcent - (ymax - ymin)/2
    end

    # plot!(pl, xlim=(xmin, xmax), ylim=(ymin, ymax))
    plot!(pl, xlim=(-3, 3), ylim=(-3, 3), aspect_ratio =  1.0, legend=false)
    return pl
end

function vertex2LibGEOS(X, Y)
    str = "POLYGON(("
    for i in 1:length(X)
        str *= string(X[i])
        str *= " "
        str *= string(Y[i])
        str *= ","
    end
    str *= string(X[1])
    str *= " "
    str *= string(Y[1])
    str *= "))"
    return readgeom(str)
end

function MYPolygon2LibGEOS(P::MYPolygon)
    str = "POLYGON(("
    for i in 1:P.n
        str *= string(P.vertexes[1, i])
        str *= " "
        str *= string(P.vertexes[2, i])
        str *= ","
    end
    str *= string(P.vertexes[1, 1])
    str *= " "
    str *= string(P.vertexes[2, 1])
    str *= "))"
    return readgeom(str)
end

function make_hole_GEOS(P::MYPolygon, H::MYPolygon)
    str = "POLYGON(("
    for i in 1:P.n
        str *= string(P.vertexes[1, i])
        str *= " "
        str *= string(P.vertexes[2, i])
        str *= ","
    end
    str *= string(P.vertexes[1, 1])
    str *= " "
    str *= string(P.vertexes[2, 1])

    str *= "),("

    for i in H.n:-1:1
        str *= string(H.vertexes[1, i])
        str *= " "
        str *= string(H.vertexes[2, i])
        str *= ","
    end
    str *= string(H.vertexes[1, end])
    str *= " "
    str *= string(H.vertexes[2, end])
    str *= "))"
    return readgeom(str)
end

# 図形と図形の共通部分
function intersectionArea(P1::MYPolygon, P2::MYPolygon)
    geP1 = MYPolygon2LibGEOS(P1)
    geP2 = MYPolygon2LibGEOS(P2)
    gePi = LibGEOS.intersection(geP1, geP2)
    return LibGEOS.area(gePi)
end

function unionArea(P1::MYPolygon, P2::MYPolygon)
    geP1 = MYPolygon2LibGEOS(P1)
    geP2 = MYPolygon2LibGEOS(P2)
    gePu = LibGEOS.union(geP1, geP2)
    return LibGEOS.area(gePu)
end

function SymmetricDifference(P1::MYPolygon, P2::MYPolygon)
    geP1 = MYPolygon2LibGEOS(P1)
    geP2 = MYPolygon2LibGEOS(P2)
    geSD = LibGEOS.symmetricDifference(geP1, geP2)

    return LibGEOS.area(geSD)
end

max_stars_repo_name: toshiki-kosen/SolvingTangramPuzzle

max_issues_repo_name: toshiki-kosen/SolvingTangramPuzzle

max_forks_repo_name: toshiki-kosen/SolvingTangramPuzzle


text:
module PolygonIO

# Gloabl Imports
using JSON3
using JSONTables
using HTTP

######### Import modules & utils ################
include("custom_structs.jl")
include("utils.jl")
include("urls.jl")
include("crypto_api.jl")
include("forex_api.jl")
include("reference_api.jl")
include("stocks_api.jl")
include("streaming_socket.jl")

######### Global export of user API  ################
    # Reference API
export PolyOpts, tickers, ticker_types, ticker_details, ticker_details_vX,
        ticker_news, markets, locales, stock_splits, stock_dividends,
        stock_financials, market_holidays, market_status, stock_exchanges,
        condition_mappings, crypto_exchanges,

    # Stock API
        stocks_trades, stocks_quotes_nbbo, stocks_last_trade_symbol, stocks_last_quote_symbol,
        stocks_daily_open_close, stocks_grouped_daily_bars, stocks_previous_close,
        stocks_aggregates_bars, stocks_snapshot_all_tickers, stocks_snapshot_ticker,
        stocks_snapshot_gainers_losers,

    # Crypto API
        last_trade_crypto_pair, crypto_daily_open_close, historic_crypto_trades,
        crypto_grouped_daily_bars, crypto_previous_close, crypto_aggregates_bars,
        crypto_snapshot_all_tickers, crypto_snapshot_ticker, crypto_snapshot_ticker_full_book,
        crypto_snapshot_gainers_losers,

    # Forex API
        historic_forex_ticks, real_time_currency_conversion, last_quote_currency_pair,
        forex_grouped_daily_bars, forex_previous_close, forex_aggregates_bars,
        forex_snapshot_ticker, forex_snapshot_all_tickers, forex_snapshot_gainers_losers

end


text:
struct OrnsteinUhlenbeck{T1,T2,T3}
  Θ::T1
  μ::T2
  σ::T3
end
# http://www.math.ku.dk/~susanne/StatDiff/Overheads1b.pdf
function (X::OrnsteinUhlenbeck)(dW,W,dt,u,p,t,rng) #dist
  if typeof(dW) <: AbstractArray
    rand_val = wiener_randn(rng,dW)
  else
    rand_val = wiener_randn(rng,typeof(dW))
  end
  drift = X.μ .+ (W[end] .- X.μ) .* exp.(-X.Θ*dt)
  diffusion = X.σ .* sqrt.((1 .- exp.(-2X.Θ*dt))./(2X.Θ))
  drift .+ rand_val .* diffusion .- W[end]
end

#=
http://www.tandfonline.com/doi/pdf/10.1080/14697688.2014.941913?needAccess=true
http://www.tandfonline.com/doi/full/10.1080/14697688.2014.941913?src=recsys
=#

#=
(qX + r) = Θ(μ-x) = Θμ - Θx
q = -Θ
r = Θμ
https://arxiv.org/pdf/1011.0067.pdf page 18
=#
function ou_bridge(dW,ou,W,W0,Wh,q,hu,p,t,rng) end
function ou_bridge!(rand_vec,ou,W,W0,Wh,q,h,u,p,t,rng) end

function OrnsteinUhlenbeckProcess(Θ,μ,σ,t0,W0,Z0=nothing;kwargs...)
  ou = OrnsteinUhlenbeck(Θ,μ,σ)
  NoiseProcess{false}(t0,W0,Z0,ou,nothing;kwargs...)
end

struct OrnsteinUhlenbeck!{T1,T2,T3}
  Θ::T1
  μ::T2
  σ::T3
end

function (X::OrnsteinUhlenbeck!)(rand_vec,W,dt,u,p,t,rng) #dist!
  wiener_randn!(rng,rand_vec)
  @.. rand_vec = X.μ+(W[end]-X.μ)*exp(-X.Θ*dt) + rand_vec*X.σ*sqrt((1-exp.(-2*X.Θ.*dt))/(2*X.Θ)) - W[end]
end
function OrnsteinUhlenbeckProcess!(Θ,μ,σ,t0,W0,Z0=nothing;kwargs...)
  ou = OrnsteinUhlenbeck!(Θ,μ,σ)
  NoiseProcess{true}(t0,W0,Z0,ou,nothing;kwargs...)
end


text:
@testset "Random spatial vector" begin
    for i = 1:10
        @test norm(rand(SpatialVector{Float64})) ≈ 1.
    end;
end;


text:
"""
Helper infrastructure to compile and sample models using `cmdstan`.

[`StanModel`](@ref) wraps a model definition (source code), while [`stan_sample`](@ref) can
be used to sample from it.

[`stan_compile`](@ref) can be used to pre-compile a model without sampling. A
[`StanModelError`](@ref) is thrown if this fails, which contains the error messages from
`stanc`.
"""
module StanBase

using Reexport

@reexport using Unicode, DelimitedFiles, Distributed
@reexport using StanDump, StanRun, StanSamples
@reexport using MCMCChains
@reexport using Parameters

using DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF

import StanRun: stan_sample, stan_cmd_and_paths, default_output_base
import StanSamples: read_samples

include("utils/findall.jl")
include("stanmodel/shared_fields.jl")
include("stanmodel/top_level_types.jl")
include("stanmodel/help_types.jl")
include("stanmodel/HelpModel.jl")
include("stanmodel/update_model_file.jl")
include("stanmodel/number_of_chains.jl")
include("stanrun/cmdline.jl")
include("stanrun/stan_sample.jl")
include("stansamples/stan_summary.jl")
include("stansamples/read_summary.jl")

"""
The directory which contains the cmdstan executables such as `bin/stanc` and
`bin/stansummary`. Inferred from the environment variable `JULIA_CMDSTAN_HOME` or `ENV["JULIA_CMDSTAN_HOME"]`
when available.

If these are not available, use `set_cmdstan_home!` to set the value of CMDSTAN_HOME.

Example: `set_cmdstan_home!(homedir() * "/Projects/Stan/cmdstan/")`

Executing `versioninfo()` will display the value of `JULIA_CMDSTAN_HOME` if defined.
"""
CMDSTAN_HOME=""

function __init__()
  global CMDSTAN_HOME = if isdefined(Main, :JULIA_CMDSTAN_HOME)
    Main.JULIA_CMDSTAN_HOME
  elseif haskey(ENV, "JULIA_CMDSTAN_HOME")
    ENV["JULIA_CMDSTAN_HOME"]
  elseif haskey(ENV, "CMDSTAN_HOME")
    ENV["CMDSTAN_HOME"]
  else
    @warn("Environment variable CMDSTAN_HOME not set. Use set_cmdstan_home!.")
    ""
  end
end

"""Set the path for the `CMDSTAN_HOME` environment variable.

Example: `set_cmdstan_home!(homedir() * "/Projects/Stan/cmdstan/")`
"""
set_cmdstan_home!(path) = global CMDSTAN_HOME = path

const src_path = @__DIR__

"""
# `rel_path_stanbase`

Relative path using the StanBase.jl src directory. This approach has been copied from
[DynamicHMCExamples.jl](https://github.com/tpapp/DynamicHMCExamples.jl)

### Example to get access to the data subdirectory
```julia
rel_path_stanbase("..", "data")
```
"""
rel_path_stanbase(parts...) = normpath(joinpath(src_path, parts...))

stan_help = stan_sample

export
  @shared_fields_stanmodels,
  CmdStanModels,
  HelpModel,
  cmdline,
  stan_help,
  stan_sample,
  read_summary,
  stan_summary,
  set_cmdstan_home!,
  findall

end # module