text:
#############################################################################
#############################################################################
#
# This file implements factorization
#
#############################################################################
#############################################################################

"""
Factors a polynomial over the field Z_p.

Returns a vector of tuples of (irreducible polynomials (mod p), multiplicity) such that their product of the list (mod p) is f. Irreducibles are fixed points on the function factor.
"""

function factor(f::PolynomialModP)::Vector{Tuple{Polynomial,Int}}
    #Cantor Zassenhaus factorization
    f.prime == 3 && return [(mod(f.terms, f.prime),1)]
    f_modp = PolynomialModP(mod(f.terms, f.prime), f.prime)
    degree(f_modp.terms) ≤ 1 && return [(f_modp.terms,1)]

    # make f primitive
    ff = PolynomialModP(prim_part(f_modp.terms)(f.prime), f.prime)
    # @show "after prim:", ff

    ff_der = PolynomialModP(derivative(prim_part(f_modp.terms)(f.prime)), f.prime)

     # make f square-free
    squares_poly = gcd(f, ff_der)
    ff = PolynomialModP(ff ÷ PolynomialModP(squares_poly, f.prime), f.prime)
    # @show "after square free:", ff

    # make f monic
    old_coeff = leading(ff.terms).coeff
    ff = ff ÷ old_coeff
    # @show "after monic:", ff

    dds = dd_factor(ff)

    ret_val = Tuple{Polynomial,Int}[]

    for (k,dd) in enumerate(dds)
        sp = dd_split(dd, k, f.prime)
        sp = map((p)->(p ÷ leading(p).coeff)(f.prime),sp) #makes the polynomials inside the list sp, monic
        for mp in sp
            push!(ret_val, (mp, multiplicity(f_modp,PolynomialModP(mp,f.prime))))
        end
    end

    #Append the leading coefficient as well
    push!(ret_val, (leading(f_modp.terms).coeff* one(Polynomial), 1) )
    # if f.prime == 3 && degree(mod(expand_factorization(ret_val),f.prime)) == 0 && mod(f.terms-expand_factorization(ret_val),f.prime) != 0
    #     ret_val = [(ret_val[1][1]*(Term(1,3)), ret_val[1][2])]
    # end
    return ret_val
end

"""
Expand a factorization.
"""
function expand_factorization(factorization::Vector{Tuple{Polynomial,Int}})::Polynomial
    length(factorization) == 1 && return first(factorization[1])^last(factorization[1])
    return *([first(tt)^last(tt) for tt in factorization]...)
end

"""
Compute the number of times g divides f
"""
function multiplicity(f::PolynomialModP, g::PolynomialModP)::Int
    degree(gcd(f, g)) == 0 && return 0
    return 1 + multiplicity(PolynomialModP((f ÷ g),g.prime), g)
end

"""
Distinct degree factorization.

Given a square free polynomial `f` returns a list, `g` such that `g[k]` is a product of irreducible polynomials of degree `k` for `k` in 1,...,degree(f) ÷ 2, such that the product of the list (mod `prime`) is equal to `f` (mod `prime`).
"""
function dd_factor(f::PolynomialModP)::Array{Polynomial}
    x = x_poly()
    w = deepcopy(x)
    g = Array{Polynomial}(undef,degree(f.terms)) #Array of polynomials indexed by degree

    #Looping over degrees
    for k in 1:degree(f.terms)
        w = rem(PolynomialModP(^(w,f.prime, f.prime), f.prime), f)
        g[k] = gcd(PolynomialModP(w - x, f.prime), f)
        f = PolynomialModP((f ÷ PolynomialModP(g[k], f.prime)), f.prime)
    end

    #edge case for final factor
    f.terms != one(Polynomial) && push!(g,f.terms)

    return g
end

"""
Distinct degree split.

Returns a list of irreducible polynomials of degree `d` so that the product of that list (mod prime) is the polynomial `f`.
"""

function dd_split(f::Polynomial, d::Int, prime::Int)::Vector{Polynomial}
    f = PolynomialModP(mod(f,prime), prime)
    degree(f.terms) == d && return [f.terms]
    degree(f.terms) == 0 && return []
    w = PolynomialModP(mod(rand(Polynomial, degree = d, monic = true), prime), prime)
    n_power = (prime^d-1) ÷ 2
    g = PolynomialModP(gcd(PolynomialModP(^(w,n_power) - one(Polynomial), prime), f), prime)
    ḡ = (f ÷ g) # g\bar + [TAB]
    return vcat(dd_split(g.terms, d, prime), dd_split(ḡ, d, prime) )
end


text:
if "./src" in LOAD_PATH
    # print("in src")
else
    dir = "./src"
    push!(LOAD_PATH,dir)
    # print("not in src")
end

module aerosolGen

using Distributions

# "The Mechanism of Breath Formation"
# DOI: 10.1089/jamp.2008.0720

# Particle size distribution for young subjects (b-1-0-f3-m-m)
# Mean = 1.015056269 μm
# Variance = 0.457766837 μm

# Determining parameters for log(X)
mean = 1.015056269e-6
var = 0.209550477e-6
std = sqrt(var)

# mean = 243900.6582
# var = 93909632285

# std = 0.457766837e-6

mu = log((mean^2)/(sqrt(var + mean^2)))
sigma = sqrt(log(1+((var)/(mean^2))))

# log_mean = exp(mean + var/2)
# log_var = (exp(var)-1)*exp(2*mean+var)

# log_mean = log(mean)-0.5*log((std/mean)^2 + 1)
# log_std = sqrt(log((std/mean)^2 + 1))

# println(log_mean)
# println("\n")
# println(log_std)
# println(mu)
# println("\n")
# println(sigma)

# println(log_mean + "\n" + log_var)

distrib0 = Distributions.LogNormal(mu, sigma)
# distrib0 = Distributions.LogNormal(mean, var)
range0 = 0.5e-10:0.5e-9:5e-8
# range0 = 0.5:0.1:4
# range0=-0.5e-6:0.2e-6:0.5e-6
discrete0 = pdf(distrib0, range0)
# discrete1 = cdf(distrib0, range0)
# discrete0=pdf(distrib0, )

# mean(distrib0)

# println(discrete0[1])
println(discrete0)
# println(discrete1)
# println(params(distrib0))
# println(meanlogx(distrib0))

end

ext: jl


text:
using DenoiseDNA
using Test

@testset "DenoiseDNA.jl" begin
    # Write your tests here.
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Blosc2_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Blosc2")
JLLWrappers.@generate_main_file("Blosc2", UUID("d43303dc-dd0e-56c6-b0a8-331f4c8c9bfb"))
end  # module Blosc2_jll


text:
## Reactive signal update performance: o, h, l, c tuple from TimeArray values-component
using TimeSeries
using Reactive

nma = 50
ohlc_BA = TimeSeries.readtimearray("../../data/OHLC_BA.csv")

ohlc_tup = (ohlc_BA.values[1,1], ohlc_BA.values[1,2], ohlc_BA.values[1,3], ohlc_BA.values[1,4])
s_ohlc = Reactive.Input(ohlc_tup)
s_close = lift(s -> s[4], s_ohlc, typ=Float64)

vsig = [s_close.value]
function vsigupd!(vsig::Vector{Float64}, nma::Int64, nmax::Int64)
    for i in 2:nmax
        ohlc_tup = (ohlc_BA.values[i,1], ohlc_BA.values[i,2], ohlc_BA.values[i,3], ohlc_BA.values[i,4])
        push!(s_ohlc, ohlc_tup)
        push!(vsig, s_close.value)
    end
end

println("OHLC updates o, h, l, c tuple:")
println("reactive signal updates timing")
expr_run = parse("@time vsigupd!(vsig, nma, length(ohlc_BA))")
for i = 1:5
    eval(expr_run)
end

### results
# elapsed time: 0.171515629 seconds (20491784 bytes allocated, 21.96% gc time)
### tuple is less efficient than vector, and requires more data restructuring anyway

Profile.clear()  # in case we have any previous profiling data
@profile vsigupd!(vsig, nma, length(ohlc_BA))
Profile.print()


text:
test_str = "16
10
15
5
1
11
7
19
6
12
4"

test_str2 = "28
33
18
42
31
14
46
20
48
47
24
23
49
45
19
38
39
11
1
32
25
35
8
17
7
9
4
2
34
10
3"

test_list = split(test_str, '\n')
test_list2 = split(test_str2, '\n')

list = readlines("day10_input.txt")

function full_connection(list)
    sorted = sort(map(x -> parse(Int, x), list))
    diff = sorted[2:end] - sorted[1:end-1]
    diff1 = count(==(1), diff)
    diff1 += sorted[1] == 1 ? 1 : 0
    diff3 = length(list) - diff1 + 1
    return diff1, diff3, diff1 * diff3
end

function num_connections(list)
    sorted = sort(map(x -> parse(Int, x), list))
    diff = vcat(1, sorted[2:end] - sorted[1:end-1], 3)
    comb = 1
    local_comb_map = [1, 1, 2, 4, 7]
    i = 1
    while i ≤ length(diff) - 1
        j = findnext(x -> x==3, diff, i)
        n = j - i
        if n in 0:4
            local_comb = local_comb_map[n+1]
        else
            local_comb =  sum([binomial(n-1,k) for k in 0:n-1]) - (n-3)(n-2)/2
        end
        i = j + 1
        comb *= local_comb
    end
    return comb
end

@show full_connection(test_list)

@show full_connection(test_list2)

@show full_connection(list)

@show num_connections(test_list)

@show num_connections(test_list2)

@show num_connections(list)

test_conn = (full_connection(test_list)[3] == 35) *
    (full_connection(test_list2)[3] == 220) *
    (full_connection(list)[3] == 2400)

test_num = (num_connections(test_list) == 8) *
    (num_connections(test_list2) == 19208) *
    (num_connections(list) == 338510590509056)

@show test_conn

@show test_num

nothing

#=
```julia
julia> @btime full_connection(list);
  7.722 μs (6 allocations: 4.95 KiB)

julia> @btime num_connections(list);
  11.199 μs (28 allocations: 6.81 KiB)
```
 =#

#=
# The number of combinations for varying n:

When there are n consecutive 1's, we can't remove the last adapter, otherwise
the step will be larger than three. But we can remove each of the first n-1
adapters. We can also change any two of the first n-1 adapters. Finally,
if n is greater than 4, we can also change any three to n-1 adapters provided
there aren't more than two consecutive removals. So, that is what we can do
    1 option of not removing anything (i.e. n-1 choose 0 options)
    n-1 choose 1 options of removing exactly one of the first n-1 adapters
    n-1 choose 2 options of removing exactly two of the first n-1 adapters
    n-1 choose 3 options of removing exactly three of the first n-1 adapters,
        discounting all the consecutive choices. The number of consecutive
        choices can be found by counting where they start. They can start at
        the first position, then at the second, up to the n-3 position, since
        we can change the last one. Hence, we can change n-3 adapters
    n-1 choose k options of removing exactly k of the first n-1 adapters,
        discounting all the n-k consecutive choices of k adapters.
        1+2+⋯+n-3, which is
        an arithmetic progression that sum up to (n-3)(1+n-3)/2 = (n-2)(n-3)/2
    ...
    n-1 choose n-1 options of removing exactly n-1 of the first n-1 adapters,
        discounting 1 consecutive choice of n-1 options, which is thus
        n = 0, 1, 2 and 3, if n-1 ≥ 3, and is zero if n ≥ 4.

Hence, the general form is

    local_comb = sum([binomial(n-1,k) for k in 0:n-1]) - sum([n-k for k in 3:n-1])
        = sum([binomial(n-1,k) for k in 0:n-1]) - (1 + 2 + 3 + ⋯ + n - 3)
        = sum([binomial(n-1,k) for k in 0:n-1]) - (n-3)(n-2)/2
        where the last term is present only if n ≥ 4

Besides, since lower values of n should be more frequent, it is faster if we
just calculate some of them before hand, e.g. for n = 0, 1, 2, 3, 4:

    n = 0 or 1 => local_comb = 1
    n = 2 => local_comb = 2
    n = 3 => local_comb = 4
    n = 4 => local_comb = 7
    n ≥ 5 => local_comb = sum([binomial(n-1,k) for k in 0:n-1]) - (n-3)(n-2)/2

That gives me

