text:
# functions to add hybridizations to a semi-directed network
# subject to topological constraints; no level restriction

"""
    addhybridedge!(net::HybridNetwork, nohybridladder::Bool, no3cycle::Bool,
                   constraints=TopologyConstraint[]::Vector{TopologyConstraint};
                   maxattempts=10::Int, fixroot=false::Bool)

Randomly choose two edges in `net` then: add hybrid edge from edge 1 to edge 2
of length 0.01. The two halves of edge 1 (and of edge 2) have equal lengths.
The hybrid partner edge (top half of edge 2, if fixroot is true) will point
towards the newly-created node on the middle of the original edge 2.

If the resulting network is a DAG, satisfies the constraint(s),
does not contain any 3-cycle (if `no3cycle=true`), and does not have
a hybrid ladder (if `nohybridladder=true`) then the proposal is successful:
`net` is modified, and the function returns the newly created hybrid node and
newly created hybrid edge.

If the resulting network is not acceptable, then a new set of edges
is proposed (using a blacklist) until one is found acceptable, or until
a maximum number of attempts have been made (`maxattempts`).
If none of the attempted proposals are successful, `nothing` is returned
(without causing an error).

After a pair of edges is picked, the "top" half of edge2 is proposed
as the partner hybrid edge with probability 0.8 if `fixroot` is false,
(to avoid changing the direction of edge2 with more than 50% chance)
and with probability 1.0 if `fixroot` is true.
If this choice does not work and if `fixroot` is false,
the other half of edge2 is proposed as the partner hybrid edge.
Note that choosing the "bottom" half of edge2 as the partner edge
requires to flip the direction of edge 2, and to move the root accordingly
(to the original child of edge2).

# examples

```jldoctest
julia> net = readTopology("((S1,(((S2,(S3)#H1),(#H1,S4)))#H2),(#H2,S5));");

julia> using Random

julia> Random.seed!(134);

julia> PhyloNetworks.addhybridedge!(net, true, true)
(PhyloNetworks.Node:
 number:9
 name:H3
 hybrid node
 attached to 3 edges, numbered: 5 16 17
, PhyloNetworks.Edge:
 number:17
 length:0.01
 minor hybrid edge with gamma=0.028819678088948808
 attached to 2 node(s) (parent first): 8 9
)

julia> writeTopology(net, round=true)
"((S1,((((#H1,S4),((S2,(S3)#H1))#H3:::0.971))#H2,#H3:0.01::0.029)),(#H2,S5));"

```
"""
function addhybridedge!(net::HybridNetwork, nohybridladder::Bool, no3cycle::Bool,
        constraints=TopologyConstraint[]::Vector{TopologyConstraint};
        maxattempts=10::Int, fixroot=false::Bool)
    all(con.type == 1 for con in constraints) || error("only type-1 constraints implemented so far")
    numedges = length(net.edge)
    blacklist = Set{Tuple{Int,Int}}()
    nmax_blacklist = numedges * (numedges-1) # all sets of edge1 -> edge2
    nattempts = 0
    while nattempts < maxattempts && length(blacklist) < nmax_blacklist
        e1 = Random.rand(1:numedges) # presumably faster than Random.randperm or Random.shuffle
        edge1 = net.edge[e1]
        e2 = Random.rand(1:(numedges-1)) # e2 must be different from e1: only numedges-1 options
        edge2 = net.edge[(e2<e1 ? e2 : e2+1)]
        (e1,e2) ∉ blacklist || # try another pair without adding to the # of attempts
            continue           # if (e1,e2) was already attempted
        nattempts += 1
        ## check that constraints are met
        p1 = getParent(edge1)
        p2 = getParent(edge2)
        constraintsmet = true
        for con in constraints
            if con.type == 1 # forbid going out of (edge1) or into (edge2) the species group
                if con.node === p1 || con.node === p2
                    push!(blacklist, (e1,e2))
                    constraintsmet = false
                    break # of constraint loop
                end
            end
        end
        constraintsmet || continue # try another pair if constraints not met
        ## check for no 3-cycle, if no3cycle is requested
        if no3cycle && hybrid3cycle(edge1, edge2) # does not depend on which partner edge is chosen
            push!(blacklist, (e1,e2))
            continue
        end
        ## check for no hybrid ladder, if requested:
         # edge2 cannot be a hybrid edge or the child of a hybrid node
        if nohybridladder && (edge2.hybrid || getParent(edge2).hybrid)
            push!(blacklist, (e1,e2))
            continue
        end
        hybridpartnernew = (fixroot ? true : rand() > 0.2) # if true: partner hybrid = new edge above edge 2
        ## check that the new network will be a DAG: no directional conflict
        if directionalconflict(net, p1, edge2, hybridpartnernew)
            if fixroot # don't try to change the direction of edge2
                push!(blacklist, (e1,e2))
                continue
            end # else: try harder: change direction of edge2 and move root
            hybridpartnernew = !hybridpartnernew # try again with opposite
            if directionalconflict(net, p1, edge2, hybridpartnernew)
                push!(blacklist, (e1,e2))
                continue
            end # else: switching hybridpartnernew worked
        end
        newgamma = rand()*0.5; # in (0,.5) to create a minor hybrid edge
        return addhybridedge!(net, edge1, edge2, hybridpartnernew, 0.01, newgamma)
    end
    # if we get here: none of the max number of attempts worked - return nothing.
    # if an attempt had worked, we would have returned something.
    return nothing
end

"""
    addhybridedge!(net::HybridNetwork, edge1::Edge, edge2::Edge, hybridpartnernew::Bool,
                   edgelength=-1.0::Float64, gamma=-1.0::Float64)

Add hybridization to `net` coming from `edge1` going into `edge2`.
2 new nodes and 3 new edges are created: `edge1` are `edge2` are both cut into 2 edges,
and a new edge is created linking the 2 new "middle" nodes, pointing from `edge1` to `edge2`.
The new node in the middle of `edge1` is a tree node.
The new node in the middle of `edge2` is a hybrid node.
Its parent edges are the newly created hybrid edge (with γ = gamma, missing by default),
and either the newly edge "above" `edge2` if `hybridpartnernew=true`,
or the old `edge2` otherwise (which would reverse the direction of `edge2` and others).

Should be called from the other method, which performs a bunch of checks.
Updates `containRoot` attributes for edges below the new hybrid node.

Output: new hybrid node (middle of the old `edge2`) and new hybrid edge.

# examples

```jldoctest
julia> net = readTopology("((S8,(((S1,(S5)#H1),(#H1,S6)))#H2),(#H2,S10));");

julia> hybnode, hybedge = PhyloNetworks.addhybridedge!(net, net.edge[13], net.edge[8], true, 0.0, 0.2)
(PhyloNetworks.Node:
 number:9
 name:H3
 hybrid node
 attached to 3 edges, numbered: 8 16 17
, PhyloNetworks.Edge:
 number:17
 length:0.0
 minor hybrid edge with gamma=0.2
 attached to 2 node(s) (parent first): 8 9
)

julia> writeTopology(net)
"((S8,(((S1,(S5)#H1),((#H1,S6))#H3:::0.8))#H2),(#H2,(S10,#H3:0.0::0.2)));"

```
"""
function addhybridedge!(net::HybridNetwork, edge1::Edge, edge2::Edge, hybridpartnernew::Bool,
                        edgelength::Float64=-1.0, gamma::Float64=-1.0)
    gamma == -1.0 || (gamma <= 1.0 && gamma >= 0.0) || error("invalid γ to add a hybrid edge")
    gbar = (gamma == -1.0 ? -1.0 : 1.0 - gamma) # 1-gamma, with γ=-1 as missing
    newnode1_tree, edgeabovee1 = breakedge!(edge1, net) # new tree node
    newnode2_hybrid, edgeabovee2 = breakedge!(edge2, net) # new hybrid node
    newnode2_hybrid.hybrid = true
    pushHybrid!(net, newnode2_hybrid) # updates net.hybrid and net.numHybrids
    # new hybrid edge, minor if γ missing (-1)
    hybrid_edge = Edge(maximum(e.number for e in net.edge) + 1, edgelength, true, gamma, gamma>0.5) # number, length, hybrid, gamma, isMajor
    # partner edge: update hybrid status, γ and direction
    if hybridpartnernew
        edgeabovee2.hybrid = true
        edgeabovee2.gamma = gbar
        if gamma>0.5
            edgeabovee2.isMajor = false
        end
    else
        c2 = getChild(edge2) # child of edge2 before we switch its direction
        i2 = findfirst(isequal(c2), net.node)
        net.root = i2 # makes c2 the new root node
        edge2.hybrid = true
        edge2.gamma = gbar
        if gamma>0.5
            edge2.isMajor = false
        end
        edge2.isChild1 =  !edge2.isChild1 # reverse the direction of edge2
        edgeabovee2.isChild1 = !edgeabovee2.isChild1
    end
    # parse hyb names to find the next available. assignhybridnames! would do them all
    rx = r"^H(\d+)$"
    hnum = net.numHybrids # to name the new hybrid, potentially
    for n in net.node
        m = match(rx, n.name)
        if m !== nothing
            hi = parse(Int, m[1])
            hnum > hi || (hnum = hi+1)
        end
    end
    newnode2_hybrid.name = "H$hnum"
    setNode!(hybrid_edge, [newnode2_hybrid, newnode1_tree]) # [child node, parent node] to match isChild1=true
    setEdge!(newnode1_tree, hybrid_edge)
    setEdge!(newnode2_hybrid, hybrid_edge)
    pushEdge!(net, hybrid_edge)
    if hybridpartnernew
        # @debug "triple-check it's a DAG" directEdges!(net)
        norootbelow!(edge2)
    else
        directEdges!(net)
        norootbelow!(edgeabovee2)
    end
    return newnode2_hybrid, hybrid_edge
end

"""
    hybrid3cycle(edge1::Edge, edge2::Edge)

Check if proposed hybrid edge from `edge1` into `edge2` would create a 3 cycle,
that is, if `edge1` and `edge2` have a node in common.
(This move cannot create a 2-cycles because new nodes would be created in the
middle of edges 1 and 2.)
"""
function hybrid3cycle(edge1::Edge, edge2::Edge)
    !isempty(findall(in(edge1.node), edge2.node))
end

"""
    directionalconflict(net::HybridNetwork, parent::Node, edge::Edge,
                        hybridpartnernew::Bool)

Check if creating a hybrid edge down of `parent` node into the middle of `edge`
would create a directed cycle in `net`, i.e. not a DAG. The proposed hybrid
would go in the direction of `edge` down its child node if `hybridpartnernew`
is true. Otherwise, both halves of `edge` would have their direction reversed,
for the hybrid to go towards the original parent node of `edge`.
Does *not* modify the network.

Output: `true` if a conflict would arise (non-DAG), `false` if no conflict.
"""
function directionalconflict(net::HybridNetwork, parent::Node, edge2::Edge, hybridpartnernew::Bool)
    if hybridpartnernew # all edges would retain their directions: use isChild1 fields
        c2 = getChild(edge2)
        return parent === c2 || isdescendant(parent, c2)
    else # after hybrid addition, edge 2 would be reversed: "up" toward its own parent
        if !edge2.containRoot || edge2.hybrid
            return true # direction of edge2 cannot be reversed
        else # net would be a DAG with reversed directions, could even be rooted on edge2
            p2 = getParent(edge2)
            return parent === p2 || isdescendant_undirected(parent, p2, edge2)
        end
    end
end


text:
include("online_stats_plots.jl")

"""
    TurkieCallback(args...; kwargs....)

## Arguments
- Option 1 :
    `model::DynamicPPL.Model, plots::Series/AbstractVector=[:histkde, Mean(Float32), Variance(Float32), AutoCov(20, Float32)]`

For each of the variables of the given model each `plot` from `plots` will be plotted
Multidimensional variable will be automatically have indices added to them
- Option 2 :
    `vars::NamedTuple/Dict`
Will plot each pair of symbol and series of plots.
Note that for multidimensional variable you should pass a Symbol as `Symbol("m[1]")` for example.
See the docs for some examples.
## Keyword arguments
- `window=1000` : Use a window for plotting the trace
- `refresh=false` : Restart the plots from scratch everytime `sample` is called again (still WIP)
"""
TurkieCallback

