text:
module BenchDMRG

using BenchmarkTools
using ITensors

suite = BenchmarkGroup()

let N = 100
  sites = siteinds("S=1", N)
  ampo = OpSum()
  for j in 1:(N - 1)
    ampo .+= ("Sz", j, "Sz", j + 1)
    ampo .+= (0.5, "S+", j, "S-", j + 1)
    ampo .+= (0.5, "S-", j, "S+", j + 1)
  end
  H = MPO(ampo, sites)
  state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
  psi = productMPS(sites, state)
  sweeps = Sweeps(5)
  maxdim!(sweeps, 10, 20, 100, 100, 200)
  cutoff!(sweeps, 1e-11)

  # Precompile
  dmrg(H, psi, sweeps; outputlevel=0)

  suite["1d_S=1_heisenberg"] = @benchmarkable dmrg($H, $psi, $sweeps; outputlevel=0)
end

let N = 100
  sites = siteinds("S=1", N; conserve_qns=true)
  ampo = OpSum()
  for j in 1:(N - 1)
    ampo .+= ("Sz", j, "Sz", j + 1)
    ampo .+= (0.5, "S+", j, "S-", j + 1)
    ampo .+= (0.5, "S-", j, "S+", j + 1)
  end
  H = MPO(ampo, sites)
  state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
  psi = productMPS(sites, state)
  sweeps = Sweeps(5)
  maxdim!(sweeps, 10, 20, 100, 100, 200)
  cutoff!(sweeps, 1E-10)

  # Precompile
  dmrg(H, psi, sweeps; outputlevel=0)

  suite["1d_S=1_heisenberg_qn"] = @benchmarkable dmrg($H, $psi, $sweeps; outputlevel=0)
end

#let Ny = 6,
#    Nx = 12
#  N = Nx * Ny
#  sites = siteinds("S=1/2", N)
#  lattice = square_lattice(Nx, Ny; yperiodic = false)
#  ampo = OpSum()
#  for b in lattice
#    ampo .+= (0.5, "S+", b.s1, "S-", b.s2)
#    ampo .+= (0.5, "S-", b.s1, "S+", b.s2)
#    ampo .+= (     "Sz", b.s1, "Sz", b.s2)
#  end
#  H = MPO(ampo, sites)
#  state = [isodd(n) ? "Up" : "Dn" for n=1:N]
#  psi = productMPS(sites, state)
#  sweeps = Sweeps(10)
#  maxdim!(sweeps, 10, 20, 100, 100, 200, 400, 800)
#  cutoff!(sweeps, 1e-8)
#
#  # Precompile
#  dmrg(H, psi, sweeps; outputlevel = 0)
#
#  suite["2d_S=1/2_heisenberg"] = @benchmarkable dmrg($H, $psi, $sweeps;
#                                                   outputlevel = 0)
#end

end

BenchDMRG.suite


text:
module TestGroupBy

using Relations
using Base.Test

name = ["Niamh", "Roger", "Genevieve", "Aiden"]
age = [27, 63, 26, 17]
eye_color = ["green", "brown", "brown", "blue"]

n = rand(10:1000)
A = rand(n)
B = rand(n)
C = rand(1:30, n)

for T in (Array, NullableArray)
    people = Relation(name = T(name), age = T(age), eye_color = T(eye_color))
    _people = copy(people)
    res = @with people groupby(eye_color, age > 26)
    @test isequal(people, _people)
    for level in res.group_levels
        for i in res.group_indices[level]
            row = res.src[i]
            for attr in fieldnames(level)
                @test isequal(getfield(row, attr), getfield(level, attr))
            end
        end
    end
    # TODO: Replace this test with something appropriate
    # @test isequal(res.groupbys, q.graph.args)

    src = Relation(A=T(A), B=T(B), C=T(C))
    _src = copy(src)
    res = @with src groupby(A > B, log(A+1) > log(A)+1, C)
    for level in res.group_levels
        for i in res.group_indices[level]
            row = res.src[i]
            for attr in fieldnames(level)
                @test isequal(getfield(row, attr), getfield(level, attr))
            end
        end
    end
end

end

ext: jl


text:
using Base.Broadcast: Broadcasted, DefaultArrayStyle

abstract type AbstractTimeSeriesStyle{N} <: Broadcast.AbstractArrayStyle{N} end

struct TimeArrayStyle{N} <: AbstractTimeSeriesStyle{N} end
TimeArrayStyle(::Val{N}) where N = TimeArrayStyle{N}()
TimeArrayStyle{M}(::Val{N}) where {N,M} = TimeArrayStyle{N}()

# Determin the output type
Base.BroadcastStyle(::Type{<:TimeArray{T,N}}) where {T,N} = TimeArrayStyle{N}()

Base.broadcastable(x::AbstractTimeSeries) = x

Base.Broadcast.instantiate(bc::Broadcasted{<:TimeArrayStyle}) =
    # skip the default axes checking
    Broadcast.flatten(bc)

function Base.copy(bc′::Broadcasted{<:TimeArrayStyle})
    tas = find_ta(bc′)

    check_column_lens(tas)

    n = length(tas)
    col′ = (n == 1) ? colnames(tas[1]) : _new_cnames.(colnames.(tas)...)
    meta′ = (n == 1) ? meta(tas[1]) : (allequal(meta.(tas)) ? meta(tas[1]) : nothing)

    # obtain shared timestamp
    tstamp_idx = overlap(timestamp.(tas)...)

    # replace TimeArray objects into Array in the Broadcasted arguments
    j = 0
    args = Any[]
    for (i, arg) ∈ enumerate(bc′.args)
        x = if arg isa TimeArray
            j += 1
            if typeof(arg).parameters[2] == 1  # 1D array
                view(values(arg), tstamp_idx[j])
            else
                view(values(arg), tstamp_idx[j], :)
            end
        else
            arg
        end
        push!(args, x)
    end

    TimeArray(view(timestamp(tas[1]), tstamp_idx[1]),
              broadcast(bc′.f, args...),
              col′,
              meta′)
end

@inline function check_column_lens(tas::Tuple)
    length(tas) <= 1 && return

    # if we have more than one TimeArray
    lens = Set(i for i ∈ map(ta -> length(colnames(ta)), tas) if i ≠ 1)
    length(lens) > 1 && throw(
        DimensionMismatch(
            "TimeArrays must have the same number of columns, " *
            "or one must be a single column"))
end

"""
find all TimeArray in the Broadcasted args and return a tuple of them.
"""
function find_ta(bc)
    ret = tuple()
    for i ∈ bc.args
        if i isa TimeArray
            ret = tuple(ret..., i)
        end
    end
    ret
end

@generated function _new_cnames(args::Vararg{Symbol, N}) where N
    expr = :(Symbol(args[1]))
    for i ∈ 2:N
        push!(expr.args, "_", :(args[$i]))
    end
    expr
end


text:
using OpenSAFT

# runtest.jl will be developed soon!


text:
#=
ev_sim_trace(arribos,demandas,salidas, potencias,C,policy)

Realiza la simulacion hasta terminar con los vehiculos de la lista.
arribos: lista de tiempos de arribo. Debe estar ordenada.
demandas: lista de demandas de carga (en energia)
salidas: lista de tiempos de salida.
potencias: lista de potencias de carga de los vehiculos
C: potencia maxima combinada
policy: una de las politicas definidas en EVSim
snapshots: tiempo de capturas
salidaReportada: vector opcional que altera los deadlines reportados.
=#
function ev_sim_trace(arribos,demandas,salidas,potencias,policy,C,snapshots; salidaReportada=nothing)

    num = length(arribos); #no. de vehiculos a procesar.
    prog = Progress(num+1, dt=0.5, desc="Simulando... ");

    eventos = 3*num+1+length(snapshots);

    #guardo parametros
    params = Dict(
        "TotalArrivals" => length(arribos),
        "AvgEnergy" => mean(demandas),
        "AvgDeadline" => mean(salidas-arribos),
        "SimTime" => salidas[end],
        "Capacity" => C,
        "Policy" => get_policy_name(policy),
        "AvgReportedDeadline" => NaN,
        "Snapshots" => length(snapshots)
    )

    if salidaReportada!=nothing
        params["AvgReportedDeadline"] = mean(salidaReportada-arribos)
    end

    #valores iniciales
    T=zeros(eventos);
    X=zeros(UInt16,length(T));   #charging vehicles
    Y=zeros(UInt16,length(T));   #already charged
    P=zeros(Float64,length(T));   #used power

    snaps = Array{Snapshot}(undef,0);

    t=0.0;
    x=0;
    y=0;
    p=0.0;

    T[1]=t;
    X[1]=x;
    Y[1]=y;
    P[1]=p;

    i=1;    #event counter
    m=1;    #snapshot counter

    charging = Array{EVinstance}(undef,0);
    alreadyCharged = Array{EVinstance}(undef,0);
    finished = Array{EVinstance}(undef,0);

    powerAllocation = Array{Float64}(undef,0);

    arrivals=0;
    expired=0;

    nextArr = arribos[arrivals+1]-t; #inicializo al primer arribo
    nextCharge = Inf;
    nextDepON = Inf;
    nextDepOFF = Inf;
    nextSnapshot = snapshots[m];

    dt,caso = findmin([nextArr;nextCharge;nextDepON;nextDepOFF;nextSnapshot])

    while dt<Inf

        t=t+dt;
        nextArr=nextArr-dt;
        nextSnapshot = nextSnapshot - dt;

        map(v->update_vehicle(v,dt),charging);
        map(v->update_vehicle(v,dt),alreadyCharged);

        if caso==1          #arribo
            arrivals=arrivals+1;
            if arrivals<num
                nextArr = arribos[arrivals+1]-t;
            else
                nextArr = Inf;
            end

            x=x+1;

            if salidaReportada==nothing
                push!(charging,EVinstance(t,salidas[arrivals],demandas[arrivals],potencias[arrivals]));
            else
                push!(charging,EVinstance(t,salidas[arrivals],salidaReportada[arrivals],demandas[arrivals],potencias[arrivals]));
            end

        elseif caso==2      #charge completed
            x=x-1;
            y=y+1;

            #guardo el auto que termina
            aux,k = findmin([ev.currentWorkload for ev in charging]); ##TODO agregar un assert de que aux==0?
            ev = charging[k];

            push!(finished,ev);
            push!(alreadyCharged,ev);
            deleteat!(charging,k)

            ev.currentPower=0.0;
            ev.departureWorkload=0.0;
            ev.completionTime=t;

        elseif caso==3      #departure without full charge
            expired=expired+1;
            x=x-1;

            #guardo el auto que termina
            aux,k = findmin([ev.currentDeadline for ev in charging]);
            ev=charging[k];

            push!(finished,ev);
            deleteat!(charging,k)

            ev.currentPower=0.0;
            ev.departureWorkload = ev.currentWorkload;
            ev.completionTime=t;

        elseif caso==4      #departure after full charge
            y=y-1;
            aux,k = findmin([ev.currentDeadline for ev in alreadyCharged]);
            deleteat!(alreadyCharged,k);

        elseif caso==5      #take snapshot

            snapshot=Snapshot(t,deepcopy(charging),deepcopy(alreadyCharged));
            push!(snaps,snapshot);

            if (m<length(snapshots))
                m=m+1;
                nextSnapshot = snapshots[m]-t;
            else
                nextSnapshot = Inf;
            end

        end

        if x>0
            #Apply the policy to all vehicles. Return total user power
            p = policy(charging,C);

            if minimum([ev.currentWorkload for ev in charging])==0
                nextCharge=0;
            else
                nextCharge = minimum([ev.currentWorkload/ev.currentPower for ev in charging]);
            end
            nextDepON = minimum([ev.currentDeadline for ev in charging]);
        else
            nextCharge = Inf;
            nextDepON = Inf;
            p=0.0;
        end

        if y>0
            nextDepOFF = minimum([ev.currentDeadline for ev in alreadyCharged]);
        else
            nextDepOFF = Inf;
        end

        i=i+1;
        T[i]=t;
        X[i]=x;
        Y[i]=y;
        P[i]=p;

        dt,caso = findmin([nextArr;nextCharge;nextDepON;nextDepOFF;nextSnapshot])

        ProgressMeter.update!(prog,arrivals);

    end

    T=T[1:i];
    X=X[1:i];
    Y=Y[1:i];
    P=P[1:i];

    pD=expired/arrivals;
    next!(prog);

    trace = TimeTrace(T,X,Y,P);
    stats = SimStatistics([],[],[],[],NaN,NaN,NaN,NaN);
    return EVSim(params,trace,finished,snaps,stats)