```julia
julia> @btime num_connections(list)
  15.115 μs (45 allocations: 8.52 KiB)
```
 =#


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7]);F = [5929*x1^2*x5^2+125*x1^2*x7^2-7084*x1*x2^2*x5+180*x1*x2^2*x7+9702*x1*x2*x3*x5+340*x1*x3*x5*x7-2002*x1*x3*x6*x7-232*x1*x4^2*x7+4158*x1*x4*x5^2-6622*x1*x4*x5*x7+248*x1*x4*x7^2-1166*x1*x5^2*x7+814*x1*x5*x6*x7+12760*x1*x5*x7^2+308*x1*x7^3+4141*x2^4-5796*x2^3*x3+3969*x2^2*x3^2+7650*x2^2*x3*x5-5220*x2^2*x4^2-2484*x2^2*x4*x5+9536*x2^2*x4*x7+306*x2^2*x7^2+3402*x2*x3*x4*x5-5418*x2*x3*x4*x7+9072*x2*x3*x7^2+7225*x3^2*x5^2+8281*x3^2*x6^2-9860*x3*x4^2*x5+10540*x3*x4*x5*x7+9646*x3*x5^2*x6-6734*x3*x5*x6^2-13832*x3*x5*x6*x7+13090*x3*x5*x7^2+3364*x4^4-7192*x4^3*x7+729*x4^2*x5^2-2322*x4^2*x5*x7-3239*x4^2*x7^2+3888*x4*x5*x7^2+3356*x4*x7^3+2809*x5^4-3922*x5^3*x6-8056*x5^3*x7+1369*x5^2*x6^2+5624*x5^2*x6*x7+5776*x5^2*x7^2+11113*x7^4-1980*x1^2*x7+16380*x1*x3*x6+9540*x1*x5^2-6660*x1*x5*x6-13680*x1*x5*x7+8100*x1^2, -14168*x1*x2*x5+360*x1*x2*x7+9702*x1*x3*x5+16564*x2^3-17388*x2^2*x3+7938*x2*x3^2+15300*x2*x3*x5-10440*x2*x4^2-4968*x2*x4*x5+19072*x2*x4*x7+612*x2*x7^2+3402*x3*x4*x5-5418*x3*x4*x7+9072*x3*x7^2, 9702*x1*x2*x5+340*x1*x5*x7-2002*x1*x6*x7-5796*x2^3+7938*x2^2*x3+7650*x2^2*x5+3402*x2*x4*x5-5418*x2*x4*x7+9072*x2*x7^2+14450*x3*x5^2+16562*x3*x6^2-9860*x4^2*x5+10540*x4*x5*x7+9646*x5^2*x6-6734*x5*x6^2-13832*x5*x6*x7+13090*x5*x7^2+16380*x1*x6, -464*x1*x4*x7+4158*x1*x5^2-6622*x1*x5*x7+248*x1*x7^2-10440*x2^2*x4-2484*x2^2*x5+9536*x2^2*x7+3402*x2*x3*x5-5418*x2*x3*x7-19720*x3*x4*x5+10540*x3*x5*x7+13456*x4^3-21576*x4^2*x7+1458*x4*x5^2-4644*x4*x5*x7-6478*x4*x7^2+3888*x5*x7^2+3356*x7^3, 11858*x1^2*x5-7084*x1*x2^2+9702*x1*x2*x3+340*x1*x3*x7+8316*x1*x4*x5-6622*x1*x4*x7-2332*x1*x5*x7+814*x1*x6*x7+12760*x1*x7^2+7650*x2^2*x3-2484*x2^2*x4+3402*x2*x3*x4+14450*x3^2*x5-9860*x3*x4^2+10540*x3*x4*x7+19292*x3*x5*x6-6734*x3*x6^2-13832*x3*x6*x7+13090*x3*x7^2+1458*x4^2*x5-2322*x4^2*x7+3888*x4*x7^2+11236*x5^3-11766*x5^2*x6-24168*x5^2*x7+2738*x5*x6^2+11248*x5*x6*x7+11552*x5*x7^2+19080*x1*x5-6660*x1*x6-13680*x1*x7, -2002*x1*x3*x7+814*x1*x5*x7+16562*x3^2*x6+9646*x3*x5^2-13468*x3*x5*x6-13832*x3*x5*x7-3922*x5^3+2738*x5^2*x6+5624*x5^2*x7+16380*x1*x3-6660*x1*x5, 250*x1^2*x7+180*x1*x2^2+340*x1*x3*x5-2002*x1*x3*x6-232*x1*x4^2-6622*x1*x4*x5+496*x1*x4*x7-1166*x1*x5^2+814*x1*x5*x6+25520*x1*x5*x7+924*x1*x7^2+9536*x2^2*x4+612*x2^2*x7-5418*x2*x3*x4+18144*x2*x3*x7+10540*x3*x4*x5-13832*x3*x5*x6+26180*x3*x5*x7-7192*x4^3-2322*x4^2*x5-6478*x4^2*x7+7776*x4*x5*x7+10068*x4*x7^2-8056*x5^3+5624*x5^2*x6+11552*x5^2*x7+44452*x7^3-1980*x1^2-13680*x1*x5];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
# add WebAssembly     # in package mode
# add IRTools

using WebAssembly: i32, f64, irfunc
using IRTools.All


text:
module IMPES

using LinearAlgebra, SparseArrays

include("types.jl")
include("flowFunctions.jl")
include("solvers.jl")
include("utilities.jl")
include("T_matrix.jl")

export RelPerms, SingleRelPerms, MixedRelPerms, Reservoir, Fluids, Grid, BoundaryConditions
export Simulation, FlowResults
export IMPES!
export CompTransport, IMPEC!
export recoveryFactor, injPoreVolumes, T_matrix
end


text:
using CLI
using Base.Test

@test parse(CLI.CLILongFlag, "--help") == CLI.CLILongFlag("help", "")
@test parse(CLI.CLILongFlag, "--help=true") == CLI.CLILongFlag("help", "true")
@test parse(CLI.CLILongFlag, "help=true") == nothing

@test parse(CLI.CLIShortFlag, "-h") == CLI.CLIShortFlag('h', "")
@test parse(CLI.CLIShortFlag, "-h=true") == CLI.CLIShortFlag('h', "true")
@test parse(CLI.CLIShortFlag, "-htrue") == CLI.CLIShortFlag('h', "true")
@test parse(CLI.CLIShortFlag, "help=true") == nothing


text:
export ExplicitOp

mutable struct ExplicitOp{T,F1<:FuncOrNothing,F2<:FuncOrNothing} <: AbstractLinearOperator{T,Function,F1,F2}
  nrow :: Int
  ncol :: Int
  symmetric :: Bool
  hermitian :: Bool
  prod :: Function
  tprod :: F1
  ctprod :: F2
end

"""
    ExplicitOp(shape::NTuple{D,Int64}, tr::Trajectory, correctionmap::Array{ComplexF64,D}
            ; echoImage::Bool=false, kargs...) where D

generates a `ExplicitOp` which explicitely evaluates the MRI Fourier signal encoding operator.

# Arguments:
* `shape::NTuple{D,Int64}`             - size of image to encode/reconstruct
* `tr::Trajectory`                     - Trajectory with the kspace nodes to sample
* `correctionmap::Array{ComplexF64,D}` - fieldmap for the correction of off-resonance effects
* `echoImage::Bool=false`              - if true sampling times will only be considered relative to the echo time
                                         this results in complex valued image even for real-valued input.
* `kargs`                              - additional keyword arguments
"""
function ExplicitOp(shape::NTuple{D,Int64}, tr::Trajectory, correctionmap::Array{ComplexF64,D}
                        ; echoImage::Bool=false, kargs...) where D

  nodes,times = kspaceNodes(tr), readoutTimes(tr)
  nrow = size(nodes,2)
  ncol = prod(shape)

  # if echo image is desired echo time is needed as an offset
  if echoImage
      echoOffset = echoTime(tr)
  else
      echoOffset = 0.0
  end

  return ExplicitOp{ComplexF64,Nothing,Function}(nrow, ncol, false, false
            , x->produ(x, shape, nodes, times, echoOffset, correctionmap)
            , nothing
            , y->ctprodu(y, shape, nodes, times, echoOffset, correctionmap))
end

function produ(x::Vector{T}, shape::NTuple{2,Int64},
                    nodes::Matrix{Float64}, times::Vector{Float64}, echoOffset::Float64,
                    correctionmap::Matrix{ComplexF64}) where T<:Union{Real,Complex}

   if isempty(correctionmap)
       disturbanceTerm = zeros(ComplexF64,shape...)
   else
       disturbanceTerm = correctionmap
   end

   x= reshape(x,shape)
   out = zeros(ComplexF64,size(nodes,2))
   for nx=1:shape[1]
       for ny=1:shape[2]
           for k=1:size(nodes,2)
               phi = (nodes[1,k]*(nx-shape[1]/2-1)+
                      nodes[2,k]*(ny-shape[2]/2-1))
               e = exp(-2*1im*pi*phi - (times[k]-echoOffset)*disturbanceTerm[nx,ny])
               out[k] += x[nx,ny] * e
           end
       end
   end
   return vec(out)
end

function produ(x::Vector{T}, shape::NTuple{3,Int64},
                    nodes::Matrix{Float64}, times::Vector{Float64}, echoOffset::Float64,
                    correctionmap::Array{ComplexF64,3}) where T<:Union{Real,Complex}

   if isempty(correctionmap)
       disturbanceTerm = zeros(ComplexF64,shape...)
   else
       disturbanceTerm = correctionmap
   end

   x = reshape(x,shape)
   out = zeros(ComplexF64,size(nodes,2))
   for nx=1:shape[1]
       for ny=1:shape[2]
           for nz=1:shape[3]
               for k=1:size(nodes,2)
                   phi = (nodes[1,k]*(nx-shape[1]/2-1)+nodes[2,k]*(ny-shape[2]/2-1)+nodes[3,k]*(nz-shape[3]/2-1))
                   e = exp(-2*1im*pi*phi - (times[k]-echoOffset)*disturbanceTerm[nx,ny,nz])
                   out[k] += x[nx,ny,nz] * e
               end
           end
       end
   end
   return vec(out)
end

function ctprodu(x::Vector{T}, shape::NTuple{2,Int64},
                        nodes::Matrix{Float64}, times::Vector{Float64}, echoOffset::Float64,
                        correctionmap::Matrix{ComplexF64}) where T<:Union{Real,Complex}

  if isempty(correctionmap)
      disturbanceTerm = zeros(ComplexF64,shape...)
  else
      disturbanceTerm = correctionmap
  end

  out = zeros(ComplexF64,shape)
  for nx=1:shape[1]
      for ny=1:shape[2]
          for k=1:size(nodes,2)
              phi = (nodes[1,k]*(nx-shape[1]/2-1)+nodes[2,k]*(ny-shape[2]/2-1))
              e = exp(-2*1im*pi*phi - (times[k]-echoOffset)*disturbanceTerm[nx,ny])
              out[nx,ny] += x[k] * conj(e)
          end
      end
  end
  return vec(out)
end

function ctprodu(x::Vector{T}, shape::NTuple{3,Int64},
                        nodes::Matrix{Float64}, times::Vector{Float64}, echoOffset::Float64,
                        correctionmap::Array{ComplexF64,3}) where T<:Union{Real,Complex}

  if isempty(correctionmap)
      disturbanceTerm = zeros(ComplexF64,shape...)
  else
      disturbanceTerm = correctionmap
  end

  out = zeros(ComplexF64,shape)
  for nx=1:shape[1]
      for ny=1:shape[2]
          for nz=1:shape[3]
              for k=1:size(nodes,2)
                 phi = (nodes[1,k]*(nx-shape[1]/2-1)+nodes[2,k]*(ny-shape[2]/2-1)+nodes[3,k]*(nz-shape[3]/2-1))
                 e = exp(-2*1im*pi*phi - (times[k]-echoOffset)*disturbanceTerm[nx,ny,nz])
                 out[nx,ny,nz] += x[k] * conj(e)
              end
          end
      end
  end
  return vec(out)
end

function adjoint(op::ExplicitOp{T}) where T
  return LinearOperator{T,Function,Nothing,Function}(op.ncol, op.nrow, op.symmetric, op.hermitian,
                        op.ctprod, nothing, op.prod)
end


text:
# Extend Distributions type hiearchy to account for non-euclidean supports
abstract type NonEuclidean <: Dist.ValueSupport end

const UnivariateFiniteSuper = Dist.Distribution{Dist.Univariate,NonEuclidean}

struct UnivariateFinite{L,U,T<:Real} <: UnivariateFiniteSuper
    decoder::CategoricalDecoder{L,U}
    prob_given_class::LittleDict{U,T}
end

"""$(MMI.UNIVARIATE_FINITE_DOCSTRING)"""
UnivariateFinite(a...) = MMI.UnivariateFinite(a...)

# Fallbacks
MMI.UnivariateFinite(::FI, d::AbstractDict) = throw(
    ArgumentError("The support of a `UnivariateFinite` can consist " *
                  "only of `CategoricalString` or `CategoricalValue` " *
                  "elements, and probabilities must be `<: Real`."))
MMI.UnivariateFinite(::FI, c::AbstractVector, p) = throw(
    ArgumentError("The `classes` must have type `AbstractVector{T}` where " *
                  "`T` is of type `CategoricalString` or `CategoricalValue` " *
                  "Perhaps you have `T=Any`?"))

# Univariate Finite from a dictionary of pairs CLASS => PROB
function MMI.UnivariateFinite(
                prob_given_class::AbstractDict{L,T}
                ) where {L<:CategoricalElement,T<:Real}
    # retrieve decoder and classes from element
    an_element     = first(keys(prob_given_class))
    parent_decoder = decoder(an_element)
    parent_classes = classes(an_element)
    given_classes  = keys(prob_given_class)
    given_probs    = values(prob_given_class)

    # check that the probabilities form a probability vector
    Dist.@check_args(UnivariateFinite, Dist.isprobvec(given_probs |> collect))

    # it's expected here that given_classes ⊆ parent_classes
    # then we form a dictionary mapping int(key) => prob
    pairs = ((int(c) => prob_given_class[c])
                for c in intersect(given_classes, parent_classes))

    return UnivariateFinite(parent_decoder, LittleDict(pairs...))
end

# Univariate Finite from a vector of classes and vector of probs
function MMI.UnivariateFinite(::FI, c::AbstractVector{C}, p::AbstractVector{P}
                              ) where C <: CategoricalElement where P <: Real
    # check that the vectors have appropriate length
    Dist.@check_args(UnivariateFinite, length(c) == length(p))
    # it's necessary to force the typing of the LittleDict otherwise it
    # may just convert to an abstractvector without keeping the 'categorical'.
    prob_given_class = LittleDict{C,P}(c[i] => p[i] for i in eachindex(c))
    return MMI.UnivariateFinite(prob_given_class)
end

# ------------------------------------------------------------------------
# utils for univariate finite objects

"""
    classes(d::UnivariateFinite)

A list of categorial elements in the common pool of classes used to
construct `d`.

    v = categorical(["yes", "maybe", "no", "yes"])
    d = UnivariateFinite(v[1:2], [0.3, 0.7])
    classes(d) # CategoricalArray{String,1,UInt32}["maybe", "no", "yes"]

"""
MMI.classes(d::UnivariateFinite) = classes(d.decoder.pool)


text:
using PowerSystems
cost = VariableCost([(1.0, 1.0), (2.0, 1.1), (3.0, 1.2)])
slopes = get_slopes(cost)
res = [1.0, 10.0, 10.0]
for (ix, v) in enumerate(slopes)
    @test isapprox(v, res[ix])
end

bps = get_breakpoint_upperbounds(cost)
res = [1.0, 0.1, 0.1]
for (ix, v) in enumerate(bps)
    @test isapprox(v, res[ix])
end

@testset "Test MarketBidCost with Polynomial Cost Timeseries with Service Forecast " begin
    initial_time = Dates.DateTime("2020-01-01")
    resolution = Dates.Hour(1)
    other_time = initial_time + resolution
    name = "test"
    horizon = 24
    service_data = Dict(initial_time => ones(horizon), other_time => ones(horizon))
    polynomial_cost = repeat([(999.0, 1.0)], 24)
    data_polynomial =
        SortedDict(initial_time => polynomial_cost, other_time => polynomial_cost)
    sys = create_rts_system()
    generators = collect(get_components(ThermalStandard, sys))
    generator = get_component(ThermalStandard, sys, get_name(generators[1]))
    market_bid = MarketBidCost(nothing)
    set_operation_cost!(generator, market_bid)
    forecast = IS.Deterministic("variable_cost", data_polynomial, resolution)
    set_variable_cost!(sys, generator, forecast)
    for s in generator.services
        forecast = IS.Deterministic(get_name(s), service_data, resolution)
        set_service_bid!(sys, generator, s, forecast)
    end

    cost_forecast = get_variable_cost(generator, market_bid, start_time = initial_time)
    @test first(TimeSeries.values(cost_forecast)).cost ==
          first(data_polynomial[initial_time])

    for s in generator.services
        service_cost = get_services_bid(generator, market_bid, s, start_time = initial_time)
        @test first(TimeSeries.values(service_cost)).cost ==
              first(service_data[initial_time])
    end