struct TurkieCallback{TN<:NamedTuple,TD<:AbstractDict}
    figure::Figure
    data::Dict{Symbol, MovingWindow}
    axis_dict::Dict
    vars::TN
    params::TD
    iter::Observable{Int}
end

function TurkieCallback(model::Model, plots::Series; kwargs...)
    return TurkieCallback(model, collect(plots.stats); kwargs...)
end

function TurkieCallback(model::Model, plots::AbstractVector = [:histkde, Mean(Float32), Variance(Float32), AutoCov(20, Float32)]; kwargs...)
    vars, vals = _params_to_array([VarInfo(model)])
    return TurkieCallback(
        (;Pair.(vars, Ref(plots))...); # Return a named Tuple
        kwargs...
    )
end

function TurkieCallback(vars::Union{Dict, NamedTuple}; kwargs...)
    return TurkieCallback((;vars...), Dict{Symbol,Any}(kwargs...))
end

function TurkieCallback(vars::NamedTuple, params::Dict)
# Create a scene and a layout
    outer_padding = 5
    resolution = get!(params, :resolution, (1200, 700))
    fig = Figure(;resolution=resolution, figure_padding=outer_padding)
    window = get!(params, :window, 1000)
    refresh = get!(params, :refresh, false)
    params[:t0] = 0
    iter = Observable(0)
    data = Dict{Symbol, MovingWindow}(:iter => MovingWindow(window, Int))
    axis_dict = Dict()
    for (i, variable) in enumerate(keys(vars))
        plots = vars[variable]
        data[variable] = MovingWindow(window, Float32)
        axis_dict[(variable, :varname)] = fig[i, 1, Left()] = Label(fig, string(variable), textsize = 30)
        axis_dict[(variable, :varname)].padding = (0, 60, 0, 0)
        onlineplot!(fig, axis_dict, plots, iter, data, variable, i)
    end
    on(iter) do i
        if i > 1 # To deal with autolimits a certain number of samples are needed
            for variable in keys(vars)
                for p in vars[variable]
                    autolimits!(axis_dict[(variable, p)])
                end
            end
        end
    end
    MakieLayout.trim!(fig.layout)
    display(fig)
    return TurkieCallback(fig, data, axis_dict, vars, params, iter)
end

function addIO!(cb::TurkieCallback, io)
    cb.params[:io] = io
end

function (cb::TurkieCallback)(rng, model, sampler, transition, state, iteration; kwargs...)
    if iteration == 1
        if cb.params[:refresh]
            refresh_plots!(cb)
        end
        cb.params[:t0] = cb.iter[]
    end
    fit!(cb.data[:iter], iteration + cb.params[:t0]) # Update the iteration value
    for (variable, val) in zip(_params_to_array([transition])...)
        if haskey(cb.data, variable) # Check if symbol should be plotted
            fit!(cb.data[variable], Float32(val)) # Update its value
        end
    end
    cb.iter[] = cb.iter[] + 1
    if haskey(cb.params, :io)
        recordframe!(cb.params[:io])
    end
end

function refresh_plots!(cb)
    #TODO
end

ext: jl


text:
#= Sigmoid.jl

a quick script to perform a sigmoid transformation using the equation:

σ(X) = 1 /(1 + ℯ ^ vX)

function:
    sigmoid(signal;v)

where signal is the input timeseries and the parameter v determines the how shallow /
steep the sigmoid curve will be. Smaller values of v flatten the curve