end

### DataFrame Compatibility

#Recibe un dataframe que tiene:
#arribos, demandas, salidas, potencias y opcionalmente salidaReportada
#se encarga de desarmar el dataframe y llamar a ev_sim_trace anterior
function ev_sim_trace(df::DataFrame,policy,C,snapshots)

    if "salidaReportada" in names(df)
        ev_sim_trace(df[!,:arribos],df[!,:demandas],df[!,:salidas],df[!,:potencias],policy,C,snapshots; salidaReportada=df[!,:salidaReportada])
    else
        ev_sim_trace(df[!,:arribos],df[!,:demandas],df[!,:salidas],df[!,:potencias],policy,C,snapshots)
    end
end


text:
"""
    EGPbetapower(σ, ξ, δ)

*EGPbetapower* corresponds to the fourth extended GP model of Naveau et al. (2016).

It is a four parameters family: δ is a threshold tuning parameter, κ is a parameter that controls the lower tail behavior, σ is a scale parameter, and ξ is a shape parameter.

```julia
EGPbetapower(σ, ξ, δ, κ)   # EGP of Naveau et al. (2016) (type 4) with parameters σ, ξ, δ, κ.

params(d)           # Get the parameters, i.e. (σ, ξ, δ, κ)
```

Reference :

* Naveau, P., Huser, R., Ribereau, P., and Hannart, A. (2016), Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection, Water Resour. Res., 52, 2753– 2769, doi:10.1002/2015WR018552.
"""
struct EGPbetapower{T<:Real} <: ContinuousUnivariateDistribution
    σ::T    # scale parameter
    ξ::T    # shape parameter
    δ::T    # threshold tuning parameter
    κ::T    # lower tail parameter

    function EGPbetapower{T}(σ::T, ξ::T, δ::T, κ::T) where {T <: Real}
        new{T}(σ, ξ, δ, κ)
    end
end

function EGPbetapower(σ::T, ξ::T, δ::T, κ::T; check_args=true) where {T <: Real}
    check_args && @check_args(EGPbetapower, σ > zero(σ) && δ > zero(δ) && κ > zero(κ))
    return EGPbetapower{T}(σ, ξ, δ, κ)
end

EGPbetapower(σ::Real, ξ::Real, δ::Real, κ::Real) = EGPbetapower(promote(σ, ξ, δ, κ)...)
EGPbetapower(σ::Integer, ξ::Integer, δ::Integer, κ::Integer) = EGPbetapower(float(σ), float(ξ), float(δ), float(κ))

minimum(d::EGPbetapower) = 0.0
maximum(d::EGPbetapower) = Inf * (d.ξ >= 0) - ( d.σ / d.ξ ) * (d.ξ < 0)
insupport(d::EGPbetapower, x::Real) = minimum(d) <= x <= maximum(d)

#### Parameters

scale(d::EGPbetapower) = d.σ
shape(d::EGPbetapower) = d.ξ    # noms à vérifier...
thresh(d::EGPbetapower) = d.δ
lowertail(d::EGPbetapower) = d.κ

params(d::EGPbetapower) = (d.σ, d.ξ, d.δ, d.κ)
partype(::EGPbetapower{T}) where {T} = T

#### Evaluation

function logpdf(d::EGPbetapower{T}, x::Real) where T<:Real
    μ = 0
    (σ, ξ, δ, κ) = params(d)

    pd = GeneralizedPareto(μ, 1, ξ)
    beta_pd = Beta(1/δ, 2)

    #g(v::Real, δ::Real, κ::Real) = (κ/2)*(1 - cdf(Beta(1/δ,2), (1-v)^δ))^(κ/2 - 1) * δ*(1-v)^(δ-1) * pdf(Beta(1/δ, 2), (1 - v)^δ)
    lg(v::Real, δ::Real, κ::Real) = log(κ/2) + (κ/2 - 1)*log(1 - Distributions.cdf(Beta(1/δ,2), (1-v)^δ)) + log(δ) + (δ-1)*log(1-v) + Distributions.logpdf(beta_pd, (1 - v)^δ)

    lf = -log(σ) + Distributions.logpdf.(pd, x/σ) + lg(Distributions.cdf(pd, x/σ), δ, κ)  # log de la vraisemblance

    return lf
end

pdf(d::EGPbetapower, x::Real) = exp(logpdf(d, x))

function logcdf(d::EGPbetapower{T}, x::Real) where T<:Real
    μ = 0
    (σ, ξ, δ, κ) = params(d)

    pd = GeneralizedPareto(μ, 1, ξ)

    G(v::Real, δ::Real, κ::Real) = (1 - Distributions.cdf(Beta(1/δ, 2), (1 - v)^δ))^(κ/2)

    lf = log(G(Distributions.cdf(pd, x/σ), δ, κ))

    return lf
end

cdf(d::EGPbetapower, x::Real) = exp(logcdf(d, x))

function quantile(d::EGPbetapower{T}, p::Real) where T<:Real
    @assert zero(p)<p<one(p) "the quantile level should be between 0 and 1."

    (σ, ξ, δ, κ) = params(d)

    invG(v::Real, δ::Real, κ::Real) = 1 - Distributions.quantile(Beta(1/δ, 2), (1 - v^(2/κ)))^(1/δ)

    x = (σ/ξ)*(((1 - invG(p,δ,κ))^(-ξ)) - 1)

    return x
end

#### Sampling

function rand(rng::AbstractRNG, d::EGPbetapower)
    # Generate a Float64 random number uniformly in (0,1].
    u = 1 - rand(rng)

    invG(v::Real, δ::Real, κ::Real) = 1 - Distributions.quantile(Beta(1/δ, 2), (1 - v^(2/κ)))^(1/δ)

    return (d.σ/d.ξ)*(((1 - invG(u,d.δ,d.κ))^(-d.ξ)) - 1)
end

sampler(d::EGPbetapower) = d


text:
using ODLCommonTools
using SummationByParts
using PumiInterface
using PdePumiInterface

dmg_name =  ".null"
smb_name = "tri8q.smb"
order = 1

sbp = TriSBP{Float64}(degree=order)

#mesh = PumiMesh2{Float64}(dmg_name, smb_name, order, sbp, dofpernode=4, shape_type=1)

arr = [3 4 5; 7 6 1]
min, max = PdePumiInterface.getMinandMax(arr)

println("min = ", min, ", max = ", max)


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ 3ef58460-c5a7-11eb-3bc6-f73e35e42d97
md"""
### Test Code for Module Tree
"""

# ╔═╡ 2e2c77ad-49e3-4e21-a119-1bee9118d848
module script include("./script.jl") end

# ╔═╡ 6e2eec65-d0e6-4e27-a2f4-11b4b820a805
Tree = script.Tree

# ╔═╡ 8279d02f-0684-4a9b-b1ab-705e74d7cd6e
module tools include("../tools.jl") end

# ╔═╡ b5ef1e5b-8078-490c-8a2e-79087653918f
JuTools = tools.JuTools

# ╔═╡ 4ed8ef3a-c7ee-48ce-baa8-84a4ecff939d
import Plots

# ╔═╡ 9db5eb6f-5822-4af5-95d8-997a9cdf3084
import Random

# ╔═╡ 360ed8fc-f1a3-4970-a3e8-a00e78734fc7
begin
X_data1, Y_data1 = JuTools.data_generate_cluster_2d(pos1=(20.0, 20.0), pos2=(50.0, 50.0),
    radius1=5.0, radius2=5.0, random_scale=8.0, data_size=1000)
size(X_data1), size(Y_data1)
end

# ╔═╡ 0eba36a8-316a-47a7-9d6c-f25eefa80192
begin
max_depth = 3
tree1 = Tree.create_decision_tree(X_data1, Y_data1, max_depth=max_depth)
end