end

@testset "Test MarketBidCost with PWL Cost Timeseries" begin
    initial_time = Dates.DateTime("2020-01-01")
    resolution = Dates.Hour(1)
    other_time = initial_time + resolution
    name = "test"
    horizon = 24
    pwl_cost = repeat([repeat([(999.0, 1.0)], 5)], 24)
    data_pwl = SortedDict(initial_time => pwl_cost, other_time => pwl_cost)
    sys = create_rts_system()
    generators = collect(get_components(ThermalStandard, sys))
    generator = get_component(ThermalStandard, sys, get_name(generators[1]))
    market_bid = MarketBidCost(nothing)
    set_operation_cost!(generator, market_bid)
    forecast = IS.Deterministic("variable_cost", data_pwl, resolution)
    set_variable_cost!(sys, generator, forecast)

    cost_forecast = get_variable_cost(generator, market_bid, start_time = initial_time)
    @test first(TimeSeries.values(cost_forecast)).cost == first(data_pwl[initial_time])
end

@testset "Test ReserveDemandCurve with Cost Timeseries" begin
    initial_time = Dates.DateTime("2020-01-01")
    resolution = Dates.Hour(1)
    other_time = initial_time + resolution
    name = "test"
    horizon = 24
    polynomial_cost = repeat([(999.0, 1.0)], 24)
    data_polynomial =
        SortedDict(initial_time => polynomial_cost, other_time => polynomial_cost)
    pwl_cost = repeat([repeat([(999.0, 1.0)], 5)], 24)
    data_pwl = SortedDict(initial_time => pwl_cost, other_time => pwl_cost)
    for d in [data_polynomial, data_pwl]
        @testset "Add deterministic from $(typeof(d)) to ReserveDemandCurve variable cost" begin
            sys = System(100)
            reserve = ReserveDemandCurve{ReserveUp}(nothing)
            add_component!(sys, reserve)
            forecast = IS.Deterministic("variable_cost", d, resolution)
            set_variable_cost!(sys, reserve, forecast)
            cost_forecast = get_variable_cost(reserve, start_time = initial_time)
            @test first(TimeSeries.values(cost_forecast)).cost == first(d[initial_time])
        end
    end
end


text:
"""

    oparse(pt::Type{T<:Parser}, c::Union{Corpus,Sentence},
           ndeps:Integer, feats=nothing; usepmap=false)

The oracle parser executes the moves that will lead to the best (least
cost) parse.  It returns the final parse(s) and, if `feats` is
specified, the input features and correct moves for training.

## Arguments:
* pt::ParserType: ArcHybridR1, ArcEagerR1, ArcHybrid13, ArcEager13
* s::Sentence or c::Corpus: input sentence(s), a single parser is returned for s, a vector of parsers for c
* ndeps::Integer: number of dependency types
* feats::Fvec: (optional) specification of features, a (p,x,y) tuple returned if specified, only p if not
* usepmap::Bool: (optional) performs parallel processing

"""
function oparse{T<:Parser}(pt::Type{T}, s::Sentence, feats=nothing; o...)
    p = pt(s)
    op_work(p, s, feats)
end

function oparse{T<:Parser}(pt::Type{T}, c::Corpus, feats=nothing; usepmap::Bool=false)
    if usepmap
        op_pmap(pt, c, feats)
    else
        # op_main returns x,y in vectors for pcat
        p = op_main(pt, c, feats)
        if isa(p[1],Vector); (p[1],p[2][1],p[3][1]); else; p; end
    end
end

function op_pmap{T<:Parser}(pt::Type{T}, c::Corpus, feats)
    d = distribute(c)           # TODO: test what happens to common s.vocab
    p = pmap(procs(d)) do x
        op_main(pt, localpart(d), feats)
    end
    return pcat(p)
end

function op_main{T<:Parser}(pt::Type{T}, c::Corpus, feats)
    pa = map(pt, c)
    op_work(pa, c, feats)
end

# With no feats specified, just parse the sentences

function op_work{T<:Parser}(pa::Vector{T}, c::Corpus, feats::Void=nothing)
    for i=1:length(c)
        op_work(pa[i], c[i], feats)
    end
    return pa
end

function op_work(p::Parser, s::Sentence, feats::Void=nothing)
    c = Array(Cost, p.nmove)
    totalcost = 0
    while anyvalidmoves(p)
        movecosts(p, s.head, s.deprel, c)
        (bestcost,bestmove) = findmin(c)
        totalcost += bestcost
        move!(p, bestmove)
    end
    if totalcost != truecost(p,s); warn("Cost mismatch"); end
    return p
end

# With feats specified, parse the sentences and compute the features for training

function op_work{T<:Parser}(pa::Vector{T}, c::Corpus, feats::Fvec, x=nothing, y=nothing, nx=0)
    (x,y) = initoparse(pa,c,feats,x,y,nx)
    for i=1:length(c)
        op_work(pa[i], c[i], feats, x, y, nx)
        nx += nmoves(pa[i], c[i])
    end
    return (pa,typeof(x)[x],typeof(y)[y])  # must return tuple of three vectors for pcat
end

function op_work(p::Parser, s::Sentence, feats::Fvec, x=nothing, y=nothing, nx=0)
    (x,y) = initoparse(p,s,feats,x,y,nx)
    c = Array(Cost, p.nmove)
    totalcost = 0; nx0 = nx
    while anyvalidmoves(p)
        movecosts(p, s.head, s.deprel, c)
        (bestcost,bestmove) = findmin(c)
        totalcost += bestcost
        features(p, s, feats, x, (nx+=1))
        y[:, nx] = zero(eltype(y))
        y[bestmove, nx] = one(eltype(y))
        move!(p, bestmove)
    end
    if nx0 + nmoves(p,s) != nx; error(); end
    if totalcost != truecost(p,s); error("Cost mismatch"); end
    return (p,x,y)
end

function initoparse(p, s, f, x, y, n)
    if f == nothing; error(); end
    (xrows,xcols) = xsize(p,s,f)
    (yrows,ycols) = ysize(p,s)
    if xcols != ycols; error("Bad size"); end
    ytype = wtype(s)
    xtype = (if isa(f,DFvec); wtype(s); elseif isa(f,SFvec); SFtype; else; error(); end)
    mincols = xcols + n
    if y==nothing; y = zeros(ytype, yrows, mincols); end
    if x == nothing; x = Array(xtype, xrows, mincols); end
    if !(isa(y,AbstractMatrix{ytype}) && (size(y,1)==yrows) && (size(y,2)>=mincols)); error("Bad y"); end
    if !(isa(x,AbstractMatrix{xtype}) && (size(x,1)==xrows) && (size(x,2)>=mincols)); error("Bad x"); end
    return (x,y)
end

function pcat(p)
    if isa(p[1], Vector)
        q = similar(p[1], 0)
        for pi in p
            append!(q, pi)
        end
        return q
    else
        pp = similar(p[1][1],0)
        xx = similar(p[1][2],0)
        yy = similar(p[1][3],0)
        for i=1:length(p)
            append!(pp, p[i][1])
            append!(xx, p[i][2])
            append!(yy, p[i][3])
        end
        return (pp, xx, yy)
    end
end

### DEAD CODE

# function oparse{T<:Parser}(pt::Type{T}, c::Corpus, ndeps::Integer, feats::Fvec)
#     pa = map(s->pt(wcnt(s), ndeps), c)
#     oparse(pa, c, ndeps, feats)
# end

# function oparse{T<:Parser}(pt::Type{T}, c::Corpus, ndeps::Integer)
#     pa = map(s->pt(wcnt(s), ndeps), c)
#     oparse(pa, c, ndeps)
# end

# function oparse{T<:Parser}(pt::Type{T}, c::Corpus, ndeps::Integer, ncpu::Integer)
#     @date Main.resetworkers(ncpu)
#     sa = distribute(c)                                  # distributed sentence array
#     pa = map(s->pt(wcnt(s), ndeps), sa)                 # distributed parser array
#     @sync for p in procs(sa)
#         @spawnat p oparse(localpart(pa), localpart(sa), ndeps)
#     end
#     pa = convert(Vector{pt}, pa)
#     @date Main.rmworkers()
#     return pa
# end

# function oparse{T<:Parser}(pt::Type{T}, c::Corpus, ndeps::Integer, ncpu::Integer, feats::DFvec)
#     # Sparse fvec not supported yet, we need SharedSparseArray
#     @date Main.resetworkers(ncpu)
#     sa = distribute(c)                                  # distributed sentence array
#     pa = map(s->pt(wcnt(s), ndeps), sa)                 # distributed parser array
#     xtype = wtype(c[1])
#     x = SharedArray(xtype, xsize(pa[1],c,feats))    # shared x array
#     y = SharedArray(xtype, ysize(pa[1],c))          # shared y array
#     fill!(y, zero(xtype))
#     nx = zeros(Int, length(c))                      # 1+nx[i] is the starting x column for i'th sentence
#     p1 = pt(1,ndeps)
#     for i=1:length(c)-1
#         nx[i+1] = nx[i] + nmoves(p1, c[i])
#     end
#     @sync for p in procs(sa)
#         @spawnat p oparse(localpart(pa), localpart(sa), ndeps, feats, x, y, nx[localindexes(sa)[1][1]])
#     end
#     pa = convert(Vector{pt}, pa)
#     @date Main.rmworkers()
#     return (pa, sdata(x), sdata(y))
# end

# function oparse{T<:Parser}(pt::Type{T}, s::Sentence, ndeps::Integer)
#     p = pt(wcnt(s), ndeps)
#     oparse(p, s, ndeps)
# end

# function oparse{T<:Parser}(pt::Type{T}, s::Sentence, ndeps::Integer, feats::Fvec)
#     p = pt(wcnt(s), ndeps)
#     oparse(p, s, ndeps, feats)
# end

# function oparse{T<:Parser}(pt::Type{T}, c::Corpus, ndeps::Integer; usepmap::Bool=false)
#     if usepmap
#         d = distribute(c)
#         p = pmap(procs(d)) do x
#             oparse(pt, localpart(d), ndeps)
#         end
#         pcat(p)
#     else
#         pa = map(s->pt(wcnt(s), ndeps), c)
#         oparse(pa, c, ndeps)
#     end
# end

# function oparse{T<:Parser}(pt::Type{T}, c::Corpus, ndeps::Integer, feats::DFvec; usepmap::Bool=false)
#     if usepmap
#         d = distribute(c)
#         p = pmap(procs(d)) do x
#             oparse(pt, localpart(d), ndeps, feats)
#         end
#         pcat(p)
#     else
#         pa = map(s->pt(wcnt(s), ndeps), c)
#         oparse(pa, c, ndeps, feats)
#     end
# end

# function op_work(p::Parser, s::Sentence, ndeps::Integer, feats::SFvec, x=nothing, y=nothing, nx=0)
#     (x,y) = initoparse(p,s,ndeps,feats,x,y,nx)
#     c = Array(Position, p.nmove)
#     totalcost = 0; nx0 = nx
#     while anyvalidmoves(p)
#         movecosts(p, s.head, s.deprel, c)
#         (bestcost,bestmove) = findmin(c)
#         totalcost += bestcost
#         # We let features directly write in x.rowval[x.colptr[nx]:x.colptr[nx+1]-1]
#         features(p, s, feats, x, (nx+=1))
#         # We do not need this check any more, x.m is set to typemax
#         # features() sorts rowval so the max can be found at the end
#         # maxrow = x.rowval[x.colptr[nx+1]-1]
#         # If max rowval exceeds matrix height, we update the height
#         # maxrow > x.m && (x.m = maxrow)
#         y[:, nx] = zero(eltype(y))
#         y[bestmove, nx] = one(eltype(y))
#         move!(p, bestmove)
#     end
#     if nx0 + nmoves(p,s) != nx; error(); end
#     if totalcost != truecost(p,s); error("Cost mismatch"); end
#     return (p,x,y)
# end


text:
using Test
using FinDiff2018
using Plots

#=
Piecewise Linear Continuous Functions

	We approximate functions by a piecewise linear continuous approximation,
	with a regular grid spacing
=#

fsinpi = samplePLC(Float64, Float64, 4, sinpi)
#fsinpi2 = samplePLC(Float64, Float64, 4, sinpi)
#@show fsinpi

#@test (fsinpi + fsinpi2) == (fsinpi2 + fsinpi)

(evaluate(fsinpi, 0.3333), sinpi(0.3333))
xs = collect(range(0, stop=1, length=1001))
#@show xs

plot(sinpi.(xs))
plot!([evaluate(fsinpi, x) for x in xs])

rm("C:/Julia/FinDiff2018/test/Plot1_sinPLC.png", force=true)
savefig("Plot1_sinPLC.png")

fsinpi2 = refine(fsinpi)
plot([evaluate(fsinpi, x) for x in xs])
plot!([evaluate(fsinpi2, x)+0.1 for x in xs])

rm("C:/Julia/FinDiff2018/test/Plot1_errorPLC.png", force=true)
savefig("Plot1_errorPLC.png")

############################################################
##			Advection Equation
############################################################
s0 = initialGaussian(18)
plot([evaluate(s0.u, x) for x in xs])

rm("C:/Julia/FinDiff2018/test/Plot2_AdStateInitGuassian.png", force=true)
savefig("Plot2_AdStateInitGuassian.png")

r0 = rhsAdvection(s0)
plot([evaluate(r0.u, x) for x in xs])

rm("C:/Julia/FinDiff2018/test/Plot2_RHSadvection.png", force=true)
savefig("Plot2_RHSadvection.png")

############################################################
##			Time evolution
############################################################
s1 = euler(rhsAdvection, 0.1, s0)

plot([evaluate(s0.u, x) for x in xs])
plot!([evaluate(s1.u, x) for x in xs])

rm("C:/Julia/FinDiff2018/test/Plot3_EulerRHSadvection.png", force=true)
savefig("Plot3_EulerRHSadvection.png")