From França et al (2018; https://doi.org/10.3389/fphys.2018.01767):
    The parameter v was chosen based on its effect on the estimated multifractal width for
    three types of time series: icEEG (NHNN1-channel 1), surrogate EEG (temporally
    shuffled values of the original time series from NHNN1-channel 1) and a simulated
    random series (with the same mean and variance), across the range v = [0.1, 2.0]
    in steps of 0.1. To find the optimal value for the parameter v, we needed to balance
    the trade-off between the three series in terms of presenting the most distinct Δα
    values (Appendix Figure E3A in Supplementary Material), while showing minimum
    distortion on the recording, or maximum correlation with the original time series

=#

function sigmoid(signal; v)
    sigmoid_transform = 1 ./(1 .+ ℯ .^(v .* signal))
    return sigmoid_transform
end


text:
# -*- mode: julia; -*-
module WAVPlay
import WAV: wavplay, wavwrite

# some standard Win32 API types and constants, see [MS-DTYP]
const BOOL = Cint
const DWORD = Culong
const TRUE = 1
# PlaySound flags from Winmm.h
const SND_SYNC      = 0x0
const SND_ASYNC     = 0x1
const SND_NODEFAULT = 0x2
const SND_MEMORY    = 0x4
const SND_FILENAME  = 0x20000

function wavplay(data, fs)
    # produce an in-memory WAV file ...
    buf=IOBuffer()
    wavwrite(data, buf, Fs=fs)
    wav = take!(buf)
    # ... and pass it to PlaySound
    success = ccall((:PlaySoundA, "Winmm.dll"), stdcall, BOOL,
                    (Ptr{Cvoid}, Ptr{Cvoid}, DWORD),
                    wav, C_NULL, SND_MEMORY | SND_SYNC | SND_NODEFAULT)
    Base.windowserror("PlaySound", success != TRUE)
end
end # module


text:
# Testing discrete univariate distributions

using Distributions
import JSON
using Base.Test
using Compat

function verify_and_test_drive(jsonfile, selected, n_tsamples::Int)
    R = JSON.parsefile(jsonfile)
    for dct in R
        ex = parse(dct["expr"])
        @assert ex.head == :call
        dsym = ex.args[1]
        dname = string(dsym)

        # test whether it is included in the selected list
        # or all are selected (when selected is empty)
        if !isempty(selected) && !(dname in selected)
            continue
        end

        # perform testing
        println("    testing $(ex)")
        dtype = eval(dsym)
        d = eval(ex)
        if dtype == TruncatedNormal
            @test isa(d, Truncated{Normal{Float64}})
        else
            @assert isa(dtype, Type) && dtype <: UnivariateDistribution
            @test isa(d, dtype)
        end

        # verification and testing
        verify_and_test(dtype, d, dct, n_tsamples)
    end
end

_parse_x(d::DiscreteUnivariateDistribution, x) = round(Int, x)
_parse_x(d::ContinuousUnivariateDistribution, x) = Float64(x)

_json_value(x::Number) = x

_json_value(x::AbstractString) =
    x == "inf" ? Inf :
    x == "-inf" ? -Inf :
    x == "nan" ? NaN :
    error("Invalid numerical value: $x")

function verify_and_test(D::Union{Type,Function}, d::UnivariateDistribution, dct::Dict, n_tsamples::Int)
    # verify properties
    #
    # Note: properties include all applicable params and stats
    #

    # D can be a function, e.g. TruncatedNormal
    if isa(D, Type)
        assert(isa(d, D))
    end

    # test various constructors for promotion, all-Integer args, etc.
    pars = params(d)

    # promotion constructor:
    float_pars = map(x -> isa(x, AbstractFloat), pars)
    if length(pars) > 1 && sum(float_pars) > 1
        mixed_pars = Any[pars...]
        first_float = findfirst(float_pars)
        mixed_pars[first_float] = Float32(mixed_pars[first_float])

        @test typeof(D(mixed_pars...)) == typeof(d)
    end

    # promote integer arguments to floats, where applicable
    if sum(float_pars) >= 1 && !any(map(isinf, pars)) && !isa(d, Geometric)
        int_pars = map(x -> ceil(Int, x), pars)
        @test typeof(D(int_pars...)) == typeof(d)
    end

    # verify properties (params & stats)
    pdct = dct["properties"]
    for (fname, val) in pdct
        expect_v = _json_value(val)
        f = eval(Symbol(fname))
        @assert isa(f, Function)
        tol_v = abs(expect_v) * 1e-8 + 1e-12
        Base.Test.test_approx_eq(f(d), expect_v, tol_v, "$fname(d)", "expect_v")
    end

    # verify logpdf and cdf at certain points
    pts = dct["points"]
    for pt in pts
        x = _parse_x(d, pt["x"])
        p = _json_value(pt["pdf"])
        lp = _json_value(pt["logpdf"])
        cf = _json_value(pt["cdf"])

        ptol = p * 1e-8 + 1e-16
        lptol = 1e-12
        cftol = 1e-12
        if isa(d, NoncentralHypergeometric)
            lptol = 1e-4
            cftol = 1e-8
        end

        Base.Test.test_approx_eq(pdf(d, x), p, ptol, "logpdf(d, $x)", "lp")
        Base.Test.test_approx_eq(logpdf(d, x), lp, lptol, "logpdf(d, $x)", "lp")

        # cdf method is not implemented for Skellam & NormalInverseGaussian
        if !isa(d, Union{Skellam, NormalInverseGaussian})
            Base.Test.test_approx_eq(cdf(d, x), cf, cftol, "cdf(d, $x)", "cf")
        end
    end

    # verify quantiles
    if !isa(d, Union{Skellam, VonMises, NormalInverseGaussian})
        qts = dct["quans"]
        for qt in qts
            q = Float64(qt["q"])
            x = Float64(qt["x"])
            @test isapprox(quantile(d, q), x, atol=1.0e-8)
        end
    end

    try
        m = mgf(d,0.0)
        @test m == 1.0
    catch e
        isa(e, MethodError) || throw(e)
    end
    try
        c = cf(d,0.0)
        @test c == 1.0
        # test some extra values: should all be well-defined
        for t in (0.1,-0.1,1.0,-1.0)
            @test !isnan(cf(d,t))
        end
    catch e
        isa(e, MethodError) || throw(e)
    end

    # generic testing
    if isa(d, Cosine)
        n_tsamples = floor(Int, n_tsamples / 10)
    elseif isa(d, NoncentralBeta) ||
           isa(d, NoncentralChisq) ||
           isa(d, NoncentralF) ||
           isa(d, NoncentralT)
        n_tsamples = min(n_tsamples, 100)
    end

    if !isa(d, Union{Skellam,
                     VonMises,
                     NoncentralHypergeometric,
                     NormalInverseGaussian})
        test_distr(d, n_tsamples)
    end
end

## main

for c in ["discrete",
          "continuous"]

    title = string(uppercase(c[1]), c[2:end])
    println("    [$title]")
    println("    ------------")
    jsonfile = joinpath(dirname(@__FILE__), "ref", "$(c)_test.ref.json")
    verify_and_test_drive(jsonfile, ARGS, 10^6)
    println()
end


text:
using Test

using QuantumCircuits
using QuantumCircuits.QCircuits.Circuit

qc1 = QCircuit(2)
qc1.x(0)
qc1.h(1)
qc1.cx(0, 1)
@test toPythonQiskitCode(qc1) == "qr1 = QuantumRegister(2)\ncr1 = ClassicalRegister(2)\nqc = QuantumCircuit(qr1, cr1)\nqc.x(0)\nqc.h(1)\nqc.cx(0, 1)\n"

qr1 = QuantumRegister(3, "q")
cr1 = ClassicalRegister(2, "c")
qc1 = QCircuit(qr1, cr1)
qc1.x(0)
qc1.h(1)
qc1.cx(0, 1)
qc1.measure([0, 1], [0, 1])
@test toPythonQiskitCode(qc1) == "qr1 = QuantumRegister(3, \"q\")\ncr1 = ClassicalRegister(2, \"c\")\nqc = QuantumCircuit(qr1, cr1)\nqc.x(0)\nqc.h(1)\nqc.cx(0, 1)\nqc.measure(0, 0)\nqc.measure(1, 1)\n"


text:
#=
Title: Discrete Simulation of 2-Dimensional Light Paths in Continuous Media
Author:
=#

#finds the magnitude of a vector in 2-space
function magnitude(vector) #vector type list
    return sqrt(vector[1]^2+vector[2]^2)
end

#finds the angle between two vectors
function get_angle(a,b)#a,b are type list
    if magnitude(a)>magnitude(b)
        hypotenuse = magnitude(a)
        adjacent = magnitude(b)
    else
        hypotenuse = magnitude(b)
        adjacent = magnitude(a)
    end
    if hypotenuse == 0
        hypotenuse = 0.00001
        break
    end
    return acos(adjacent / hypotenuse)
end

#finds the lowest angle between the norm and the velocity
function choose_angle(normal,velocity)
    positive = get_angle(normal,velocity)
    negative = get_angle(-normal,velocity)
    if abs(positive)<=abs(negative)
        return positive
    else
        return negative
    end
end

#rotates a matrix by some angle in radians
function instantaneous_rotation(current_velocity,current_angle)
    x = velocity[1]*cos(angle)-velocity[2]*sin(angle)
    y = velocity[1]*sin(angle)+velocity[2]*cos(angle)
    return [x,y]
end

#the continuous refractive index at any point
function refractive_index(point) #point type list (x,y)
    return (point[1]^2)*(point[2]^2)+1
end

#Gradient of level plane W(x,y,n) is normal to surface n
#the particular normal at some point
function projected_normal(point) #point type list (x,y)
    x = 2*point[1]*point[2]^2
    y = 2*point[2]*point[1]^2
    return [x,y] #we only need the projection onto 2-space
end

#generates a random initial velocity
function generate_velocity()
    candidate = [rand(-c:c)-rand(),0]
    candidate[2] = sqrt(c^2 - candidate[1]^2)
    return candidate
end

#generates a random initial position
function generate_position()
    return [rand(-10:10)+rand(),rand(-10:10)+rand()]
end

#=
Given n1 = c/v1 and n2 = c/v2, n1/n2 = (c/v1)(v2/c) = (v2/v1)(c/c)= (v2/v1)
so (v1 n1)/n2 = v2. Or in this case, n1 is the last n(x,z), v1 is the current
velocity, and n2 is the current n(x,z), and v2 is the velocity we are trying
to calculate.
=#

#Main Body of Code

#initial conditions of the photon:
c = 2.99e8
current_position,current_velocity = generate_position(),generate_velocity()
current_n = c / magnitude(current_velocity)

running = true
while running
    #have to figure out how to update current_n
    incident_normal = projected_normal(current_position)
    incident_angle = choose_angle(incident_normal,current_velocity)
    #the refractive angle between the refractive velocity and the refractive normal
    refractive_angle = asin((current_n*sin(incident_angle))/refractive_index(current_position))
    #=
    given you want to start with some current velocity, this function gives
    you the factor by which to multiply the current velocity, in order to get the
    desired refractive velocity
    =#
    scaling_factor = c/(magnitude(current_velocity)*refractive_index(current_position))
    next_velocity = scaling_factor*current_velocity#temporary... not the real next_velocity. still needs to be rotated

    current_velocity,current_position = next_velocity,next_position
    current_n = refractive_index(current_position)
end

#=
If p in R^2 is a function of x and y, then we can parametrize
x and y so that x = x(t) and y = y(t). Then, if v is current_velocity
in R^2, and v is parametrized with t as well, so that v[1] = x'{t}
and v[2] = y'(t). Then x'(t) is defined as the limit as h->0 of
(x(t+h)-x(t))/h. If instead of 0, we allow h->epsilon, then
x'(t) = (x(t+epsilon)-x(t))/epsilon. We can similarly say for y'(t)
that y'(t) = (y(t+epsilon)-y(t))/epsilon.
Then epsilon*x'(t) = x(t+epsilon) - x(t) and
epsilon*y'(t) = y(t+epsilon) - y(t).
Then epsilon*x'(t) + x(t) = x(t+epsilon)
and epsilon*y'(t) + y(t) = y(t+epsilon).
=#

#=
General sketch:
*photon = [current_velocity,current_position]
*use current_velocity and current_position to find next_velocity.
*use next_velocity to find next_position
*assign: current_velocity,current_position = next_velocity,next_position

currents = [current_position,current_velocity,current_n,current_angle]
nexts = [next_position,next_velocity,refractive_index(current_position),
next_angle]

refractive_index(current_position) = c / magnitude(next_velocity)
magnitude(next_velocity) = c / refractive_index(current_position)

next_velocity is proportional to current_velocity. Therefore
magnitude(next_velocity) = k*magnitude(current_velocity)
=#

max_stars_repo_name: AdityaSChanda/lightpaths

max_issues_repo_name: AdityaSChanda/lightpaths

max_forks_repo_name: AdityaSChanda/lightpaths


text:
@inline UnsafeAtomics.load(x) = UnsafeAtomics.load(x, seq_cst)
@inline UnsafeAtomics.store!(x, v) = UnsafeAtomics.store!(x, v, seq_cst)
@inline UnsafeAtomics.cas!(x, cmp, new) = UnsafeAtomics.cas!(x, cmp, new, seq_cst, seq_cst)
@inline UnsafeAtomics.modify!(ptr, op, x) = UnsafeAtomics.modify!(ptr, op, x, seq_cst)

#! format: off
# https://github.com/JuliaLang/julia/blob/v1.6.3/base/atomics.jl#L23-L30
if Sys.ARCH == :i686 || startswith(string(Sys.ARCH), "arm") ||
   Sys.ARCH === :powerpc64le || Sys.ARCH === :ppc64le
    const inttypes = (Int8, Int16, Int32, Int64,
                      UInt8, UInt16, UInt32, UInt64)
else
    const inttypes = (Int8, Int16, Int32, Int64, Int128,
                      UInt8, UInt16, UInt32, UInt64, UInt128)
end

# https://github.com/JuliaLang/julia/blob/v1.6.3/base/atomics.jl#L331-L341
const llvmtypes = IdDict{Any,String}(
    Bool => "i8",  # julia represents bools with 8-bits for now. # TODO: is this okay?
    Int8 => "i8", UInt8 => "i8",
    Int16 => "i16", UInt16 => "i16",
    Int32 => "i32", UInt32 => "i32",
    Int64 => "i64", UInt64 => "i64",
    Int128 => "i128", UInt128 => "i128",
    Float16 => "half",
    Float32 => "float",
    Float64 => "double",
)
#! format: on

const OP_RMW_TABLE = [
    (+) => :add,
    (-) => :sub,
    right => :xchg,
    (&) => :and,
    (⊼) => :nand,
    (|) => :or,
    (⊻) => xor,
    max => :max,
    min => :min,
]

for (op, rmwop) in OP_RMW_TABLE
    fn = Symbol(rmwop, "!")
    @eval @inline UnsafeAtomics.$fn(x, v) = UnsafeAtomics.$fn(x, v, seq_cst)
    @eval @inline UnsafeAtomics.$fn(ptr, x, ord) =
        first(UnsafeAtomics.modify!(ptr, $op, x, ord))
end

# Based on: https://github.com/JuliaLang/julia/blob/v1.6.3/base/atomics.jl
for typ in inttypes
    lt = llvmtypes[typ]
    rt = "$lt, $lt*"

    for ord in orderings
        ord in (release, acq_rel) && continue

        @eval function UnsafeAtomics.load(x::Ptr{$typ}, ::$(typeof(ord)))
            return llvmcall(
                $("""
                %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                %rv = load atomic $rt %ptr $ord, align $(sizeof(typ))
                ret $lt %rv
                """),
                $typ,
                Tuple{Ptr{$typ}},
                x,
            )
        end
    end

    for ord in orderings
        ord in (acquire, acq_rel) && continue

        @eval function UnsafeAtomics.store!(x::Ptr{$typ}, v::$typ, ::$(typeof(ord)))
            return llvmcall(
                $("""
                %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                store atomic $lt %1, $lt* %ptr $ord, align $(sizeof(typ))
                ret void
                """),
                Cvoid,
                Tuple{Ptr{$typ},$typ},
                x,
                v,
            )
        end
    end

    for success_ordering in (monotonic, acquire, release, acq_rel, seq_cst),
        failure_ordering in (monotonic, acquire, seq_cst)

        @eval function UnsafeAtomics.cas!(
            x::Ptr{$typ},
            cmp::$typ,
            new::$typ,
            ::$(typeof(success_ordering)),
            ::$(typeof(failure_ordering)),
        )
            success = Ref{Int8}()
            GC.@preserve success begin
                old = llvmcall(
                    $(
                        """
                        %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                        %rs = cmpxchg $lt* %ptr, $lt %1, $lt %2 $success_ordering $failure_ordering
                        %rv = extractvalue { $lt, i1 } %rs, 0
                        %s1 = extractvalue { $lt, i1 } %rs, 1
                        %s8 = zext i1 %s1 to i8
                        %sptr = inttoptr i$WORD_SIZE %3 to i8*
                        store i8 %s8, i8* %sptr
                        ret $lt %rv
                        """
                    ),
                    $typ,
                    Tuple{Ptr{$typ},$typ,$typ,Ptr{Int8}},
                    x,
                    cmp,
                    new,
                    Ptr{Int8}(pointer_from_objref(success)),
                )
            end
            return (old = old, success = !iszero(success[]))
        end
    end

    for (op, rmwop) in OP_RMW_TABLE
        rmw = string(rmwop)
        fn = Symbol(rmw, "!")
        if (rmw == "max" || rmw == "min") && typ <: Unsigned
            # LLVM distinguishes signedness in the operation, not the integer type.
            rmw = "u" * rmw
        end
        for ord in orderings
            @eval function UnsafeAtomics.modify!(
                x::Ptr{$typ},
                ::typeof($op),
                v::$typ,
                ::$(typeof(ord)),
            )
                old = llvmcall(
                    $("""
                    %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                    %rv = atomicrmw $rmw $lt* %ptr, $lt %1 $ord
                    ret $lt %rv
                    """),
                    $typ,
                    Tuple{Ptr{$typ},$typ},
                    x,
                    v,
                )
                return (old, $op(old, v))
            end
        end
    end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Git"
version = v"2.34.1"

# Collection of sources required to build Git
sources = [
    ArchiveSource("https://mirrors.edge.kernel.org/pub/software/scm/git/git-$(version).tar.xz",
                  "3a0755dd1cfab71a24dd96df3498c29cd0acd13b04f3d08bf933e81286db802c"),
    ArchiveSource("https://github.com/git-for-windows/git/releases/download/v$(version).windows.1/Git-$(version)-32-bit.tar.bz2",
                  "109d69b92e5383d40765064bcd25183af05b12c64a95597419a703004a7c8521"; unpack_target = "i686-w64-mingw32"),
    ArchiveSource("https://github.com/git-for-windows/git/releases/download/v$(version).windows.1/Git-$(version)-64-bit.tar.bz2",
                  "0d962f5894b94b93a966a2c12f77330e5f68aacae775059fb30dd61f2c08ef00"; unpack_target = "x86_64-w64-mingw32"),
]

# Bash recipe for building across all platforms
script = raw"""
install_license ${WORKSPACE}/srcdir/git-*/COPYING

if [[ "${target}" == *-ming* ]]; then
    # Fast path for Windows: just copy the content of the tarball to the prefix
    cp -r ${WORKSPACE}/srcdir/${target}/mingw${nbits}/* ${prefix}
    exit
fi

cd $WORKSPACE/srcdir/git-*/

# We need a native "tclsh" to cross-compile
apk add tcl

CACHE_VALUES=()
if [[ "${target}" != *86*-linux-* ]]; then
    # Cache values of some tests, let's hope to have got them right
    CACHE_VALUES+=(
        ac_cv_iconv_omits_bom=no
        ac_cv_fread_reads_directories=yes
        ac_cv_snprintf_returns_bogus=no
    )
else
    # Explain Git that we aren't cross-compiling on these platforms
    sed -i 's/cross_compiling=yes/cross_compiling=no/' configure
fi

if [[ "${target}" == *-apple-* ]]; then
    LDFLAGS="-lgettextlib"
elif [[ ${target} == *-freebsd* ]]; then
    LDFLAGS="-L${prefix}/lib -lcharset"
fi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --with-curl \
    --with-expat \
    --with-openssl \
    --with-iconv=${prefix} \
    --with-libpcre2 \
    --with-zlib=${prefix} \
    "${CACHE_VALUES[@]}" \
    CPPFLAGS="-I${prefix}/include" \
    LDFLAGS="${LDFLAGS}"
make -j${nproc}
make install INSTALL_SYMLINKS="yes, please"
"""

platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("git", :git),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # Need a host gettext for msgfmt
    HostBuildDependency("Gettext_jll"),
    Dependency("LibCURL_jll"),
    Dependency("Expat_jll"; compat="2.2.10"),
    Dependency("OpenSSL_jll"; compat="1.1.10"),
    # I believe Gettext is needed for macOS (and probably only here)
    Dependency("Gettext_jll"; compat="=0.21.0"),
    Dependency("Libiconv_jll"),
    Dependency("PCRE2_jll", v"10.35.0"),
    Dependency("Zlib_jll"),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Load the Julia package manager, which is itself a Julia package
using Pkg

try
    # Load DECAES
    @eval using DECAES
catch e
    # DECAES not found; install DECAES and try loading again
    Pkg.add("DECAES")
    @eval using DECAES
end

# Call the command line interface entrypoint function
main()


text:
using Revise
using gtpmip

graph = solve_hpd("test/hpd/domain.hpd", "test/hpd/problem.hpd"; max_levels=200)
1


text:
println("==== Test module JTetra ====")

include("testPersistency.jl")
include("testDeRhamComplex.jl")
include("testDerivatives.jl")

println("==== JTetra: All tests passed! ====")

max_stars_repo_name: JuliaInv/JTetra

max_issues_repo_name: JuliaInv/JTetra

max_forks_repo_name: JuliaInv/JTetra


text:
function solve()
    n, m = [parse(Int, x) for x in split(readline())]
    edges = [Set{Int}() for _ in 1:n]
    for _ in 1:m
        a, b = [parse(Int, x) for x in split(readline())]
        push!(edges[a], b)
        push!(edges[b], a)
    end

    frontier = [1]
    parent = zeros(Int, n)
    parent[1] = 1
    while length(frontier) > 0
        v = popfirst!(frontier)
        for e in edges[v]
            if parent[e] == 0
                parent[e] = v
                push!(frontier, e)
            end
        end
    end

    println("Yes")
    for p in parent[2:end]
        println(p)
    end
end

solve()

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
# schema and helper methods
mutable struct Schema
    schema::Vector{SchemaElement}
    name_lookup::Dict{AbstractString,SchemaElement}

    function Schema(elems::Vector{SchemaElement})
        name_lookup = Dict{AbstractString,SchemaElement}()
        name_stack = AbstractString[]
        nchildren_stack = Int[]

        for idx in 1:length(elems)
            sch = elems[idx]
            # TODO: may be better to have fully qualified names here to avoid chashes
            nested_name = (idx == 1) ? sch.name : join([name_stack; sch.name], '.')
            name_lookup[nested_name] = sch

            if !isempty(nchildren_stack)
                @debug("$(sch.name) is a child. remaining $(nchildren_stack[end]-1) children")
                if nchildren_stack[end] == 1
                    pop!(nchildren_stack)
                    pop!(name_stack)
                else
                    nchildren_stack[end] -= 1
                end
            end

            if (idx > 1) && (Thrift.isfilled(sch, :num_children) && (sch.num_children > 0))
                @debug("$(sch.name) has $(sch.num_children) children")
                push!(nchildren_stack, sch.num_children)
                push!(name_stack, sch.name)
            end
        end
        new(elems, name_lookup)
    end
end

leafname(schname::AbstractString) = istoplevel(schname) ? schname : leafname(split(schname, '.'))
leafname(schname::Vector) = schname[end]

parentname(schname::Vector) = istoplevel(schname) ? schname : schname[1:(end-1)]
#parentname(schname::AbstractString) = join(parentname(split(schname, '.')), '.')
parentname(schname::AbstractString) = schname

#istoplevel(schname::AbstractString) = !('.' in schname)
istoplevel(schname::AbstractString) = true
istoplevel(schname::Vector) = !(length(schname) > 1)

elem(sch::Schema, schname::AbstractString) = sch.name_lookup[schname]
elem(sch::Schema, schname::Vector) = elem(sch, join(schname, '.'))

function isrequired(sch::Schema, schname)
    schelem = elem(sch, schname)
    Thrift.isfilled(schelem, :repetition_type) && (schelem.repetition_type == FieldRepetitionType.REQUIRED)
end

function isrepeated(sch::Schema, schname)
    schelem = elem(sch, schname)
    Thrift.isfilled(schelem, :repetition_type) && (schelem.repetition_type == FieldRepetitionType.REPEATED)
end

function bit_or_byte_length(sch::Schema, schname)
    schelem = elem(sch, schname)
    Thrift.isfilled(schelem, :type_length) ? schelem.type_length : 0
end

function max_repetition_level(sch::Schema, schname)
    lev = isrepeated(sch, schname) ? 1 : 0
    istoplevel(schname) ? lev : (lev + max_repetition_level(sch, parentname(schname)))
end

function max_definition_level(sch::Schema, schname)
    lev = isrequired(sch, schname) ? 0 : 1
    istoplevel(schname) ? lev : (lev + max_definition_level(sch, parentname(schname)))
end

abstract type SchemaConverter end

# parquet schema to Protobuf schema converter
mutable struct ProtoConverter <: SchemaConverter
    to::IO
end

schema(conv::ProtoConverter, sch::Schema, schema_name::Symbol) = schema(conv, sch.schema, schema_name)
schema(conv::ProtoConverter, sch::Vector{SchemaElement}, schema_name::Symbol) = schema_to_proto_schema(conv.to, sch, schema_name)

function schema_to_proto_schema(io::IO, sch::Vector{SchemaElement}, schema_name::Symbol)
    nchildren = Int[]
    lev0 = IOBuffer()
    ios = IO[lev0]
    println(lev0, "message ", schema_name, " {")
    for schemaelem in sch
        _sch_to_proto(schemaelem, ios, nchildren)
    end
    lev0 = ios[end]
    println(lev0, "}")
    write(io, String(take!(lev0)))
    nothing
end

function _sch_to_proto(sch::SchemaElement, ios::Vector{IO}, nchildren::Vector{Int}=Int[])
    lchildren = length(nchildren)

    if isfilled(sch, :_type)
        jtypestr = PLAIN_PROTOBUF_TYPES[sch._type+1]
    else
        # this is a composite type
        jtypestr = sch.name * "Type"
    end
    # we are not looking at converted types yet

    repntype = ""
    if isfilled(sch, :repetition_type)
        (sch.repetition_type == FieldRepetitionType.REQUIRED) && (repntype = "required")
        (sch.repetition_type == FieldRepetitionType.OPTIONAL) && (repntype = "optional")
        (sch.repetition_type == FieldRepetitionType.REPEATED) && (repntype = "repeated")
    end

    fldid = isfilled(sch, :field_id) ? (" = " * string(sch.field_id)) : ""

    if lchildren > 0
        println(ios[end], "    ", repntype, " ", jtypestr, " ", sch.name, fldid, ";")
    end

    if isfilled(sch, :num_children)
        if lchildren > 0
            lvlio = IOBuffer()
            println(lvlio, "struct ", jtypestr, " {")
            push!(ios, lvlio)
        end
        push!(nchildren, sch.num_children)
    elseif lchildren > 0
        nchildren[lchildren] -= 1
        if nchildren[lchildren] == 0
            pop!(nchildren)
            lvlio = pop!(ios)
            if !isempty(ios)
                println(lvlio, "}")
                prevlvlio = pop!(ios)
                println(lvlio, "")
                print(lvlio, String(take!(prevlvlio)))
            end
            push!(ios, lvlio)
        end
    end
end

# parquet schema to Thrift schema converter
mutable struct ThriftConverter <: SchemaConverter
    to::IO
end

schema(conv::ThriftConverter, sch::Schema, schema_name::Symbol) = schema(conv, sch.schema, schema_name)
schema(conv::ThriftConverter, sch::Vector{SchemaElement}, schema_name::Symbol) = schema_to_thrift_schema(conv.to, sch, schema_name)

function schema_to_thrift_schema(io::IO, sch::Vector{SchemaElement}, schema_name::Symbol)
    nchildren = Int[]
    lev0 = IOBuffer()
    ios = IO[lev0]
    println(lev0, "struct ", schema_name, " {")
    for schemaelem in sch
        _sch_to_thrift(schemaelem, ios, nchildren)
    end
    lev0 = ios[end]
    println(lev0, "}")
    write(io, String(take!(lev0)))
    nothing
end

function _sch_to_thrift(sch::SchemaElement, ios::Vector{IO}, nchildren::Vector{Int}=Int[])
    lchildren = length(nchildren)

    if isfilled(sch, :_type)
        jtypestr = PLAIN_THRIFT_TYPES[sch._type+1]
    else
        # this is a composite type
        jtypestr = sch.name * "Type"
    end
    # we are not looking at converted types yet

    if isfilled(sch, :repetition_type) && (sch.repetition_type == FieldRepetitionType.REPEATED)  # array type
        jtypestr = "list<" * jtypestr * ">"
    end

    repntype = ""
    if isfilled(sch, :repetition_type)
        (sch.repetition_type == FieldRepetitionType.REQUIRED) && (repntype = "required")
        (sch.repetition_type == FieldRepetitionType.OPTIONAL) && (repntype = "optional")
    end

    fldid = isfilled(sch, :field_id) ? (string(sch.field_id) * ":") : ""

    if lchildren > 0
        println(ios[end], "    ", fldid, " ", repntype, " ", jtypestr, " ", sch.name, (nchildren[lchildren] == 0) ? "." : "")
    end

    if isfilled(sch, :num_children)
        if lchildren > 0
            lvlio = IOBuffer()
            println(lvlio, "struct ", jtypestr, " {")
            push!(ios, lvlio)
        end
        push!(nchildren, sch.num_children)
    elseif lchildren > 0
        nchildren[lchildren] -= 1
        if nchildren[lchildren] == 0
            pop!(nchildren)
            lvlio = pop!(ios)
            if !isempty(ios)
                println(lvlio, "}")
                prevlvlio = pop!(ios)
                println(lvlio, "")
                print(lvlio, String(take!(prevlvlio)))
            end
            push!(ios, lvlio)
        end
    end
end

# parquet schema to Julia types
mutable struct JuliaConverter <: SchemaConverter
    to::Union{Module,IO}
end

schema(conv::JuliaConverter, sch::Union{Schema, Vector{SchemaElement}}, schema_name::Symbol) = schema_to_julia_types(conv.to, sch, schema_name)

function schema_to_julia_types(mod::Module, sch::Schema, schema_name::Symbol)
    io = IOBuffer()
    schema_to_julia_types(io, sch, schema_name)
    typestr = "begin\n" * String(take!(io)) * "\nend"
    parsedtypes = Meta.parse(typestr)
    Core.eval(mod, parsedtypes)
end

schema_to_julia_types(io::IO, sch::Schema, schema_name::Symbol) = schema_to_julia_types(io, sch.schema, schema_name)

function schema_to_julia_types(io::IO, sch::Vector{SchemaElement}, schema_name::Symbol)
    nchildren = Int[]
    lev0 = IOBuffer()
    ios = IO[lev0]
    println(lev0, "mutable struct ", schema_name)
    println(lev0, "    ", schema_name, "() = new()")
    for schemaelem in sch
        _sch_to_julia(schemaelem, ios, nchildren)
    end
    lev0 = ios[end]
    println(lev0, "end")
    write(io, String(take!(lev0)))
    nothing
end

function _sch_to_julia(sch::SchemaElement, ios::Vector{IO}, nchildren::Vector{Int}=Int[])
    lchildren = length(nchildren)

    if isfilled(sch, :_type)
        jtype = PLAIN_JTYPES[sch._type+1]
        jtypestr = string(jtype)
    else
        # this is a composite type
        jtypestr = "var\"" * replace(sch.name, "\""=>"\\\"") * "Type\""
    end
    # we are not looking at converted types yet

    if (isfilled(sch, :_type) && (sch._type == _Type.BYTE_ARRAY || sch._type == _Type.FIXED_LEN_BYTE_ARRAY)) ||
       (isfilled(sch, :repetition_type) && (sch.repetition_type == FieldRepetitionType.REPEATED))  # array type
        jtypestr = "Vector{" * jtypestr * "}"
    end

    if lchildren > 0
        println(ios[end], "    var\"", replace(sch.name, "\""=>"\\\""), "\"::", jtypestr)
    end

    if isfilled(sch, :num_children)
        if lchildren > 0
            lvlio = IOBuffer()
            println(lvlio, "mutable struct ", jtypestr)
            println(lvlio, "    ", jtypestr, "() = new()")
            push!(ios, lvlio)
        end
        push!(nchildren, sch.num_children)
    elseif lchildren > 0
        nchildren[lchildren] -= 1
        if nchildren[lchildren] == 0
            pop!(nchildren)
            lvlio = pop!(ios)
            if !isempty(ios)
                println(lvlio, "end")
                prevlvlio = pop!(ios)
                println(lvlio, "")
                print(lvlio, String(take!(prevlvlio)))
            end
            push!(ios, lvlio)
        end
    end
end


text:
# object like json schema types
# https://json-schema.org/draft/2020-12/json-schema-core.html
# https://json-schema.org/understanding-json-schema/index.html
# https://opis.io/json-schema
using StructTypes
using StructTypes: OrderedStruct, UnorderedStruct, Mutable, DictType
using StructTypes: ArrayType
using StructTypes: StringType, NumberType, BoolType, NullType
using NamedTupleTools
using Dates
using UUIDs
using Base.Docs

const array_types = [ArrayType()]
const object_types = [OrderedStruct(), UnorderedStruct(), Mutable(), DictType()]
const primative_types = [StringType(), NumberType(), BoolType(), NullType()]

function union_types(t)
    l = []
    if t.a isa Union
        push!(l, union_types(t.a)...)
    else
        push!(l, t.a)
    end

    if t.b isa Union
        push!(l, union_types(t.b)...)
    else
        push!(l, t.b)
    end
    l
end

# whats is the signature for compound types
# bellow doesn't functions as intented
# json_type(::Type{Vector})  = array

json_type(::Type{Bool}) = "boolean"
json_type(::Type{<:Integer}) = "integer"
json_type(::Type{<:Real}) = "number"
json_type(::Type{Nothing}) = "null"
json_type(::Type{Missing}) = "null"
json_type(::Type{String}) = "string"
json_type(::Type{<:Enum}) = "string"

array_type(::Type{<:Vector{T}}) where {T} = T
array_type(::Type{<:Array{T}}) where {T} = T

json_schema(::Type{Headers{T}}, d=Dict{Symbol,Any}()) where {T} = json_schema(T, d)
json_schema(::Type{Body{T}}, d=Dict{Symbol,Any}()) where {T} = json_schema(T, d)

# maybe it would be better to define a trait for this
function doc_str(t::Type{T}) where {T}
    # Docs.getdocs could be usefull here

    return nothing

    md = Docs.doc(T)
    s = repr(md)
    if startswith(s, "No documentation found.")
        return nothing
    else
        return s
    end
end

function json_schema(::Type{T}, d=Dict{Symbol,Any}()) where {T}

    if T isa Union
        return JSONSchema(
            oneOf=unique(json_schema.(union_types(T)))
        )
    end

    sT = StructTypes.StructType(T)

    # http://json-schema.org/understanding-json-schema/reference/object.html
    if sT in object_types
        d[:type] = "object"
        d[:description] = doc_str(T)
        if !(sT == DictType())
            properties = Dict{String,Any}()
            StructTypes.foreachfield(T) do i, field, field_type
                field_schema = json_schema(field_type)
                push!(properties, string(field) => field_schema)
            end
            d[:properties] = properties
        end
        # http://json-schema.org/understanding-json-schema/reference/array.html
    elseif sT in array_types
        d[:type] = "array"
        if T <: Tuple
            d[:prefixItems] = map(x -> json_schema(x), T.types)
        end
        if T <: Array
            aT = array_type(T)
            if aT isa Union
                error("union types in arrays are not supported")
            end
            d[:items] = json_schema(aT)
        end
    elseif sT in primative_types
        if T <: Enum
            d[:description] = doc_str(T)
            d[:enum] = Base.Enums.namemap(T) |> values |> collect .|> String
        else
            d[:type] = json_type(T)
        end
    else
        error("unable to determin json_schema for $T")
    end

    schema = JSONSchema(; d...)

    if isnothing(schema.type) &&
       isnothing(schema.anyOf) &&
       isnothing(schema.enum) &&
       isnothing(schema.oneOf) &&
       isnothing(schema.not)
        @warn "No data type in schema"
    end

    return schema
end

# @enum JSONSchemaFormat

# str enum regex object boolean null array integer number
JSONSchemaType = String

function json_schema_format(x::Type{T}) where {T}

    # renaming valid formats
    # "hostname"
    # "idn-hostname"
    # "ipv4"
    # "ipv6"
    # "uri-reference"
    # "iri"
    # "ir-reference"
    # "uri-template"
    # "json-pointer"

    return if x <: DateTime
        "date-time"
    elseif x <: Time
        "time"
    elseif x <: Date
        "date"
    elseif x <: TimePeriod
        "duration"
    elseif x <: AbstractString
        "email"
        "idn-email"
    elseif x <: UUID
        "uuid"
    elseif x <: Regex
        "regex"
    elseif x <: URI
        "uri"
    else
        error("unsupported type $x")
    end
end

Base.@kwdef struct JSONSchema
    # can it really be and array of types?
    type::Union{Nothing,JSONSchemaType} = nothing

    enum::Union{Nothing,Vector{String}} = nothing

    # schema composition
    allOf::Union{Nothing,Vector{JSONSchema}} = nothing
    anyOf::Union{Vector{JSONSchema},Nothing} = nothing
    oneOf::Union{Nothing,Vector{JSONSchema}} = nothing
    not::Union{Nothing,Vector{JSONSchema}} = nothing

    # string related fields #
    maxLength::Union{Int,Nothing} = nothing
    minLength::Union{Int,Nothing} = nothing
    pattern::Union{Regex,Nothing} = nothing
    format::Union{String,Nothing} = nothing

    # number releated fields #
    multipleOf::Union{Int,Nothing} = nothing
    minimum::Union{Int,Nothing} = nothing
    exclusiveMinimum::Union{Int,Nothing} = nothing
    maximum::Union{Int,Nothing} = nothing
    exclusiveMaximum::Union{Int,Nothing} = nothing

    # object related fields #
    properties::Union{Dict{String,Any},Nothing} = nothing
    parternProperties::Union{Dict{Regex,Any},Nothing} = nothing
    additionalProperties::Union{Union{Bool,Any},Nothing} = nothing
    required::Union{Array{String},Nothing} = nothing
    propertiesNames::Union{Pair{String,String},Nothing} = nothing
    minProperties::Union{Int,Nothing} = nothing
    maxProperties::Union{Int,Nothing} = nothing
    const_::Union{Any,Nothing} = nothing

    # array releated fields #
    items::Union{JSONSchema,Nothing,Bool} = nothing
    prefixItems::Union{Vector{JSONSchema},Nothing} = nothing
    additionalItem::Union{Bool,Nothing} = nothing
    contains::Union{Vector{JSONSchema},Nothing} = nothing
    unqiueItems::Union{Bool,Nothing} = nothing
    minItems::Union{Int,Nothing} = nothing
    maxItems::Union{Int,Nothing} = nothing

    # generic fields #
    title::Union{String,Nothing} = nothing
    description::Union{String,Nothing} = nothing
    examples::Union{Array{Any},Nothing} = nothing
    deprecated::Union{Bool,Nothing} = nothing
    readOnly::Union{Bool,Nothing} = nothing
    writeOnly::Union{Bool,Nothing} = nothing
end

StructTypes.StructType(::Type{JSONSchema}) = StructTypes.Struct()
StructTypes.omitempties(::Type{JSONSchema}) = true

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
module TestProject

using StaticArrays

function dot(x)
    v = SVector(x...)
    return v'v
end

end # module


text:
# # Deep Continuous Bag of Words (DeepCBOW) Model

include("bow.jl")

# This time after mapping the words to embeddings and adding up these embeddings we pass
# them through a couple linear + tanh layers.

function pred(words)
    h = sum(W_emb[:,words], dims=2)
    for (w,b) in zip(W_h, b_h)
        h = tanh.(w * h .+ b)
    end
    return W_sm * h .+ b_sm
end

# Again, we only need to change the prediction function and the parameters:

EMB_SIZE = 64
HID_SIZE = 64
HID_LAY = 2
W_emb = param(EMB_SIZE, nwords)
W_h = [ param(HID_SIZE, lay == 1 ? EMB_SIZE : HID_SIZE) for lay in 1:HID_LAY ]
b_h = [ param(HID_SIZE) for lay in 1:HID_LAY ]
W_sm = param(ntags, HID_SIZE)
b_sm = param(ntags)

# Bonus question: can this model represent a function that the CBOW model cannot represent?

train()


text:
module Exam

using DataFrames
using CSV
using YAML
using Plots
using DataStructures

export read_exam
export read_grading
export evaluate

 # for half points Float64 is exact!
const MyFixedPointNumber = Float64

include("bonus.jl")
include("plotting.jl")

"""
    get_grades(exam)

Get a list of all grades
"""
function get_grades(exam)
    grades = [elem["grade"] for elem in exam["grades"]]
end

"""
    read_exam(path)

Read the `.yml` file for the exam.
"""
function read_exam(path)
    YAML.load_file(path)
end

"""

Read the `.csv` file with the results of the exam.
"""
function read_grading(path)
    CSV.read(path, DataFrame)#, types=Dict(1 =>String, 2=> Int, 3=> Bool, 4 => Float16))
end

function calculate_total_points!(results, exam)
    ntasks = length(exam["tasks"])
    total_points = sum(eachcol(select(results, ncol(results) - ntasks + 1:ncol(results))))
    insertcols!(results, :total_points => total_points)
    return results
end

function calculate_grades!(results, exam)
    # sum up the task points
    total_points = results[:, :total_points]

    # calculate a grade according to the grading scheme
    grades = map(total_points) do p
        grade = MyFixedPointNumber(last(exam["grades"])["grade"])
        best_points = last(exam["grades"])["min_points"]

        for g in exam["grades"]
            if p ≥ g["min_points"] && g["min_points"] ≥ best_points
                grade = g["grade"]
                best_points = g["min_points"]
            end
        end
        grade
    end

    if exam["bonus"]["available"]
        if exam["bonus"]["system"] == "increment_grade"
            grades = increment_grade(grades, results[:, :bonus], exam)
        end
    end

    insertcols!(results, :grade => grades)
    return results
end

"""
    evaluate(results_path, exam_path; <kwargs>)

# Keywords and default values
* `histogram=true`: export grade histogram
* `anonymize=false`: remove names and only print id

"""
function evaluate(results_path, exam_path;
                    histogram=true,
                    anonymize=false,
                    dump_results_CSV=true,
                )

    exam = read_exam(exam_path)
    results = read_grading(results_path)

    for i in 1:length(exam["grades"])
        exam["grades"][i]["grade"] = MyFixedPointNumber(exam["grades"][i]["grade"])
    end

    # calculate points and insert into dataframe
    calculate_total_points!(results, exam)
    # insert grades into table
    calculate_grades!(results, exam)

    # plot histogram
    if histogram
        grade_histogram(results, exam)
    end

    if anonymize
        select!(results, Not(:name))
    end

    if dump_results_CSV
        CSV.write("results_debug.csv", results)
    end

    return results
end

end # module


text:
module Science_ABM

using Agents
using Distributions
using Statistics
using YAML

include("types.jl")
include("step_functions.jl")
include("create_model.jl")
include("data_collection.jl")

end # module


text:
using EVQueues, PGFPlots

lambda=120.0;
mu=1.0;
gamma=0.5;
#C=80.0;
C=60.0;

Tfinal=100.0;

edf = ev_edf(lambda,mu,gamma,Tfinal,C,snapshots=[Tfinal])
llf = ev_llf(lambda,mu,gamma,Tfinal,C,snapshots=[Tfinal])
llr = ev_llr(lambda,mu,gamma,Tfinal,C,snapshots=[Tfinal])

#compute_statistics!(sim)

fig = Axis([
                Plots.Linear([ev.requestedEnergy for ev in edf.EVs][6000:10:end],[ev.requestedEnergy-ev.departureWorkload for ev in edf.EVs][6000:10:end], style="solid,only marks=true,blue", legendentry="EDF"),
                Plots.Linear([ev.requestedEnergy for ev in llf.EVs][6000:10:end],[ev.requestedEnergy-ev.departureWorkload for ev in llf.EVs][6000:10:end], style="solid,only marks=true,red", legendentry="LLF"),
                Plots.Linear([ev.requestedEnergy for ev in llr.EVs][6000:10:end],[ev.requestedEnergy-ev.departureWorkload for ev in llr.EVs][6000:10:end], style="solid,only marks=true,green", legendentry="LLR"),
           ],
           legendPos="north west", xlabel="Requested Service (\$S\$)", ylabel="Attained Service (\$S_r\$)", xmin=0, xmax=3, ymin=0, ymax=3, width="0.7\\columnwidth", height="0.4\\columnwidth"
    );

save("/tmp/sigma_compare.tex",fig,include_preamble=false)


text:
using Base.Test, Logging
Logging.configure(level = Logging.WARNING) # set to DEBUG to see error tables

@time @testset "DeformableRegistration" begin

@testset "submodules" begin
    include("ImageProcessing.jl")
    include("Transformation.jl")
    include("Interpolation.jl")
    include("Distance.jl")
    include("Regularizer.jl")
end

# test full parametric and nonparametric image registration
@testset "example registration cases" begin
    include("ParametricRegistration.jl")
    include("NonparametricRegistration.jl")
end

end


text:
@testset "Second order" begin
    @testset "Bits $i" for i in 3:5
        x = 0:1 << i
        y_approx = fpcos.(Int16.(x), Val(i))
        y = cos.(x / 1 << i * π / 2) .* 1 << i
        @test all(abs.(y_approx .- y) .< 2.0)
    end
end


text:
# AtmosLESDefaultPerturbations
#
# Computes perturbations from the horizontal averages for various
# fields.

using ..Atmos
using ..Mesh.Topologies
using ..Mesh.Grids
using Thermodynamics

"""
    setup_atmos_default_perturbations(
        ::AtmosLESConfigType,
        interval::String,
        out_prefix::String;
        writer::AbstractWriter,
        interpol = nothing,
    )

Create the "AtmosLESDefaultPerturbations" `DiagnosticsGroup` which
contains the following diagnostic variables, all of which are
perturbations from the horizontal average of that variable.

- u_prime
- v_prime
- w_prime
- avg_rho_prime
- temp_prime
- pres_prime
- thd_prime
- et_prime
- ei_prime
- ht_prime
- hi_prime

When an `EquilMoist` moisture model is used, the following additional
diagnostic variables are also output:

- qt_prime
- ql_prime
- qv_prime
- thv_prime
- thl_prime

The perturbations are output on `x`, `y`, `z` dimensions of an
interpolated grid (`interpol` _must_ be specified) and a (unlimited)
`time` dimension at the specified `interval`.
"""
function setup_atmos_default_perturbations(
    ::AtmosLESConfigType,
    interval::String,
    out_prefix::String;
    writer = NetCDFWriter(),
    interpol = nothing,
)
    # TODO: remove this
    @assert !isnothing(interpol)

    return DiagnosticsGroup(
        "AtmosLESDefaultPerturbations",
        Diagnostics.atmos_les_default_perturbations_init,
        Diagnostics.atmos_les_default_perturbations_fini,
        Diagnostics.atmos_les_default_perturbations_collect,
        interval,
        out_prefix,
        writer,
        interpol,
    )
end

# Compute sums for density-averaged horizontal averages
#
# These are trimmed down versions of `atmos_les_default_simple_sums!()`:
# - operate on interpolated grid (no `MH` scaling)
# - skip `w_ht_sgs` and `w_qt_sgs`
function atmos_les_default_perturbations_sums!(
    atmos::AtmosModel,
    state,
    thermo,
    sums,
)
    sums.u += state.ρu[1]
    sums.v += state.ρu[2]
    sums.w += state.ρu[3]
    sums.avg_rho += state.ρ
    sums.rho += state.ρ * state.ρ
    sums.temp += thermo.temp * state.ρ
    sums.pres += thermo.pres * state.ρ
    sums.thd += thermo.θ_dry * state.ρ
    sums.et += state.energy.ρe
    sums.ei += thermo.e_int * state.ρ
    sums.ht += thermo.h_tot * state.ρ
    sums.hi += thermo.h_int * state.ρ

    atmos_les_default_perturbations_sums!(
        moisture_model(atmos),
        state,
        thermo,
        sums,
    )

    return nothing
end
function atmos_les_default_perturbations_sums!(
    ::AbstractMoistureModel,
    state,
    thermo,
    sums,
)
    return nothing
end
function atmos_les_default_perturbations_sums!(
    moist::EquilMoist,
    state,
    thermo,
    sums,
)
    sums.moisture.qt += state.moisture.ρq_tot
    sums.moisture.ql += thermo.moisture.q_liq * state.ρ
    sums.moisture.qv += thermo.moisture.q_vap * state.ρ
    sums.moisture.thv += thermo.moisture.θ_vir * state.ρ
    sums.moisture.thl += thermo.moisture.θ_liq_ice * state.ρ

    return nothing
end

# Perturbations from horizontal averages
function vars_atmos_les_default_perturbations(m::AtmosModel, FT)
    @vars begin
        u_prime::FT
        v_prime::FT
        w_prime::FT
        avg_rho_prime::FT             # ρ
        temp_prime::FT
        pres_prime::FT
        thd_prime::FT                 # θ_dry
        et_prime::FT                  # e_tot
        ei_prime::FT                  # e_int
        ht_prime::FT
        hi_prime::FT

        moisture::vars_atmos_les_default_perturbations(moisture_model(m), FT)
    end
end
vars_atmos_les_default_perturbations(::AbstractMoistureModel, FT) = @vars()
function vars_atmos_les_default_perturbations(m::EquilMoist, FT)
    @vars begin
        qt_prime::FT                  # q_tot
        ql_prime::FT                  # q_liq
        qv_prime::FT                  # q_vap
        thv_prime::FT                 # θ_vir
        thl_prime::FT                 # θ_liq
    end
end
num_atmos_les_default_perturbation_vars(m, FT) =
    varsize(vars_atmos_les_default_perturbations(m, FT))
atmos_les_default_perturbation_vars(m, array) =
    Vars{vars_atmos_les_default_perturbations(m, eltype(array))}(array)

# Compute the perturbations from horizontal averages
function atmos_les_default_perturbations!(
    atmos::AtmosModel,
    state,
    thermo,
    ha,
    vars,
)
    u = state.ρu[1] / state.ρ
    vars.u_prime = u - ha.u
    v = state.ρu[2] / state.ρ
    vars.v_prime = v - ha.v
    w = state.ρu[3] / state.ρ
    vars.w_prime = w - ha.w
    vars.avg_rho_prime = state.ρ - ha.avg_rho
    vars.temp_prime = thermo.temp - ha.temp
    vars.pres_prime = thermo.pres - ha.pres
    vars.thd_prime = thermo.θ_dry - ha.thd
    et = state.energy.ρe / state.ρ
    vars.et_prime = et - ha.et
    vars.ei_prime = thermo.e_int - ha.ei
    vars.ht_prime = thermo.h_tot - ha.ht
    vars.hi_prime = thermo.h_int - ha.hi

    atmos_les_default_perturbations!(
        moisture_model(atmos),
        atmos,
        state,
        thermo,
        ha,
        vars,
    )

    return nothing
end
function atmos_les_default_perturbations!(
    ::AbstractMoistureModel,
    ::AtmosModel,
    state,
    thermo,
    ha,
    vars,
)
    return nothing
end
function atmos_les_default_perturbations!(
    m::EquilMoist,
    atmos::AtmosModel,
    state,
    thermo,
    ha,
    vars,
)
    qt = state.moisture.ρq_tot / state.ρ
    vars.moisture.qt_prime = qt - ha.moisture.qt
    vars.moisture.ql_prime = thermo.moisture.q_liq - ha.moisture.ql
    vars.moisture.qv_prime = thermo.moisture.q_vap - ha.moisture.qv
    vars.moisture.thv_prime = thermo.moisture.θ_vir - ha.moisture.thv
    vars.moisture.thl_prime = thermo.moisture.θ_liq_ice - ha.moisture.thl

    return nothing
end

function atmos_les_default_perturbations_init(
    dgngrp::DiagnosticsGroup,
    currtime,
)
    FT = eltype(Settings.Q)
    atmos = Settings.dg.balance_law
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)

    if !(dgngrp.interpol isa InterpolationBrick)
        @warn """
            Diagnostics $(dgngrp.name): requires `InterpolationBrick`!
            """
        return nothing
    end

    if mpirank == 0
        # get dimensions for the interpolated grid
        dims = dimensions(dgngrp.interpol)

        # set up the variables we're going to be writing
        vars = OrderedDict()
        varnames = map(
            s -> startswith(s, "moisture.") ? s[10:end] : s,
            flattenednames(vars_atmos_les_default_perturbations(atmos, FT)),
        )
        for varname in varnames
            vars[varname] = (tuple(collect(keys(dims))...), FT, OrderedDict())
        end

        # create the output file
        dprefix = @sprintf("%s_%s", dgngrp.out_prefix, dgngrp.name)
        dfilename = joinpath(Settings.output_dir, dprefix)
        noov = Settings.no_overwrite
        init_data(dgngrp.writer, dfilename, noov, dims, vars)
    end

    return nothing