# ╔═╡ 9284844b-a29c-49c2-bc78-2aa66bfa6460
let
plot_X1 = [X_data1[i, 1] for i in 1:(size(X_data1)[1]) if Y_data1[i] == 0.0]
plot_Y1 = [X_data1[i, 2] for i in 1:(size(X_data1)[1]) if Y_data1[i] == 0.0]
plot_X2 = [X_data1[i, 1] for i in 1:(size(X_data1)[1]) if Y_data1[i] == 1.0]
plot_Y2 = [X_data1[i, 2] for i in 1:(size(X_data1)[1]) if Y_data1[i] == 1.0]
X_test_1 = 0:70
X_test_2 = 0:70
X_test_contour = hcat(vec(repeat(reshape(X_test_1, :, 1), 1, 71)'), vec(repeat(X_test_2, 71, 1)))
X_test_contour = reshape(Tree.predict(X_test_contour, tree1), 71, :)
Plots.gr()
Plots.contour(X_test_1, X_test_2, X_test_contour, fill=true, c=Plots.cgrad(:matter, rev=false, categorical=false),
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X1, plot_Y1, leg=false, c="red")
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.xlims!((0.0, 70.0))
Plots.ylims!((0.0, 70.0))
Plots.title!("Decision Tree C4.5, max_depth=$max_depth")
end

# ╔═╡ 52e82e96-9a5e-4294-b516-e4bcce028ef2
JuTools.compute_accuracy(Tree.predict(X_data1, tree1), Y_data1)

# ╔═╡ ad8c0261-e9fd-414d-8ed5-357b61162216
tree2 = Tree.create_decision_tree(X_data1, Y_data1, max_depth=10)

# ╔═╡ cbbfed54-c98d-4f4b-a347-4e437625e90b
let
max_depth = 10
plot_X1 = [X_data1[i, 1] for i in 1:(size(X_data1)[1]) if Y_data1[i] == 0.0]
plot_Y1 = [X_data1[i, 2] for i in 1:(size(X_data1)[1]) if Y_data1[i] == 0.0]
plot_X2 = [X_data1[i, 1] for i in 1:(size(X_data1)[1]) if Y_data1[i] == 1.0]
plot_Y2 = [X_data1[i, 2] for i in 1:(size(X_data1)[1]) if Y_data1[i] == 1.0]
X_test_1 = 0:70
X_test_2 = 0:70
X_test_contour = hcat(vec(repeat(reshape(X_test_1, :, 1), 1, 71)'), vec(repeat(X_test_2, 71, 1)))
X_test_contour = reshape(Tree.predict(X_test_contour, tree2), 71, :)
Plots.gr()
Plots.contour(X_test_1, X_test_2, X_test_contour, fill=true, c=Plots.cgrad(:matter, rev=false, categorical=false),
    linewidth=0, background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X1, plot_Y1, leg=false, c="red")
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.xlims!((0.0, 70.0))
Plots.ylims!((0.0, 70.0))
Plots.title!("Decision Tree C4.5, max_depth=$max_depth")
end

# ╔═╡ 7027f706-328e-4181-8d66-6267e55f8982
JuTools.compute_accuracy(Tree.predict(X_data1, tree2), Y_data1)

# ╔═╡ Cell order:
# ╟─3ef58460-c5a7-11eb-3bc6-f73e35e42d97
# ╠═2e2c77ad-49e3-4e21-a119-1bee9118d848
# ╠═6e2eec65-d0e6-4e27-a2f4-11b4b820a805
# ╠═8279d02f-0684-4a9b-b1ab-705e74d7cd6e
# ╠═b5ef1e5b-8078-490c-8a2e-79087653918f
# ╠═4ed8ef3a-c7ee-48ce-baa8-84a4ecff939d
# ╠═9db5eb6f-5822-4af5-95d8-997a9cdf3084
# ╠═360ed8fc-f1a3-4970-a3e8-a00e78734fc7
# ╠═0eba36a8-316a-47a7-9d6c-f25eefa80192
# ╠═9284844b-a29c-49c2-bc78-2aa66bfa6460
# ╠═52e82e96-9a5e-4294-b516-e4bcce028ef2
# ╠═ad8c0261-e9fd-414d-8ed5-357b61162216
# ╠═cbbfed54-c98d-4f4b-a347-4e437625e90b
# ╠═7027f706-328e-4181-8d66-6267e55f8982

max_stars_repo_name: teamclouday/JuliaLearn

max_issues_repo_name: teamclouday/JuliaLearn

max_forks_repo_name: teamclouday/JuliaLearn


text:
export modelimage

"""
    $(TYPEDEF)
Container for non-analytic model that contains a image cache which will hold the image,
a Fourier transform cache C, which usually an instance of a <: FourierCache which usually
holds a interpolator that allows you to compute visibilities.

# Notes
This is an internal implementation detail that shouldn't usually be called directly.
Instead the user should use the exported function `modelimage`, for example

```julia
using Comrade
m = ExtendedRing(20.0, 5.0)

# This creates an version where the image is dynamically specified according to the
# radial extent of the image
mimg = modelimage(m) # you can also optionally pass the number of pixels nx and ny

# Or you can create an IntensityMap
img = intensitymap(m, 100.0, 100.0, 512, 512)
mimg = modelimage(m, img)
```
"""
struct ModelImage{M,I,C} <: AbstractModelImage{M}
    model::M
    image::I
    cache::C
end
@inline visanalytic(::Type{<:ModelImage{M}}) where {M} = visanalytic(M)
@inline imanalytic(::Type{<:ModelImage{M}}) where {M} = imanalytic(M)
@inline isprimitive(::Type{<:ModelImage{M}}) where {M} = isprimitive(M)

function Base.show(io::IO, mi::ModelImage)
   #io = IOContext(io, :compact=>true)
   #s = summary(mi)
   ci = first(split(summary(mi.cache), "{"))
   println(io, "ModelImage")
   println(io, "\tmodel: ", summary(mi.model))
   println(io, "\timage: ", summary(mi.image))
   println(io, "\tcache: ", ci)
end

model(m::AbstractModelImage) = m.model
flux(mimg::ModelImage) = flux(intensitymap!(mimg.image, mimg.model))

# function intensitymap(mimg::ModelImage)
#     intensitymap!(mimg.image, mimg.model)
#     mimg.image
# end

radialextent(m::ModelImage) = fov(m.image)[1]/2

#@inline visibility_point(m::AbstractModelImage, u, v) = visibility_point(model(m), u, v)

@inline intensity_point(m::AbstractModelImage, x, y) = intensity_point(model(m), x, y)

"""
    $(SIGNATURES)
Construct a `ModelImage` from a `model`, `image` and the optionally
specified visibility algorithm `alg`

# Notes
For analytic models this is a no-op and just return the model.
For non-analytic models this wraps the model in a object with an image
and precomputes the fourier transform using `alg`.
"""
@inline function modelimage(model::M, image::ComradeBase.AbstractIntensityMap; alg=FFT()) where {M}
    return modelimage(visanalytic(M), model, image; alg)
end

@inline function modelimage(::IsAnalytic, model, image; alg=FFT())
    return model
end

@inline function modelimage(::NotAnalytic, model, image; alg=FFT())
    intensitymap!(image, model)
    cache = create_cache(alg, image)
    return ModelImage(model, image, cache)
end

"""
    modelimage(m)
Construct a `ModelImage` where just the model `m` is specified

# Notes
If m `IsAnalytic()` is the visibility domain this is a no-op and just returns the model itself.
Otherwise `modelimage` will *guess* a reasonable field of view based on the `radialextent`
function. One can optionally pass the number of pixels nx and ny in each direction.
"""
function modelimage(m::M;
                    fovx=2*radialextent(m),
                    fovy=2*radialextent(m),
                    nx=512,
                    ny=512,
                    alg=FFT()) where {M}
    if visanalytic(M) == IsAnalytic()
        return m
    else
        T = typeof(intensity_point(m, 0.0, 0.0))
        img = IntensityMap(zeros(T,ny,nx), fovx, fovy)
        modelimage(m, img; alg)
        #throw(ArgumentError("$m is not an analytic model a image must be specified"))
    end
end


text:
function ACQR_em(y,Ai,Ci,Qi,Ri,m1i,P1i,max_iter::Int,tol::Float64,txo::Bool)
	#
	# estimate A, C, Q, R, m1, P1 using the EM algorithm for model
	#
	# x_{t+1} = A*x_{t} + w_{t}
	# y_{t}   = C*x_{t} + v_{t}
	#
	# cov(w_{t},v_{t}) = [Q 0;0 R]
	# x1 -> N(m1,P1)
	#
	# javier.cara@upm.es, 2106-02
	#

	# data by rows
	y,ny,nt = byrow(y)

	# number of states
	nx = size(Ai,1)

	# initial values
	A = Ai
	C = Ci
	Q = Qi
	R = Ri
	m1 = m1i
	P1 = P1i

	# log-likelihood values
	loglikv = zeros(max_iter)

	# Syy does not depend on the iterations
	Syy = zeros(ny,ny)
	for t in 1:nt
		Syy = Syy + y[:,t]*y[:,t]'
	end

	tol1 = 1.0
	iter = 1
	while (iter <= max_iter) && (tol1 > tol)
		time1 = time()

		# E-step
		# ---------------------------------------------------------------------------------
		# Kalmanfilter
		(xtt,Ptt,xtt1,Ptt1,et,St,Kt,loglik) = ACQR_kfilter(y,A,C,Q,R,m1,P1)
		(xtN,PtN,Pt1tN) = ACQR_ksmoother(A,xtt,Ptt,xtt1,Ptt1)

		loglikv[iter] = loglik
		if iter > 1
			tol1 = abs( (loglikv[iter] - loglikv[iter-1])/loglikv[iter-1] )
		end

		# initial values
		Sxx = zeros(nx,nx)
		Sx1x = zeros(nx,nx)
		Syx = zeros(ny,nx)

		# matrices Sxx, Sx1x, Syx, Sx1x1
		for t in 1:nt
			Sxx = Sxx + PtN[:,:,t] + xtN[:,t]*xtN[:,t]'
			Sx1x = Sx1x + Pt1tN[:,:,t] + xtN[:,t+1]*xtN[:,t]'
			Syx = Syx + y[:,t]*xtN[:,t]'
		end
		Sx1x1 = Sxx - (PtN[:,:,1] + xtN[:,1]*xtN[:,1]') + (PtN[:,:,nt+1] + xtN[:,nt+1]*xtN[:,nt+1]')

		# M-step
		# -------------------------------------------------------------------------------------
		# Matrices m1 y P1
		m1 = xtN[:,1]
		P1 = PtN[:,:,1]

		# Matrix A
		A = Sx1x/Sxx

		# Matrix Q
		M1 = Sx1x*A'
		Q = Sx1x1 - M1 - M1' + A*Sxx*A'
		Q = 1/nt*Q
		Q = (Q + Q')/2 # to make sure it's a symmetric matrix

		# Matrix C
		C = Syx/Sxx

		# Matrix R
		M1 = Syx*C'
		R = Syy - M1' - M1 + C*Sxx*C'
		R = 1/nt*R
		R = (R + R')/2 # to make sure it's a symmetric matrix

		etime = time() - time1
		if txo
			println( "Iter " * @sprintf("%3d",iter) * ",   @time = " * @sprintf("%.2E",etime) * ",   logLik = " * @sprintf("%.6E",loglik) * ",   tol = ", @sprintf("%.2E",tol1) )
		end

		iter += 1

	end

	loglikv = loglikv[1:(iter-1)]

	# Akaike Information Criterion
	P = ny*nx + nx*nx + nx*(nx+1)/2 + ny*(ny+1)/2
	aic = -2*loglikv[end] + 2*P

	# output
	return A,C,Q,R,m1,P1,loglikv,aic

end


text:
function positional_kmer_kernel( str1, str2 )
    s = 0
    prev = 0
    for i = 1:min( Base.length(str1), Base.length(str2) )
        if str1[i] == str2[i]
            prev += 1
            s += prev
        else
            prev = 0
        end
    end
    return s
end

function positional_kmer_kernel( n, str1, str2 )
    s = 0
    prev= 0
    for i = 1:min( Base.length(str1), Base.length(str2) )
        if str1[i] == str2[i]
            prev = min( n, prev + 1 )
            s += prev
        else
            prev = 0
        end
    end
    return s
end

function positional_kmer_kernel( n )
    kernel(a,b) = positional_kmer_kernel(n,a,b)
    return kernel
end

max_stars_repo_name: tbepler/ets-analysis-scripts

max_issues_repo_name: tbepler/ets-analysis-scripts

max_forks_repo_name: tbepler/ets-analysis-scripts


text:
ang = pi/3
c = DCM(ang,'x',false)
b = [6;2;6]
d = c*c*b
print(d)

Earth_a = 149600000
Mars_a = 227920000
Sun_mu = 1.3271*10^11
e = circHohmann(Earth_a, Mars_a, Sun_mu)


text:
function L1Frame(arg0::CelestialBody, arg1::CelestialBody)
    return L1Frame((CelestialBody, CelestialBody), arg0, arg1)
end


text:
include("handler/router.jl")
include("handler/base.jl")
include("handler/index_page.jl")
include("handler/state.jl")
include("handler/handlers.jl")


text:
abstract type AbstractSimulation end
"""
    Simulation

Contains all data needed by mesher and em-solver to run a simulation.

The simulation object shouldn't include any information about space or time discretization
or any information related to the solver or numerical algorithms used.
Its only purpose is to hold all the data needed by the solver and mesher to run a simulation.

The workflow is as follows:
- Construct all simulation componenets
- Call the mesh function to get a mesh object
- Discretize the simulation data
- Initalize fields for storing the solution
- Call the solve function to get the solution

"""
@with_kw struct Simulation <: AbstractSimulation
    background
    structures
    sources
    boundaryconditions
    # symmetries = NoSymmetry()
    # monitors = []
end


text:
#---------------------------------------------------------
# # [ImagePhantoms overview](@id 01-overview)
#---------------------------------------------------------

#=
This page explains the Julia package
[`ImagePhantoms`](https://github.com/JuliaImageRecon/ImagePhantoms.jl).

This page was generated from a single Julia file:
[01-overview.jl](@__REPO_ROOT_URL__/01-overview.jl).
=#

#md # In any such Julia documentation,
#md # you can access the source code
#md # using the "Edit on GitHub" link in the top right.

#md # The corresponding notebook can be viewed in
#md # [nbviewer](http://nbviewer.jupyter.org/) here:
#md # [`01-overview.ipynb`](@__NBVIEWER_ROOT_URL__/01-overview.ipynb),
#md # and opened in [binder](https://mybinder.org/) here:
#md # [`01-overview.ipynb`](@__BINDER_ROOT_URL__/01-overview.ipynb).

# ### Setup

# Packages needed here.

using ImagePhantoms
using ImageGeoms: ImageGeom, axesf
using MIRTjim: jim, prompt
using UnitfulRecipes
using Unitful: mm
using InteractiveUtils: versioninfo

# The following line is helpful when running this example.jl file as a script;
# this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

# ### Overview

#=
When developing image reconstruction methods,
it can be helpful to simulate data (e.g., sinograms)
using software-defined images called phantoms.

The simplest method here is to make a Shepp-Logan phantom image
similar its use in other packages.
=#

image = shepp_logan(256) # CT version by default
jim(image, "SheppLogan", clim=(0.9, 1.1), yflip=false)

# ### Sinograms and spectra

#=
Often for image reconstruction algorithm development,
we need not only the phantom image, but also its sinogram and spectrum.
We start with the vector of ellipses that defines the phantom,
using a typical field of view (FOV) of 200mm for a head:
=#

objects = shepp_logan(SheppLoganToft(); fovs=(200mm,200mm))

#=
From that collection we can compute images, sinograms and spectra.
It is convenient (but not required)
to use `ImageGeoms` to help with the sampling.
=#

ig = ImageGeom(dims=(200,256), deltas=(1mm,1mm))
image = phantom(axes(ig)..., objects)
jim(axes(ig)..., image, xlabel="x", ylabel="y", title="SheppLoganToft")

# Here is the sinogram corresponding to this phantom,
# computed analytically from the ellipse parameters:

r = LinRange(-100mm,100mm,401)
ϕ = deg2rad.(0:180)
sino = radon(r, ϕ, objects)
jim(r, ϕ, sino, title="Sinogram", yflip=false, aspect_ratio=:none)

# Here is the 2D spectrum (Fourier transform) of this phantom,
# computed analytically from the ellipse parameters:

kspace = spectrum(axesf(ig)..., objects)
jim(axesf(ig)..., log10.(abs.(kspace/(1mm)^2)), xlabel="ν₁", ylabel="ν₂", title="log10|Spectrum|")

# ## Reproducibility

# This page was generated with the following version of Julia:

io = IOBuffer()
versioninfo(io)
split(String(take!(io)), '\n')

# And with the following package versions

import Pkg; Pkg.status()


text:
load("utils/req.jl")
req("utils/utils.jl")

# == split_fdef ===============================================================

function split_fdef(fdef::Expr)
    @expect (fdef.head == :function) || (fdef.head == :(=))
    @expect length(fdef.args) == 2
    signature, body = tuple(fdef.args...)
    @expect is_expr(signature, :call)
    @expect length(signature.args) >= 1
    (signature, body)
end
split_fdef(f::Any) = error("split_fdef: expected function definition, got\n$f")

# == @failexpect ==============================================================

const failexp_handlers = Dict{Function,Function}()

# add failed expect message generator, e g
# @failexpect is_expr(ex, head) = "expected expr(:$head,...), got $ex"
macro failexpect(fdef)
    signature, body = split_fdef(fdef)
    f, args = signature.args[1], signature.args[2:end]

    quote; let
        const handler = get_failexp_handler(($f), ($string(f)))
        handler($args...) = $body
    end; end
end
function get_failexp_handler(f::Function, fname::String)
    if !has(failexp_handlers, f)
        hname = gensym("failexpect_handler_for__"*fname)
        failexp_handlers[f] = @eval begin
            ($hname)(args...) = nothing
            ($hname)
        end
    end
    failexp_handlers[f]
end

function fail_expect(predexpr, pred::Function, args::Tuple)
    if has(failexp_handlers, pred); msg = failexp_handlers[pred](args...)
    else                            msg = nothing; end

    is(msg, nothing) ? fail_expect(predexpr) : error(msg)
end

function code_checkexpect_call(predexpr::Expr)
    @assert predexpr.head == :call
    pred, args = predexpr.args[1], predexpr.args[2:end]
    qpredexpr =  quoted_expr(predexpr)

    quote; let args = ($quoted_tuple(args))
        ($pred)(args...) ? nothing : fail_expect(($qpredexpr), ($pred), args)
    end; end
end
# tie into @expect
function code_checkexpect(ex::Expr)
    if is_expr(ex, :call); code_checkexpect_call(ex)
    else;                  default_code_checkexpect(ex);   end
end

@failexpect is_expr(ex, head::Symbol) = "expected expr(:$head,...), got $ex"

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
path = "docs/src/assets/"

using ISA, Plots
t = 0.0:6.25e-5:3.0
𝐶₀ = AMFMtriplet(t->0.2+0.8cos(11t), t->200.0, 0.0)
𝐶₁ = AMFMtriplet(t->exp(-abs(t/3)), t->100t, 0.1)
𝐶₂ = AMFMtriplet(t->𝒩ᵤ(t; μ=1.5, σ=1.0), t->150 + 125sin(5t), π)
𝐶₃ = AMFMtriplet(t->u(t-1.725)-u(t-2.475), t->50, 0.0)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂,𝐶₃])

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

function getFrames(cameraStart,cameraEnd,nSteps)
    a = cameraStart[1]:(cameraEnd[1].-cameraStart[1])./(nSteps-1):cameraEnd[1]
    b = cameraStart[2]:(cameraEnd[2].-cameraStart[2])./(nSteps-1):cameraEnd[2]
    viewAngles = [(a[i],b[i]) for i∈1:nSteps]
    return [plot(𝑆, timeaxis=t, camera=cam) for cam∈viewAngles]
end

function holdFrame(camera,nTimes)
    return [plot(𝑆, timeaxis=t, camera=camera) for i∈1:nTimes]
end

view3D = (20,80)
viewTF = (0,90)
viewRT = (0,0)
viewRF = (90,0)

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

nTimes = 30
camera = view3D
frames = holdFrame(camera,nTimes)

nSteps = 20
cameraStart = view3D
cameraEnd = viewTF
frames = cat(frames,getFrames(cameraStart,cameraEnd,nSteps),dims=1)

nTimes = 30
camera = viewTF
frames = cat(frames,holdFrame(camera,nTimes),dims=1)

nSteps = 20
cameraStart = viewTF
cameraEnd = view3D
frames = cat(frames,getFrames(cameraStart,cameraEnd,nSteps),dims=1)

nTimes = 30
camera = view3D
frames = cat(frames,holdFrame(camera,nTimes),dims=1)

nSteps = 80
cameraStart = view3D
cameraEnd = viewRT
frames = cat(frames,getFrames(cameraStart,cameraEnd,nSteps),dims=1)

nTimes = 30
camera = viewRT
frames = cat(frames,holdFrame(camera,nTimes),dims=1)

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

anim = Animation()
for i∈1:length(frames)
    plot(frames[i])
    frame(anim) #frame pngs are saved in temp directory anim.dir
end
mp4(anim, path*"movieTest.mp4", fps=30)


text:
using PhomericsTest
using Test

@testset "Operators.jl" begin
    @test Matrix(DX(1.0, 2, 2)) == [0.0 0.5 0.0 0.0; -0.5 0.0 0.0 0.0; 0.0 0.0 0.0 0.5; 0.0 0.0 -0.5 0.0]
    @test Matrix(DXX(1.0, 2, 2)) == [-2.0 1.0 0.0 0.0; 1.0 -2.0 0.0 0.0; 0.0 0.0 -2.0 1.0; 0.0 0.0 1.0 -2.0]
    @test Matrix(DY(1.0, 2, 2)) == [0.0 0.0 0.5 0.0; 0.0 0.0 0.0 0.5; -0.5 0.0 0.0 0.0; 0.0 -0.5 0.0 0.0]
    @test Matrix(DYY(1.0, 2, 2)) == [-2.0 0.0 1.0 0.0; 0.0 -2.0 0.0 1.0; 1.0 0.0 -2.0 0.0; 0.0 1.0 0.0 -2.0]
end


text:
module UCIMLRepo

using HTTP
using DataArrays
using DataFrames
using PyCall

@pyimport lxml.html as lh

export ucirepodata,ucirepoinfo,ucirepolist

function download_http(url::String)
    println("fetching from the following url : ",url)
    request = HTTP.request("GET", url)
    if request.status != 200
        error("download failed")
    end
#converting received data into ASCII format and returning the value
    rcvd_data = String(request.body)
end

#This function converts the input string into a dataframe
function parse_rcvd_data(rcvd_data::String)
    rows = String[]
    rows = split(rcvd_data,'\n',keep = false)
    df = DataFrame()

    for i = 1:size(rows)[1]
        df[i] = split(rows[i],',',keep =true)
    end
    df
end

#TODO: Use readtable("file.csv") instead of splitting the array.

function parse_rcvd_info(rcvd_info::String)
    rows = String[]
    rows = split(rcvd_info,'\n',keep =false)

    for i = 1:size(rows)[1]
        println(rows[i])
    end
end
#TODO: used the parsed info to create the header line for the dataframe

# Exported Function. Loads data into a dataframe from the repository as specified by the parameter
# usage: df = ucirepodata("iris")
# you can alternatively specify the exact url to *any* .data file that you want to load as so:
# df = ucirepodata("http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data",false)

function ucirepodata(DsetNameOrUrl::String, DorU::Bool=true)
    if(DorU)
#creating the url to be used for download
        url = string("http://archive.ics.uci.edu/ml/machine-learning-databases/", DsetNameOrUrl,"/",DsetNameOrUrl,".data")

    else
        url = DsetNameOrUrl
    end

#loading received data into a string type variable
    rcvd_data = download_http(url)

    df = DataFrame()
#parsing received data and loading into an appropriate dataframe
    df = parse_rcvd_data(rcvd_data)
end

#Exported Function. Displays all relevant information about the dataset
#Usage: ucirepoinfo("iris")
# you can alternatively specify the exact url to *any* .data file that you want to load and view, as so:
# ucirepoinfo("http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.names",false)

function ucirepoinfo(DsetNameOrUrl::String, DorU::Bool=true)
    if(DorU)
	#creating the url to be used for download
        url = string("http://archive.ics.uci.edu/ml/machine-learning-databases/", DsetNameOrUrl,"/",DsetNameOrUrl,".names")

    else
        url = DsetNameOrUrl
    end

#loading received data into a string type variable
    rcvd_info = download_http(url)

#parsing recieved info and displaying on STDOUT
    parse_rcvd_info(rcvd_info)
end

#displays a list of all the available datasets.
#work in progress. Function not exported yet.
#need to work on parsing the html data rcvd to extract a list of available datasets

function ucirepolist()
    doc = lh.parse("http://archive.ics.uci.edu/ml/datasets.html")
    text = doc[:xpath]("/html//table[2]//table[2]")
    rows = text[2][:xpath]("/html//table[2]//table[2]//tr[1]")
    i =2;

    print("name|default task\n")
    while i<= 284

         name 			= rows[1][:xpath](string("""/html//table[2]//td[2]/table[2]//tr[""",i,"""]//td[2]/p//text()"""))[1]
         dtype 			= rows[1][:xpath](string("""/html//table[2]//td[2]/table[2]//tr[""",i,"""]//td[2]/p//text()"""))[2]
         #dtask[i-1] 		= rows[1][:xpath](string("""/html//table[2]//td[2]/table[2]//tr[""",i,"""]//td[3]/p//text()"""))
         #atypes[i-1] 		= rows[1][:xpath](string("""/html//table[2]//td[2]/table[2]//tr[""",i,"""]//td[4]/p//text()"""))
         #instances[i-1] 	= rows[1][:xpath](string("""/html//table[2]//td[2]/table[2]//tr[""",i,"""]//td[5]/p//text()"""))
         #attributes[i-1] 	= rows[1][:xpath](string("""/html//table[2]//td[2]/table[2]//tr[""",i,"""]//td[6]/p//text()"""))
         #year[i-1] 		= rows[1][:xpath](string("""/html//table[2]//td[2]/table[2]//tr[""",i,"""]//td[7]/p//text()"""))

	 print(name," | ", dtype,"\n")
	 i += 1
    end

end

end # module


text:
import TakingBroadcastSeriously: @unfuse, broadcast_
using Base.Test

# A fake custom array type

struct FooArray{T,N} <: AbstractArray{T,N}
  data::Array{T,N}
end

Base.size(xs::FooArray, a...) = size(xs.data, a...)
Base.getindex(xs::FooArray, a...) = getindex(xs.data, a...)

# Unfuse stuff

@unfuse FooArray

blist = []

function broadcast_(::typeof(sin), xs::FooArray)
  push!(blist, sin)
  FooArray(sin.(xs.data))
end

function broadcast_(::typeof(cos), xs::FooArray)
  push!(blist, cos)
  FooArray(cos.(xs.data))
end

@testset "Seriously" begin
    xs = rand(5)
    xs′ = FooArray(xs)

    @test cos.(sin.(xs)) == cos.(sin.(xs′))

    @test blist == [sin, cos]

    @test xs.^3.0 == xs′.^3.0
    @test 50 .* xs == 50 .* xs′
    @test xs.^3 == xs′.^3
end


text:
const CONTEXT_CACHE = Dict{Module, Context}()
const FLAGS = CXTranslationUnit_DetailedPreprocessingRecord | CXTranslationUnit_SkipFunctionBodies | CXTranslationUnit_PrecompiledPreamble #| CXTranslationUnit_SingleFileParse
const HEADER = joinpath(@__DIR__, "$(nameof(@__MODULE__)).h")

Base.String(ctx::Context) = String(typeof(ctx))
header(ctx::Context) = header(typeof(ctx))
header(::Type{C}) where {C<:Context}  = HEADER
language(ctx::Context) = language(typeof(ctx))
language(::Type{Context{lang}}) where {lang} = lang

getcontext(mod::Module) = get(CONTEXT_CACHE, mod, nothing)
Base.close(ctx::Context) = delete!(CONTEXT_CACHE, ctx.mod)

getjl(ctx::Context, args...; kwargs...) = getjl(typeof(ctx), args...; annotate = getblock(ctx).flags.annot, kwargs...)
getjl(::Type{C}, str::String; kwargs...) where {C<:Context} = getjl(C, Symbol(str); kwargs...)
getjl(::Type{C}, sym::Symbol; annotate::Bool = false) where {C<:Context} = esc(Symbol(replace(String(sym), r"^(struct|union|enum)\s+" => (annotate ? s"\1_" : ""))))

getname(ctx::Context, args...; kwargs...) = getname(typeof(ctx), args...; kwargs...)
getname(::Type{C}, str::String; kwargs...) where {C<:Context} = getname(C, Symbol(str); kwargs...)
function getname(::Type{C}, sym::Symbol; isbitstype::Bool = false) where {C<:Context}
	sym = Symbol("$(language(C))\"$(sym)\"")
	return isbitstype ? Symbol("($(sym))") : sym
end

gettype(ctx::Context, args...; kwargs...) = gettype(typeof(ctx), args...; kwargs...)
gettype(::Type{C}, str::String; kwargs...) where {C<:Context} = gettype(C, Symbol(str); kwargs...)
gettype(::Type{C}, sym::Symbol; kwargs...) where {C<:Context} = esc(getname(C, sym; kwargs...))

# TODO: prevent new compiler contexts within a module from overwriting libraries of prior contexts
function getlib(ctx::Context)
	dir = get_scratch!(@__MODULE__, string(Scratch.find_uuid(ctx.mod)), ctx.mod)
	file = "libcbinding-$(length(ctx.blocks)).$(dlext)"
	return QuoteNode(Symbol(joinpath(dir, file)))
end

function getlib(ctx::Context, sym)
	for (handle, lib) in ctx.libs
		hdl = dlsym(handle, sym, throw_error = false)
		isnothing(hdl) || hdl == C_NULL || return lib isa Symbol ? QuoteNode(lib) : lib
	end
	getblock(ctx).flags.quiet || @warn "Failed to find `$(sym)` in:\n"*join(map(((handle, lib),) -> "  $(lib isa Symbol ? lib : "or the Julia process")", collect(ctx.libs)), '\n')
	return Nothing
end

function getcode(ctx::Context, cursor::CXCursor)
	# WARNING:  clang_getCursorPrettyPrinted has bugs:
	# struct S22 {
	# 	union { int i; float f; struct { int j; }; } (*(*f)(struct S22 x, int i))(int, float);
	# };
	policy = clang_getCursorPrintingPolicy(cursor)
	if policy != C_NULL
		code = _string(clang_getCursorPrettyPrinted, cursor, policy)
		isempty(code) || return code
	end

	# if easy way didn't produce anything, we need to do it the messy way
	ismacro = cursor.kind == CXCursor_MacroDefinition
	tokens = tokenize(ctx.tu[], cursor)
	tokens = ismacro ? tokens[2:end] : tokens
	tokens = join(map(token -> clang_getTokenKind(token) == CXToken_Comment ? "" : string(ctx.tu[], token), tokens), ' ')
	tokens = replace(tokens, "{" => "{\n")
	tokens = replace(tokens, " ;" => ";\n")
	tokens = replace(tokens, " }" => "\n}")
	tokens = replace(tokens, " ," => ",")
	tokens = replace(tokens, ": " => ":")
	tokens = replace(tokens, " :" => ":")
	tokens = replace(tokens, "( " => "(")
	tokens = replace(tokens, " )" => ")")
	tokens = replace(tokens, "[ " => "[")
	tokens = replace(tokens, " [" => "[")
	tokens = replace(tokens, " ]" => "]")
	tokens = replace(tokens, " (" => "(")
	tokens = replace(tokens, r"(\w)(\(\*)\s+" => s"\1 \2")
	tokens = replace(tokens, "  " => " ")
	tokens = replace(tokens, r"\n\s+" => "\n")

	indents = 0
	code = ""
	isParens = false
	wasComma = false
	wasEOL = false
	for c in tokens
		if c == '}'
			indents -= 1
		elseif c == '{'
			indents += 1
		elseif c == '('
			isParens = true
		elseif c == ')'
			isParens = false
		end

		if wasEOL
			code *= repeat("  ", indents)
		end
		if !(!isParens && wasComma)
			code *= c
		end

		wasEOL = c == '\n'
		wasComma = c == ','
		if !isParens && wasComma
			code *= '\n'
			wasEOL = true
		end
	end
	code = strip(code)

	if ismacro
		code = "#define $(string(cursor)) $(code)"
	end

	return code
end

getexprs(ctx::Context) = getexprs(ctx, clang_getTranslationUnitCursor(ctx.tu[]))

function getexprs(ctx::Context, syms, blocks...)
	expr = filter(!isnothing, collect(blocks))

	blk = getblock(ctx)
	if blk.flags.jlsyms
		for (sym, jlsym, docs) in syms
			isnothing(jlsym) || push!(expr, :(const $(jlsym) = $(sym)))
		end
	end

	if !blk.flags.priv
		push!(expr, :(export $(map(((sym, jlsym, docs),) -> sym, syms)...)))
		blk.flags.jlsyms && push!(expr, :(export $(filter(!isnothing, map(((sym, jlsym, docs),) -> jlsym, syms))...)))
	end

	if !blk.flags.nodocs
		for (sym, jlsym, docs) in syms
			isnothing(docs) && continue
			push!(expr, :(@doc $(docs) $(startswith(String(sym.args[1]), '@') ? QuoteNode(sym.args[1]) : sym.args[1])))
			isnothing(jlsym) && continue
			blk.flags.jlsyms && push!(expr, :(@doc $(docs) $(startswith(String(jlsym.args[1]), '@') ? QuoteNode(jlsym.args[1]) : jlsym.args[1])))
		end
	end

	push!(expr, first(first(syms)))
	return expr
end

function getexprs_tu(ctx::Context, cursor::CXCursor)
	exprs = []

	for child in children(cursor)
		if !getblock(ctx).flags.libc
			range = getlocation(child)
			isnothing(range) && continue  # ignore built-ins and such

			first(range).file == header(ctx) || haskey(ctx.hdrs, first(range).file) || continue
			first(range).file != header(ctx) || first(range).line > ctx.line || continue
		end

		append!(exprs, filter(!isnothing, getexprs(ctx, child)))
	end

	return exprs
end

function getexprs_include(ctx::Context, cursor::CXCursor)
	exprs = []

	file = clang_getIncludedFile(cursor)
	file = _string(clang_getFileName, file)

	if !isnothing(file) && !haskey(ctx.hdrs, file)
		crng = getlocation(cursor)
		basedir = get(ctx.hdrs, first(crng).file, "")
		isimplicit = !isempty(basedir) && startswith(file, basedir)
		isexplicit = first(crng).file == header(ctx)

		if isexplicit || isimplicit || getblock(ctx).flags.libc
			push!(ctx.hdrs, file => isexplicit && getblock(ctx).flags.implic ? "$(dirname(file))/" : basedir)
			push!(exprs, :(include_dependency($(file))))
		end
	end

	return exprs
end

getblock(ctx::Context) = isempty(ctx.blocks) ? nothing : last(ctx.blocks)
getblock(ctx::Context, ind::Integer) = ind in eachindex(ctx.blocks) ? ctx.blocks[ind] : nothing
function getblock(ctx::Context, loc::CodeLocation)
	loc.file == header(ctx) || return nothing

	ind = searchsortedfirst(ctx.blocks, loc.line, lt = (a::CodeBlock, b::Integer) -> a.lines.stop < b)
	return getblock(ctx, ind)
end

function getlines(ctx::Context, blk::CodeBlock)
	lines = []

	io = IOBuffer(ctx.src.data[1:ctx.src.size])
	for (ind, line) in enumerate(eachline(io))
		ind in blk.lines && push!(lines, line)
	end

	return lines
end

function Base.pop!(ctx::Context)
	blk = getblock(ctx)
	if !isnothing(blk)
		pop!(ctx.blocks)
		io = IOBuffer(take!(ctx.src))
		ctx.src = IOBuffer()
		for (ind, line) in enumerate(eachline(io))
			ind in blk.lines || println(ctx.src, line)
		end
	end

	return ctx
end

# WARNING:  risky hack to trick Julia dlopen into ccalling the C dlopen with NULL (meaning to dlopen current process rather than a library)
struct _NullCString <: AbstractString end
Base.cconvert(::Type{Cstring}, ::_NullCString) = Cstring(C_NULL)

function configure!(ctx::Context)
	# find packaged include
	includedir = joinpath(dirname(dirname(libclang_path())), "lib", "clang")
	for entry in readdir(includedir)
		m = match(r"^(\d+\.\d+\.\d+)$", entry)
		if !isnothing(m)
			includedir = joinpath(includedir, m.captures[1], "include")
			break
		end
	end

	pushfirst!(ctx.args,
		"-isystem", includedir,  # these -isystem args are needed since Clang_jll packaging seems to have clang's InstallDir not set/found correctly
		"-isystem", joinpath(dirname(dirname(libclang_path())), "include"),
		# NOTE: in wrap!, ctx.args[5:7] is used separately from ctx.args[8:end] in the compile command
		"-x", lowercase(String(ctx)),
		"-Wno-empty-translation-unit",
	)

	libs = []
	libpaths = []
	for ind in eachindex(ctx.args)
		if ctx.args[ind] == "-L"
			if ind < length(ctx.args)
				push!(libpaths, ctx.args[ind+1])
				ctx.args[ind] = ""
				ctx.args[ind+1] = ""
			end
		elseif ctx.args[ind] == "-l"
			if ind < length(ctx.args)
				push!(libs, ctx.args[ind+1])
				ctx.args[ind] = ""
				ctx.args[ind+1] = ""
		end
		elseif startswith(ctx.args[ind], "-L")
			push!(libpaths, ctx.args[ind][3:end])
			ctx.args[ind] = ""
		elseif startswith(ctx.args[ind], "-l")
			push!(libs, ctx.args[ind][3:end])
			ctx.args[ind] = ""
		end
	end

	for lib in libs
		l = find_library(lib, libpaths)
		l = isempty(l) ? find_library("lib"*lib, libpaths) : l
		isempty(l) && error("Failed to find library $(lib) in $(join(libpaths, ", "))")

		lib = dlopen(l)
		lib == C_NULL && error("Failed to dlopen library $(l)")
		push!(ctx.libs, lib => Symbol(l))
	end
	lib = dlopen(_NullCString())
	lib == C_NULL && error("Failed to dlopen Julia process")
	push!(ctx.libs, lib => Nothing)

	ctx.ind = clang_createIndex(false, true)
	ctx.ind == C_NULL && error("Failed to create index")
end

function Base.push!(ctx::Context, loc::LineNumberNode, flags::NamedTuple, str::String)
	str = endswith(str, '\n') ? str : str*'\n'  # ensure new line here so countlines is consistent
	print(ctx.src, str)

	blk = getblock(ctx)
	start = isnothing(blk) ? 1 : blk.lines.stop+1
	stop  = start + countlines(IOBuffer(str))-1
	push!(ctx.blocks, CodeBlock(
		loc,
		start:stop,
		flags,
		String[],
	))
end

function parse!(ctx::Context)
	empty!(ctx.hdrs)
	empty!(ctx.macros)

	unsaved = [CXUnsavedFile(
		Filename = pointer(header(ctx)),
		Contents = pointer(ctx.src.data),
		Length   = ctx.src.size,
	)]

	if ctx.tu[] == C_NULL
		ptrs = map(pointer, ctx.args)
		err = clang_parseTranslationUnit2(ctx.ind, header(ctx), ptrs, length(ptrs), unsaved, length(unsaved), FLAGS, ctx.tu)
	else
		err = clang_reparseTranslationUnit(ctx.tu[], length(unsaved), unsaved, FLAGS)
	end
	err == CXError_Success || error("Failed to parse translation unit $(err)")

	isfatal = false
	for ind in 1:clang_getNumDiagnostics(ctx.tu[])
		diag = nothing
		try
			diag = clang_getDiagnostic(ctx.tu[], ind-1)

			sev = clang_getDiagnosticSeverity(diag)
			loc = clang_getDiagnosticLocation(diag)
			loc = getlocation(loc)
			isnothing(loc) && continue
			(loc.file == header(ctx) && loc.line > ctx.line) || sev in (CXDiagnostic_Error, CXDiagnostic_Fatal) || continue

			blk = getblock(ctx, loc)
			isquiet = !isnothing(blk) && blk.flags.quiet
			if !isquiet
				msg = _string(clang_formatDiagnostic, diag, clang_defaultDiagnosticDisplayOptions())
				msg = replace(msg, (header(ctx)*r"(?:\:\d+)*") => (!isnothing(blk) ? "$(blk.loc.file):$(blk.loc.line)" : @__FILE__))

				if !isnothing(blk)
					lines = getlines(ctx, blk)
					line = (loc.line - blk.lines.start)+1

					pre = length(lines[line]) - length(lstrip(lines[line]))
					here = lines[line][1:pre]*repeat(" ", loc.column-pre-1)*"^~~~~~ here"
					insert!(lines, line+1, here)

					msg *= "\n"*join(lines[max(1, line-2):min(length(lines), line+3)], '\n')
				end
			end

			if sev in (CXDiagnostic_Error, CXDiagnostic_Fatal)
				isfatal = true
				isquiet || @error msg  (_module = !isnothing(blk) ? ctx.mod : @__MODULE__)  (_file = !isnothing(blk) ? String(blk.loc.file) : @__FILE__)  (_line = !isnothing(blk) ? blk.loc.line : @__LINE__)
			elseif sev == CXDiagnostic_Warning
				isquiet || @warn  msg  (_module = !isnothing(blk) ? ctx.mod : @__MODULE__)  (_file = !isnothing(blk) ? String(blk.loc.file) : @__FILE__)  (_line = !isnothing(blk) ? blk.loc.line : @__LINE__)
			else
				isquiet || @debug msg  (_module = !isnothing(blk) ? ctx.mod : @__MODULE__)  (_file = !isnothing(blk) ? String(blk.loc.file) : @__FILE__)  (_line = !isnothing(blk) ? blk.loc.line : @__LINE__)
			end
		finally
			clang_disposeDiagnostic(diag)
		end
	end

	isfatal && error("Errors parsing $(String(ctx)) code")
end

function wrap!(ctx::Context)
	getblock(ctx).flags.wrap || return
	isempty(getblock(ctx).inlines) && return

	mktemp() do path, file
		blk = 1
		io = IOBuffer(ctx.src.data[1:ctx.src.size])
		for (ind, line) in enumerate(eachline(io))
			println(file, line)
			if ind == ctx.blocks[blk].lines.stop
				foreach(inline ->  println(file, inline), ctx.blocks[blk].inlines)
				blk += 1
			end
		end
		close(file)

		libclang.Clang_jll.clang() do bin
			run(`$(bin) $(ctx.args[5:7]) -w -O2 -fPIC -shared -o $(getlib(ctx).value) $(path) $(ctx.args[8:end])`)  # TODO: add -rpath for all ctx.libs?
		end
	end
end

function advance!(ctx::Context)
	getblock(ctx).flags.defer && return

	stop = ctx.line
	root = clang_getTranslationUnitCursor(ctx.tu[])
	for cursor in children(root)
		range = getlocation(cursor)
		isnothing(range) && continue

		first(range).file == header(ctx) || continue
		if first(range).line > ctx.line && last(range).line > stop
			stop = last(range).line
		end
	end
	ctx.line = stop
end

function clang_cmd(mod::Module, loc::LineNumberNode, lang::Symbol, str::String)
	cmd = :(``)
	cmd.args[2] = loc
	cmd.args[3] = str

	return quote
		CONTEXT_CACHE[$(mod)] = Context{$(QuoteNode(lang))}($(mod), $(esc(cmd)).exec...)
		nothing
	end
end

function clang_str(mod::Module, loc::LineNumberNode, lang::Symbol, str::String, opts::String = "")
	# TODO: allow interpolation in str?

	flags = (;
		libc    = 'c' in opts,
		defer   = 'd' in opts,
		nofunc  = 'f' in opts,
		implic  = 'i' in opts,
		jlsyms  = 'J' in opts || 'j' in opts,
		annot   = 'J' in opts,
		nomacro = 'm' in opts,
		notify  = 'n' in opts,
		priv    = 'p' in opts,
		quiet   = 'q' in opts,
		ref     = 'r' in opts,
		skip    = 's' in opts,
		notype  = 't' in opts,
		nodocs  = 'u' in opts,
		novar   = 'v' in opts,
		wrap    = 'w' in opts,
	)

	ref = getref(Context{lang}, str)
	isnothing(ref) || return ref
	flags.ref && error("Failed to parse $(String(Context{lang})) expression $(language(Context{lang}))$(repr(str))")

	haskey(CONTEXT_CACHE, mod) || error("Compiler context ($(language(Context{lang}))`...`) not created before declaring $(String(Context{lang})) bindings ($(language(Context{lang}))\"...\") in $(mod)")
	ctx = CONTEXT_CACHE[mod]

	exprs = []
	try
		push!(ctx, loc, flags, str)
		parse!(ctx)
		exprs = getexprs(ctx)
		wrap!(ctx)
		advance!(ctx)
	catch
		pop!(ctx)
		rethrow()
	end

	return quote
		$(exprs...)
	end
end


text:
baremodule Try

export @?, Ok, Err

module InternalPrelude
abstract type AbstractResult{T,E} end
function _IsOkError end
end  # module InternalPrelude

struct Ok{T} <: InternalPrelude.AbstractResult{T,Union{}}
    value::T
end

struct Err{E} <: InternalPrelude.AbstractResult{Union{},E}
    value::E
    backtrace::Union{Nothing,typeof(InternalPrelude.backtrace())}
end

function unwrap end
function unwrap_err end

function oktype end
function errtype end
function isok end
function iserr end

function enable_errortrace end
function disable_errortrace end

# Core exceptions
struct IsOkError <: InternalPrelude.Exception
    ok::InternalPrelude.AbstractResult

    InternalPrelude._IsOkError(ok) = new(ok)
end

function map end

macro and_return end
function var"@?" end
function var"@return" end
macro and end
macro or end

function and end
function and_then end
function or end
function or_else end
function unwrap_or_else end

module Internal

import ..Try: @and_return, @?, @and, @or
using ..Try: Err, Ok, Try
using ..Try.InternalPrelude: AbstractResult, _IsOkError

include("ExternalDocstrings.jl")
using .ExternalDocstrings: @define_docstrings

include("core.jl")
include("errortrace.jl")
include("tools.jl")
include("branch.jl")
include("show.jl")

end  # module Internal

Internal.@define_docstrings

end  # baremodule Try


text:
using MyWorkflow
using Documenter

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

makedocs(;
    modules=[MyWorkflow],
    authors="sakurairihito <sakurairihito@gmail.com> and contributors",
    repo="https://github.com/Sakurairihito/MyWorkflow.jl/blob/{commit}{path}#{line}",
    sitename="MyWorkflow.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Sakurairihito.github.io/MyWorkflow.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Sakurairihito/MyWorkflow.jl",
    devbranch="main",
)

max_stars_repo_name: sakurairihito/playground

max_issues_repo_name: sakurairihito/playground

max_forks_repo_name: sakurairihito/playground


text:
import Base.Iterators: flatten

function puzzle07(path::String=joinpath(@__DIR__, "..", "data/07.txt"))
    input = readlines(path)
    programnames = []
    weights = []
    for line in input
        append!(programnames, [[m.match for m in eachmatch(r"[a-zA-Z]+", line)]])
        for m in eachmatch(r"[0-9]+", line)
            w = parse(Int, m.match)
            append!(weights, [w])
        end
    end
    weights = [first.(programnames) weights]

    graph = Array{Any}
    for i in 1:size(programnames, 1)
        program = programnames[i]
        parent = first(program)
        children = String.(program[1:end .!= 1])
        for child in children
            weight = weights[weights[:,1] .== child, 2]
            graph = vcat(graph, [parent child weight])
        end
    end

    for item in graph[:, 1]
        if !in(item, graph[:, 2])
            global partone
            partone = item
            break
        end
    end

    # Walk the tree
    parent = partone
    children = graph[graph[:,1] .== parent, 2]

    return [partone, nothing]
end

function walktree(node)
    children = graph[graph[:,1] .== node, 2]
    if size(x, 1) == 0
        weight = first(graph[graph[:, 2] .== node, 3])
        return weight
    else
        weights = walktree.(children)
        return sum(weights)
    end
end


text:
function JLD.save(solver::Union{DeepQLearningSolver, DeepRecurrentQLearningSolver},
                  policy::Union{DQNPolicy, LSTMPolicy};
                  weights_file::String = "weights.jld",
                  problem_file::String = "problem.jld")
    saver = tf.train.Saver()
    save(problem_file, "solver", solver, "env", policy.env)
    train.save(saver, policy.sess, weights_file)
end

function restore(;problem_file::String="problem.jld", weights_file::String="weights.jld", graph=Graph())
    problem = load(problem_file)
    solver = problem["solver"]
    env = problem["env"]
    train_graph = build_graph(solver, env, graph)
    policy = restore_policy(env, solver, train_graph, weights_file)
    return policy
end

function restore_policy(env::AbstractEnvironment, solver::DeepQLearningSolver, train_graph::TrainGraph, weights_file::String)
    saver = train.Saver()
    tf.train.restore(saver, train_graph.sess, weights_file)
    policy = DQNPolicy(train_graph.q, train_graph.s, env, train_graph.sess)
end

function restore_policy(env::AbstractEnvironment, solver::DeepRecurrentQLearningSolver, train_graph::RecurrentTrainGraph, weights_file::String)
    saver = train.Saver()
    tf.train.restore(saver, train_graph.sess, weights_file)
    policy = LSTMPolicy(train_graph.sess, env, solver.arch, solver.dueling)
end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
################################################################################

# load sample
x = [
  10362.170502 2234.273811 2244.926224 647.264278 655.394653 571.437164 662.987797 486.684425 395.458192 392.330780 349.113530 403.472491 409.692116 500.905411 547.848330 789.415714;
  9246.340912 1466.947083 1906.052313 1209.676670 606.815348 520.541983 571.563523 557.675089 473.436043 402.273785 345.846071 505.062222 376.109373 439.998122 406.717417 398.954038;
  8926.864639 2024.670303 1371.395947 1277.803127 978.470009 456.044172 580.524843 546.832397 342.843889 706.637840 408.624054 324.005280 337.768628 331.042826 444.800913 449.033676;
  10628.121246 1498.386951 1187.115562 925.775265 694.869560 257.954501 592.700003 492.023813 452.887017 563.844071 473.371099 296.826092 479.452673 288.918503 369.519865 430.178178;
  11443.098831 1500.468246 1641.784470 929.626314 695.049788 473.568182 568.370424 542.348984 449.226664 351.930719 440.918041 346.611481 539.186420 325.506171 361.940255 451.275470;
  6898.439972 1674.825584 2142.319328 923.535397 562.254679 542.897036 544.597461 447.280304 338.127405 249.294769 331.028492 335.807143 421.090866 236.912484 241.648853 423.595156;
  9221.546188 1511.961166 1357.589783 1006.436895 532.545324 605.403816 468.057093 388.927434 439.460726 281.329487 469.301092 309.418703 340.099140 381.352773 227.883408 402.339500;
  11321.235947 1362.566475 944.705898 696.605307 411.700334 517.204907 383.636882 489.580448 301.713282 371.661125 476.757751 266.741815 154.978138 412.402929 279.845093 382.712366;
  8597.829269 1833.756462 1628.936111 1163.014809 697.425591 718.742764 431.683627 524.160545 522.474073 404.779799 286.859573 291.428293 371.560964 297.910174 426.403369 289.336908;
  9112.851349 1626.732288 1602.905491 1432.479465 663.499470 804.251534 739.028986 512.910912 397.432220 389.500046 407.156647 517.501736 347.876340 241.200329 416.171814 341.232077;
]

# setup hidden Markov model
hmm = setup(x)

################################################################################

max_stars_repo_name: DanielRivasMD/HiddenMarkovModelReaders

max_issues_repo_name: DanielRivasMD/HiddenMarkovModelReaders


text:
abstract type Copula{T} end
value_eltype(::Copula{T}) where T = T

"""
    simulate_copula(t::Int, copula::Copula; rng::AbstractRNG = Random.GLOBAL_RNG)

Returns t realizations from a given copula

```jldoctest
julia> using Random

julia> Random.seed!(43);

julia> simulate_copula(2, Gumbel_cop(3, 1.5))
2×3 Array{Real,2}:
 0.740038  0.918928  0.950674
 0.637826  0.483514  0.123949

julia> Random.seed!(43);

julia> simulate_copula(3, Clayton_cop(2, 1.))
3×2 Array{Float64,2}:
 0.562482  0.896247
 0.968953  0.731239
 0.749178  0.38015

 julia> Random.seed!(43);

 julia> simulate_copula(2, Clayton_cop(2, -.5))
 2×2 Array{Float64,2}:
  0.180975  0.818017
  0.888934  0.863358

julia> Random.seed!(43);

julia> simulate_copula(2, Gumbel_cop_rev(3, 1.5))
2×3 Array{Real,2}:
 0.259962  0.081072  0.0493259
 0.362174  0.516486  0.876051

  julia> Random.seed!(43);

 julia> simulate_copula(2, Clayton_cop_rev(2, -0.5))
 2×2 Array{Float64,2}:
   0.819025  0.181983
   0.111066  0.136642

julia> Random.seed!(43);

julia> simulate_copula(4, AMH_cop(2, 0.5))
4×2 Array{Float64,2}:
 0.483939  0.883911
 0.962064  0.665769
 0.707543  0.25042
 0.915491  0.494523

julia> Random.seed!(43);

julia> simulate_copula(4, AMH_cop(2, -0.5))
4×2 Array{Float64,2}:
 0.180975  0.820073
 0.888934  0.886169
 0.408278  0.919572
 0.828727  0.335864

julia> Random.seed!(43);

julia> simulate_copula(4, AMH_cop_rev(2, 0.5))
4×2 Array{Float64,2}:
 0.516061   0.116089
 0.0379356  0.334231
 0.292457   0.74958
 0.0845089  0.505477

julia> simulate_copula(4, AMH_cop_rev(2, -0.5))
4×2 Array{Float64,2}:
 0.819025  0.179927
 0.111066  0.113831
 0.591722  0.0804284
 0.171273  0.664136

julia> Random.seed!(43);

julia> simulate_copula(4, Frank_cop(2, 3.5))
4×2 Array{Float64,2}:
 0.650276  0.910212
 0.973726  0.789701
 0.690966  0.358523
 0.747862  0.29333

julia> Random.seed!(43);

julia> simulate_copula(4, Frank_cop(2, 0.2, SpearmanCorrelation))
4×2 Array{Float64,2}:
 0.504123  0.887296
 0.962936  0.678791
 0.718628  0.271543
 0.917759  0.51439

julia> Random.seed!(43);

julia> c = Chain_of_Archimedeans([4., 11.], "frank")
Chain_of_Archimedeans(3, [4.0, 11.0], ["frank", "frank"])

julia> simulate_copula(1, c)
1×3 Array{Float64,2}:
 0.180975  0.492923  0.679345

julia> c = Chain_of_Archimedeans([.5, .7], ["frank", "clayton"], KendallCorrelation)
Chain_of_Archimedeans(3, [5.736282707019972, 4.666666666666666], ["frank", "clayton"])

julia> Random.seed!(43);

julia> simulate_copula(1, c)
1×3 Array{Float64,2}:
 0.180975  0.408582  0.646887

julia> Random.seed!(43)

julia> simulate_copula(10, Chain_of_Frechet([0.6, 0.4], [0.3, 0.5]))
  10×3 Array{Float64,2}:
  0.996764  0.996764  0.996764
  0.204033  0.795967  0.204033
  0.979901  0.979901  0.0200985
  0.120669  0.879331  0.120669
  0.453027  0.453027  0.453027
  0.800909  0.199091  0.800909
  0.54892   0.54892   0.54892
  0.933832  0.933832  0.0661679
  0.396943  0.396943  0.396943
  0.804096  0.851275  0.955881

julia> Random.seed!(43);

julia> simulate_copula(10, Gaussian_cop([1. 0.5; 0.5 1.]))
10×2 Array{Real,2}:
 0.589188  0.815308
 0.708285  0.924962
 0.747341  0.156994
 0.227634  0.183116
 0.227575  0.957376
 0.271558  0.364803
 0.445691  0.52792
 0.585362  0.23135
 0.498593  0.48266
 0.190283  0.594451

julia> Random.seed!(43);

julia> simulate_copula(10, Student_cop([1. 0.5; 0.5 1.], 10))
10×2 Array{Real,2}:
 0.658199  0.937148
 0.718244  0.92602
 0.809521  0.0980325
 0.263068  0.222589
 0.187187  0.971109
 0.245373  0.346428
 0.452336  0.524498
 0.57113   0.272525
 0.498443  0.48082
 0.113788  0.633349

julia> Random.seed!(43);

julia> f = Frechet_cop(3, 0.5)
Frechet_cop(3, 0.5, 0.0)

julia> simulate_copula(1, f)
1×3 Array{Real,2}:
0.180975  0.775377  0.888934

julia> Random.seed!(43);

julia> f = Frechet_cop(2, 0.5, 0.2)
Frechet_cop(2, 0.5, 0.2)

julia> simulate_copula(1, f)
1×2 Array{Real,2}:
0.180975  0.775377

julia> Random.seed!(43);

julia> cop = Marshall_Olkin_cop([1.,2.,3.])
Marshall_Olkin_cop(2, [1.0, 2.0, 3.0])

julia> simulate_copula(1, cop)
1×2 Array{Float64,2}:
  0.854724  0.821831

julia> Random.seed!(43);

julia> c = Hierarchical_Gumbel_cop([5., 4., 3.])
Hierarchical_Gumbel_cop(4, [5.0, 4.0, 3.0])

julia> simulate_copula(3, c)
3×4 Array{Float64,2}:
 0.100353  0.207903  0.0988337  0.0431565
 0.347417  0.217052  0.223734   0.042903
 0.73617   0.347349  0.168348   0.410963

julia> Random.seed!(43);

julia> c1 = Clayton_cop(2, 2.)
Clayton_cop(2, 2.0)

julia> c2 = Clayton_cop(2, 3.)
Clayton_cop(2, 3.0)

julia> cp = Nested_Clayton_cop([c1, c2], 1, 1.1)
Nested_Clayton_cop(Clayton_cop[Clayton_cop(2, 2.0), Clayton_cop(2, 3.0)], 1, 1.1)

julia> simulate_copula(4, cp)
4×5 Array{Float64,2}:
0.514118  0.84089   0.870106  0.906233  0.739349
0.588245  0.85816   0.935308  0.944444  0.709009
0.59625   0.665947  0.483649  0.603074  0.153501
0.200051  0.304099  0.242572  0.177836  0.0851603

julia> c1 = AMH_cop(2, .7)
AMH_cop(2, 0.7)

julia> c2 = AMH_cop(2, .8)
AMH_cop(2, 0.8)

julia> cp = Nested_AMH_cop([c1, c2], 1, 0.2)
Nested_AMH_cop(AMH_cop[AMH_cop(2, 0.7), AMH_cop(2, 0.8)], 1, 0.2)

julia> Random.seed!(43);

julia> simulate_copula(4, cp)
4×5 Array{Float64,2}:
 0.557393  0.902767  0.909853  0.938522  0.586068
 0.184204  0.866664  0.699134  0.226744  0.102932
 0.268634  0.383355  0.179023  0.533749  0.995958
 0.578143  0.840169  0.743728  0.963226  0.576695

julia> c1 = Frank_cop(2, 4.)
Frank_cop(2, 4.0)

julia> c2 = Frank_cop(2, 5.)
Frank_cop(2, 5.0)

julia> c = Nested_Frank_cop([c1, c2],1, 2.0)
Nested_Frank_cop(Frank_cop[Frank_cop(2, 4.0), Frank_cop(2, 5.0)], 1, 2.0)

julia> Random.seed!(43);

julia> simulate_copula(1, c)
1×5 Array{Float64,2}:
 0.642765  0.901183  0.969422  0.9792  0.74155

julia> Random.seed!(43);

julia> c1 = Gumbel_cop(2, 2.)
Gumbel_cop(2, 2.0)

julia> c2 = Gumbel_cop(2, 3.)
Gumbel_cop(2, 3.0)

julia> cp = Nested_Gumbel_cop([c1, c2], 1, 1.1)
Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 2.0), Gumbel_cop(2, 3.0)], 1, 1.1)

julia> simulate_copula(4, cp)
4×5 Array{Float64,2}:
 0.387085   0.693399   0.94718   0.953776  0.583379
 0.0646972  0.0865914  0.990691  0.991127  0.718803
 0.966896   0.709233   0.788019  0.855622  0.755476
 0.272487   0.106996   0.756052  0.834068  0.661432

julia> a = Gumbel_cop(2, 5.)
Gumbel_cop(2, 5.0)

julia> b = Gumbel_cop(2, 6.)
Gumbel_cop(2, 6.0)

julia> c = Gumbel_cop(2, 5.5)
Gumbel_cop(2, 5.5)

julia> p1 = Nested_Gumbel_cop([a,b], 1, 2.)
Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.0), Gumbel_cop(2, 6.0)], 1, 2.0)

julia> p2 = Nested_Gumbel_cop([c], 2, 2.1)
Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.5)], 2, 2.1)

julia> copula = Double_Nested_Gumbel_cop([p1, p2], 1.5)
Double_Nested_Gumbel_cop(Nested_Gumbel_cop[Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.0), Gumbel_cop(2, 6.0)], 1, 2.0), Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.5)], 2, 2.1)], 1.5)

julia> Random.seed!(43);

julia> simulate_copula(5, copula)
5×9 Array{Float64,2}:
 0.598555   0.671584   0.8403     0.846844  0.634609  0.686927  0.693906  0.651968    0.670812
 0.0518892  0.191236   0.0803859  0.104325  0.410727  0.529354  0.557387  0.370518    0.592302
 0.367914   0.276196   0.382616   0.470171  0.264135  0.144503  0.13097   0.00687015  0.01417
 0.632727   0.596879   0.244176   0.338809  0.58771   0.147539  0.219099  0.287937    0.0569943
 0.310365   0.0483216  0.119312   0.107155  0.336619  0.279602  0.262756  0.438432    0.403061
```
"""
function simulate_copula(t, copula; rng = Random.GLOBAL_RNG)
    U = zeros(value_eltype(copula), t, copula.n)
    simulate_copula!(U, copula; rng = rng)
    return U
end


text:
module UBI
#=
Might as well do this while we're at it.
A very basic basic income; based roughly on

Painter, Anthony, Jamie Cooke, and Ahmed Aima. 2019. ‘A Basic Income for Scotland’. Royal Society for the encouragement of Arts, Manufactures and Commerce. https://www.thersa.org/globalassets/pdfs/rsa-a-basic-income-for-scotland.pdf.
=#

using ScottishTaxBenefitModel
using .ModelHousehold:
    BenefitUnit,
    Household,
    Person,
    get_benefit_units

using .STBParameters:
    UBISys

using .STBIncomes
using .Definitions

using .Results: BenefitUnitResult, HouseholdResult

export calc_UBI!

function calc_UBI!(
    benefit_unit_result :: BenefitUnitResult,
    benefit_unit        :: BenefitUnit,
    ubisys              :: UBISys )
    scp = 0.0
    bu = benefit_unit
    bur = benefit_unit_result # shortcuts
    for (pid,pers) in bu.people
        if pers.age < ubisys.adult_age
            bi = ubisys.child_amount
        elseif pers.age < ubisys.retirement_age
            bi = ubisys.adult_amount
        else
            bi = ubisys.universal_pension
        end
        bur.pers[pid].income[BASIC_INCOME] = bi
     end
end

function calc_UBI!(
    household_result :: HouseholdResult,
    hh               :: Household,
    ubisys           :: UBISys )
    if ubisys.abolished
        return
    end
    bus = get_benefit_units( hh )
    for bn in eachindex( bus )
        calc_UBI!(
            household_result.bus[bn],
            bus[bn],
            ubisys
        )
    end
end

end # module


text:
let
    name = "Turing"
    uuid = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
    url = "https://github.com/TuringLang/Turing.jl"
    julia = true
    license = "MIT"
    is_wrapper = false
    metadata_pkg(TuringModel; name, uuid, url, julia, license, is_wrapper)
end

let
    input = Table(Continuous)
    target = AbstractVector{Continuous}
    descr = "Bayesian model defined via Turing"
    path = "$PKG.TuringModel"
    metadata_model(TuringModel; input, target, descr, path)
end


text:
@testset "Client" begin
    @testset "Token" begin
        # Tokens should always be prepended by "Bot ".
        c = Client("token")
        @test c.token == "Bot token"
        # Unless it's already there.
        c = Client("Bot token")
        @test c.token == "Bot token"
    end

    @testset "logkws" begin
        # We don't include the connection version or shard unless they're relevant.
        kws = Dict(logkws(c))
        @test haskey(kws, :time)
        @test !haskey(kws, :shard)
        @test !haskey(kws, :conn)
        c.shards += 1
        c.conn.io = IOBuffer()
        c.conn.v = 1
        kws = Dict(logkws(c))
        @test get(kws, :shard, nothing) == 0
        @test get(kws, :conn, nothing) == 1

        # We can add our own keywords.
        kws = Dict(logkws(c; x=1))
        @test get(kws, :x, nothing) == 1

        # We can also overwrite defaults.
        kws = logkws(c; time=0)
        @test length(count(p -> p.first === :time, kws)) == 1
        @test get(Dict(kws), :time, nothing) == 0

        # And we can hide defaults with undef.
        kws = Dict(logkws(c; time=undef))
        @test !haskey(kws, :time)
    end

    @testset "Utils" begin
        # We haven't connected to the gateway yet.
        @test me(c) === nothing

        disable_cache!(c)
        @test !c.use_cache
        enable_cache!(c) do
            @test c.use_cache
        end
        @test !c.use_cache
        enable_cache!(c)
        @test c.use_cache
        disable_cache!(c) do
            @test !c.use_cache
        end
        @test c.use_cache
    end
end


text:
module AVRDUDE

const version = "0.0.1"

bin_dir = Pkg.dir("AVRDUDE", "deps", "avrdude", "v"*version)

const bin_location = abspath(joinpath(bin_dir, "avrdude", @winows ? ".exe" : ""))

end # module


text:
using Thermodynamics
using Thermodynamics.TemperatureProfiles
using Thermodynamics.TestedProfiles
using UnPack
using CLIMAParameters
using CLIMAParameters.Planet
using Plots
import Thermodynamics
Thermodynamics.print_warning() = false

struct EarthParameterSet <: AbstractEarthParameterSet end;
const param_set = EarthParameterSet();
FT = Float64;
thermo_dir = dirname(dirname(pathof(Thermodynamics)));
include(joinpath(thermo_dir, "docs", "plot_helpers.jl"));
profiles = TestedProfiles.PhaseEquilProfiles(param_set, Array{FT});
@unpack ρ, e_int, q_tot = profiles

dims = (10, 10, 10);
ρ = range(min(ρ...), stop = max(ρ...), length = dims[1]);
e_int = range(min(e_int...), stop = max(e_int...), length = dims[2]);
q_tot = range(min(q_tot...), stop = max(q_tot...), length = dims[3]);

ρ_all = Array{FT}(undef, prod(dims));
e_int_all = Array{FT}(undef, prod(dims));
q_tot_all = Array{FT}(undef, prod(dims));

linear_indices = LinearIndices((1:dims[1], 1:dims[2], 1:dims[3]));
TS = Array{Union{ThermodynamicState, Nothing}}(undef, prod(dims));
TS_no_err = Array{ThermodynamicState}(undef, prod(dims));

@inbounds for i in linear_indices.indices[1]
    @inbounds for j in linear_indices.indices[2]
        @inbounds for k in linear_indices.indices[3]
            p = linear_indices[i, j, k]
            ρ_all[p] = ρ[i]
            e_int_all[p] = e_int[j]
            q_tot_all[p] = q_tot[k]

            Thermodynamics.error_on_non_convergence() = false
            TS_no_err[p] = PhaseEquil_ρeq(param_set, ρ[i], e_int[j], q_tot[k])
            Thermodynamics.error_on_non_convergence() = true
            # @show p/prod(linear_indices.indices)*100
            try
                TS[p] = PhaseEquil_ρeq(param_set, ρ[i], e_int[j], q_tot[k])
            catch
                TS[p] = nothing
            end
        end
    end
end

# Full 3D scatter plot
function save_masked_plot3D(TS_no_err, mask, title, filename)
    ρ_mask = ρ_all[mask]
    T_mask = air_temperature.(TS_no_err[mask])
    q_tot_mask = q_tot_all[mask]
    pts = (ρ_mask, T_mask, q_tot_mask)
    Plots.plot(pts..., seriestype = :scatter, markersize = 7)
    plot!(
        xlabel = "Density",
        ylabel = "Temperature",
        zlabel = "Total specific humidity",
        title = "$title",
    )
    savefig(filename)
end;

save_masked_plot3D(TS_no_err, TS .== nothing, "NC", "Scatter3DNonConverged.svg");
save_masked_plot3D(TS_no_err, TS .!= nothing, "C", "Scatter3DConverged.svg");

# 2D binned scatter plots
function save_masked_plot2D_slices(TS_no_err, mask, title, filename)
    ρ_mask = ρ_all[mask]
    T_mask = air_temperature.(TS_no_err[mask])
    q_tot_mask = q_tot_all[mask]
    save_binned_surface_plots(ρ_mask, T_mask, q_tot_mask, title, filename)
end;

save_masked_plot2D_slices(
    TS_no_err,
    TS .== nothing,
    "NC",
    "Slices2DNonConverged.svg",
);
save_masked_plot2D_slices(
    TS_no_err,
    TS .!= nothing,
    "C",
    "Slices2DConverged.svg",
);

@warn "Note that the temperature axis for the non-converged
plot is not necessarily accurate, since the temperatures are
the result of a non-converged saturation adjustment"


text:
module RedBlackGraph

include("avus.jl")
include("shortestpaths/floyd-warshall.jl")
include("relationship.jl")
include("util/relationship_reader.jl")

export
# Avus number system
red_one, AInteger, AInt8, AInt16, AInt32, AInt64, AInt128, MSB,

# Shortest Path / Transitive Closure
floyd_warshall_transitive_closure, floyd_warshall_transitive_closure!,

# Relationship calculation
lookup_relationship, calculate_relationship, Relationship,

# Utilities
read_graph

end