sol1 = solveAdvection(1.0, 100, 0.5);
sol2 = solveAdvection(1.0, 200, 0.5);

function plotSolution(sol, every)
    plot!([evaluate(sol1.states[1].u, x) for x in xs])
    for state in sol1.states[every:every:end]
        plot!([evaluate(state.u, x) for x in xs])
    end
    plot!()
end

plot()
plotSolution(sol1, 10)

rm("C:/Julia/FinDiff2018/test/Plot3_solveAdvection.png", force=true)
savefig("Plot3_solveAdvection.png")

############################################################
##			Discretization Error, Convergence
############################################################

sol1 = solveAdvection(1.0, 100, 0.5);
sol2 = solveAdvection(1.0, 200, 0.5);
sol3 = solveAdvection(1.0, 400, 0.5);

@show (sol1.states[2].time, sol2.states[3].time, sol3.states[5].time)

ue1 = refine(sol1.states[2].u) + (-1.0) * sol2.states[3].u
ue2 = refine(sol2.states[3].u) + (-1.0) * sol3.states[5].u;

plot([evaluate(ue1, x) for x in xs])
plot!([evaluate(ue2, x) for x in xs])

rm("C:/Julia/FinDiff2018/test/Plot4_descreteError_Convergence.png", force=true)
savefig("Plot4_descreteError_Convergence.png")


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

@doc raw"""NodeMetricsList is a list of NodeMetrics.

    IoK8sApiMetricsV1beta1NodeMetricsList(;
        metadata=nothing,
        items=nothing,
    )

    - metadata::IoK8sApimachineryPkgApisMetaV1ListMeta : Standard list metadata. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - items::Vector{IoK8sApiMetricsV1beta1NodeMetrics} : List of node metrics.
"""
mutable struct IoK8sApiMetricsV1beta1NodeMetricsList <: SwaggerModel
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ListMeta } # spec name: metadata
    items::Any # spec type: Union{ Nothing, Vector{IoK8sApiMetricsV1beta1NodeMetrics} } # spec name: items

    function IoK8sApiMetricsV1beta1NodeMetricsList(;metadata=nothing, items=nothing)
        o = new()
        validate_property(IoK8sApiMetricsV1beta1NodeMetricsList, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiMetricsV1beta1NodeMetricsList, Symbol("items"), items)
        setfield!(o, Symbol("items"), items)
        o
    end
end # type IoK8sApiMetricsV1beta1NodeMetricsList

const _property_map_IoK8sApiMetricsV1beta1NodeMetricsList = Dict{Symbol,Symbol}(Symbol("metadata")=>Symbol("metadata"), Symbol("items")=>Symbol("items"))
const _property_types_IoK8sApiMetricsV1beta1NodeMetricsList = Dict{Symbol,String}(Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ListMeta", Symbol("items")=>"Vector{IoK8sApiMetricsV1beta1NodeMetrics}")
Base.propertynames(::Type{ IoK8sApiMetricsV1beta1NodeMetricsList }) = collect(keys(_property_map_IoK8sApiMetricsV1beta1NodeMetricsList))
Swagger.property_type(::Type{ IoK8sApiMetricsV1beta1NodeMetricsList }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiMetricsV1beta1NodeMetricsList[name]))}
Swagger.field_name(::Type{ IoK8sApiMetricsV1beta1NodeMetricsList }, property_name::Symbol) =  _property_map_IoK8sApiMetricsV1beta1NodeMetricsList[property_name]

function check_required(o::IoK8sApiMetricsV1beta1NodeMetricsList)
    true
end

function validate_property(::Type{ IoK8sApiMetricsV1beta1NodeMetricsList }, name::Symbol, val)
end


text:
# This is an attempt to include the required for n-body simulations fields into structures
# Indeed, tools of ConcreteAbstractions.jl seem to be more suitable for the fields inheriting
DiffEqBase.@def position_velocity_mass begin
    r::SVector{3,cType}
    v::SVector{3,cType}
    m::mType
end

abstract type Body
end

struct MassBody{cType <: Real,mType <: Real} <: Body
    @position_velocity_mass
end

struct ChargedParticle{cType <: Real,mType <: Real,qType <: Real} <: Body
    @position_velocity_mass
    q::qType
end

struct MagneticParticle{cType <: Real,mType <: Real,mmType <: Real} <: Body
    @position_velocity_mass
    mm::SVector{3,mmType}
end

struct WaterMolecule <: Body
    O::MassBody
    H1::MassBody
    H2::MassBody
end

function generate_bodies_in_cell_nodes(n::Int, m::Real, v_dev::Real, L::Real; rng=MersenneTwister(n))
    velocities = v_dev * randn(rng, Float64, (3, n))
    bodies = MassBody[]

    count = 1
    dL = L / (ceil(n^(1 / 3)))
    for x = dL/2:dL:L, y = dL/2:dL:L, z = dL/2:dL:L
        if count > n
            break
        end
        r = SVector(x, y, z)
        v = SVector{3}(velocities[:,count])
        body = MassBody(r, v, m)
        push!(bodies, body)
        count += 1
    end
    return bodies
end


text:
function relevant_guide_key(attr)
    if RecipesBase.is_explicit(attr, :letter)
        letter = attr[:letter]
        Symbol(letter, :guide)
    else
        :guide
    end
end

function accumulate_unit_info(vui, attr, u, serno; ax = nothing)
    if length(vui) == 0
        # Store pre-existing guide strings in "series 0"
        gx = get(attr, :xguide, nothing)
        gy = get(attr, :yguide, nothing)
        gz = get(attr, :zguid, nothing)
        addinfo = SeriesUnitInfo(0, gx, nothing, nothing, :x, nothing)
        push!(vui, addinfo)
        addinfo = SeriesUnitInfo(0, gy, nothing, nothing, :y, nothing)
        push!(vui, addinfo)
        addinfo = SeriesUnitInfo(0, gz, nothing, nothing, :z, nothing)
        push!(vui, addinfo)
    end

    letter = get(attr, :letter, nothing)
    if isnothing(letter)
        letter = ax
    end
    relevant_axis_guide = relevant_guide_key(attr)
    label = get(attr, :label, nothing)
    sertyp = get(attr, :seriestype, nothing)
    if sertyp == :postfixseries
        sertyp = nothing
    end
    if typeof(label) <: AbstractArray
        label = label[serno]
    end
    if typeof(sertyp) <: AbstractArray && sertyp[serno] != :postfixseries
        sertyp = sertyp[serno]
    end
    addinfo = SeriesUnitInfo(serno, u, label, relevant_axis_guide, letter, sertyp)
    push!(vui, addinfo)
    vui
end

units_from_string(st) = Set{String}([ma.captures[1] for ma in eachmatch(r"\[(.*?)\]", st)])

function axis_guide_with_units(letter, attrdic, unitinfo::Vector{SeriesUnitInfo})
    existing_guide_string = string(get(attrdic, Symbol(letter, :guide), ""))
    # Collect units for the axis, one instance only
    existing_units = units_from_string(existing_guide_string )
    new_u = Set{String}()
    for uinf in unitinfo[4:end]
        uinf.letter == letter && push!(new_u, string(uinf.unit_foo))
    end
    add_units = setdiff(new_u,  existing_units)
    add_string = if length(add_units) == 0
        ""
    else
        "[" * join(collect(add_units), "], [") * "]"
    end
    if length(add_string) == 0
        existing_guide_string
    else
        if length(existing_units) > 0
            existing_guide_string * ", " * add_string
        else
            existing_guide_string * " " * add_string
        end
    end
end

function axis_unit_series_no(axisletter, unitinfo::Vector{SeriesUnitInfo})
    for uinf in reverse(unitinfo)
        uinf.letter == axisletter && return uinf.unit_foo, uinf.serno
    end
    error("Could not determine axis unit and series number")
end

function sertype(index, unitinfo::Vector{SeriesUnitInfo})
    for uinf in reverse(unitinfo)
        uinf.serno == index && !isnothing(uinf.sertyp) && return uinf.sertyp
    end
    default(:seriestype)
end

_modified_ribbon(ribbon, u) = ribbon / u
_modified_ribbon(ribbon::Function, u) = x-> ribbon(x) / (1u)

function modified_ribbon(ribbon::T, serno, letter, u) where T<:Union{Missing,<:Quantity}
    if letter == :y
        ribbon / u
    else
        ribbon
    end
end

function modified_ribbon(ribbon::T, serno, letter, u) where T<:Tuple
    if letter ==:y
        map(ribbon) do r
            modified_ribbon(r, serno, letter, u)
        end
    else
        ribbon
    end
end

function modified_ribbon(ribbon::T, serno, letter, u) where T<:AbstractRange
    if letter ==:y
        broadcast(se-> _modified_ribbon(se, u), ribbon)
    else
        ribbon
    end
end

function modified_ribbon(ribbon::T, serno, letter, u) where T<:Matrix
    _debug_recipes[1] && printstyled(color=:yellow, "\n    MechGluePlots modified_ribbon($ribbon<:Matrix, $serno, $letter, $u)\n")
    if letter ==:y
        if size(ribbon, 2) > 1
            # Branch to modify
            serib = ribbon[:, serno]
            # Modified branch
            modrib_se = broadcast(se-> _modified_ribbon(se, u), serib)
            # We are only modifying the branch of this tree relevant to the series, but we need to
            # return all of the tree. Most likely, an Array{Any}
            hcat(ribbon[:, 1:serno-1], modrib_se, ribbon[:, serno + 1:end])
        else
            ribbon / u
        end
    else
        ribbon
    end
end

# Never hit during tests
function modified_ribbon(ribbon::T, serno, letter, u) where T<:Matrix{Function}
    _debug_recipes[1] && printstyled(color=:green, "\n    MechGluePlots modified_ribbon($ribbon<:Matrix{Function}, $serno, $letter, $u)\n")
    if size(ribbon, 2) > 1
        # Branch to modify
        serib = ribbon[:, serno]
        # Modified branch
        modrib_se = broadcast(se-> modified_ribbon(se, serno, letter, u), serib)
        # We are only modifying the branch of this tree relevant to the series, but we need to
        # return all of the tree.
        hcat(ribbon[:, 1:serno-1], modrib_se, ribbon[:, serno + 1:end])
    else
        error("This form of the ribbon argument is not implemented with quantity plots")
    end
end
# Never hit during tests
function modified_ribbon(ribbon::T, serno, letter, u) where T<:Vector{Function}
    _debug_recipes[1] && printstyled(color=:green, "\n    MechGluePlots modified_ribbon($ribbon<:Vector{Function}, $serno, $letter, $u)\n")
    if size(ribbon, 1) > 1
        # Branch to modify
        se = ribbon[serno]
        # Modified branch
        modrib_se = modified_ribbon(se, serno, letter, u)
        # We are only modifying the branch of this tree relevant to the series, but we need to
        # return all of the tree.
        vcat(ribbon[1:serno-1], modrib_se, ribbon[serno + 1:end])
    else
        error("This form of the ribbon argument is not implemented with quantity plots")
    end
end

function modified_ribbon(ribbon::Function , serno, letter, u)
    _debug_recipes[1] && printstyled(color=:yellow, "\n    MechGluePlots modified_ribbon($ribbon::Function, $serno, $letter, $u)\n")
    if letter == :y
        # u is the defined output unit
        x-> ribbon(x) / (1u)
    else
        # u is the defined input unit
        x -> ribbon(x∙u)
    end
end

function modified_ribbon(ribbon::T , serno, letter, u) where T<:Vector
    if letter == :x
        ribbon
    else
        _modified_ribbon(ribbon, u)
    end
end
# Never hit during tests
function modified_ribbon(ribbon::T , serno, letter, u) where T
    @info T
    @info ribbon
    error("This form of the ribbon argument is not implemented with quantity plots")
end

"""
    extract_series(argobj, serno)
Extract the part of 'argobj' structure that belongs to 'series' no,
without evaluating - argobj may be a function, a tuple, etc.

This mimicks the way the 'argobj' argument is interpreted by 'Plots.jl',
but this leaves evaluating for other parts in the pipeline.
"""
function extract_series(argobj::Tuple{S, T}, serno) where {S, T}
    (extract_series(argobj[1], serno),
    extract_series(argobj[2], serno))
end
extract_series(argobj, serno) = argobj[serno]
extract_series(argobj::Matrix, serno) = argobj[:, serno]
extract_series(argobj::Nothing, serno) = argobj
extract_series(argobj::Number, serno) = argobj
extract_series(argobj::Function, serno) = argobj
extract_series(argobj::AbstractRange, serno) = argobj

##########################
#  Print debug functions
#  To use, call:
#  RecipesBase.debug(true)
##########################

function print_prettyln(v::T) where T <: Vector{SeriesUnitInfo}
    headings =  string.(fieldnames(SeriesUnitInfo))
    lens = map(s -> length(s), headings )
    heading = "               " * join(map(headings, lens) do h, l
        cepad(h, l + 2)
    end)
    printstyled(stderr, heading, "\n", color=:yellow)
    for si in v
        fields = fieldnames(SeriesUnitInfo)
        vs = [map(x-> string(getfield(si, x)), fields)..., "]"]
        vspad = map((h, l) -> cepad(h,l), vs, lens)
        line = " SeriesUnitInfo(" * join(vspad, ',') * ")"
        println(stderr, line)
    end
end

function print_prettyln(attr::Dict{Symbol, Any})
    for (ke, va) in attr
        if ke == :unitinfo
            printstyled(lpad("unitinfo", 18), "    =>\n", color = :green)
            print_prettyln(va)
        else
            printstyled(lpad(ke, 18), color=:green)
            print(IOContext(stderr, :limit=>true), "    => ", va, "\n")
        end
    end
    println(" ")
end
function print_prettyln(attr)
    println("{", typeof(attr), "}")
    println(attr)
end
function print_prettyln(attr::Plots.RecipesPipeline.DefaultsDict)
    println("    Explicit")
    print_prettyln(attr.explicit)
    if length(keys(attr.defaults)) < 10
        println("    Defaults")
        print_prettyln(attr.defaults)
    else
        println("    Defaults (not shown, length $(length(keys(attr.defaults))))")
    end
end
"""
    cepad(s, n; p = ' ')
centerpad string s within length n
"""
cepad(s, n; p = ' ') = rpad(lpad(s,div(n + textwidth(string(s)), 2), p), n, p)


text:
# Utils
function module_dir()
    return joinpath(@__DIR__, "..")
end