end

"""
    atmos_les_default_perturbations_collect(dgngrp, currtime)

Perform a global grid traversal to compute various diagnostics.
"""
function atmos_les_default_perturbations_collect(
    dgngrp::DiagnosticsGroup,
    currtime,
)
    interpol = dgngrp.interpol
    if !(interpol isa InterpolationBrick)
        @warn """
            Diagnostics $(dgngrp.name): requires `InterpolationBrick`!
            """
        return nothing
    end

    mpicomm = Settings.mpicomm
    dg = Settings.dg
    Q = Settings.Q
    mpirank = MPI.Comm_rank(mpicomm)
    atmos = dg.balance_law
    grid = dg.grid
    grid_info = basic_grid_info(dg)
    topl_info = basic_topology_info(grid.topology)
    Nqk = grid_info.Nqk
    Nqh = grid_info.Nqh
    npoints = prod(grid_info.Nq)
    nrealelem = topl_info.nrealelem
    nvertelem = topl_info.nvertelem
    nhorzelem = topl_info.nhorzrealelem

    # get needed arrays onto the CPU
    if array_device(Q) isa CPU
        ArrayType = Array
        state_data = Q.realdata
        aux_data = dg.state_auxiliary.realdata
    else
        ArrayType = CuArray
        state_data = Array(Q.realdata)
        aux_data = Array(dg.state_auxiliary.realdata)
    end
    FT = eltype(state_data)

    # Compute thermo variables
    thermo_array = Array{FT}(undef, npoints, num_thermo(atmos, FT), nrealelem)
    @traverse_dg_grid grid_info topl_info begin
        state = extract_state(dg, state_data, ijk, e, Prognostic())
        aux = extract_state(dg, aux_data, ijk, e, Auxiliary())

        thermo = thermo_vars(atmos, view(thermo_array, ijk, :, e))
        compute_thermo!(atmos, state, aux, thermo)
    end

    # Interpolate the state and thermo variables.
    interpol = dgngrp.interpol
    istate =
        ArrayType{FT}(undef, interpol.Npl, number_states(atmos, Prognostic()))
    interpolate_local!(interpol, Q.realdata, istate)

    ithermo = ArrayType{FT}(undef, interpol.Npl, num_thermo(atmos, FT))
    interpolate_local!(interpol, ArrayType(thermo_array), ithermo)

    # FIXME: accumulating to rank 0 is not scalable
    all_state_data = accumulate_interpolated_data(mpicomm, interpol, istate)
    all_thermo_data = accumulate_interpolated_data(mpicomm, interpol, ithermo)

    if mpirank == 0
        # get dimensions for the interpolated grid
        dims = dimensions(interpol)

        nx = length(dims["x"][1])
        ny = length(dims["y"][1])
        nz = length(dims["z"][1])

        # collect horizontal sums
        simple_sums = [
            zeros(FT, num_atmos_les_default_simple_vars(atmos, FT))
            for _ in 1:nz
        ]
        @traverse_interpolated_grid nx ny nz begin
            statei = Vars{vars_state(atmos, Prognostic(), FT)}(view(
                all_state_data,
                lo,
                la,
                le,
                :,
            ))
            thermoi = thermo_vars(atmos, view(all_thermo_data, lo, la, le, :))
            simple = atmos_les_default_simple_vars(atmos, simple_sums[le])
            atmos_les_default_perturbations_sums!(
                atmos,
                statei,
                thermoi,
                simple,
            )
        end

        # compute horizontal averages
        simple_avgs = [
            zeros(FT, num_atmos_les_default_simple_vars(atmos, FT))
            for _ in 1:nz
        ]
        for le in 1:nz
            simple_avgs[le] .= simple_sums[le] ./ (nx * ny)
        end

        # complete density averaging
        simple_varnames = map(
            s -> startswith(s, "moisture.") ? s[10:end] : s,
            flattenednames(vars_atmos_les_default_simple(atmos, FT)),
        )
        for vari in 1:length(simple_varnames)
            for le in 1:nz
                ha = atmos_les_default_simple_vars(atmos, simple_avgs[le])
                avg_rho = ha.avg_rho
                if simple_varnames[vari] != "avg_rho"
                    simple_avgs[le][vari] /= avg_rho
                end
            end
        end

        # now compute the perturbations from the horizontal averages
        perturbations_array = Array{FT}(
            undef,
            nx,
            ny,
            nz,
            num_atmos_les_default_perturbation_vars(atmos, FT),
        )
        @traverse_interpolated_grid nx ny nz begin
            statei = Vars{vars_state(atmos, Prognostic(), FT)}(view(
                all_state_data,
                lo,
                la,
                le,
                :,
            ))
            thermoi = thermo_vars(atmos, view(all_thermo_data, lo, la, le, :))
            ha = atmos_les_default_simple_vars(atmos, simple_avgs[le])
            perturbations = atmos_les_default_perturbation_vars(
                atmos,
                view(perturbations_array, lo, la, le, :),
            )
            atmos_les_default_perturbations!(
                atmos,
                statei,
                thermoi,
                ha,
                perturbations,
            )
        end

        # prepare and write out the perturbations
        varvals = OrderedDict()
        varnames = map(
            s -> startswith(s, "moisture.") ? s[10:end] : s,
            flattenednames(vars_atmos_les_default_perturbations(atmos, FT)),
        )
        for (vari, varname) in enumerate(varnames)
            varvals[varname] = perturbations_array[:, :, :, vari]
        end

        # write output
        append_data(dgngrp.writer, varvals, currtime)
    end

    MPI.Barrier(mpicomm)
    return nothing