# Activate package
using Pkg
Pkg.activate(module_dir())

using LinearAlgebra
using LyceumBase
using LyceumMuJoCo
using LyceumMuJoCoViz
using MuJoCo
using Plots
using JLD2
using Random

include("methods.jl")

################################################################################
# Demo
################################################################################

timestep=0.01
jm, jd, mjsim = mj_model("humanoid.xml", timestep=0.01)
tsim = 1.0
N = Int(floor(tsim/timestep))

jd.qpos[3] += 1.0
Random.seed!(100)
jd.qvel .= 1.0 * (rand(27) .- 0.5)

traj = zeros(83, 0)
E = zeros(N)
Plin = zeros(N)
Pang = zeros(N)
tcompute = 0.0
for i = 1:N
    jd.ctrl .= 1e-5* 0.5
    tcompute += @elapsed mj_step(jm, jd);
    traj = hcat(traj, deepcopy(getstate(mjsim)))
    E[i] = energy(jm, jd)
	Plin[i] = norm(linear_momentum(jm, jd))
	# Pang[i] = angular_momentum(jm, jd)
end
tsim / t
plot(E)
plot(Plin)
plot(Pang)

LyceumMuJoCoViz.visualize(mjsim, trajectories=traj)

jd.qpos
jd.qvel


text:
param(:g,'n',:k,3)'

Pr"ω_i,j,3"

σx(:i)*σy(:i)
σz(3)*σy(2)

a(:m)*adag(:n)*adag(3)*a(3)*adag(3)*a(:n)*adag(:m)

σy(1)'
σy(:m)'
(adag(1)*a(:n))' == adag(:n)*a(1)
a(1) * (σy(1) * a(1))'

let
    tmp = ∑(:i,adag(:i)*a(:i))
    adag(:n)*tmp
    a(:n)*tmp
    σz(:n)*tmp
    tmp*2
    2*tmp
    tmp*param(:g,'n',:n)
    tmp*a(:n)
    tmp*adag(:n)
    tmp*(adag(:n)*adag(:n))
    tmp*(adag(:n)*a(:n))
end
∑(:i,σx(:i))*σy(:j)
∑(:i,σx(:i))*σy(:a)
∑(:i,σx(:i)*σy(:a))

comm(σx(5),σy(3))
comm(σx(5),σx(5))
comm(σx(1),σz(1))
2*comm(σx(5),σy(5))
comm(σy(5),comm(σx(5),σy(5)))

comm((2//5)*Pc"h"*σx(5),3Pc"g"*σy(5))

σx(1)*σy(1)
σx(3)*σy(2)*σy(2)

σx(5)*σy(5)*σz(5)
σx(5)*σy(5)*σx(5)

Pc"g_2"*Pc"g_1"*(Pc"g_3")'*Pr"b"*Pc"a"*Pr"d_3"

normal_form(a(2)*adag(2)*adag(3)*σy(1)*adag(2)*adag(4)*σz(1)*a(2)*adag(3))

adag(2)*a(2)*adag(2)*σy(1) + 2*param(:g,'c',1)

-1//2 + σp(1) * σm(1)

σp(2)
σm(2)
σp(1)*σm(1)

comm(σp(1),σp(1))
comm(σm(:n),σp(:n))
comm(a(2),adag(:m))

comm(a(1),adag(1)*a(1))
comm(a(1),adag(1)*a(2)*a(1))
comm(a(1),a(2)*adag(1)*a(1)*adag(2))

expval_as_corrs(a(2))
expval_as_corrs(3*a(2))

expval_as_corrs(a(2)*a(2))
expval_as_corrs(2*a(1)*a(2)*a(3))
expval_as_corrs(a(1)*a(2)*a(3)*a(4))
expval_as_corrs(adag(2)*a(1)*σz(1))

expval_as_corrs(param(:g,'r')*adag(3)*a(2)*a(2))
expval_as_corrs(-1*param(:g,'r',1)*σz(1))
expval_as_corrs(adag(:n)*σy(:k)*a(:m)*adag(:m))

let
    H = ∑(:i,param(:ω,'r',:i)*adag(:i)*a(:i)) + ∑(:j,(1//2)*param(:ωe,'r',:j)*σz(:j)) +
        ∑(:i,∑(:j,param(:g,'r',(:i,:j))*(adag(:i)+a(:i))*σx(:j)))
    comm(adag(:n)*a(:m),H)
    comm(adag(:n)*a(:n),H)
    for op in [a(:n),adag(:n),σx(:k),σy(:k),σz(:k)]
        latex(op),latex(comm(op,H))
    end
    for op in [adag(:n)*σz(:k),a(:n)*σz(:k)]
        latex(op),latex(comm(op,H))
    end

    for op in [adag(:n)*σz(:k),a(:n)*σz(:k)]
        for x = 1:5
            op = comm(op,H)
        end
    end

    expval(H)
    #expval_as_corrs(H)
end

vacA(Avac(σp(1)*σm(1)))
vacA(Avac(σm(1)*σp(1)))

vacExpVal(adag(),a())

for i2 = (:n,:m)
    stateop = (1/√2)*adag(:n)*adag(:n)*adag(:n) + (1/√2)*adag(i2)
    for A in [adag(:n)*a(:n),adag(:n)*adag(:n)*a(:n)*a(:n)]
        latex(A),latex(vacExpVal(A,stateop))
    end
end

function dotest(H,A,n=5)
    for _ = 1:n
        A = normal_form(comm(H,A))
    end
    A
end

H = ∑(:n,∑(:m,∑(:k,Pr"ω_n,m"*adag(:n,:k)*a(:m,:k)))) + ∑(:i,1//2*Pr"ν_i"*σz(:i)) + ∑(:n,∑(:k,∑(:i,Pr"g_i,n,k"*σx(:i)*(adag(:n,:k)+a(:n,:k)))))
dotest(H,σz(:j))

allops = (a(1),f(:α),σz(:i),a(2),f(:β),σz(:j),a(3),f(:γ),σz(:k))
for N = 1:9
    ops = allops[1:N]
    A = normal_form(8*Pr"g_i"*prod(ops))
    expval_as_corrs(A)
end


text:
"""
    ArbReal(x)
    ArbReal(x; [bits])
    ArbReal(x; [digits [, base=10]])

- `ArbReal`(x) == `ArbReal`(x, bits=workingprecision(ArbReal))
- `ArbReal`(x, digits=digits, base=2) == `ArbReal`(x, bits=digits)
- `ArbReal`(x, digits=digits, base=10) == `ArbReal`(x, digits=digits)

ArbReals are arbitrary precision binary floating-point regions with some specific precision established at construction.

The precision of an `ArbReal` must be >= 24 bits or >= 8 digits.

Internally, an `ArbReal` represents an interval given by its midpoint and radius (half-width of the interval).
This representation serves as a _ball_ over the real numbers.  Calculations with `ArbReals` generate results
that are assured to contain the true mathematical result.  There is no a priori assurance on the width of
the interval that results.  Good practice is to set the precision substantively greater than the precision
of the resultant width required, and to check the radius of results.

See also: [`ArbFloat`](@ref), [`ArbComplex`](@ref), [`ball`](@ref), [`setball`](@ref), [`midpoint`](@ref), [`radius`](@ref)
""" ArbReal

# ArbReal structs hold balls given as midpoint, radius

mutable struct ArbReal{P} <: Real     # P is the precision in bits
                    #      midpoint
    mid_exp::Int    # fmpz         exponent of 2 (2^exp)
    mid_size::UInt  # mp_size_t    nwords and sign (lsb holds sign of significand)
    mid_d1::UInt    # significand  unsigned, immediate value or the initial span
    mid_d2::UInt    #   (d1, d2)   the final part indicating the significand, or 0
                    #      radius
    rad_exp::Int    # fmpz       exponent of 2 (2^exp)
    rad_man::UInt   # mp_limb_t  radius, unsigned by definition

    function ArbReal{P}() where {P}
        z = new{P}(0,0,0,0,0,0)
        ccall(@libarb(arb_init), Cvoid, (Ref{ArbReal},), z)
        finalizer(arb_clear, z)
        return z
    end
end

# for use within structs, e.g. ArbRealMatrix
const PtrToArbReal = Ptr{ArbReal} # arb_ptr
const PtrToPtrToArbReal = Ptr{Ptr{ArbReal}} # arb_ptr*

arb_clear(x::ArbReal{P}) where {P} = ccall(@libarb(arb_clear), Cvoid, (Ref{ArbReal},), x)

ArbReal{P}(x::ArbReal{P}) where {P} = x
ArbReal(x::ArbReal{P}) where {P} = x

ArbReal{P}(x::Missing) where {P} = missing
ArbReal(x::Missing) = missing

@inline sign_bit(x::ArbReal{P}) where {P} = isodd(x.mid_size)

function ArbReal(x::T; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Number}
    bits > 0 && bits < MINIMUM_PRECISION_BASE2 && throw(DomainError("bit precision $bits < $MINIMUM_PRECISION_BASE2"))
    digits > 0 && digits < MINIMUM_PRECISION_BASE10 && throw(DomainError("digit precision $digits < $MINIMUM_PRECISION_BASE10"))
    if base === 10
        bits = digits > 0 ? bits4digits(digits)+extrabits() : (bits > 0 ? bits+extrabits() : DEFAULT_PRECISION.x)
    elseif base === 2
        bits = bits > 0 ? bits+extrabits() : (digits > 0 ? digits+extrabits() : DEFAULT_PRECISION.x)
    else
        throw(ErrorException("base expects 2 or 10"))
    end
    ArbReal{bits}(x)
end

swap(x::ArbReal{P}, y::ArbReal{P}) where {P} = ccall(@libarb(arb_swap), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), x, y)

function copy(x::ArbReal{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_set), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), z, x)
    return z
end

function ArbReal{P}(x::Int32) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_set_si), Cvoid, (Ref{ArbReal}, Clong), z, x)
    return z
end
ArbReal{P}(x::T) where {P, T<:Union{Int8, Int16}} = ArbReal{P}(Int32(x))
ArbReal{P}(x::T) where {P, T<:Union{Int64, Int128}} = ArbReal{P}(BigInt(x))

function ArbReal{P}(x::UInt32) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_set_ui), Cvoid, (Ref{ArbReal}, Culong), z, x)
    return z
end
ArbReal{P}(x::T) where {P, T<:Union{UInt8, UInt16}} = ArbReal{P}(UInt32(x))
ArbReal{P}(x::T) where {P, T<:Union{UInt64, UInt128}} = ArbReal{P}(BigInt(x))

function ArbReal{P}(x::Float64) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_set_d), Cvoid, (Ref{ArbReal}, Cdouble), z, x)
    return z
end
ArbReal{P}(x::T) where {P, T<:Union{Float16, Float32}} = ArbReal{P}(Float64(x))

function ArbReal{P}(x::BigFloat) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_set_interval_mpfr), Cvoid, (Ref{ArbReal}, Ref{BigFloat}, Ref{BigFloat}, Clong), z, x, x, P)
    return z
end
ArbReal{P}(x::BigInt) where {P} = ArbReal{P}(BigFloat(x))
ArbReal{P}(x::Rational{T}) where {P, T<:Signed} = ArbReal{P}(BigFloat(x))

BigFloat(x::ArbReal{P}) where {P} = BigFloat(ArbFloat{P}(x))
BigInt(x::ArbReal{P}) where {P} = BigInt(trunc(BigFloat(ArbFloat{P}(x))))

function Base.Integer(x::ArbReal{P}) where {P}
    if isinteger(x)
       abs(x) <= typemax(Int64) ? Int64(x) : BigInt(x)
    else
       throw(InexactError("$x"))
    end
end

function ArbReal{P}(x::ArbFloat{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), z, x)
    return z
end

for (F,A) in ((:floor, :arf_floor), (:ceil, :arf_ceil))
    @eval begin
        function $F(x::ArbReal{P}) where {P}
            z = ArbFloat{P}()
            ccall(@libarb($A), Cvoid, (Ref{ArbFloat}, Ref{ArbFloat}), z, x)
            return ArbReal{P}(z)
        end
        $F(::Type{T}, x::ArbReal{P}) where {P,T<:Integer} = T($F(x))
    end
end

function trunc(x::ArbReal{P}) where {P}
    newlo = trunc(ArbFloat(lowerbound(x)))
    newhi = trunc(ArbFloat(upperbound(x)))
    result = setinterval(newlo, newhi)
    return result
end

trunc(::Type{T}, x::ArbReal{P}) where {P, T<:Integer} = T(trunc(x))

function modf(x::ArbReal{P}) where {P}
    ipart = trunc(x)
    fpart = x - ipart
    return (fpart, ipart)
end

fmod(fpartipart::Tuple{ArbReal{P}, ArbReal{P}}) where {P} =
    fpartipart[1] + fpartipart[2]
fmod(fpart::ArbReal{P}, ipart::ArbReal{P}) where {P}=
    fpart + ipart

div(x::ArbReal{P}, y::ArbReal{P}) where {P} =
    trunc(x / y)

rem(x::ArbReal{P}, y::ArbReal{P}) where {P} =
    x - (div(x,y) * y)

function divrem(x::ArbReal{P}, y::ArbReal{P}) where {P}
    dv = div(x,y)
    rm = x - (dv * y)
    return (dv, rm)
end

fld(x::ArbReal{P}, y::ArbReal{P}) where {P} =
    floor(x / y)

mod(x::ArbReal{P}, y::ArbReal{P}) where {P} =
    x - (fld(x,y) * y)

function fldmod(x::ArbReal{P}, y::ArbReal{P}) where {P}
    fd = fld(x,y)
    md = x - (fd * y)
    return (fd, md)
end

cld(x::ArbReal{P}, y::ArbReal{P}) where {P} =
    ceil(x / y)

function midpoint(x::ArbReal{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_get_mid_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), z, x)
    return z
end

function radius(x::ArbReal{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_get_rad_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), z, x)
    return z
end

function midpoint_byref(x::ArbReal{P}) where {P}
    w = ArbReal{P}()
    ccall(@libarb(arb_get_mid_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), w, x)
    z = ArbFloat{P}(w)
    return z
end

function radius_byref(x::ArbReal{P}) where {P}
    w = ArbReal{P}()
    ccall(@libarb(arb_get_rad_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), w, x)
    z = ArbFloat{P}(w)
    return z
end

midpoint(x::ArbReal{P}, ::Type{ArbFloat{P}}) where {P} = midpoint_byref(x)
radius(x::ArbReal{P}, ::Type{Mag}) where {P} = radius_byref(x)