end # function collect

function atmos_les_default_perturbations_fini(
    dgngrp::DiagnosticsGroup,
    currtime,
) end


text:
@sumProductRule(:node_type     => Multiplication,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Nothing, Message{Gaussian}, Message{PointMass}),
                :name          => SPMultiplicationOutNGP)

@sumProductRule(:node_type     => Multiplication,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Nothing, Message{PointMass}, Message{Gaussian}),
                :name          => SPMultiplicationOutNPG)

@sumProductRule(:node_type     => Multiplication,
                :outbound_type => Message{PointMass},
                :inbound_types => (Nothing, Message{PointMass}, Message{PointMass}),
                :name          => SPMultiplicationOutNPP)

@sumProductRule(:node_type     => Multiplication,
                :outbound_type => Message{GaussianWeightedMeanPrecision},
                :inbound_types => (Message{Gaussian}, Nothing, Message{PointMass}),
                :name          => SPMultiplicationIn1GNP)

@sumProductRule(:node_type     => Multiplication,
                :outbound_type => Message{PointMass},
                :inbound_types => (Message{PointMass}, Nothing, Message{PointMass}),
                :name          => SPMultiplicationIn1PNP)

@sumProductRule(:node_type     => Multiplication,
                :outbound_type => Message{GaussianWeightedMeanPrecision},
                :inbound_types => (Message{Gaussian}, Message{PointMass}, Nothing),
                :name          => SPMultiplicationAGPN)