function radius(x::ArbReal{P}, ::Type{ArbFloat{P}}) where {P}
    z = ArbFloat{P}()
    mag = radius_byref(x)
    ccall(@libarb(arf_set_mag), Cvoid, (Ref{ArbFloat}, Ref{Mag}), z, mag)
    return z
end

# a type specific hash function helps the type to 'just work'
const hash_arbreal_lo = (UInt === UInt64) ? 0x37e642589da3416a : 0x5d46a6b4
const hash_0_arbreal_lo = hash(zero(UInt), hash_arbreal_lo)
# two values of the same precision
#    with identical midpoint significands and identical radus_exponentOf2s hash equal
# they are the same value, one is less accurate yet centered about the other
Base.hash(z::ArbReal{P}, h::UInt) where {P} =
    hash(z.d1 ⊻ z.exp,
         (h ⊻ hash(z.d2 ⊻ (~reinterpret(UInt,P)), hash_arbreal_lo)
            ⊻ hash_0_arbreal_lo))


text:
function SplittablesBase.halve(x::AbstractReadOnlyContainers)
    left, right = halve(parent(x))
    return (constructorof(typeof(x))(left), constructorof(typeof(x))(right))
end

SplittablesBase.amount(x::AbstractReadOnlyContainers) = amount(parent(x))


text:
using ExistingProcessManagers
using Test
using Distributed

@testset "Unit tests" begin
    @testset "constructors.jl" begin
        hosts_and_ports = [
            ("127.0.0.1", 9601), # the host is "127.0.0.1", and the port number is 9601
            ("127.0.0.1", 9602), # the host is "127.0.0.1", and the port number is 9602
        ]
        @test ExistingProcessManager(hosts_and_ports) isa ExistingProcessManager
    end
end


text:
"""Number of pixels in full CAR ring of this WCS."""
fullringsize(wcs::AbstractWCSTransform) =
    round(Int64, abs(2π / (getunit(wcs) * getcdelt(wcs)[1])))

"""Number of rings in a fullsky version of this WCS."""
fullringnum(wcs::AbstractWCSTransform) =
    1 + round(Int64, abs(π / (getunit(wcs) * getcdelt(wcs)[2])))

"""Get index of ring given angle"""
function first_last_rings_in_fullsky(shape, wcs::AbstractWCSTransform)
    Δδ = getcdelt(wcs)[end] * getunit(wcs)
    Δθ = abs(Δδ)
    δ₁ = pix2sky(shape, wcs, 1, 1)[2]
    δ₂ = pix2sky(shape, wcs, 1, shape[2])[2]
    θ₁ = π/2 - δ₁
    θ₂ = π/2 - δ₂
    i1 = round(Int64, θ₁ / Δθ) + 1
    i2 = round(Int64, θ₂ / Δθ) + 1

    return i1:sign(i2-i1):i2
end

# libsharp wants ascending colatitude (θ) and increasing right ascension (ϕ)
function get_flip_slices(shape, wcs)
    Δα, Δδ = getcdelt(wcs) .* getunit(wcs)
    flipx_slice = (Δα ≥ 0) ? (1:shape[1]) : (shape[1]:-1:1)
    flipy_slice = (Δδ ≤ 0) ? (1:shape[2]) : (shape[2]:-1:1)
    return flipx_slice, flipy_slice
end

"""Clenshaw-Curtis quadrature weights, make a Libsharp GeomInfo"""
function make_cc_geom_info(shape, wcs₀::AbstractWCSTransform)

    flipx_slice, flipy_slice = get_flip_slices(shape, wcs₀)
    _, wcs = slice_geometry(shape, wcs₀, flipx_slice, flipy_slice, shape[3:end])
    subinds = first_last_rings_in_fullsky(shape, wcs)
    @assert first(subinds) ≤ last(subinds)  # vertical angle must be increasing
    nringstot = fullringnum(wcs)
    ppring = fullringsize(wcs)
    phi0 = pix2sky(shape, wcs, 1, 2)[1]  # since RA = ϕ

    geom_info_ptr = Ref{Ptr{Cvoid}}()
    μ = chebyshevjacobimoments1(Float64, nringstot, 0.0, 0.0)
    weights = Cdouble.(clenshawcurtisweights(μ)[subinds] .* 2π ./ ppring)
    theta = Cdouble.(range(0, π, length=nringstot)[subinds])

    # constant stuff
    nsubrings = length(subinds)
    phi0s = fill(Cdouble(phi0), nsubrings)
    nph = fill(Cint(ppring), nsubrings)
    stride = ones(Cint, nsubrings)
    offsets = Cptrdiff_t.(ppring .* (0:(nsubrings-1)))

    ccall(
        (:sharp_make_geom_info, Libsharp.libsharp2),
        Cvoid,
        (Cint,       Ref{Cint}, Ref{Cptrdiff_t}, Ref{Cint}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Ptr{Cvoid}}),
        Cint(nsubrings),nph,  offsets,           stride,    phi0s,         theta,        weights,      geom_info_ptr
    )

    return GeomInfo(geom_info_ptr[])
end

"""Complete rings for spherical harmonic transforms. Also converts to AbstractArray{Float64,2}."""
function create_sht_band(m::Enmap)
    shape = size(m)
    flipx_slice, flipy_slice = get_flip_slices(shape, getwcs(m))
    _, wcs = slice_geometry(shape, getwcs(m), flipx_slice, flipy_slice, shape[3:end])
    target_size = (fullringsize(wcs), size(m,2), shape[3:end]...)
    converted = Enmap(zeros(Float64, target_size), wcs)
    other_slices = axes(m)[3:end]

    # extend the ring by adding more columns to the end
    converted.data[1:shape[1], :, other_slices...] .= view(m.data, flipx_slice, flipy_slice, other_slices...)
    return converted
end

function create_sht_band(shape, wcs::AbstractWCSTransform)
    target_size = (fullringsize(wcs), shape[2:end]...)
    return Enmap(zeros(Float64, target_size), wcs)
end

"""Generate an estimate of the Nyquist frequency for a map"""
getlmax(wcs) = fullringsize(wcs) ÷ 2

"""perform forward SHT of an intensity map"""
function map2alm(input_map::Enmap{T,2}; lmax=nothing, mmax=lmax) where T
    if isnothing(lmax)
        lmax = getlmax(getwcs(input_map))
        mmax = lmax
    end
    band = create_sht_band(input_map)
    alm_info = make_triangular_alm_info(lmax, mmax, 1)
    geom_info = make_cc_geom_info(size(band), getwcs(band))
    nalms = alm_count(alm_info)
    npix = map_size(geom_info)
    alm = Alm(lmax, mmax, zeros(ComplexF64, nalms))
    flat_map = reshape(band.data, npix)

    sharp_execute!(
        SHARP_MAP2ALM, 0,
        [alm.alm],
        [flat_map],
        geom_info, alm_info, SHARP_DP
    )
    return alm
end

# spin 2 transforms ----

# 2-tuple of scalar enmaps means QU transform
function map2alm(maps::NTuple{2,EM}; lmax=nothing, mmax=lmax) where {T, EM<:Enmap{T,2}}
    if isnothing(lmax)
        lmax = getlmax(getwcs(first(maps)))
        mmax = lmax
    end
    band_maps = create_sht_band.(maps)
    geom_info = make_cc_geom_info(size(band_maps[1]), getwcs(band_maps[1]))
    npix = map_size(geom_info)
    flat_maps = [reshape(m.data, npix) for m in band_maps]

    alm_info = make_triangular_alm_info(lmax, mmax, 1)
    nalm = alm_count(alm_info)
    alms = [Alm(lmax, mmax, zeros(ComplexF64, nalm)) for m in band_maps]

    sharp_execute!(
        SHARP_MAP2ALM, 2,
        [a.alm for a in alms],
        flat_maps,
        geom_info, alm_info, SHARP_DP)

    return alms[1], alms[2]
end

# 3-tuple of scalar enmaps means IQU transform
function map2alm(maps::NTuple{3,EM}; lmax=nothing, mmax=lmax) where {T, EM<:Enmap{T,2}}
    if isnothing(lmax)
        lmax = getlmax(getwcs(first(maps)))
        mmax = lmax
    end
    alm_i = map2alm(maps[1]; lmax=lmax, mmax=mmax)
    alm_e, alm_b = map2alm(maps[2:3]; lmax=lmax, mmax=mmax)
    return alm_i, alm_e, alm_b
end

# handle enmaps that are stacks of I, QU, and IQU in the third dim
function map2alm(input_map::Enmap{T,3}; lmax=nothing, mmax=lmax) where T
    if isnothing(lmax)
        lmax = getlmax(getwcs(input_map))
        mmax = lmax
    end
    if size(input_map,3) == 1  # this I map for some reason has a single third dim
        sliced_1d_map = @view input_map[:,:,1]
        return map2alm(sliced_1d_map; lmax=lmax, mmax=mmax)
    elseif size(input_map,3) == 2  # QU map
        return _map2almQU(input_map; lmax=lmax, mmax=mmax)
    elseif size(input_map,3) == 3  # IQU map
        alm_t = map2alm(@view(input_map[:,:,1]); lmax=lmax, mmax=mmax)
        alm_e, alm_b = _map2almQU(@view(input_map[:,:,2:3]); lmax=lmax, mmax=mmax)
        return alm_t, alm_e, alm_b
    end
    throw(ArgumentError("SHTs require shape (nx,ny,ncomp) with 1 ≤ ncomp ≤ 3, for I, QU, and IQU."))
end

# special case QU: have to manually ccall with the columns when
# Enmap is a stack (nx, ny, 2)
# (internal) perform forward SHT of an QU map
function _map2almQU(input_map::Enmap{T,3}; lmax, mmax) where T

    @assert size(input_map,3) == 2  # check that we have Q and U channels

    band_map = create_sht_band(input_map)
    geom_info = make_cc_geom_info(size(band_map), getwcs(band_map))
    npix = map_size(geom_info)
    flat_map = reshape(band_map.data, npix * 2)

    alm_info = make_triangular_alm_info(lmax, mmax, 1)
    nalm = alm_count(alm_info)
    alm_e = Alm(lmax, mmax, zeros(ComplexF64, nalm))
    alm_b = Alm(lmax, mmax, zeros(ComplexF64, nalm))

    # have to manually ccall here since I didn't add the pointers directly
    GC.@preserve alm_e alm_b flat_map ccall(
        (:sharp_execute, Libsharp.libsharp2),
        Cvoid,
        (Cint, Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Cint,
            Ref{Cdouble}, Ref{Culonglong}),
        SHARP_MAP2ALM, 2,
        [pointer(alm_e.alm), pointer(alm_b.alm)],
        [pointer(flat_map, 1), pointer(flat_map, 1+npix)],
        geom_info.ptr, alm_info.ptr, SHARP_DP,
        Ptr{Cdouble}(C_NULL), Ptr{Culonglong}(C_NULL))

    return alm_e, alm_b
end

## inverse

"""perform inverse SHT of an intensity map"""
function alm2map(alm::Alm, shape, wcs::AbstractWCSTransform)

    band = create_sht_band(shape[1:2], wcs)
    geom_info = make_cc_geom_info(size(band), getwcs(band))
    npix = map_size(geom_info)
    flat_map = reshape(band.data, npix)
    alm_info = make_triangular_alm_info(alm.lmax, alm.mmax, 1)

    sharp_execute!(
        SHARP_ALM2MAP, 0,
        [alm.alm],
        [flat_map],
        geom_info, alm_info, SHARP_DP)

    Δα, Δδ = getcdelt(wcs)
    band_width = size(band,1)
    xslice = (Δα ≥ 0) ? (1:shape[1]) : (band_width:-1:(band_width-shape[1]+1))
    yslice = (Δδ ≤ 0) ? (1:shape[2]) : (shape[2]:-1:1)

    return Enmap(band.data[xslice, yslice], wcs)
end

function alm2map(alms::NTuple{2,Alm}, shape, wcs::AbstractWCSTransform)
    lmax, mmax = first(alms).lmax, first(alms).mmax
    alm_info = make_triangular_alm_info(first(alms).lmax, first(alms).mmax, 1)

    band_maps = (create_sht_band(shape[1:2], wcs), create_sht_band(shape[1:2], wcs))
    geom_info = make_cc_geom_info(size(first(band_maps)), getwcs(first(band_maps)))
    npix = map_size(geom_info)
    flat_maps = [reshape(band_maps[1].data, npix), reshape(band_maps[2].data, npix)]

    npix = map_size(geom_info)
    alm_info = make_triangular_alm_info(lmax, mmax, 1)

    sharp_execute!(
        SHARP_ALM2MAP, 2,
        [a.alm for a in alms],
        flat_maps,
        geom_info, alm_info, SHARP_DP)

    Δα, Δδ = getcdelt(wcs)
    band_width = size(first(band_maps),1)
    xslice = (Δα ≥ 0) ? (1:shape[1]) : (band_width:-1:(band_width-shape[1]+1))
    yslice = (Δδ ≤ 0) ? (1:shape[2]) : (shape[2]:-1:1)

    return [Enmap(b.data[xslice, yslice], wcs) for b in band_maps]
end

alm2map(alms::NTuple{3,Alm}, shape, wcs::AbstractWCSTransform) =
    (alm2map(alms[1], shape, wcs), alm2map(alms[2:3], shape, wcs)...)

function alm2map(alms::Vector{A}, shape, wcs::AbstractWCSTransform) where {A <: Alm}
    if length(alms) == 1
        return alm2map(alms[1], shape, wcs)
    elseif length(alms) == 2
        return alm2map((alms[1], alms[2]), shape, wcs)
    elseif length(alms) == 3
        return alm2map((alms[1], alms[2], alms[3]), shape, wcs)
    end
end