@sumProductRule(:node_type     => Multiplication,
                :outbound_type => Message{PointMass},
                :inbound_types => (Message{PointMass}, Message{PointMass}, Nothing),
                :name          => SPMultiplicationAPPN)

@sumProductRule(:node_type     => Multiplication,
                :outbound_type => Message{Gamma},
                :inbound_types => (Nothing, Message{Union{Gamma, Wishart}}, Message{PointMass}),
                :name          => SPMultiplicationOutNΓP)

@sumProductRule(:node_type     => Multiplication,
                :outbound_type => Message{Gamma},
                :inbound_types => (Nothing, Message{PointMass}, Message{Union{Gamma, Wishart}}),
                :name          => SPMultiplicationOutNPΓ)

@sumProductRule(:node_type     => Multiplication,
                :outbound_type => Message{Gamma},
                :inbound_types => (Message{Union{Gamma, Wishart}}, Nothing, Message{PointMass}),
                :name          => SPMultiplicationIn1ΓNP)

@sumProductRule(:node_type     => Multiplication,
                :outbound_type => Message{Gamma},
                :inbound_types => (Message{Union{Gamma, Wishart}}, Message{PointMass}, Nothing),
                :name          => SPMultiplicationAΓPN)


text:
"""
The data assimilation system (Optimal Interpolation)
Load:
    x_a_init.txt
    y_o.txt
Save:
    x_b_oi_nmc.txt
    x_a_oi_nmc.txt
"""