text:
using Libdl
# Check for an iconv implementation with the GNU (non-POSIX) behavior:
# EILSEQ is returned when a sequence cannot be converted to target encoding,
# instead of succeeding and only returning the number of invalid conversions
# This non-standard behavior is required to allow replacing invalid sequences
# with a user-defined character.
# Implementations with this behavior include glibc, GNU libiconv (on which Mac
# OS X's is based) and win_iconv.
function validate_iconv(lib, iconv_open, iconv_close, iconv)
    h = dlopen_e(lib)
    h == C_NULL && return false
    # Needed to check libc
    f = dlsym_e(h, iconv_open)
    f == C_NULL && return false

    cd = ccall(f, Ptr{Nothing}, (Cstring, Cstring), "ASCII", "UTF-8")
    cd == Ptr{Nothing}(-1) && return false

    s = "café"
    a = Vector{UInt8}(undef, sizeof(s))
    inbufptr = Ref{Ptr{UInt8}}(pointer(s))
    inbytesleft = Ref{Csize_t}(sizeof(s))
    outbufptr = Ref{Ptr{UInt8}}(pointer(a))
    outbytesleft = Ref{Csize_t}(length(a))
    ret = ccall(dlsym_e(h, iconv), Csize_t,
                (Ptr{Nothing}, Ptr{Ptr{UInt8}}, Ref{Csize_t}, Ptr{Ptr{UInt8}}, Ref{Csize_t}),
                cd, inbufptr, inbytesleft, outbufptr, outbytesleft)
    ccall(dlsym_e(h, iconv_close), Nothing, (Ptr{Nothing},), cd) == -1 && return false

    return ret == -1 % Csize_t && Libc.errno() == Libc.EILSEQ
end

using BinaryProvider

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libiconv"], :libiconv)
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaStrings/IConvBuilder/releases/download/v1.15+build.3"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/IConv.aarch64-linux-gnu.tar.gz", "7a8d86f866ebee3218d590f74887d464c9980125a2e0c590acdcc0cfc5c84dd5"),
    Linux(:aarch64, :musl) => ("$bin_prefix/IConv.aarch64-linux-musl.tar.gz", "a37fb558319c37970e70bdc786d6bfe84f95ca27a72f6d334a1411d9c305ee2c"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/IConv.arm-linux-gnueabihf.tar.gz", "c3e7eaf54fe3722153e29e29f1e08f7acdbe432fb1bcf8e0d72f15e710ff0c03"),
    Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/IConv.arm-linux-musleabihf.tar.gz", "200db91435875f57ea827ef69b31b6d4c11f089aacc695846699155488245b1d"),
    Linux(:i686, :glibc) => ("$bin_prefix/IConv.i686-linux-gnu.tar.gz", "c0dea5b3c6ebfc0ae0d3b3c1affe25a93bb11c5da8dd34c2a40ffb132e61e4ee"),
    Linux(:i686, :musl) => ("$bin_prefix/IConv.i686-linux-musl.tar.gz", "96cccd8a539c3ad4a70ddade4396f91015248d6fad799adcf0f79f3940d1c672"),
    Windows(:i686) => ("$bin_prefix/IConv.i686-w64-mingw32.tar.gz", "7d8b600337c4153f9a4ce2115f071aabc9a8c5565cb95cad4c6cc3dd6df17fb3"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/IConv.powerpc64le-linux-gnu.tar.gz", "8e7edc2adbd331cb1f0a944f7eb26ed56ac884830d4d65785671c0e620be84e6"),
    MacOS(:x86_64) => ("$bin_prefix/IConv.x86_64-apple-darwin14.tar.gz", "5b585a372c272aed10823f5ea4473849396ee928f9f6b7643411138549349bf2"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/IConv.x86_64-linux-gnu.tar.gz", "6f58cd9722f492aa0e519c9d91082c0b0776bd53cbff17f643172fec2c4d4347"),
    Linux(:x86_64, :musl) => ("$bin_prefix/IConv.x86_64-linux-musl.tar.gz", "c72aaabc45465fdb2b7d74df2bd4a0282cf1d2a46e55868b129003547396b69f"),
    FreeBSD(:x86_64) => ("$bin_prefix/IConv.x86_64-unknown-freebsd11.1.tar.gz", "69a5585871a61547e71cd06d32269e972f2d7ab488f82ffb852b3b5beaa83d8f"),
    Windows(:x86_64) => ("$bin_prefix/IConv.x86_64-w64-mingw32.tar.gz", "c553d5d36ad831c13572b1c06230b41471b1dea560fbc03d5213cf4304de74cf"),
)

# Detect already present libc iconv or libiconv
# (notably for Linux, Mac OS and other Unixes)
global found_iconv = false
for lib in ("libc", "libc-bin", "libiconv", "iconv")
    if lib in ("libc", "libc-bin")
        global iconv_open = :iconv_open
        global iconv_close = :iconv_close
        global iconv = :iconv
    else
        global iconv_open =  :libiconv_open
        global iconv_close = :libiconv_close
        global iconv = :libiconv
    end
    if validate_iconv(lib, iconv_open, iconv_close, iconv)
        global found_iconv = true
        write(joinpath(@__DIR__, "deps.jl"),
              """
              ## This file autogenerated by Pkg.build(\\"StringEncodings\\").
              ## Do not edit.
              const libiconv = "$lib"
              function check_deps()
                  global libiconv
                  if Libdl.dlopen_e(libiconv) == C_NULL
                      error("\$(libiconv) cannot be opened. Please re-run Pkg.build(\\"StringEncodings\\"), and restart Julia.")
                  end

              end
              """)
        break
    end
end

if !found_iconv
    iconv_open =  :libiconv_open
    iconv_close = :libiconv_close
    iconv = :libiconv
    # Install unsatisfied or updated dependencies:
    unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
    if haskey(download_info, platform_key())
        url, tarball_hash = download_info[platform_key()]
        if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
            # Download and install binaries
            install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
        end
    elseif unsatisfied
        # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
        # Alternatively, you could attempt to install from a separate provider,
        # build from source or something even more ambitious here.
        error("Your platform $(triplet(platform_key())) is not supported by this package! Try installing libiconv manually.")
    end

    # Write out a deps.jl file that will contain mappings for our products
    write_deps_file(joinpath(@__DIR__, "deps.jl"), products)
end

open(joinpath(@__DIR__, "deps.jl"), "a") do io
    write(io,
          """
          const iconv_open_s = :$iconv_open
          const iconv_close_s = :$iconv_close
          const iconv_s = :$iconv
          """)
end


text:
__precompile__()

module OrdinalGWAS

using LinearAlgebra
using CSV, DataFrames, Distributions, Reexport
using SnpArrays, VCFTools, VariantCallFormat, BGEN
using MathOptInterface
const MOI = MathOptInterface
@reexport using OrdinalMultinomialModels

export ordinalgwas

include("gwas.jl")

end

max_stars_repo_name: Hua-Zhou/PolrGWAS

max_issues_repo_name: Hua-Zhou/PolrGWAS


text:
module HighDimMixedModels

# Write your package code here.

#using MixedModels: modelmatrix
using Base: Float64, AbstractVecOrTuple
using DataFrames: Dict
#using LinearAlgebra: AbstractMatrix, include
using StatsModels
using LinearAlgebra
using DataFrames
using MixedModels
using NLopt
using CategoricalArrays

import Base: *
import MixedModels: fit, fit!

# abstract type MixedModel{T} <: StatsModels.RegressionModel end # model with fixed and random effects

include("bricks.jl")
#include("bricks_Number.jl")

export highDimMixedModel,
    modelmatrix,
    fit,
    fit!,
    highDimMat,
    XMat,
    ReMat

end

max_stars_repo_name: solislemuslab/HighDimMixedModels

max_issues_repo_name: solislemuslab/HighDimMixedModels

max_forks_repo_name: solislemuslab/HighDimMixedModels


text:
module MoreEntities

using ..Ahorn, Maple

glassblockcodename = "glassBlock"
glassblockbgcodename = "BGGlassBlock"

Comment(x::Integer, y::Integer, comment::String) = Maple.Entity("comment", x=x, y=y, comment=comment)

GlassBlock(x::Integer, y::Integer, width::Integer=16, height::Integer=16, sinks::Bool=false) = Maple.Entity(glassblockcodename, x=x, y=y, width=width, height=height, sinks=sinks)
GlassBlockBG(x::Integer, y::Integer, width::Integer=16, height::Integer=16, sinks::Bool=false) = Maple.Entity(glassblockbgcodename, x=x, y=y, width=width, height=height, sinks=sinks)

RefillUpdated(x::Integer, y::Integer, twoDash::Bool, oneUse::Bool) = Maple.Entity("refill", x=x, y=y, twoDash=twoDash, oneUse=oneUse)
KeyUpdated(x::Integer, y::Integer) = Maple.Entity("key", x=x, y=y)
DarkChaserUpdated(x::Integer, y::Integer, canChangeMusic::Bool=false) = Entity("darkChaser", x=x, y=y, canChangeMusic=canChangeMusic)
function BadelineBossUpdated(x::Integer, y::Integer, nodes::Array{Tuple{T, T}, 1}=Tuple{Integer, Integer}[], patternIndex::Integer=1, startHit::Bool=false, cameraPastY::Number=120.0, lockCamera::Bool=true, canChangeMusic::Bool=false) where {T <: Integer}
    return Entity("finalBoss", x=x, y=y, startHit=startHit, nodes=nodes, patternIndex=patternIndex, cameraPastY=cameraPastY, lockCamera=lockCamera, canChangeMusic=canChangeMusic)
end

BigWaterfall(x::Integer, y::Integer, width::Integer=16, height::Integer=32, layer::String="FG") = Entity("bigWaterfall", x=x, y=y, width=width, height=height, layer=layer)

placements = Dict{String, Ahorn.EntityPlacement}(
    "Intro Car" => Ahorn.EntityPlacement(
        Maple.IntroCar
    ),
    "Starjump Block" => Ahorn.EntityPlacement(
        Maple.StarJumpBlock,
	    "rectangle",
	    Dict{String, Any}(
		    "sinks" => true
	    )
    ),
    "Trapdoor" => Ahorn.EntityPlacement(
        Maple.TrapDoor
    ),
    "Heart Gate" => Ahorn.EntityPlacement(
	    (x, y) -> Maple.HeartDoor(x, y, 16, 96, 4),
	    "rectangle",
	    Dict{String, Any}(
		    "requires" => 4
	    )
    ),
    "Sound Source" => Ahorn.EntityPlacement(
        Maple.SoundSource
    ),
    "Memorial (Everest)" => Ahorn.EntityPlacement(
        (x, y) -> Maple.EverestMemorial(x, y, false, "MEMORIAL", "scenery/memorial/memorial")
    ),
    "Core Text (Everest)" => Ahorn.EntityPlacement(
        Maple.EverestCoreMessage
    ),
    "Core Text" => Ahorn.EntityPlacement(
        Maple.CoreMessage
    ),
    "Glass Block" => Ahorn.EntityPlacement(
        (x, y) -> GlassBlock(x, y, 16, 16, false),
        "rectangle",
        Dict{String,Any}(
            "width" => 16,
            "height" => 16
        )
    ),
    "Glass Block (BG)" => Ahorn.EntityPlacement(
        (x, y) -> GlassBlockBG(x, y, 16, 16, false),
        "rectangle",
        Dict{String,Any}(
            "width" => 16,
            "height" => 16
        )
    ),
    "Dream Mirror" => Ahorn.EntityPlacement(
        Maple.DreamMirror
    ),
    "Lightbeam" => Ahorn.EntityPlacement(
        # note to self: lightbeams are RECTANGLES you fuckin dumb eediot
        (x, y) -> Maple.Lightbeam(x, y, 16, 8, 0),
	    "rectangle"
    ),
    "Floating Debris" => Ahorn.EntityPlacement(
        Maple.FloatingDebris,
	    "point"
    ),
    "Foreground Debris" => Ahorn.EntityPlacement(
        Maple.ForegroundDebris,
	    "point"
    ),
    "Waterfall (Big)" => Ahorn.EntityPlacement(
        BigWaterfall,
        "rectangle"
    ),
    "Refill (1.2.5.0)" => Ahorn.EntityPlacement(
	    (x, y) -> RefillUpdated(x, y, false, false)
    ),
    "Key (1.2.6.0)" => Ahorn.EntityPlacement(
        KeyUpdated
    ),
    "Comment (Ahorn)" => Ahorn.EntityPlacement(
        (x, y) -> Comment(x, y, "Insert your comment here..."),
        "point",
        Dict{String,Any}(
            "comment" => "Insert your comment here..."
        )
    ),
    "Badeline Chaser (Everest)" => Ahorn.EntityPlacement(
	    (x, y) -> DarkChaserUpdated(x, y, false)
    ),
    "Badeline Boss (Everest)" => Ahorn.EntityPlacement(
	    BadelineBossUpdated
    ),
)

function nodeLimits(entity::Maple.Entity)
    if entity.name == "key"
        return true, 0, -1
    end
end

function editingOptions(entity::Maple.Entity)
    if entity.name == "bigWaterfall"
        return true, Dict{String, Any}(
            "layer" => String["FG", "BG"]
        )
    end
end

function minimumSize(entity::Maple.Entity)
    if entity.name == "starJumpBlock" ||
        entity.name == glassblockcodename ||
        entity.name == glassblockbgcodename ||
        entity.name == "bigWaterfall"
        return true, 8, 8
    elseif entity.name == "heartGemDoor"
        return true, 16, 1
    elseif entity.name == "lightbeam"
        return true, 6, 6
    end
end

function resizable(entity::Maple.Entity)
    if entity.name == "starJumpBlock" ||
        entity.name == "heartGemDoor" ||
        entity.name == glassblockcodename ||
        entity.name == glassblockbgcodename ||
        entity.name == "bigWaterfall" ||
        entity.name == "lightbeam"
        return true, true, true
    end
end

function renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Entity)
    if entity.name == "key"
        x, y = Ahorn.entityTranslation(entity)
        nodes = get(entity.data, "nodes", ())

        if !isempty(nodes)
            nx, ny = Int.(nodes[1])
            Ahorn.drawSprite(ctx, "collectables/key/idle00.png", nx, ny)
            Ahorn.drawArrow(ctx, x, y, nx, ny, Ahorn.colors.selection_selected_fc, headLength=6)
        end
    end
end

function selection(entity::Maple.Entity)
    if entity.name == "darkChaser"
        x, y = Ahorn.entityTranslation(entity)

        return true, Ahorn.Rectangle(x - 2, y - 16, 12, 16)
    elseif entity.name == "floatingDebris"
        x, y = Ahorn.entityTranslation(entity)
        return true, Ahorn.Rectangle(x - 6, y - 6, 12, 12)
    elseif entity.name == "foregroundDebris"
        x, y = Ahorn.entityTranslation(entity)
        return true, Ahorn.Rectangle(x - 24, y - 1, 24, 24)
    elseif entity.name == "refill"
        x, y = Ahorn.entityTranslation(entity)
        return true, Ahorn.Rectangle(x - 6, y - 6, 12, 12)
    elseif entity.name == "key"
        x, y = Ahorn.entityTranslation(entity)
        nodes = get(entity.data, "nodes", ())

        if isempty(nodes)
            return true, Ahorn.Rectangle(x - 8, y - 8, 16, 16)
        else
            nx, ny = Int.(nodes[1])
            return true, [Ahorn.Rectangle(x - 8, y - 8, 16, 16), Ahorn.Rectangle(nx - 8, ny - 8, 16, 16)]
        end
    elseif entity.name == "introCar"
        x, y = Ahorn.entityTranslation(entity)
        return true, Ahorn.Rectangle(x - 22, y - 18, 47, 18)
    elseif entity.name == "trapdoor"
        x, y = Ahorn.entityTranslation(entity)
        return true, Ahorn.Rectangle(x, y + 4, 24, 8)
    elseif entity.name == "soundsource"
        x, y = Ahorn.entityTranslation(entity)
        return true, Ahorn.Rectangle(x, y, 8, 8)
    elseif entity.name == "everest/coreMessage" || entity.name == "coreMessage"
        x, y = Ahorn.entityTranslation(entity)
        return true, Ahorn.Rectangle(x, y, 8, 8)
    elseif entity.name == "heartGemDoor"
        x, y = Ahorn.entityTranslation(entity)
	    width = Int(get(entity.data, "width", 8))
	    height = Int(get(entity.data, "height", 8))
        return true, Ahorn.Rectangle(x, y-216, width, 432)
    elseif entity.name == "everest/memorial"
	    x, y = Ahorn.entityTranslation(entity)
        return true, Ahorn.Rectangle(x - 20, y - 60, 40, 60)
    elseif entity.name == "dreammirror"
        x, y = Ahorn.entityTranslation(entity)
        return true, Ahorn.Rectangle(x - 32, y - 32, 64, 32)
    elseif entity.name == "lightbeam"
        x, y = Ahorn.entityTranslation(entity)
        width = Int(get(entity.data, "width", 8))
        height = Int(get(entity.data, "height", 8))
        return true, Ahorn.Rectangle(x-width/2, y, width, height)
    elseif entity.name == "comment"
        x, y = Ahorn.entityTranslation(entity)
        return true, Ahorn.Rectangle(x - 8, y - 8, 16, 16)
    # generic rectangle selection
    elseif entity.name == glassblockcodename ||
        entity.name == glassblockbgcodename ||
        entity.name == "starJumpBlock" ||
        entity.name == "bigWaterfall"
        x, y = Ahorn.entityTranslation(entity)
        width = Int(get(entity.data, "width", 8))
        height = Int(get(entity.data, "height", 8))
        return true, Ahorn.Rectangle(x, y, width, height)
    end
end

bigWaterfallColor = (135, 206, 250, 1) ./ (255, 255, 255, 1.0) .* (0.4, 0.6, 0.7, 0.7)

function radToDegree(rad::Number)
    return rad * (180/π)
end

function render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Entity, room::Maple.Room)
    if entity.name == "comment"
        Ahorn.drawSprite(ctx, "objects/ahornrender/comment.png", 0, 0)
        return true
    elseif entity.name == "darkChaser"
        Ahorn.drawSprite(ctx, "characters/badeline/sleep00.png", 4, -16)

        return true
    elseif entity.name == "key"
        Ahorn.drawSprite(ctx, "collectables/key/idle00.png", 0 - 8, 0 - 8)
        return true
    elseif entity.name == "refill"
        twoDash = Bool(get(entity.data, "twoDash", false))
        if twodash
            Ahorn.drawSprite(ctx, "objects/refillTwo/idle00.png", 0, 0)
        else
            Ahorn.drawSprite(ctx, "objects/refill/idle00.png", 0, 0)
        end

        return true
    elseif entity.name == "floatingDebris"
        x, y = Ahorn.entityTranslation(entity)
        Ahorn.drawSprite(ctx, "debris/b.png", 0, 0)
        return true
    elseif entity.name == "foregroundDebris"
        x, y = Ahorn.entityTranslation(entity)
        Ahorn.drawSprite(ctx, "scenery/fgdebris/rock_a00.png", 0, 0)
        return true
    elseif entity.name == "introCar"
	    x, y = Ahorn.entityTranslation(entity)
	    Ahorn.drawSprite(ctx, "scenery/car/wheels.png", 0, -9)
        Ahorn.drawSprite(ctx, "scenery/car/body.png", 0, -9)

        return true
    elseif entity.name == "trapdoor"
	    x, y = Ahorn.entityTranslation(entity)
        Ahorn.drawSprite(ctx, "objects/door/trap01.png", 12, 16)

        return true
    elseif entity.name == "soundsource"
	    x, y = Ahorn.entityTranslation(entity)
        Ahorn.drawSprite(ctx, "objects/ahornrender/sound.png", 4, 4)

        return true
    elseif entity.name == "everest/coreMessage" || entity.name == "coreMessage"
        x, y = Ahorn.entityTranslation(entity)
        Ahorn.drawSprite(ctx, "objects/ahornrender/textsample.png", 0, 0)
        return true
    elseif entity.name == "everest/memorial"
	    sprite = get(entity.data, "sprite", "scenery/memorial/memorial")
        Ahorn.drawSprite(ctx, sprite, 0, -32)

        return true
    elseif entity.name == "starJumpBlock"
        x = Int(get(entity.data, "x", 0))
        y = Int(get(entity.data, "y", 0))

        width = Int(get(entity.data, "width", 32))
        height = Int(get(entity.data, "height", 32))

	    Ahorn.drawRectangle(ctx, 0, 0, width, height, (0, 0, 0, 0) ./ 255, (255, 255, 255, 255) ./ 255)

        return true
    elseif entity.name == glassblockcodename || entity.name == glassblockbgcodename
        x = Int(get(entity.data, "x", 0))
        y = Int(get(entity.data, "y", 0))

        width = Int(get(entity.data, "width", 32))
        height = Int(get(entity.data, "height", 32))

	    Ahorn.drawRectangle(ctx, 0, 0, width, height, (0, 0, 0, 150) ./ 255, (255, 255, 255, 255) ./ 255)

        return true
    elseif entity.name == "heartGemDoor"
        x, y = Ahorn.entityTranslation(entity)
	    width = Int(get(entity.data, "width", 8))
	    height = Int(get(entity.data, "height", 8))
	    rwidth, rheight = room.size

	    Ahorn.drawRectangle(ctx, 0, 32, width, rheight, (76, 168, 214, 102) ./ 255, (108, 214, 235, 255) ./ 255)
	    Ahorn.drawRectangle(ctx, 0, -32-rheight, width, rheight, (76, 168, 214, 102) ./ 255, (108, 214, 235, 255) ./ 255)
        Ahorn.drawRectangle(ctx, 0, -32, width, 64, (76, 168, 214, 102) ./ 255, (108, 214, 235, 0) ./ 255)

        return true
    elseif entity.name == "dreammirror"
	    x, y = Ahorn.entityTranslation(entity)
	    Ahorn.drawSprite(ctx, "objects/mirror/glassbreak00.png", 0, -14)
        Ahorn.drawSprite(ctx, "objects/mirror/frame.png", 0, -16)

        return true
    elseif entity.name == "lightbeam"
        x, y = Ahorn.entityTranslation(entity)
	    width = Int(get(entity.data, "width", 8))
	    height = Int(get(entity.data, "height", 8))
        ang = Int(get(entity.data, "rotation", 0))

        # no idea how to make this work, please help

        # nx = cos(ang) * 32
        # ny = sin(radToDegree(ang)) * 32

        Ahorn.drawRectangle(ctx, 0-width/2, 0, width, height, (255, 255, 255, 150) ./ 255, (255, 255, 255, 255) ./ 255)
        # Ahorn.drawArrow(ctx, 0, 0, nx, ny, (1.0, 1.0, 1.0, 1.0), headLength=6)
        return true
    elseif entity.name == "bigWaterfall"
        x = Int(get(entity.data, "x", 0))
        y = Int(get(entity.data, "y", 0))
        width = Int(get(entity.data, "width", 0))
        height = Int(get(entity.data, "height", 0))

        Ahorn.drawRectangle(ctx, 0, 0, width, height, bigWaterfallColor, (0.0, 0.0, 0.0, 0.0))
        return true
    else
        return false
    end
end

end

max_stars_repo_name: thefoxbot/MoreEntities

max_issues_repo_name: thefoxbot/MoreEntities

max_forks_repo_name: thefoxbot/MoreEntities


text:
using Test
using ProbabilisticCircuits

include("../helper/pc_equals.jl")

@testset "Load and save a small SPN" begin

    function test_my_circuit(pc)

        @test pc isa ProbCircuit

        # Testing number of nodes and parameters
        @test  9 == num_parameters(pc)
        @test 20 == num_nodes(pc)

        # Testing Read Parameters
        EPS = 1e-7
        or1 = children(children(pc)[1])[2]
        @test abs(or1.log_probs[1] - (-1.6094379124341003)) < EPS
        @test abs(or1.log_probs[2] - (-1.2039728043259361)) < EPS
        @test abs(or1.log_probs[3] - (-0.916290731874155))  < EPS
        @test abs(or1.log_probs[4] - (-2.3025850929940455)) < EPS

        or2 = children(children(pc)[1])[1]
        @test abs(or2.log_probs[1] - (-2.3025850929940455))  < EPS
        @test abs(or2.log_probs[2] - (-2.3025850929940455))  < EPS
        @test abs(or2.log_probs[3] - (-2.3025850929940455))  < EPS
        @test abs(or2.log_probs[4] - (-0.35667494393873245)) < EPS

        @test abs(pc.log_probs[1] - (0.0)) < EPS
        @test check_parameter_integrity(pc)
    end

    # first load circuit from PSDD file
    pc1 = zoo_psdd("little_4var.psdd")

    mktempdir() do tmp

        spn_path = "$tmp/example.spn"
        write(spn_path, pc1)

        pc2 = read(spn_path, ProbCircuit)

        test_my_circuit(pc2)
        test_pc_equals(pc1, pc2)

    end

end

@testset "Can save SPNs with nonbinary multiplications" begin

    pc1 = fully_factorized_circuit(ProbCircuit, 10)

    mktempdir() do tmp

        spn_path = "$tmp/example.spn"
        write(spn_path, pc1)

        pc2 = read(spn_path, ProbCircuit)

        test_pc_equals(pc1, pc2)
    end

end


text:
include("splinesParallel.jl")
include("splinesSerial.jl")

knots = [0,0,0,1,1,2,2,3,3,4,4,4]
p = sqrt(2)/2.0
controls = [[-1,0,1], [-p,p,p], [0,1,1], [p,p,p],[1,0,1], [p,-p,p],
         [0,-1,1], [-p,-p,p], [-1,0,1]]
pc1 = plarBezier(S1, [[-1,0,0,1],[-1,0,1,1]])
pc2 = plarBezier(S1, [[-p,p,0,p],[-p,p,p,p]])
pc3 = plarBezier(S1, [[0,1,0,1],[0,1,1,1]])
pc4 = plarBezier(S1, [[p,p,0,p],[p,p,p,p]])
pc5 = plarBezier(S1, [[1,0,0,1],[1,0,1,1]])
pc6 = plarBezier(S1, [[p,-p,0,p],[p,-p,p,p]])
pc7 = plarBezier(S1, [[0,-1,0,1],[0,-1,1,1]])
pc8 = plarBezier(S1, [[-p,-p,0,p],[-p,-p,p,p]])
pc9 = plarBezier(S1, [[-1,0,0,1],[-1,0,1,1]])
controls = [pc1,pc2,pc3,pc4,pc5,pc6,pc7,pc8,pc9]

ptnurbs = pt_nurbs(S2,2,knots,controls)
point = [0.7, 0.2]
result = ptnurbs(point)


text:
using NCDatasets
using DataStructures
using ArgParse, JSON
using DataStructures

println("""
This program adjust the forcing file such that HMXL and TEMP/SALT are conform.
That is, if HMXL occupies n grids, then TEMP/SALT in top n grids will be set
to the top grid value.
""")

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

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

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin

        "--zdomain-file"
            help = "Vertical domain file containing z_w_top, z_w_bot in meters with z=0 is the surface."
            arg_type = String
            required = true

        "--topo-file"
            help = "File containing Nz_bot."
            arg_type = String
            required = true

        "--forcing-file"
            help = "Forcing file that contains HMXL, SALT and TEMP"
            arg_type = String
            required = true

        "--output-file"
            help = "The output file."
            arg_type = String
            required = true

    end

    return parse_args(s)
end

parsed = DataStructures.OrderedDict(parse_commandline())
JSON.print(parsed,4)

pleaseRun(`julia $(@__DIR__)/flatten_surface.jl
    --zdomain-file $(parsed["zdomain-file"])
    --topo-file $(parsed["topo-file"])
    --file-HMXL $(parsed["forcing-file"])
    --file-TRACER $(parsed["forcing-file"])
    --varname-TRACER TEMP
    --output-file "_tmp_flattened_TEMP.nc"
    --output-dimnames nlon nlat z_t
`)

pleaseRun(`cp $(parsed["forcing-file"]) $(parsed["output-file"])`)
pleaseRun(`ncks -A -v TEMP _tmp_flattened_TEMP.nc $(parsed["output-file"])`)

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# Required packages

using Vann
using PyPlot
using DataFrames

# Load data

path_inputs = Pkg.dir("Vann", "data/atnasjo")

date, tair, prec, q_obs, frac = load_data(path_inputs)

# Compute potential evapotranspiration

epot = epot_zero(date)

# Model choices

snow_choice = TinStandard
hydro_choice = Gr4j

tstep = 24.0

time = date[1]

st_snow = eval(Expr(:call, snow_choice, tstep, time, frac))
st_hydro = eval(Expr(:call, hydro_choice, tstep, time))

# Run calibration

param_snow, param_hydro = run_model_calib(st_snow, st_hydro, tair, prec, epot, q_obs)

println(param_snow)
println(param_hydro)

# Reinitilize model

st_snow = eval(Expr(:call, snow_choice, tstep, time, param_snow, frac))
st_hydro = eval(Expr(:call, hydro_choice, tstep, time, param_hydro))

# Run model with best parameter set

q_sim, snow_out, hydro_out = run_model(st_snow, st_hydro, tair, prec, epot; return_all = true)

# Plot results

plot_sim(snow_out)
plot_sim(hydro_out)

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
using DiffEqBenchmarks
DiffEqBenchmarks.weave_file(".","test.jmd")