include("setting.jl");
using DelimitedFiles;
using LinearAlgebra;
using OrdinaryDiffEq;

# load initial condition
x_a_init = readdlm("data/x_a_init.txt");

# load observations
y_o_save = readdlm("data/y_o.txt");
y_o_time = y_o_save[:, 1];
y_o_data = y_o_save[:, 2:end];

# merge into sync
y_o_save = fill(NaN, (nT + 1, N))
for j = 1:N
    ke = 0
    for i = 1:nT+1
        # find nearest time
        ks = ke + 1
        Ta = (i - 1) * dT
        Ts = Ta - dT / 2
        Te = Ta + dT / 2
        ke = ks
        for k = ks:length(y_o_time)
            if y_o_time[k] > Te
                break
            end
            ke = k
        end
        min_dt = dT
        for k = ks:ke
            dt = abs(y_o_time[k] - Ta)
            if !isnan(y_o_data[k, j]) && dt < min_dt
                y_o_save[i, j] = y_o_data[k, j]
                min_dt = dt
            end
        end
    end
end

# initial x_b: no values at the initial time (assign NaN)
x_b_save = Matrix{Float64}(undef, nT + 1, N);
x_b_save[1, :] .= NaN;

# initial x_a: from x_a_init
x_a_save = similar(x_b_save);
x_a_save[1, :] = x_a_init';

# initial background error variance
σ = 0.2
R = Diagonal(repeat([σ^2], N))
B = R

α = 2.5e-1
stable = false
tol = 1e-6
step = 1

# observation operator
function observation_operators(y_o::Matrix{Float64})::Vector{Matrix{Float64}}
    len = size(y_o, 1)
    hs = Vector{Matrix{Float64}}()
    for i = 1:len
        row = map(x -> isnan(x) ? 0.0 : 1.0, y_o[i, :])
        push!(hs, Diagonal(row))
    end
    hs
end
H = observation_operators(y_o_save)

while !stable
    @show step
    # OI
    for tt = 1:nT
        tts = tt - 1
        tti = tt + 1
        Ts = tts * dT
        Ta = tt * dT
        # println("Cycle = ", tt, ", Ts = ", round(Ts; digits = 3), ", Ta = ", round(Ta; digits = 3))

        #--------------
        # forecast step
        #--------------
        prob = ODEProblem(lorenz96, x_a_save[tt, :], (Ts, Ta))
        sol = solve(prob, Tsit5())
        x_b_save[tti, :] = sol[end]

        #--------------
        # analysis step
        #--------------
        # background
        x_b = sol[end]

        # observation
        y_o = replace(y -> isnan(y) ? 0.0 : y, y_o_save[tti, :])

        # innovation
        d = y_o - H[tti] * x_b

        # analysis scheme
        K = B * H[tti]' / (H[tti] * B * H[tti]' + R)
        x_a = x_b + K * d

        x_a_save[tti, :] = x_a
    end # end tt loop

    # NMC
    x_a = x_a_save
    cycles = size(x_a, 1)

    ti_start = 4 + 1
    t1 = 4.0 * dT
    t2 = 8.0 * dT
    Pb = ones(Float64, N, N)

    for d = 0:20
        coeffs = Vector()

        for i in ti_start:cycles
            # integrate time period t1, end in t
            prob = ODEProblem(lorenz96, x_a[i, :], (0.0, t1))
            sol = solve(prob, Tsit5())
            xf_t1 = sol[end]

            # integrate time period t2, end in t
            prob = ODEProblem(lorenz96, x_a[i-4, :], (0.0, t2))
            sol = solve(prob, Tsit5())
            xf_t2 = sol[end]

            δf = xf_t2 - xf_t1
            pb = δf * δf'
            for j in 1:N
                pos = (j + d + 40 - 1) % 40 + 1
                neg = (j - d + 40 - 1) % 40 + 1
                push!(coeffs, pb[j, pos])
                push!(coeffs, pb[j, neg])
            end
        end

        filter!(x -> !isnan(x), coeffs)
        E = sum(coeffs) / length(coeffs)
        for i in 1:N
            pos = (i + d + 40 - 1) % 40 + 1
            neg = (i - d + 40 - 1) % 40 + 1
            global Pb[i, pos] = α * E
            global Pb[i, neg] = α * E
        end
    end

    Δ = sqrt(sum((B - Pb) .^ 2))
    global stable = Δ < tol
    global B = Pb
    global step += 1
    @show Δ
    if step == 30
        println("not converged")
        break
    end
end

open("data/x_a_oi_nmc.txt", "w") do xa
    writedlm(xa, x_a_save)
end

open("data/x_b_oi_nmc.txt", "w") do xb
    writedlm(xb, x_b_save)
end

open("data/b_nmc.txt", "w") do nmc
    writedlm(nmc, B)
end

max_stars_repo_name: Typas/Data-Assimilation-Project

max_issues_repo_name: Typas/Data-Assimilation-Project

max_forks_repo_name: Typas/Data-Assimilation-Project


text:
using AbstractPlotting
using MakieGallery
using Test

const MINIMAL = get(ENV, "ABSTRACTPLOTTING_MINIMAL", "false")

# does this machine have a OPENGL?
const OPENGL = haskey(ENV, "OPENGL") || haskey(ENV, "GITLAB_CI") || try success(pipeline(`glxinfo`, `grep version`)) catch; false end # if it's Gitlab, it must be JuliaGPU
OPENGL || (MINIMAL = "true")

@show OPENGL MINIMAL

OPENGL && begin @info "OpenGL detected"; using GLMakie end
OPENGL || @warn "No OpenGL detected!  Software tests only."

include("conversions.jl")
include("quaternions.jl")
include("projection_math.jl")

@testset "basic functionality" begin
    scene = scatter(rand(4))
    @test scene[Axis].ticks.title_gap[] == 3
    scene[Axis].ticks.title_gap = 4
    @test scene[Axis].ticks.title_gap[] == 4
    @test scene[Axis].tickmarks.length[] == (3, 3)
end

if MINIMAL == "false"

    # Load all entries in the database

    database = MakieGallery.load_database()

    ## Exceptions are made on basis of title and not index,
    # because index may change as MakieGallery changes.

    # All these require FFMpeg and need to save,
    # and are therefore ignored
    ffmpeg_exs  = [
        "Animation", "Lots of Heatmaps","Animated Scatter",
        "Chess Game", "Record Video",  "Animated surface and wireframe",
        "Moire", "pong", "pulsing marker", "Travelling wave",
        "Type recipe for molecule simulation", "Cobweb plot"
    ]

    # All these require GLMakie and so are ignored
    glmakie_exs = [
        "Textured Mesh", "Load Mesh", "Wireframe of a Mesh",
        "FEM mesh 3D", "Normals of a Cat", "Line GIF"
    ]

    # Requires GDAL (a GL package) so ignored
    gdal_exs = [
        "WorldClim visualization"
    ]

    # Requires GLMakie and ModernGL, so ignored
    moderngl_exs = [
        "Explicit frame rendering"
    ]

    # use Stepper plots, which save, so ignored
    save_exs = [
        "Axis theming", "Labels", "Color Legend",
        "Stepper demo"
    ]

    # hopefullly fixed by next tag of MakieGallery (already in AbstractPlotting
    color_exs = ["colormaps"]

    # curl fails with this for some reason, so it has been ignored.
    curl_exs = ["Earth & Ships"]
    stepper_exclude = ["Tutorial plot transformation", "Legend"]
    if OPENGL
        excluded_examples = Set()
    else
        # combine all exceptions into a single Set
        excluded_examples = excluded_examples = union(
            ffmpeg_exs, glmakie_exs, gdal_exs, moderngl_exs,
            save_exs, color_exs, curl_exs, stepper_exclude
        )
    end

else
    # Load only short tests
    database = MakieGallery.load_tests()
    excluded_examples = Set(["Comparing contours, image, surfaces and heatmaps"])
end

if !OPENGL # run software only tests...
    # Make sure we don't include Makie in the usings
    empty!(MakieGallery.plotting_backends)
    push!(MakieGallery.plotting_backends, "AbstractPlotting")
    @test AbstractPlotting.current_backend[] isa Missing
    @info "Starting minimal software tests"
    filter!(database) do example
        !(example.title in excluded_examples) &&
        !("record" in example.tags)
    end
    @testset "Gallery short tests" begin
        # iterate over database
        @testset "$(database[i].title) (#$i)" for i in 1:length(database)
            @debug("Running " * database[i].title * "\n(index $i)\n")
            # evaluate the entry
            try
                MakieGallery.eval_example(database[i])
                @test true
            catch e
                # THis is ok, since we try to record something, which we can't
                # without backend
                if e isa LoadError && e.error isa MethodError && (
                        e.error.f == AbstractPlotting.backend_display ||
                        e.error.f == AbstractPlotting.format2mime
                    )
                    @test true
                else
                    @test rethrow(e)
                end
            end
            @assert AbstractPlotting.current_backend[] isa Missing
       end
    end
else # full MakieGallery comparisons here
    using GLMakie
    @info("Running full tests - artifacts will be stored!")
    for exc in excluded_examples
        printstyled("Excluded ", color = :yellow, bold = true)
        println(exc)
    end
    filter!(entry -> !(entry.title in excluded_examples), database)
    tested_diff_path = joinpath(@__DIR__, "tested_different")
    test_record_path = joinpath(@__DIR__, "test_recordings")
    rm(tested_diff_path, force = true, recursive = true)
    mkpath(tested_diff_path)
    rm(test_record_path, force = true, recursive = true)
    mkpath(test_record_path)
    examples = MakieGallery.record_examples(test_record_path)
    MakieGallery.run_comparison(test_record_path, tested_diff_path)
end


text:
function preconditioner(
    controlODE,
    precondition;
    θ=initial_params(controlODE.controller),
    reg_coeff=1.0f0,
    saveat=(),
    progressbar=true,
    control_range_scaling=nothing,
    plot_progress=false,
    plot_final=true,
    optimizer=LBFGS(; linesearch=BackTracking()),
    maxiters=50,
    allow_f_increases=true,
    integrator=INTEGRATOR,
    sensealg=SENSEALG,
    adtype=GalacticOptim.AutoForwardDiff(),
    kwargs...,
)
    @info "Preconditioning..."
    @infiltrate
    fixed_dudt!(du, u, p, t) = controlODE.system!(du, u, p, t, precondition; input=:time)

    prog = Progress(
        length(controlODE.tsteps[2:end]);
        desc="Pretraining in subintervals...",
        dt=0.5,
        showspeed=true,
        enabled=progressbar,
    )
    for partial_time in controlODE.tsteps[2:end]
        partial_tspan = (controlODE.tspan[1], partial_time)
        fixed_prob = ODEProblem(fixed_dudt!, controlODE.u0, partial_tspan)
        fixed_sol = solve(fixed_prob, integrator; saveat, sensealg)

        function precondition_loss(params; plot=nothing)
            plot_arrays = Dict(:reference => [], :control => [])
            sum_squares = 0.0f0
            # mean_squares = 0.0f0

            # for (time, state) in zip(fixed_sol.t, fixed_sol.u)  # Zygote error
            for (i, state) in enumerate(eachcol(Array(fixed_sol)))
                reference = precondition(fixed_sol.t[i], nothing)  # precondition(time, params)
                control = controlODE.controller(state, params)
                diff_square = (control - reference) .^ 2
                if !isnothing(control_range_scaling)
                    diff_square ./ control_range_scaling
                end
                sum_squares += sum(diff_square)
                # mean_squares += mean(diff_square)
                Zygote.ignore() do
                    if !isnothing(plot)
                        push!(plot_arrays[:reference], reference)
                        push!(plot_arrays[:control], control)
                    end
                end
            end
            Zygote.ignore() do
                if !isnothing(plot)
                    @argcheck plot in [:unicode, :pyplot]
                    reference = reduce(hcat, plot_arrays[:reference])
                    control = reduce(hcat, plot_arrays[:control])

                    if plot == :unicode
                        for r in 1:size(reference, 1)
                            p = lineplot(reference[r, :]; name="fixed")
                            lineplot!(p, control[r, :]; name="neural")
                            display(p)
                        end
                    elseif plot == :pyplot
                        cmap = ColorMap("tab20")  # binary comparisons
                        nrows = size(reference, 1)
                        rows = 1:nrows
                        fig, axs = plt.subplots(
                            nrows,
                            1;
                            sharex="col",
                            squeeze=false,
                            constrained_layout=true,
                            #tight_layout=false,
                        )
                        for r in rows
                            ax = axs[r, 1]
                            ax.plot(
                                fixed_sol.t,
                                reference[r, :],
                                "o";
                                label="fixed",
                                color=cmap(2r - 2),
                            )
                            ax.plot(
                                fixed_sol.t,
                                control[r, :];
                                label="neural",
                                color=cmap(2r - 1),
                            )
                            ax.legend()
                        end
                        fig.supxlabel("time")
                        fig.suptitle("Preconditioning")
                        fig.show()
                    end
                end
            end
            regularization = reg_coeff * mean(abs2, params)
            return sum_squares + regularization
            # return sum_squares / mean_squares + regularization
        end

        optimization = sciml_train(
            precondition_loss, θ, optimizer, adtype; maxiters, allow_f_increases, kwargs...
        )
        θ = optimization.minimizer
        pvar = plot_progress ? :unicode : nothing
        ploss = precondition_loss(θ; plot=pvar)
        next!(prog; showvalues=[(:loss, ploss)])
        if partial_time == controlODE.tsteps[end] && plot_final
            precondition_loss(θ; plot=:pyplot)
        end
    end
    return θ
end

function constrained_training(
    controlODE,
    losses;
    αs,
    δs,
    starting_params=initial_params(controlODE.controller),
    show_progressbar=false,
    plots_callback=nothing,
    datadir=nothing,
    metadata=Dict(),  # metadata is added to this dict always
    optimizer=NADAM(),  # LBFGS(; linesearch=BackTracking()),
    maxiters=100,
    allow_f_increases=true,
    sensealg=SENSEALG,
    adtype=GalacticOptim.AutoForwardDiff(),
    kwargs...,
)
    @argcheck length(αs) == length(δs)

    @info "Training with constraints..."

    prog = Progress(
        length(αs); desc="Fiacco-McCormick barrier iterations", enabled=show_progressbar
    )

    θ = starting_params
    for (α, δ) in zip(αs, δs)

        # closure to comply with optimization interface
        loss(params) = reduce(+, losses(controlODE, params; α, δ, sensealg, kwargs...))

        # if !isnothing(plots_callback)
        #     plots_callback(controlODE, θ)
        # end

        # function print_callback(params, loss)
        #     println(loss)
        #     return false
        # end

        result = sciml_train(
            loss, θ, optimizer, adtype; maxiters, allow_f_increases, kwargs...
        )

        θ =
            result.minimizer +
            1.0f-1 * std(result.minimizer) * randn(length(result.minimizer))

        objective, state_penalty, control_penalty, regularization = losses(
            controlODE, θ; α, δ, kwargs...
        )

        local_metadata = Dict(
            :objective => objective,
            :state_penalty => state_penalty,
            :control_penalty => control_penalty,
            :regularization_cost => regularization,
            :parameters => θ,
            :num_params => length(initial_params(controlODE.controller)),
            :layers => controller_shape(controlODE.controller),
            :penalty_relaxations => δs,
            :penalty_coefficients => αs,
            :tspan => controlODE.tspan,
            :tsteps => controlODE.tsteps,
        )
        metadata = merge(metadata, local_metadata)
        store_simulation("delta_$(round(δ, digits=2))", controlODE, θ; metadata, datadir)
        next!(
            prog;
            showvalues=[
                (:α, α),
                (:δ, δ),
                (:objective, objective),
                (:state_penalty, state_penalty),
                (:control_penalty, control_penalty),
                (:regularization, regularization),
            ],
        )
    end
    return θ
end


text:
#=
Given an array of numbers, find the maximum sum of any contiguous subarray of the array.

For example, given the array [34, -50, 42, 14, -5, 86], the maximum sum would be 137, since we would take elements 42, 14, -5, and 86.

Given the array [-5, -1, -8, -9], the maximum sum would be 0, since we would not take any elements.

Do this in O(N) time.
=#

function maximum_contiguous_sum(X::Array{Int,1})
    max_sum = 0
    current_subarray = []
    current_sum = 0
    for index in eachindex(X)
        if current_sum + X[index] > 0
            push!(current_subarray, X[index])
            current_sum = current_sum + X[index]
            if current_sum > max_sum
                max_sum = current_sum
            end
        else
            current_sum = 0
            current_subarray = []
        end
    end
    return max_sum
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
"""
idwi(args...)  -->Plot

iwm is a Plot Recipe in which makes a spatial interpolation of data provided applying Inverse Weight Mapping.

General equation for spatial interpolation is:
``{ Z }^{ * }=\\sum _{ i=1 }^{ n }{ { \\lambda  }_{ i }{ z }_{ i } } ``

where ``\\lambda`` is the weight Parameter

Applying the Inverse‐Distance Weighted Interpolation

``{ Z }\\left( x \\right) =\\frac { \\sum _{ i=1 }^{ n }{ { w }_{ i }{ z }_{ i } }  }{ \\sum _{ i=1 }^{ n }{ { z }_{ i } }  } \\\\ \\\\ { w }_{ i }=\\frac { 1 }{ { d }_{ i }^{ n } } ``

<br> The next table show the list of variables allowed:

|PropertyName|Args|Default|Input|Description
|---|---|---|---|---|
|x|Mandatory|--|Array{Number,1}| X values|
|y|Mandatory|--|Array{Number,1}| Y values|
|z|Mandatory|--|Array{Number,1}| Z values|
|m|Optional|m=100|m=number| Number of grids in y axis|
|n|Optional|n=100|n=number| Number of grids in x axis|
|p|Optional|p=2|p=number| IDWI exponent. The greater more influenced by near points|
"""
@userplot idwi

@recipe function f(h::idwi; m=100, n=100, p=2, surf=false)

    ylabel --> "North [m]"
    xlabel --> "East [m]"
    seriescolor --> :Spectral

x, y, z = h.args

# Build the grid of prediction points

    ## X & Y range
  Xrange=range(minimum(x),stop=maximum(x),length=n)
  Yrange=range(minimum(y),stop=maximum(y),length=m)

    ## X & Y grids
  X=repeat(Xrange',length(Yrange),1)
  Y=repeat(Yrange,1,length(Xrange))

    # X & Y grids converted to vectors
    Xvec=vec(X)
    Yvec=vec(Y)

    #Number of Sample Points and Prediction Points
    Nsp=length(x)                           #Numbers of Sample Points
    Npp=length(Xrange)*length(Yrange)       #Number of Prediction Points

    #DISTANCE BETWEEN PREDICTION AND SAMPLE POINTS---------------------------
    #Matrix of distance   I x J    Sample Points x Prediction Points

    d=zeros(Nsp,Npp)
    for i=1:Nsp
        for j=1:Npp
            d[i,j]=sqrt((x[i]-Xvec[j]).^2+(y[i]-Yvec[j]).^2)
        end
    end

    #Estimate the matrix of weights

    W=map(dis->(1 ./(dis.^p)),d)

    #Calculate value in prediction point in a vector

    Zvec=zeros(Npp)
    for i=1:Npp
        Zvec[i]=W[:,i]'*z ./ sum(W[:,i])
    end

    Z=reshape(Zvec,m,n)

#Coutour series
if surf==true
    @series begin
        seriestype := :surface
        zflip := true
        Xrange, Yrange, Z
        end
else
@series begin
    seriestype := :contour
    clabels := true
    fill --> true

    Xrange, Yrange, Z
end
end
end