text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8, x9, x10)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9", "x10"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8, x9, x10]);F = [6561*x1^2*x2^2+2601*x1^2*x3^2+528*x1^2*x3*x7-144*x1^2*x3*x9+170*x1^2*x7^2-66*x1^2*x7*x9+9*x1^2*x9^2-98*x1*x10^2*x7+5184*x1*x2*x4^2+9882*x1*x2*x4*x7-11826*x1*x2*x5*x7+2916*x1*x2*x5*x8-266*x1*x2*x6*x7-4590*x1*x3^3+2340*x1*x3^2*x5-4140*x1*x3^2*x6+2784*x1*x3*x5*x8+3780*x1*x3*x6^2+912*x1*x3*x6*x8+5400*x1*x3*x7^2-1190*x1*x3*x7*x9-742*x1*x4^2*x7+1276*x1*x5*x7*x8-348*x1*x5*x8*x9+418*x1*x6*x7*x8-114*x1*x6*x8*x9+49*x10^4+266*x10^2*x2*x6+1190*x10^2*x3*x9+742*x10^2*x4^2+361*x2^2*x6^2+3230*x2*x3*x6*x9+2014*x2*x4^2*x6+2601*x3^4-2652*x3^3*x5+4692*x3^3*x6+676*x3^2*x5^2-2392*x3^2*x5*x6-2168*x3^2*x6^2-6120*x3^2*x7^2+7225*x3^2*x9^2+9010*x3*x4^2*x9+2184*x3*x5*x6^2+3120*x3*x5*x7^2-3864*x3*x6^3-5520*x3*x6*x7^2+3833*x4^4+3904*x4^3*x7-4672*x4^2*x5*x7+1152*x4^2*x5*x8+3721*x4^2*x7^2-8906*x4*x5*x7^2+2196*x4*x5*x7*x8+5329*x5^2*x7^2-2628*x5^2*x7*x8+3688*x5^2*x8^2+2204*x5*x6*x8^2+1764*x6^4+5040*x6^2*x7^2+361*x6^2*x8^2+3600*x7^4+1728*x1*x10*x3+792*x1*x10*x7-216*x1*x10*x9+5022*x1*x2^2-602*x1*x5*x7+602*x10^2*x5+4176*x10*x5*x8+1368*x10*x6*x8+1984*x2*x4^2+3782*x2*x4*x7+1634*x2*x5*x6-4526*x2*x5*x7+1116*x2*x5*x8+7310*x3*x5*x9+4558*x4^2*x5+1296*x10^2+961*x2^2+1849*x5^2, 13122*x1^2*x2+5184*x1*x4^2+9882*x1*x4*x7-11826*x1*x5*x7+2916*x1*x5*x8-266*x1*x6*x7+266*x10^2*x6+722*x2*x6^2+3230*x3*x6*x9+2014*x4^2*x6+10044*x1*x2+1984*x4^2+3782*x4*x7+1634*x5*x6-4526*x5*x7+1116*x5*x8+1922*x2, 5202*x1^2*x3+528*x1^2*x7-144*x1^2*x9-13770*x1*x3^2+4680*x1*x3*x5-8280*x1*x3*x6+2784*x1*x5*x8+3780*x1*x6^2+912*x1*x6*x8+5400*x1*x7^2-1190*x1*x7*x9+1190*x10^2*x9+3230*x2*x6*x9+10404*x3^3-7956*x3^2*x5+14076*x3^2*x6+1352*x3*x5^2-4784*x3*x5*x6-4336*x3*x6^2-12240*x3*x7^2+14450*x3*x9^2+9010*x4^2*x9+2184*x5*x6^2+3120*x5*x7^2-3864*x6^3-5520*x6*x7^2+1728*x1*x10+7310*x5*x9, 10368*x1*x2*x4+9882*x1*x2*x7-1484*x1*x4*x7+1484*x10^2*x4+4028*x2*x4*x6+18020*x3*x4*x9+15332*x4^3+11712*x4^2*x7-9344*x4*x5*x7+2304*x4*x5*x8+7442*x4*x7^2-8906*x5*x7^2+2196*x5*x7*x8+3968*x2*x4+3782*x2*x7+9116*x4*x5, -11826*x1*x2*x7+2916*x1*x2*x8+2340*x1*x3^2+2784*x1*x3*x8+1276*x1*x7*x8-348*x1*x8*x9-2652*x3^3+1352*x3^2*x5-2392*x3^2*x6+2184*x3*x6^2+3120*x3*x7^2-4672*x4^2*x7+1152*x4^2*x8-8906*x4*x7^2+2196*x4*x7*x8+10658*x5*x7^2-5256*x5*x7*x8+7376*x5*x8^2+2204*x6*x8^2-602*x1*x7+602*x10^2+4176*x10*x8+1634*x2*x6-4526*x2*x7+1116*x2*x8+7310*x3*x9+4558*x4^2+3698*x5, -266*x1*x2*x7-4140*x1*x3^2+7560*x1*x3*x6+912*x1*x3*x8+418*x1*x7*x8-114*x1*x8*x9+266*x10^2*x2+722*x2^2*x6+3230*x2*x3*x9+2014*x2*x4^2+4692*x3^3-2392*x3^2*x5-4336*x3^2*x6+4368*x3*x5*x6-11592*x3*x6^2-5520*x3*x7^2+2204*x5*x8^2+7056*x6^3+10080*x6*x7^2+722*x6*x8^2+1368*x10*x8+1634*x2*x5, 528*x1^2*x3+340*x1^2*x7-66*x1^2*x9-98*x1*x10^2+9882*x1*x2*x4-11826*x1*x2*x5-266*x1*x2*x6+10800*x1*x3*x7-1190*x1*x3*x9-742*x1*x4^2+1276*x1*x5*x8+418*x1*x6*x8-12240*x3^2*x7+6240*x3*x5*x7-11040*x3*x6*x7+3904*x4^3-4672*x4^2*x5+7442*x4^2*x7-17812*x4*x5*x7+2196*x4*x5*x8+10658*x5^2*x7-2628*x5^2*x8+10080*x6^2*x7+14400*x7^3+792*x1*x10-602*x1*x5+3782*x2*x4-4526*x2*x5, 2916*x1*x2*x5+2784*x1*x3*x5+912*x1*x3*x6+1276*x1*x5*x7-348*x1*x5*x9+418*x1*x6*x7-114*x1*x6*x9+1152*x4^2*x5+2196*x4*x5*x7-2628*x5^2*x7+7376*x5^2*x8+4408*x5*x6*x8+722*x6^2*x8+4176*x10*x5+1368*x10*x6+1116*x2*x5, -144*x1^2*x3-66*x1^2*x7+18*x1^2*x9-1190*x1*x3*x7-348*x1*x5*x8-114*x1*x6*x8+1190*x10^2*x3+3230*x2*x3*x6+14450*x3^2*x9+9010*x3*x4^2-216*x1*x10+7310*x3*x5, -196*x1*x10*x7+196*x10^3+532*x10*x2*x6+2380*x10*x3*x9+1484*x10*x4^2+1728*x1*x3+792*x1*x7-216*x1*x9+1204*x10*x5+4176*x5*x8+1368*x6*x8+2592*x10];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
const DROPOUT       = 0
const BATCH_SIZE    = 32
const SEQ_LENGTH    = 32
const DIM_HIDDEN    = 256
const DIM_EMBED     = 256
const LSTM_N_LAYER  = 2
const N_EPOCH       = 21
const BASE_LR       = 0.01
const WEIGHT_DECAY  = 0.00001
const CLIP_GRADIENT = 1
const NAME          = :ptb
const N_GPU         = 4
const USE_GPU       = true
const DATA_TR_RATIO = 0.9
const CKPOINT_PREFIX = joinpath(dirname(@__FILE__), "checkpoints/$NAME")

const BATCH_SIZE_SMP= 10
const SAMPLE_LENGTH = 100
const SAMPLE_START  = 'a'

const UNKNOWN_CHAR  = Char(0)
const INPUT_FILE    = joinpath(dirname(@__FILE__), "input.txt")
const VOCAB_FILE    = joinpath(dirname(@__FILE__), "vocab.dat")

ext: jl


text:
using Test: @test
import Trixi

# Use a macro to avoid world age issues when defining new initial conditions etc.
# inside an elixir.
"""
    @test_trixi_include(elixir; l2=nothing, linf=nothing,
                                atol=10*eps(), rtol=0.001,
                                parameters...)

Test Trixi by calling `trixi_include(elixir; parameters...)`.
By default, only the absence of error output is checked.
If `l2` or `linf` are specified, in addition the resulting L2/Linf errors
are compared approximately against these reference values, using `atol, rtol`
as absolute/relative tolerance.
"""
macro test_trixi_include(elixir, args...)

  local l2   = get_kwarg(args, :l2, nothing)
  local linf = get_kwarg(args, :linf, nothing)
  local atol = get_kwarg(args, :atol, 500*eps())
  local rtol = get_kwarg(args, :rtol, sqrt(eps()))
  local kwargs = Pair{Symbol, Any}[]
  for arg in args
    if arg.head == :(=) && !(arg.args[1] in (:l2, :linf, :atol, :rtol))
      push!(kwargs, Pair(arg.args...))
    end
  end

  quote
    Trixi.mpi_isroot() && println("═"^100)
    Trixi.mpi_isroot() && println($elixir)

    # evaluate examples in the scope of the module they're called from
    @test_nowarn_mod trixi_include(@__MODULE__, $elixir; $kwargs...)

    # if present, compare l2 and linf errors against reference values
    if !isnothing($l2) || !isnothing($linf)
      l2_measured, linf_measured = analysis_callback(sol)

      if Trixi.mpi_isroot() && !isnothing($l2)
        @test length($l2) == length(l2_measured)
        for (l2_expected, l2_actual) in zip($l2, l2_measured)
          @test isapprox(l2_expected, l2_actual, atol=$atol, rtol=$rtol)
        end
      end

      if Trixi.mpi_isroot() && !isnothing($linf)
        @test length($linf) == length(linf_measured)
        for (linf_expected, linf_actual) in zip($linf, linf_measured)
          @test isapprox(linf_expected, linf_actual, atol=$atol, rtol=$rtol)
        end
      end
    end

    Trixi.mpi_isroot() && println("═"^100)
    Trixi.mpi_isroot() && println("\n\n")
  end
end

# Get the first value assigned to `keyword` in `args` and return `default_value`
# if there are no assignments to `keyword` in `args`.
function get_kwarg(args, keyword, default_value)
  val = default_value
  for arg in args
    if arg.head == :(=) && arg.args[1] == keyword
      val = arg.args[2]
      break
    end
  end
  return val
end

# Modified version of `@test_nowarn` that prints the content of `stderr` when
# it is not empty. This is useful for debugging failing tests.
macro test_nowarn_debug(expr)
  quote
    let fname = tempname()
      try
        ret = open(fname, "w") do f
          redirect_stderr(f) do
            $(esc(expr))
          end
        end
        stderr_content = read(fname, String)
        if !isempty(stderr_content)
          println("Content of `stderr`:\n", stderr_content)
        end
        @test isempty(stderr_content)
        ret
      finally
        rm(fname, force=true)
      end
    end
  end
end

# Modified version of `@test_nowarn` that prints the content of `stderr` when
# it is not empty and ignnores module replacements.
macro test_nowarn_mod(expr)
  quote
    let fname = tempname()
      try
        ret = open(fname, "w") do f
          redirect_stderr(f) do
            $(esc(expr))
          end
        end
        stderr_content = read(fname, String)
        if !isempty(stderr_content)
          println("Content of `stderr`:\n", stderr_content)
        end
        # We also ignore simple module redefinitions for convenience. Thus, we
        # check whether every line of `stderr_content` is of the form of a
        # module replacement warning.
        @test occursin(r"^(WARNING: replacing module .+\.\n)*$", stderr_content)
        ret
      finally
        rm(fname, force=true)
      end
    end
  end
end


text:
"""
    NS_best(B::BnBModel)

Selects node with the lowest upper lower bound. Returns `(IntvBox,LBD,UBD,id,pos)`
where `Intv` is the interval box, `LBD` is the lower bound of the node, `UBD` is the
upper bound of the node, `id` is the id number of the node, and `pos` is the position
of the node in the BnB tree.
"""
function NS_best(B::BnBModel)
  tL,ind = findmin(B.LBD)
  return splice!(B.box,ind),splice!(B.LBD,ind),splice!(B.UBD,ind),splice!(B.id,ind),splice!(B.pos,ind)
end

"""
    NS_depth_breadth(B::BnBModel)

Selects node on the top of the stack. Returns `(IntvBox,LBD,UBD,id,pos)`
where `Intv` is the intervalbox, `LBD` is the lower bound of the node, `UBD` is the
upper bound of the node, `id` is the id number of the node, and `pos` is the position
of the node in the BnB tree.
"""
function NS_depth_breadth(B::BnBModel)
  return pop!(B.box),pop!(B.LBD),pop!(B.UBD),pop!(B.id),pop!(B.pos)
end


text:
using CollegeStratCollege, Test

mdl = CollegeStratCollege;

function study_time_grid_test(perCourse)
	@testset "Study Time grid" begin
        stg = mdl.make_test_study_time_grid(perCourse);
        @test mdl.validate_stg(stg);

        @test stime_grid_per_course(stg) == perCourse;
        n = n_study_times(stg);

        nTried = 3;
        sTimeGridV = study_time_grid(stg, nTried);
        @test all(sTimeGridV .> 0.0);
        @test length(sTimeGridV) == n;
        @test max_study_time(stg, n) ≈ maximum(sTimeGridV);
        sTimeGrid2V = study_time_from_grid(stg, stg.studyTimeV, nTried);
        @test isapprox(sTimeGridV, sTimeGrid2V);
	end
end

@testset "StudyTimeGrid" begin
    for perCourse in (true, false)
        study_time_grid_test(perCourse);
    end
end

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


text:
module LogGMM

import Distributions
using DocStringExtensions
import SpecialFunctions

import ..GMM

"""
$(TYPEDSIGNATURES)

Formula copied from wikipedia
https://en.wikipedia.org/wiki/Log-normal_distribution
"""
function lognormal_lpdf(
    x::Matrix{Float64}, mu::Vector{Float64}, sigma::Vector{Float64}
)::Matrix{Float64}
    if all(sigma .< 0)
        throw(DimensionMismatch("sigma: $(sigma) is below 0"))
    end

    sigma = transpose(max.(sigma, eps(Float64)))

    Z = x * sigma .* sqrt(2 * pi)
    if size(Z) != (length(x), length(sigma))
        throw(DimensionMismatch("size(Z): $(size(Z)) does not equal $((length(x), length(sigma)))"))
    end
    E = 1//2 .* ((log.(x) .- transpose(mu)) ./ sigma) .^ 2
    result = -E .- log.(Z)
    if size(result) != (length(x), length(sigma))
        throw(DimensionMismatch("size(result): $(size(result)) does not equal$((length(x), length(sigma)))"))
    end

    return result
end

"""
$(TYPEDSIGNATURES)

LGMM1 with low, high
"""
function LGMM1(
    weights::Vector{Float64}, mus::Vector{Float64}, sigmas::Vector{Float64},
    low::Float64, high::Float64, sample_size::Int
)::Matrix{Float64}
    if low > high
        throw(ArgumentError(string(
            "low (", low, ") should not be greater than high ", high
        )))
    end
    if !(length(weights) == length(mus) == length(sigmas))
        throw(DimensionMismatch(string(
            "length(weights): ", length(weights),
            " doesn't equal length(mus): ", length(mus),
            " nor length(sigmas): ", length(sigmas),
        )))
    end

    sample_size_t = (sample_size, 1)

    n_samples = Int(prod(sample_size_t))
    # -- draw from truncated components
    samples = Float64[]
    while length(samples) < n_samples
        active_real::Real = argmax(rand(Distributions.Multinomial(1, weights)))
        draw = rand(Distributions.Normal(mus[active_real], sigmas[active_real]))
        if low <= draw < high
            push!(samples, exp(draw))
        end
    end

    return reshape(samples, sample_size_t)
end
"""
$(TYPEDSIGNATURES)

LGMM1 without low, high
"""
function LGMM1(
    # weights must be Floats or Multinomial will fail
    weights::Vector{Float64}, mus::Vector{Float64}, sigmas::Vector{Float64},
    sample_size::Int
)::Matrix{Float64}
    if !(length(weights) == length(mus) == length(sigmas))
        throw(DimensionMismatch(string(
            "length(weights): ", length(weights),
            " doesn't equal length(mus): ", length(mus),
            " nor length(sigmas): ", length(sigmas),
        )))
    end

    sample_size_t = (sample_size, 1)

    n_samples = Int(prod(sample_size_t))
    active_arr::Array = getindex.(
        permutedims(argmax(rand(Distributions.Multinomial(1, weights), n_samples); dims = 1)),
        1
    )
    if length(active_arr) != n_samples
        throw(DimensionMismatch(string(
            "length(active_arr): ", length(active_arr),
            " doesn't equal n_samples: ", n_samples
        )))
    end
    samples = exp.(rand.(Distributions.Normal.(mus[active_arr], sigmas[active_arr])))

    return reshape(samples, sample_size_t)
end
"""
$(TYPEDSIGNATURES)

LGMM1 with low, high and q
"""
function LGMM1(
    weights::Vector{Float64}, mus::Vector{Float64}, sigmas::Vector{Float64},
    low::Float64, high::Float64, q::Float64, sample_size::Int
)::Matrix{Float64}
    samples = LGMM1(weights, mus, sigmas, low, high, sample_size)

    return round.(samples ./ q) .* q
end
"""
$(TYPEDSIGNATURES)

LGMM1 without low or high
"""
function LGMM1(
    weights::Vector{Float64}, mus::Vector{Float64}, sigmas::Vector{Float64}, q::Real,
    sample_size::Int
)::Matrix{Float64}
    samples = LGMM1(weights, mus, sigmas, sample_size)

    return round.(samples ./ q) .* q
end

# LGMM1_lpdf

"""
$(TYPEDSIGNATURES)

Implements CDF as given by `0.5 + 0.5 erf( log(x) - mu / sqrt(2 sigma^2))`

The maximum is used to move negative values and `0` up to a point where they do not cause
`NaN` or `Inf`, but also don't contribute much to the CDF.
"""
function lognormal_cdf(
    x::Matrix{<: Real}, mu::Vector{Float64}, sigma::Vector{Float64}
)::Matrix{Float64}
    if length(x) != length(mu) && length(x) != 1 && length(mu) != 1
        throw(DimensionMismatch("x: $(x) and mu $(mu) are not the same length or 1"))
    end

    top = log.(max.(x, eps(Float64))) .- mu
    bottom = max.(sqrt(2) .* sigma, eps(Float64))
    z = top ./ bottom
    return 0.5 .+ (0.5 .* SpecialFunctions.erf.(z))
end

function logprob(
    samples::Matrix{Float64}, weights::Vector{Float64}, mus::Vector{Float64},
    sigmas::Vector{Float64},
    low::Real, high::Real, q::Real, p_accept::Real
)::Matrix{Float64}
    prob = zeros(size(samples))
    for (w, mu, sigma) in zip(weights, mus, sigmas)
        ubound = min.(samples .+ (q / 2.0), exp(high))
        lbound = max.(samples .- (q / 2.0), exp(low))
        lbound = max.(0, lbound)
        # -- two-stage addition is slightly more numerically accurate
        inc_amt = w .* lognormal_cdf(ubound, [mu], [sigma])
        inc_amt .-= w .* lognormal_cdf(lbound, [mu], [sigma])
        prob += inc_amt
    end

    return log.(prob) .- log(p_accept)
end
function logprob(
    samples::Matrix{Float64}, weights::Vector{Float64}, mus::Vector{Float64},
    sigmas::Vector{Float64},
    q::Real, p_accept::Real
)::Matrix{Float64}
    prob = zeros(size(samples))
    for (w, mu, sigma) in zip(weights, mus, sigmas)
        ubound = samples .+ (q / 2.0)
        lbound = samples .- (q / 2.0)
        lbound = max.(0, lbound)
        # -- two-stage addition is slightly more numerically accurate
        inc_amt = w .* lognormal_cdf(ubound, [mu], [sigma])
        inc_amt .-= w .* lognormal_cdf(lbound, [mu], [sigma])
        prob += inc_amt
    end

    return log.(prob) .- log(p_accept)
end

"""
$(TYPEDSIGNATURES)

LGMM1_lpdf without low, high
(LGMM1_lpdf with low, high was identical)
"""
function LGMM1_lpdf(
    samples::Matrix{Float64}, weights::Vector{Float64}, mus::Vector{Float64},
    sigmas::Vector{Float64}
)::Vector{Float64}
    if !(length(weights) == length(mus) == length(sigmas))
        throw(DimensionMismatch(string(
            "length(weights): ", length(weights),
            " doesn't equal length(mus): ", length(mus),
            " nor length(sigmas): ", length(sigmas),
        )))
    end

    # compute the lpdf of each sample under each component
    lpdfs = lognormal_lpdf(samples, mus, sigmas)
    rval = GMM.logsum_rows(lpdfs .+ transpose(log.(weights)))
    if length(rval) != size(samples, 1)
        throw(DimensionMismatch(string(
            "Length of rval (", length(rval), ") does not match size of samples (",
            size(samples), ")"
        )))
    end

    return rval
end

function LGMM1_lpdf(
    samples::Matrix{Float64}, weights::Vector{Float64}, mus::Vector{Float64},
    sigmas::Vector{Float64}, low::Float64, high::Float64
)::Vector{Float64}
    if !(length(weights) == length(mus) == length(sigmas))
        throw(DimensionMismatch(string(
            "length(weights): ", length(weights),
            " doesn't equal length(mus): ", length(mus),
            " nor length(sigmas): ", length(sigmas),
        )))
    end

    # compute the lpdf of each sample under each component
    lpdfs = lognormal_lpdf(samples, mus, sigmas)
    rval = GMM.logsum_rows(lpdfs .+ transpose(log.(weights)))
    if length(rval) != size(samples, 1)
        throw(DimensionMismatch(string(
            "Length of rval (", length(rval), ") does not match size of samples (",
            size(samples), ")"
        )))
    end

    return rval
end
function LGMM1_lpdf(
    samples::Matrix{Float64}, weights::Vector{Float64}, mus::Vector{Float64},
    sigmas::Vector{Float64}, q::Float64
)::Matrix{Float64}
    if !(length(weights) == length(mus) == length(sigmas))
        throw(DimensionMismatch(string(
            "length(weights): ", length(weights),
            " doesn't equal length(mus): ", length(mus),
            " nor length(sigmas): ", length(sigmas),
        )))
    end

    p_accept = 1
    return logprob(samples, weights, mus, sigmas, q, p_accept)
end
function LGMM1_lpdf(
    samples::Matrix{Float64}, weights::Vector{Float64}, mus::Vector{Float64},
    sigmas::Vector{Float64}, low::Float64, high::Float64, q::Float64
)::Matrix{Float64}
    if !(length(weights) == length(mus) == length(sigmas))
        throw(DimensionMismatch(string(
            "length(weights): ", length(weights),
            " doesn't equal length(mus): ", length(mus),
            " nor length(sigmas): ", length(sigmas),
        )))
    end

    p_accept = sum(weights .* (GMM.normal_cdf([high], mus, sigmas) - GMM.normal_cdf([low], mus, sigmas)))
    return logprob(samples, weights, mus, sigmas, low, high, q, p_accept)
end

end # module LogGMM


text:
function iterate(camera::Camera, state=nothing)
    try
        return take!(camera), camera
    catch e
        if isa(e, InvalidStateException)
            return nothing
        else
            throw(e)
        end
    end
end

IteratorSize(camera::Camera) = Base.IsInfinite()


text:
module CausalModels

include("zennascausalmodel/CausalModels.jl")

end


text:
using Hecke

boundexp = ARGS[1]
startfield = ARGS[2]
number = ARGS[3]

sprint_formatted(fmt, args...) = @eval @sprintf($fmt, $(args...))

if length(ARGS) == 4
  basename = ARGS[4]
else
  basename = ""
end

boundexp = parse(Int, boundexp)
startfield = parse(Int, startfield)
number = parse(Int, number)

@assert iseven(boundexp)

boundquad = div(boundexp, 2)

lall = Hecke.quadratic_extensions(10^boundquad)

l = Hecke.quadratic_extensions(10^boundquad, u = startfield:(startfield + number - 1))

width = length(string(length(lall)))

for (i, K) in enumerate(l)
  z = Hecke.D5_extensions(fmpz(10)^boundexp, [K])
  if length(z) != 0
    fname = basename * sprint_formatted("%0$(width)d", startfield + i - 1)
    @show fname
    file = open(fname, "a")
    write(file, "# $(K.pol)\n")
    for (p, D) in z
      write(file, "($p, $D)\n")
    end
    close(file)
  end
end


text:
module NaiveNASflux

using Reexport
@reexport using NaiveNASlib
@reexport using Flux
import Flux.Tracker: hook
using Statistics
using Setfield
using LinearAlgebra
import InteractiveUtils: subtypes
import JuMP: @variable, @constraint

export FluxLayer, FluxParLayer, FluxDense, FluxRecurrent, FluxRnn, FluxGru, FluxLstm, FluxConvolutional, FluxConv, FluxConvTranspose, FluxDepthwiseConv, FluxCrossCor, FluxTransparentLayer, FluxParInvLayer, FluxDiagonal, FluxLayerNorm, FluxParNorm, FluxBatchNorm, FluxInstanceNorm, FluxGroupNorm, FluxNoParLayer

export mutable, concat, AbstractMutableComp, MutableLayer, LazyMutable, NoParams

export named, validated, logged

export ActivationContribution, neuron_value, Ewma

export indim, outdim, actdim, layer, layertype

export idmapping, idmapping_nowarn

export KernelSizeAligned, mutate_weights

include("types.jl")
include("util.jl")
include("select.jl")
include("mutable.jl")
include("vertex.jl")
include("pruning.jl")
include("weightinit.jl")

# Modified version of Flux.treelike for mutable structs which are mutable mainly because they are intended to be wrapped in MutationVertices which are not easy to create in the manner which Flux.mapchildren is designed.
function treelike_mutable(m::Module, T, fs = treelike_fields(T))
  @eval m begin
    Flux.children(x::$T) = ($([:(x.$f) for f in fs]...),)
    function Flux.mapchildren(f, x::$T)
        $([:(x.$fn = f(x.$fn)) for fn in fs]...)
        return x
    end
  end
end

treelike_fields(T) = fieldnames(T)

macro treelike_mutable(T, fs = nothing)
  fs == nothing || Meta.isexpr(fs, :tuple) || error("@treelike_mutable T (a, b)")
  fs = fs == nothing ? [] : [:($(map(QuoteNode, fs.args)...),)]
  :(treelike_mutable(@__MODULE__, $(esc(T)), $(fs...)))
end

for subtype in subtypes(AbstractMutableComp)
  @treelike_mutable subtype
end

end # module


text:
# Author: Lukasz Bratos

# Funkcja zwraca wartosc wyrazenia sin(x) + cos(3x) dla danego x
function f(x :: Float64)
    return sin(x) + cos(3 * x)
end

# Funkcja zwraca wartosc wyrazenia cos(x) - 3sin(3x) dla danego x (co jest pochodna wyrazenia sin(x) + cos(3x))
function derivativeOfF(x :: Float64)
    return cos(x) - 3 * sin(3 * x)
end

# Funkcja wylicza wartosc pochodnej funkcji sin(x) + cos(3x) dla danego x
function derivative(func, x0, h)
    return (func(x0 + h) - func(x0)) / h
end

# Wypisywanie wynikow w formacie pod LaTeXa
for i in 0:54
    println("\$ 2^{-$i} \$ & ", 1.0 + 2.0 ^ (-i), " & ", derivative(f, 1.0, 2.0 ^ (-i)), " & ", abs(derivativeOfF(1.0) - derivative(f, 1.0, 2.0 ^ (-i))), " \\\\\n\\hline")
end

println("f'(1.0) = ", derivativeOfF(1.0))


text:
# K=funLtoK(L)
# The function funLtoK(L) converts a kernel L matrix into a (normalized)
# kernel K matrix. The K matrix has to be semi-positive definite.
#
# Code available here:
# Keeler, 2018, https://github.com/hpaulkeeler/DetPoisson_Julia
#
# #TEST
# B=[3, 2, 1; 4, 5,6; 9, 8,7];
# L=B'*B
# L =
#
#   106    98    90
#    98    93    88
#    90    88    86
# K=funLtoK(L)
# K =
#
#     0.7602    0.3348   -0.0906
#     0.3348    0.3320    0.3293
#    -0.0906    0.3293    0.7492

function funLtoK(L)
    eigenVectLK=eigvecs(L); #eigen decomposition -- vectors
    eigenValL=(eigvals(L)); #eigen decomposition -- values
    eigenValK = eigenValL./(1 .+eigenValL); #eigenvalues of K
    eigenValK=Diagonal(eigenValK); #eigenvalues of L as diagonal matrix
    K=eigenVectLK*eigenValK*(eigenVectLK'); #recombine from eigen components
    K=real(K); #make sure all values are real
    return K;
end


text:
"""
    Ground(; name)

Ground node with the potential of zero and connector `g`. Every circuit must have one ground
node.

# Connectors:
- `g`
"""
function Ground(;name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    ODESystem(eqs, t, [], []; systems=[g], name=name)
end

"""
    Resistor(; name, R)

Creates an ideal Resistor following Ohm's Law.

# States:
See [OnePort](@ref)

# Connectors:
- `p` Positive pin
- `n` Negative pin

# Parameters:
- `R`: [`Ω`] Resistance
"""
function Resistor(;name, R)
    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters R=R
    eqs = [
        v ~ i * R
    ]
    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
    Conductor(;name, G)

Creates an ideal conductor.

# States:
See [OnePort](@ref)

# Connectors:
- `p` Positive pin
- `n` Negative pin

# Parameters:
- `G`: [`S`] Conductance
"""
function Conductor(;name, G)
    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters G=G
    eqs = [
        i ~ v * G
    ]
    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
    Capacitor(; name, C)

Creates an ideal capacitor.

# States:
- `v(t)`: [`V`]
  The voltage across the capacitor, given by `D(v) ~ p.i / C`

# Connectors:
- `p` Positive pin
- `n` Negative pin

# Parameters:
- `C`: [`F`] Capacitance
- `v_start`: [`V`] Initial voltage of capacitor
"""
function Capacitor(;name, C, v_start=0.0)
    @named oneport = OnePort(;v_start=v_start)
    @unpack v, i = oneport
    pars = @parameters C=C
    eqs = [
        D(v) ~ i / C
    ]
    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
    Inductor(; name, L)

Creates an ideal Inductor.

# States:
See [OnePort](@ref)

# Connectors:
- `p` Positive pin
- `n` Negative pin

# Parameters:
- `L`: [`H`] Inductance
- `i_start`: [`A`] Initial current through inductor
"""
function Inductor(;name, L, i_start=0.0)
    @named oneport = OnePort(;i_start=i_start)
    @unpack v, i = oneport
    pars = @parameters L=L
    eqs = [
        D(i) ~ 1 / L * v
    ]
    extend(ODESystem(eqs, t, [], pars; name=name), oneport)
end

"""
    IdealOpAmp(; name)

Ideal operational amplifier (norator-nullator pair).
The ideal OpAmp is a two-port. The left port is fixed to `v1 = 0` and `i1 = 0` (nullator).
At the right port both any voltage `v2` and any current `i2` are possible (norator).

# States:
See [TwoPort](@ref)

# Connectors:
- `p1` Positive pin (left port)
- `p2` Positive pin (right port)
- `n1` Negative pin (left port)
- `n2` Negative pin (right port)
"""
function IdealOpAmp(;name)
    @named twoport = TwoPort()
    @unpack v1, v2, i1, i2 = twoport

    eqs = [
        v1 ~ 0
        i1 ~ 0
    ]
    extend(ODESystem(eqs, t, [], [], name=name), twoport)
end

"""
    Short(; name)

Short is a simple short cut branch. That means the voltage drop between both pins is zero.

# States:
See [OnePort](@ref)

# Connectors:
- `p` Positive pin
- `n` Negative pin
"""
function Short(; name)
    @named oneport = OnePort()
    @unpack v, i = oneport
    eqs = [v ~ 0]
    extend(ODESystem(eqs, t, [], []; name=name), oneport)
end


text:
using OrdinaryDiffEq, DiffEqDevTools, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear

prob = prob_ode_linear
sol =solve(prob,Rosenbrock32())
dt₀ = sol.t[2]

prob = prob_ode_2Dlinear
sol =solve(prob,ExplicitRK(tableau=constructBogakiShampine3()))
dt₀ = sol.t[2]

@test  1e-7 < dt₀ < .1
@test_throws ErrorException sol = solve(prob,Euler())
#dt₀ = sol.t[2]

sol3 =solve(prob,ExplicitRK(tableau=constructDormandPrince8_64bit()))
dt₀ = sol3.t[2]

@test 1e-7 < dt₀ < .3

T = Float32
u0 = T.([1.0;0.0;0.0])

tspan = T.((0,100))
prob = remake(prob, u0=u0, tspan=tspan)
@test_nowarn solve(prob, Euler(); dt=T(0.0001))

tspan = T.((2000,2100))
prob = remake(prob, tspan=tspan)
@test_throws ArgumentError solve(prob, Euler(); dt=T(0.0001)) # Loops forever

function rober(du,u,p,t)
  y₁,y₂,y₃ = u
  k₁,k₂,k₃ = p
  du[1] = -k₁*y₁+k₃*y₂*y₃
  du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
  du[3] =  k₂*y₂^2
  nothing
end
u0 = Float32[1.0,0.0,0.0]
tspan = (0f0, 1f5)
params = (4f-2,3f7,1f4)
prob = ODEProblem(rober,u0,tspan,params)
sol = solve(prob, Rosenbrock23())

# https://github.com/SciML/DifferentialEquations.jl/issues/743

using LinearAlgebra
function f(du, u, p, t)
    du[1] = -p[1]*u[1] + p[2]*u[2]*u[3]
    du[2] = p[1]*u[1] - p[2]*u[2]*u[3] - p[3]*u[2]*u[2]
    du[3] = u[1] + u[2] + u[3] - 1.
end
M = Diagonal([1,1,0])
p = [0.04, 10^4, 3e7]
u0 = [1.,0.,0.]
tspan = (0., 1e6)
prob = ODEProblem(ODEFunction(f, mass_matrix = M), u0, tspan, p)
sol = solve(prob, Rodas5())
@test sol.t[end] == 1e6


text:
#############
# Constants #
#############

"""
Defines valid time systems usable as inputs in the time module, and in functions
pertaining to the `Epoch` class.

Valid systems are: `"GPS"`, `"TAI"`, `"TT"`, `"UTC"`, and `"UT1"`
"""
VALID_TIME_SYSTEMS = ["GPS", "TAI", "TT", "UTC", "UT1"]

####################
# Helper Functions #
####################

string_to_nanoseconds(str::AbstractString) = parse(Float64, str)*10^(9-length(str))

function align_epoch_data(days::Integer, seconds::Integer, nanoseconds::Real)
    # Align nanoseconds to [0, 1.0e9)
    while nanoseconds < 0.0
        nanoseconds += 1.0e9
        seconds     -= 1
    end

    while nanoseconds >= 1.0e9
        nanoseconds -= 1.0e9
        seconds     += 1
    end

    # Align seconds to [0, 86400)
    while seconds < 0
        seconds += 86400
        days    -= 1
    end

    while seconds >= 86400
        seconds -= 86400
        days    += 1
    end

    return days, seconds, nanoseconds
end

################
# Time Methods #
################

export caldate_to_mjd
"""
Convert a Gregorian calendar date to the equivalent Modified Julian Date representation of that time instant.

# Aguments:
- `year::Integer` Year
- `year::Integer` Month
- `year::Integer` Day
- `hour::Integer` Hour
- `minute::Integer` Minute
- `second::Real` Seconds
- `nanoseconds::Real` Nanoseconds

Returns:
- `mjd::Float64` Modified Julian Date of Epoch
"""
function caldate_to_mjd(year::Integer, month::Integer, day::Integer, hour::Integer=0, minute::Integer=0, second::Real=0.0, nanoseconds::Real=0.0)
    status, jd, fd = iauDtf2d("TAI", year, month, day, hour, minute, second + nanoseconds/1.0e9)

    mjd = (jd - MJD_ZERO) + fd

    return mjd
end

export mjd_to_caldate
"""
Convert a Modified Julian Date to the equivalent Gregorian calendar date representation of the same instant in time.

# Aguments:
- `mjd::Real`: Modified Julian Date of Epoch

Returns:
- `year::Int`: Year
- `year::Int`: Month
- `year::Int`: Day
- `hour::Int`: Hour
- `minute::Int`: Minute
- `second::Float64`: Seconds
- `nanoseconds::Float64`: Nanoseconds
"""
function mjd_to_caldate(mjd::Real)
    status, iy, im, id, ihmsf = iauD2dtf("TAI", 9, MJD_ZERO, mjd)

    return iy, im, id, ihmsf[1], ihmsf[2], ihmsf[3], ihmsf[4]/1.0e9
end

export caldate_to_jd
"""
Convert a Gregorian calendar date to the equivalent Julian Date representation of that time instant.

# Aguments:
- `year::Integer`: Year
- `year::Integer`: Month
- `year::Integer`: Day
- `hour::Integer`: Hour
- `minute::Integer`: Minute
- `second::Real`: Seconds
- `nanoseconds::Real`: Nanoseconds

Returns:
- `mjd::Float64`: Julian Date of Epoch
"""
function caldate_to_jd(year::Integer, month::Integer, day::Integer, hour::Integer=0, minute::Integer=0, second::Real=0.0, nanoseconds::Real=0.0)
    status, jd, fd = iauDtf2d("TAI", year, month, day, hour, minute, second + nanoseconds/1.0e9)

    jd = jd + fd

    return jd
end

export jd_to_caldate
"""
Convert a Julian Date to the equivalent Gregorian calendar date representation of the same instant in time.

# Aguments:
- `jd::Real`: Julian Date of Epoch

Returns:
- `year::Int`: Year
- `year::Int`: Month
- `year::Int`: Day
- `hour::Int`: Hour
- `minute::Int`: Minute
- `second::Float64`: Seconds
- `microsecond::Float64`: Nanoseconds
"""
function jd_to_caldate(jd::Real)
    status, iy, im, id, ihmsf = iauD2dtf("TAI", 9, jd, 0)

    return iy, im, id, ihmsf[1], ihmsf[2], ihmsf[3], ihmsf[4]/1.0e9
end

export elapsed_from_epoch
"""
Compute the number of elapsed seconds since a given Epoch from the day number. Can be used to compute the elapsed time since a given Julian or Modified Julian Date.

Arguments:
- `day_number::Real`: Day number, can contain fractional days. Asummes all days are a uniform 86400.0 seconds in length.
- `day_epoch::Real`: Modified Julian Date of Epoch

Returns:
- `t::Float`: Number of elapsed seconds between the Provided Modified
    Julian date and the epoch.
"""
function elapsed_from_epoch(day_number::Real, day_epoch=0::Real)
    return (day_number - day_epoch)*86400.0
end

export days_from_elapsed
"""
Computes the day number in a given time scale given the elapsed time since epoch and the epoch itself.

Assumes all days are counted using a uniform 86400.0 seconds over the time span.

Arguments:
- `t::Real`: Elapsed seconds since the `day_epoch`.
- `day_epoch::Real`: Day number of the epoch. Common values are `SatelliteDynamics.MJD_ZERO` (to get the Julian Day number) or `SatelliteDynamics.MJD2000` (to get Modified Julian Days if reckoning time from January 1, 2000 0H)

Returns:
- `days::Float`: Number of elapsed days in the time scale.
"""
function days_from_elapsed(t::Real, day_epoch=0::Real)
    return (t/86400.0 + day_epoch)
end

###############
# Epoch Class #
###############

"""
Returns whether a symbol is a valid time system.

Valid systems are: `"GPS"`, `"TAI"`, `"TT"`, `"UTC"`, and `"UT1"`

Arguments:
- `tsys::String`: time system symbol to check for validity.

Returns:
- `valid::Bool`: Return `true` if tsys is a valid time system.
"""
function valid_time_system(tsys::String)
    return tsys in VALID_TIME_SYSTEMS
end

export Epoch
"""
The `Epoch` type represents a single instant in time. It is used throughout the
SatelliteDynamics module. It is meant to provide a clear definition of moments
in time and provide a convenient interface display time in various representations
as well as in differrent time systems. The internal data members are also chosen
such that the representation maintains nanosecond-precision in reprersenation
of time and doesn't accumulate floating-point arithmetic errors larger than
nanoseconds even after centuries.

Supports `+`, `+=`, `-`, and `-=` operators. Two Epoch's can be differenced to
return the time difference between two Epochs. If adding a `Real` number it is
interpreted as an offset in seconds to add to the Epoch.

The class also supports all arithmetic operators: `==`, `!=`, `<`, `<=`, `>`, `>=`

Arguments:
- `year::Int` Year
- `year::Int` Month
- `year::Int` Day
- `hour::Int` Hour (optional)
- `minute::Int` Minute (optional)
- `second::Real` Seconds (optional)
- `nanoseconds::Real` Nanoseconds (optional)
- `tsys::String`: Time system of the epoch at initialization

The Epoch class can be also be initialized from a string. Examples of Valid String constructors are:

```julia
epc = Epoch("2018-12-20")
epc = Epoch("2018-12-20T16:22:19.0Z")
epc = Epoch("2018-12-20T16:22:19.123Z")
epc = Epoch("2018-12-20T16:22:19.123456789Z")
epc = Epoch("2018-12-20T16:22:19Z")
epc = Epoch("20181220T162219Z")
epc = Epoch("2018-12-01 16:22:19 GPS")
epc = Epoch("2018-12-01 16:22:19.0 GPS")
epc = Epoch("2018-12-01 16:22:19.123 GPS")
epc = Epoch("2018-12-01 16:22:19.123456789 GPS")
```
"""
struct Epoch
    # All days, seconds, and nanoseconds are stored internally in the TAI time scale, conversion to from TAI is done on intpu/output interacting.
    days::Int # Total days [0, ∞)
    seconds::Int # Integer seconds [0, 86400)
    nanoseconds::Float64 # Fractional seconds [0, 1)
    tsys::String # Time system of epoch
end

function Epoch(year::Real, month::Real, day::Real, hour::Real=0, minute::Real=0,
               second::Real=0, nanosecond::Real=0.0; tsys::String="TAI")
    if !valid_time_system(tsys)
        error("Invalid time system $(String(tsys))")
    end

    # Convert date into days and fractional days in desired time system
    status, jd, fd = iauDtf2d("TAI", year, month, day, hour, minute, 0)

    # Log initial output
    # @debug "iauDtf2d returned ($status, $jd, $fd)"

    if status != 0
        error("Non-zero SOFA status returned initializing Epoch.")
    end

    # Get time system offset based on days and fractional days using SOFA
    tsys_offset     = time_system_offset(jd, fd, tsys, "TAI")
    foffset, offset = modf(tsys_offset)

    # Ensure days is an integer number, add entire fractional component to the
    # fractional days variable
    fdays, days = modf(jd)
    fd          = fd + fdays

    # Convert fractional days into total seconds still retaining fractional part
    seconds = fd * 86400.0
    fsecs, secs = modf(seconds)

    # @debug "seconds: $seconds"
    # @debug "fsecs: $fsecs, secs: $secs"

    # Now trip second of the fractional part
    fsecond, second = modf(second)
    # @debug "fsecond: $fsecond, second: $second"

    # Add the integer parts together
    seconds = secs + second + offset
    # @debug "seconds: $seconds"

    # Convert the fractional parts to nanoseconds
    nanoseconds = nanosecond + fsecs*1e9 + fsecond*1e9 + foffset*1e9
    nanoseconds = Float64(nanoseconds)
    # @debug "nanoseconds: $nanoseconds"

    # @debug "Initializing Epoch: days: $days, seconds: $seconds, nanoseconds: $nanoseconds"

    # Ensure type consistency of input:
    days        = Int(days)
    seconds     = Int(seconds)
    nanoseconds = Float64(nanoseconds)

    days, seconds, nanoseconds = align_epoch_data(days, seconds, nanoseconds)
    # @debug "Initializing Epoch: days: $days, seconds: $seconds, nanoseconds: $nanoseconds"

    # Parse inputs into Epoch internal storage
    Epoch(days, seconds, nanoseconds, tsys)
end

# All REGEX used for epoch initialization MUST BE AN EXACT MATCH.
# Otherwise the regex parse improperly match on a partial fit and fail.
VALID_EPOCH_REGEX = [
    r"^(\d{4})\-(\d{2})\-(\d{2})$",
    # r"^(\d{4})\-(\d{2})\-(\d{2})[T](\d{2})\:(\d{2})\:(\d{2})([+-])(\d{2})\:(\d{2})$",
    r"^(\d{4})\-(\d{2})\-(\d{2})[T](\d{2})\:(\d{2})\:(\d{2})[Z]$",
    r"^(\d{4})\-(\d{2})\-(\d{2})[T](\d{2})\:(\d{2})\:(\d{2})[.](\d*)[Z]$",
    r"^(\d{4})(\d{2})(\d{2})[T](\d{2})(\d{2})(\d{2})[Z]$",
    r"^(\d{4})\-(\d{2})\-(\d{2})\s(\d{2})\:(\d{2})\:(\d{2})\.*\s*(\d*)\s*([A-Z]*)$"
]

function Epoch(str::String)
    year        = 0
    month       = 0
    day         = 0
    hour        = 0
    minute      = 0
    second      = 0.0
    nanoseconds = 0.0
    tsys        = "UTC"

    m = nothing
    # Iterate through valid regex string
    for regex in VALID_EPOCH_REGEX
        m = match(regex, str)
        if !(m === nothing)
            # @debug m
            # Parse date (common to all)
            year  = parse(Int, m[1])
            month = parse(Int, m[2])
            day   = parse(Int, m[3])

            # Parse time (most have this)
            if length(m.captures) >= 6
                hour   = parse(Int, m[4])
                minute = parse(Int, m[5])
                second = parse(Float64, m[6])
            end

            # Parse additional types
            if length(m.captures) == 7
                nanoseconds = string_to_nanoseconds(m[7])
            elseif length(m.captures) == 8
                if m[7] != ""
                    nanoseconds = string_to_nanoseconds(m[7])
                end
                tsys = string(m[8])

                if !(tsys in VALID_TIME_SYSTEMS)
                    throw(ErrorException("Parsed invalid time system: \"$tsys\""))
                end
            end

            # Exit early since a match has been found
            break
        end
    end

    # No valid match found throw error
    if m === nothing
        error("Invalid Epoch string. Must be iso8061 compliant.")
    end

    # @debug "regex parsed: $year $month $day $hour $minute $second $nanoseconds $(string(tsys))"

    return Epoch(year, month, day, hour, minute, second, nanoseconds, tsys=tsys)
end

####################
# Epoch Arithmetic #
####################

function Base.:+(epc::Epoch, t::Real)
    # Immidiately separate seconds and fractional seconds
    fseconds, seconds = modf(t)
    seconds = Int(seconds) # Conver seconds to integer seconds

    # Compute time delta aligned
    dt_days        = div(seconds, 86400)
    dt_seconds     = seconds % 86400
    dt_nanoseconds = fseconds*1e9

    # Perform additon to get new epoch
    days        = epc.days + dt_days
    seconds     = epc.seconds + dt_seconds
    nanoseconds = epc.nanoseconds + dt_nanoseconds

    # Align to proper ranges in reverse order
    days, seconds, nanoseconds = align_epoch_data(days, seconds, nanoseconds)

    return Epoch(days, seconds, nanoseconds, epc.tsys)
end

function Base.:+(t::Real, epc::Epoch)
    return epc + t
end

function Base.:-(epc::Epoch, t::Real)
    return epc + (-t)
end

function Base.:-(epc_left::Epoch, epc_right::Epoch)
    # Return difference in time between two epochs
    return (epc_left.days - epc_right.days)*86400.0 + (epc_left.seconds - epc_right.seconds) + (epc_left.nanoseconds - epc_right.nanoseconds)/1e9
end

##########
# Output #
##########

function Base.show(io::IO, epc::Epoch)
    year, month, day, hour, minute, second, nanodseconds = caldate(epc, tsys="UTC")

    s = @sprintf "Epoch(%02d-%02d-%02dT%02d:%02d:%02d.%03dZ)" year month day hour minute second floor(nanodseconds/1.0e6);

    print(io, s)
end

export epoch_to_jdfd
"""
Compute the two-part date format used by SOFA.jl functions forr a given Epoch.

Arguments:
- `epc::Epoch`: Epoch
- `tsys::String`: Time system to return output in

Returns:
- `d1::Real`: First part of two part date. [days]
- `d2::Real`: Second part of two part date. [days]
"""
function epoch_to_jdfd(epc::Epoch; tsys::String=epc.tsys)
    offset = time_system_offset(epc, "TAI", tsys)

    return epc.days, (epc.seconds + offset + epc.nanoseconds/1.0e9)/86400.0
end

export caldate
"""
Return the Gregorian calendar date for a specific

Arguments:
- `epc::Epoch`: Input epoch
- `tsys::String`: Time system to compute output in.

Returns:
- `year::Int`: Year of epoch
- `month::Int`: Month of epoch
- `day::Int`: Day of epoch
- `hour::Int`: Hour of epoch
- `minute::Int`: Minute of epoch
- `second::Int`: Second of epoch
- `nanoseconds::Int`: Year of epoch
"""
function caldate(epc::Epoch; tsys::String=epc.tsys)
    offset = time_system_offset(epc, "TAI", tsys)

    status, iy, im, id, ihmsf = iauD2dtf("TAI", 9,
                                    epc.days,
                                    (epc.seconds + offset + epc.nanoseconds/1.0e9)/86400.0)

    return iy, im, id, ihmsf[1], ihmsf[2], ihmsf[3], ihmsf[4]
end

export jd
"""
Compute the Julian Date for a specific epoch

Arguments:
- `epc::Epoch`: Epoch
- `tsys::String`: Time system to return output in

Returns:
- `jd::Real`: Julian date of the epoch in the requested time system
"""
function jd(epc::Epoch; tsys::String=epc.tsys)
    offset = time_system_offset(epc, "TAI", tsys)

    return (epc.days + (epc.seconds + epc.nanoseconds/1.0e9 + offset)/86400.0)
end

export mjd
"""
Compute the Modified Julian Date for a specific epoch

Arguments:
- `epc::Epoch`: Epoch
- `tsys::String`: Time system to return output in

Returns:
- `mjd::Real`: Julian date of the epoch in the requested time system
"""
function mjd(epc::Epoch; tsys::String=epc.tsys)
    offset = time_system_offset(epc, "TAI", tsys)

    return (epc.days + (epc.seconds + epc.nanoseconds/1.0e9 + offset)/86400.0) - MJD_ZERO
end

export day_of_year
"""
Return the day-of-year number for a given `Epoch`.

January 1 0h of each year will return 1.

Arguments:
- `epc::Epoch`: Epoch
- `tsys::String`: Time system to return output in

Returns:
- `doy::Real`: Day of year number.
"""
function day_of_year(epc::Epoch; tsys::String=epc.tsys)
    # Compute MJD of first day of yearr
    year, month, day, hour, minute, second, = caldate(epc, tsys=tsys)
    mjd0 = caldate_to_mjd(year, 1, 1)

    # Compute MJD of current day
    offset = time_system_offset(epc, "TAI", tsys)
    mjd = (epc.days + (epc.seconds + epc.nanoseconds/1.0e9 + offset)/86400.0) - MJD_ZERO

    # Get day of year is the difference
    doy = mjd - mjd0 + 1.0

    return doy
end

export gmst
"""
Compute the Greenwich Mean Sidereal Time for the given Epoch.

Arguments:
- `epc::Epoch`: Epoch
- `use_degrees::Bool`: Return output in degrees (Default: false)

Returns:
- `gmst::Real`: Greenwich Mean Sidereal Time [rad/deg]
"""
function gmst(epc::Epoch; use_degrees::Bool=false)
    uta, utb = epoch_to_jdfd(epc, tsys="UT1")
    tta, ttb = epoch_to_jdfd(epc, tsys="TT")

    gmst = iauGmst06(uta, utb, tta, ttb)

    return use_degrees ? gmst*180.0/pi : gmst
end

export gast
"""
Compute the Greenwich Mean Sidereal Time for the given Epoch.

Arguments:
- `epc::Epoch`: Epoch
- `use_degrees::Bool`: Return output in degrees (Default: false)

Returns:
- `gast::Real`: Greenwich Apparent Sidereal Time [rad/deg]
"""
function gast(epc::Epoch; use_degrees::Bool=false)
    uta, utb = epoch_to_jdfd(epc, tsys="UT1")
    tta, ttb = epoch_to_jdfd(epc, tsys="TT")

    gast = iauGst06a(uta, utb, tta, ttb)

    return use_degrees ? gast*180.0/pi : gast
end

# Hash function
function Base.hash(epc::Epoch, h::UInt)
    return hash(epc.days, hash(epc.seconds, hash(epc.nanoseconds, hash(:Epoch, h))))
    # return hash(epc.days, hash(epc.seconds, hash(:Epoch, h)))
end

# Comparison operators
function Base.:(==)(epc_left::Epoch, epc_right::Epoch)
    return ((epc_left.days == epc_right.days) &&
            (epc_left.seconds == epc_right.seconds) &&
            isapprox(epc_left.nanoseconds, epc_right.nanoseconds, atol=1e-3))
end

function Base.:(!=)(epc_left::Epoch, epc_right::Epoch)
    return ((epc_left.days != epc_right.days) ||
            (epc_left.seconds != epc_right.seconds) ||
            !isapprox(epc_left.nanoseconds, epc_right.nanoseconds, atol=1e-3))
end

function Base.:(<)(epc_left::Epoch, epc_right::Epoch)
    return ((epc_left.days < epc_right.days) ||
            ((epc_left.days == epc_right.days) &&
             (epc_left.seconds < epc_right.seconds)) ||
             ((epc_left.days == epc_right.days) &&
             (epc_left.seconds == epc_right.seconds) &&
             (epc_left.nanoseconds < epc_right.nanoseconds)))
end

Base.isless(epc_left::Epoch, epc_right::Epoch) = Base.:(<)(epc_left::Epoch, epc_right::Epoch)

function Base.:(<=)(epc_left::Epoch, epc_right::Epoch)
    return ((epc_left.days < epc_right.days) ||
            ((epc_left.days == epc_right.days) &&
            (epc_left.seconds < epc_right.seconds)) ||
            ((epc_left.days == epc_right.days) &&
            (epc_left.seconds == epc_right.seconds) &&
            (epc_left.nanoseconds <= epc_right.nanoseconds)))
end

function Base.:(>)(epc_left::Epoch, epc_right::Epoch)
    return ((epc_left.days > epc_right.days) ||
            ((epc_left.days == epc_right.days) &&
             (epc_left.seconds > epc_right.seconds)) ||
             ((epc_left.days == epc_right.days) &&
             (epc_left.seconds == epc_right.seconds) &&
             (epc_left.nanoseconds > epc_right.nanoseconds)))

end

Base.isgreater(epc_left::Epoch, epc_right::Epoch) = Base.:(>)(epc_left::Epoch, epc_right::Epoch)

function Base.:(>=)(epc_left::Epoch, epc_right::Epoch)
    return ((epc_left.days > epc_right.days) ||
            ((epc_left.days == epc_right.days) &&
            (epc_left.seconds > epc_right.seconds)) ||
            ((epc_left.days == epc_right.days) &&
            (epc_left.seconds == epc_right.seconds) &&
            (epc_left.nanoseconds >= epc_right.nanoseconds)))
end

#######################
# Time System Offsets #
#######################

export time_system_offset
"""
Compute the offset between two time systems at a given Epoch.

The offset (in seconds) is computed as:

    time_system_offset = tsys_dest - tsys_src

The value returned is the number of seconds that musted be added to the source time system given the input epoch, to get the equivalent epoch.

Conversions are accomplished using SOFA C library calls.
Epoch.

Arguments:
- `jd::Real`: Part 1 of two-part date (Julian days)
- `fd::Real`: Part 2 of two-part date (Fractional days)
- `tsys_src::String`: Base time system
- `tsys_dest::String`: Destination time system

Returns:
- `offset::Float`: Offset between soruce and destination time systems in seconds.
"""
function time_system_offset(jd, fd, tsys_src::String, tsys_dest::String)
    # If no transformation is needed needed return early
    if tsys_src == tsys_dest
        return 0.0
    end

    offset = 0.0

    # Convert To TAI
    if tsys_src == "GPS"
        offset += TAI_GPS
    elseif tsys_src == "TT"
        offset += TAI_TT
    elseif tsys_src == "UTC"
        status, iy, im, id, ihmsf = iauD2dtf("UTC", 6, jd, fd) # Returns TAI-UTC
        status, dutc = iauDat(iy, im, id, (ihmsf[1]*3600 + ihmsf[2]*60 + ihmsf[3] + ihmsf[4]/1e6)/86400.0)
        offset += dutc
    elseif tsys_src == "UT1"
        # Convert UT1 -> UTC
        offset -= UT1_UTC((jd - MJD_ZERO) + fd)

        # Convert UTC -> TAI
        status, iy, im, id, ihmsf = iauD2dtf("UTC", 6, jd, fd + offset) # Returns TAI-UTC
        status, dutc = iauDat(iy, im, id, (ihmsf[1]*3600 + ihmsf[2]*60 + ihmsf[3] + ihmsf[4]/1e6)/86400.0)
        offset += dutc
    elseif tsys_src == "TAI"
        # Do nothing in this case
    end

    # Covert from TAI to source
    if tsys_dest == "GPS"
        offset += GPS_TAI
    elseif tsys_dest == "TT"
        offset += TT_TAI
    elseif tsys_dest == "UTC"
        # Initial UTC guess
        u1, u2 = jd, fd + offset/86400.0

        # Iterate to get the UTC time
        for i in 1:3
            status, d1, d2 = iauUtctai(u1, u2)

            # Adjust UTC guess
            u1 += jd - d1
            u2 += fd + offset/86400.0 - d2
        end

        # Compute Caldate from two-part date
        status, iy, im, id, ihmsf = iauD2dtf("UTC", 6, u1, u2)

        status, dutc = iauDat(iy, im, id, (ihmsf[1]*3600 + ihmsf[2]*60 + ihmsf[3] + ihmsf[4]/1e6)/86400.0)
        offset -= dutc
    elseif tsys_dest == "UT1"
        # Initial UTC guess
        u1, u2 = jd, fd + offset/86400.0

        # Iterate to get the UTC time
        for i in 1:3
            status, d1, d2 = iauUtctai(u1, u2)

            # Adjust UTC guess
            u1 += jd - d1
            u2 += fd + offset/86400.0 - d2
        end

        # Compute Caldate from two-part date
        status, iy, im, id, ihmsf = iauD2dtf("UTC", 6, u1, u2)

        status, dutc = iauDat(iy, im, id, (ihmsf[1]*3600 + ihmsf[2]*60 + ihmsf[3] + ihmsf[4]/1e6)/86400.0)
        offset -= dutc

        # Convert UTC to UT1
        offset += UT1_UTC(u1 + u2 + offset/86400.0 - MJD_ZERO)
    elseif tsys_dest == "TAI"
        # Do nothing in this case
    end

    return offset
end

function time_system_offset(epc::Epoch, tsys_src::String, tsys_dest::String)
    jd = epc.days
    fd = (epc.seconds + epc.nanoseconds/1.0e9)/86400.0
    return time_system_offset(jd, fd, tsys_src, tsys_dest)
end


text:
using Test

import ActinRingsMC as armc

function main()
    energy_test()
end

function energy_test()
    ks = 7.82e-9
    kd = 7.82e-9/5
    T = 300
    delta = 36e-9
    Xc = 12e-9
    EI = 3.6e-26
    lf = 84
    Nfil = 6
    Nsca = 4

    Lf = armc.calc_Lf(lf, delta)

    # Test energy functions
    max_height = armc.calc_max_lattice_height(Nsca, lf)
    min_height = armc.calc_min_lattice_height(Nsca, lf)
    overlap = 21
    start_height = max_height - Nsca*overlap
    start_radius = armc.calc_radius(delta, start_height)

    lattice = armc.Lattice(start_height, max_height, min_height)
    filaments = armc.generate_starting_config(lattice, Nfil, Nsca, lf, overlap)
    armc.update_occupancies!(filaments, lattice)
    sysparms = armc.SystemParams(ks, kd, T, delta, Xc, EI, Lf, lf, Nfil, Nsca)
    system = armc.System(sysparms, filaments, start_radius)

    @test isapprox(armc.overlap_energy(system, lattice), -5.729146253430745e-19)
    @test isapprox(armc.total_bending_energy(system), 1.617547027233688e-19)
    @test isapprox(armc.total_energy(system, lattice), -4.111599226197057e-19)
end

ks = 7.82e-9
kd = 7.82e-9/5
T = 300
delta = 36e-9
Xc = 12e-9
EI = 3.6e-26
lf = 84
Nfil = 6
Nsca = 4

Lf = armc.calc_Lf(lf, delta)

max_height = armc.calc_max_lattice_height(Nsca, lf)
min_height = armc.calc_min_lattice_height(Nsca, lf)
overlap = 1
start_height = max_height - Nsca*(overlap - 1)
start_radius = armc.calc_radius(delta, start_height)

lattice = armc.Lattice(start_height, max_height, min_height)
filaments = armc.generate_starting_config(lattice, Nfil, Nsca, lf, overlap)
armc.update_occupancies!(filaments, lattice)
sysparms = armc.SystemParams(ks, kd, T, delta, Xc, EI, Lf, lf, Nfil, Nsca)
system = armc.System(sysparms, filaments, start_radius)

mkpath("outs")
file = armc.prepare_vtf_file("outs/test.vtf", system)
armc.write_vtf(system, file)
close(file)

iters = 2
steps = 1e4
max_bias_diff = 5
write_interval = 1e3
radius_move_freq = 0.5
filebase = "outs/test"
analytical_biases = true
read_biases = false
biases_filename = ""
restart_iter = 0
binwidth = 1

simparms = armc.SimulationParams(
    iters,
    steps,
    max_bias_diff,
    write_interval,
    radius_move_freq,
    filebase,
    analytical_biases,
    read_biases,
    biases_filename,
    restart_iter,
    binwidth
)

armc.run!(system, lattice, simparms)
armc.run_us!(system, lattice, simparms)

filebase = "outs/test_run-1"
read_biases = true
biases_filename = "outs/test.biases"
restart_iter = 2
simparms = armc.SimulationParams(
    iters,
    steps,
    max_bias_diff,
    write_interval,
    radius_move_freq,
    filebase,
    analytical_biases,
    read_biases,
    biases_filename,
    restart_iter,
    binwidth
)

armc.run_us!(system, lattice, simparms)

main()
rm("outs", recursive=true)

max_stars_repo_name: cumberworth/ActinRingsMC

max_issues_repo_name: cumberworth/ActinRingsMC

max_forks_repo_name: cumberworth/ActinRingsMC


text:
using Revise
using Makie, CairoMakie, ZigZagBoomerang, SparseArrays, LinearAlgebra
using Random

function zz_sticky_events(xs0, i)
    first = true
    y = []
    for x in xs0
        if x[i] == 0.0  && first == true
            push!(y, x)
            first = false
        elseif x[i] == 0.0 && first == false
            continue
        else
            first = true
        end
    end
    return y
end

    Random.seed!(14)
    function ϕ(x, i, μ)
        x[i] - μ[i]
    end
    κ = 2.0
    n = 2
    μ = zeros(n)
    x0 = randn(n)
    θ0 = rand((-1.0,1.0), n)
    T = 10.0
    c = 0.1*ones(n)
    @time trace0, _ = ZigZagBoomerang.sspdmp(ϕ, 0.0, x0, θ0, T, c, κ, ZigZag(sparse(1.0I,n,n), zeros(n)), μ)
    ts0, xs0 = splitpairs(trace0)
    p1  = Scene()
        lines!(getindex.(xs0,1), getindex.(xs0,2), color=(:black, 0.4), markersize=0.1)
        y1 = zz_sticky_events(xs0, 1)
        y2 = zz_sticky_events(xs0, 2)
        scatter!(getindex.(y1,1), getindex.(y1,2),
            color = (:red, 0.3), strokecolor = (:red, 0.3), markersize = 15, marker = 'x')
        scatter!(getindex.(y2,1), getindex.(y2,2),
            color = (:blue, 0.3), strokecolor = (:blue, 0.3), markersize = 15, marker = 'x')
        scatter!([xs0[1][1]], [xs0[1][2]], color = (:black, 1.0), strokecolor = (:black, 1.0), markersize = 2.0)
        text!("x(0),y(0)",textsize = 0.15, position = (xs0[1][1]- 0.05, xs0[1][2] + 0.05))

p2 = lines(ts0, getindex.(xs0,1), color = (:red, 0.5), label = "x1", leg = true)
    lines!(p2, ts0, getindex.(xs0,2), color = (:blue, 0.5), label = "x2")
    xlabel!(p2, "t")
    ylabel!(p2, "")
p = vbox(p1,p2)

save("./figures/2d_sticky_samplers.png", title(p, "2d sticky Zig-Zag"))


text:
using ControlSystems

"""
    rlocus(G)

Draw the root locus plot
"""
function ControlSystems.rlocus(G::FOTF)
    G1 = fotf2cotf(G)
    ControlSystems.rlocus(G1)
end


text:
"

 author: Atsushi Sakai
"

function main()
    println(PROGRAM_FILE," start!!")

    println(PROGRAM_FILE," Done!!")
end

if abspath(PROGRAM_FILE) == @__FILE__
    @time main()
end

max_stars_repo_name: AtsushiSakai/dotfiles

max_issues_repo_name: AtsushiSakai/dotfiles

max_forks_repo_name: AtsushiSakai/dotfiles


text:
module TestBasics

using Transducers: _unzip, arguments
include("preamble.jl")

@testset "_unzip" begin
    @test _unzip(((Val(0), Val(1), Val(2)), (Val(5), Val(6), Val(7)))) ===
          ((Val(0), Val(5)), (Val(1), Val(6)), (Val(2), Val(7)))

    @test _unzip((
        (Val(0), Val(1), Val(2), Val(3), Val(4)),
        (Val(5), Val(6), Val(7), Val(8), Val(9)),
    )) === (
        (Val(0), Val(5)),
        (Val(1), Val(6)),
        (Val(2), Val(7)),
        (Val(3), Val(8)),
        (Val(4), Val(9)),
    )

    @test _unzip((ntuple(identity, 3), ntuple(x -> 10x, 3), ntuple(x -> 100x, 3))) ===
          ((1, 10, 100), (2, 20, 200), (3, 30, 300))

    if VERSION >= v"1.1"
        @test_inferred _unzip(((Val(0), Val(1), Val(2)), (Val(5), Val(6), Val(7))))
        @test_inferred _unzip((ntuple(Val, 10), ntuple(x -> Val(x + 10), 10)))
    end
end

@testset "arguments(zip(...))" begin
    @testset for args in Any[
        (1:1,),
        (1:1, 2:2),
        (1:1, 2:2, 3:3),
        (1:1, 2:2, 3:3, 4:4),
        (1:1, 2:2, 3:3, 4:4, 5:5),
    ]
        @test arguments(zip(args...)) === args
    end
end

end  # module


text:
module InvolutionMCMC

include("../src/Jaynes.jl")
using .Jaynes
using Gen

jmodel = @jaynes function model()
    z = ({:z} ~ Bernoulli(0.5))
    if z
        m1 = ({:m1} ~ Gamma(1, 1))
        m2 = ({:m2} ~ Gamma(1, 1))
    else
        m = ({:m} ~ Gamma(1, 1))
        (m1, m2) = (m, m)
    end
    {:y1} ~ Normal(m1, 0.1)
    {:y2} ~ Normal(m2, 0.1)
end

jprop = @jaynes function fixed_structure_proposal(tr)
    if tr[:z] == true
        {:m1} ~ Normal(tr[:m1], 0.1)
        {:m2} ~ Normal(tr[:m2], 0.1)
    else
        {:m} ~ Normal(tr[:m], 0.1)
    end
end

select_mh_structure_kernel(tr) = mh(tr, select([(:z, )]))[1]
fixed_structure_kernel(tr) = mh(tr, jprop, ())[1]

# ------------ Simple MH ------------ #

test = () -> begin
    (y1, y2) = (1.0, 1.3)
    obs = choicemap(Pair{Tuple, Any}[(:y1, ) => y1,
                                     (:y2, ) => y2,
                                     (:z, ) => false,
                                     (:m, ) => 1.2])
    tr, _ = generate(jmodel, (), obs)
    for iter=1:100
        tr = select_mh_structure_kernel(tr)
        tr = fixed_structure_kernel(tr)
    end
end
test()

# ------------ Involution DSL ------------ #

function merge_means(m1, m2)
    m = sqrt(m1 * m2)
    dof = m1 / (m1 + m2)
    (m, dof)
end

function split_mean(m, dof)
    m1 = m * sqrt((dof / (1 - dof)))
    m2 = m * sqrt(((1 - dof) / dof))
    (m1, m2)
end

sm_prop = @jaynes function split_merge_proposal(tr)
    if tr[:z]
        # currently two segments, switch to one
    else
        # currently one segment, switch to two
        {:dof} ~ Uniform(0, 1)
    end
end

@involution function split_merge_involution(model_args, proposal_args, proposal_retval)

    if @read_discrete_from_model(:z)

        # currently two segments, switch to one
        @write_discrete_to_model(:z, false)
        m1 = @read_continuous_from_model(:m1)
        m2 = @read_continuous_from_model(:m2)
        (m, dof) = merge_means(m1, m2)
        @write_continuous_to_model(:m, m)
        @write_continuous_to_proposal(:dof, dof)

    else

        # currently one segments, switch to two
        @write_discrete_to_model(:z, true)
        m = @read_continuous_from_model(:m)
        dof = @read_continuous_from_proposal(:dof)
        (m1, m2) = split_mean(m, dof)
        @write_continuous_to_model(:m1, m1)
        @write_continuous_to_model(:m2, m2)
    end
end

split_merge_kernel(tr) = mh(tr, sm_prop, (), split_merge_involution)[1]

test_involution = () -> begin
    (y1, y2) = (1.0, 1.2)
    obs = choicemap(Pair{Tuple, Any}[(:y1, ) => y1,
                                     (:y2, ) => y2,])
    tr, = generate(jmodel, (), obs)
    trs = []
    for iter=1:300
        tr = split_merge_kernel(tr)
        tr = fixed_structure_kernel(tr)
        push!(trs, tr)
    end
    trs
end
trs = test_involution()

display(lineplot(map(trs) do tr
                     tr[:z]
                 end))

end # module


text:
using CartesianDDM
using Printf
using LinearAlgebra
using SparseArrays
using IterativeSolvers
using CartesianCutCell

dofs = (8, 12)
procs = (2, 3)
overs = (4, 2)

dims = CartesianDDMDimension.(dofs, procs, overs)
context = CartesianDDMContext(dims)

A = laplacian(dofs)
dA = decompose(context, A)

P = Diagonal(A)
dP = CartesianDDMRAS(dA)

rhs = rand(prod(dofs))
drhs = decompose(context, rhs)

u = zeros(prod(dofs))
du = decompose(context, u)

res = similar(rhs)
dres = similar(drhs)

io = stdout

str = "%12s%2s" * "%12s" ^ 2 * "\n"
fmt = Printf.Format(str)

Printf.format(io, fmt, "Step", "  ", "Diagonal", "Block diag.")

str = "%12d%2s" * "%12.5E" ^ 2 * "\n"
fmt = Printf.Format(str)

itmax = 50

for it in 1:itmax
    res .= rhs - A * u
    dres .= drhs - dA * du

    u .+= P \ res
    du .+= dP \ dres

    Printf.format(io, fmt, it, "  ", norm(res), norm(dres))
end


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

name = "Edlib"
version = v"1.2.6"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/Martinsos/edlib.git", "ba4272ba68fcdbe31cbc10853de1841701e4e60a")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/edlib/
# these programs do not link for some reason, but we don't need them
mv CMakeLists.txt CMakeLists.txt.orig
egrep -v "helloWorld|runTests|aligner" <CMakeLists.txt.orig >CMakeLists.txt
cd build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DBUILD_SHARED_LIBS=ON ..
make -j$nproc
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_cxxstring_abis(platforms)

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

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
# gcc 6.1.0: there is a bug in stdlib on *-linux-musl in earlier versions (e.g. see https://patchwork.ozlabs.org/patch/544393/)
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"6.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    subsumes_string(u::String, v::String[, u_indexing])

Check whether stringified variable name `v` describes a sub-range of stringified variable `u`.

This is a very restricted version `subumes(u::VarName, v::VarName)` only really supporting:
- Scalar: `x` subsumes `x[1, 2]`, `x[1, 2]` subsumes `x[1, 2][3]`, etc.

## Note
- To get same matching capabilities as `AbstractPPL.subumes(u::VarName, v::VarName)`
  for strings, one can always do `eval(varname(Meta.parse(u))` to get `VarName` of `u`,
  and similarly to `v`. But this is slow.
"""
function subsumes_string(u::String, v::String, u_indexing=u * "[")
    return u == v || startswith(v, u_indexing)
end

"""
    inargnames(varname::VarName, model::Model)

Statically check whether the variable of name `varname` is an argument of the `model`.

Possibly existing indices of `varname` are neglected.
"""
@generated function inargnames(::VarName{s}, ::Model{_F, argnames}) where {s, argnames, _F}
    return s in argnames
end

"""
    inmissings(varname::VarName, model::Model)

Statically check whether the variable of name `varname` is a statically declared unobserved variable
of the `model`.

Possibly existing indices of `varname` are neglected.
"""
@generated function inmissings(::VarName{s}, ::Model{_F, _a, _T, missings}) where {s, missings, _F, _a, _T}
    return s in missings
end


text:
using BenchmarkTools

SUITE = BenchmarkGroup()
BCST_A = parse(Int, ENV["BCST_A"])
for i in [0, 1]
    k1 = i * BCST_A
    SUITE["k1=$k1"] = @benchmarkable nothing
end


text:
module BacktrackRTS

export Node, runCurry, trailNode, makeChoice, setChoice, makeFree, makeNode,
       setRHS, setConst, setForwardNode, bindVarNode, voidNode, rew_nf

##############################################################################
# Definition of graph nodes
mutable struct Node
  tag    :: UInt8       # 0: operation
                        # 1: forward
                        # 2: choice
                        # 3: fail
                        # 4: free
                        # 5: integer
                        # 6: float
                        # 7: char
                        # 8: partial function call
                        # 9: partial constructor call
                        # >=10: constructor
  value  :: Int64       # demanded argument (0: nothing demanded) for functions,
                        # choice identifier for choices,
                        # index of constructors or integer/char value
                        # missing arguments for partial calls
  otsk   :: Int64       # identifier to remember creation time (choice point)
  fcode                 # function code
  args   :: Array{Node} # arguments or forward node
  symbol :: Union{String,Float64} # string representation of node symbol
                                  # or value of a float
end

##############################################################################
# Auxiliaries for implementing backtracking

# Elements stored on the backtrack stack
mutable struct ChoicePoint
  choicenum :: Int64        # number of this choice point
  stack     :: Vector{Node} # the stack contents at the choice point
  trail     :: Vector{Tuple{Node,Node}} # the trail containing all node changes
                                        # to be undone on backtracking
end

# The global backtrack stack. Initially, it contains a void choice point.
backtrackstack = [ChoicePoint(0, [], [])]

# This global variable contains the value of `last(backtrackstack).choicenum`.
# Used for faster access since its value is used at various places.
currentchoicenum = 0

# Counter for trail events (just for statistics)
trailcounter = 0

# Push a node together with a copy of its contents on the trail,
# but only if it has been created before the current choice point.
# The latter condition avoids unnecessary trailing of deterministic
# computations.
function trailNode(n :: Node)
  global trailcounter
  global currentchoicenum

  if n.otsk < currentchoicenum
    #print("TRAIL NODE: "); printNode(n); println()
    trailcounter += 1
    cn = Node(n.tag, n.value, n.otsk, n.fcode, copy(n.args), n.symbol)
    currentcp = last(backtrackstack)
    push!(currentcp.trail,(n,cn))
  end
end

##############################################################################
# Evaluate all (non-deterministic) results of a main operation.
function runCurry(mainname::String, maincode, withnf::Bool, interactive::Bool,
                  withbfs::Bool, onevalue::Bool, quiet::Bool)
  global currentchoicenum = 0
  global trailcounter     = 0
  global backtrackstack   = [ChoicePoint(0, [], [])]

  stack ::Vector{Node} = []
  numchoices = 0 # number of created choice points
  nd = Node(0, 0, 0, maincode, [], mainname)
  mainnd = withnf ? makeNormalForm(nd) : nd
  finished = false
  if !quiet ; print("EVALUATE: "); printNode(mainnd); println() ; end
  cn = mainnd
  while ! (finished || (isempty(stack) && isempty(backtrackstack)))
    #traceEval(isempty(stack) ? cn : stack[1])
    if cn.tag >= 4              # free, constructor, integer,...
      if isempty(stack)
        if !quiet ; print("RESULT: ") end
        printNode(cn); println()
        finished = onevalue
        if !finished && interactive
          answer = "?"
          while answer != "" && answer != "y" && answer != "n" && answer != "a"
            print("More values? [Y(es)/n(o)/a(ll)] ")
            answer = chomp(readline())
            if answer == "a"
              interactive = false
            elseif answer == "n"
              finished = true
            end
          end
        end
        cn.tag = 3 # initiate backtracking
      else
        cn = pop!(stack)
        cn.fcode(cn)              # apply rewrite rule
      end

    elseif cn.tag == 1  # forward node?
      #println("FORWARD")
      cn = cn.args[1]
      if !isempty(stack)
        topfun = last(stack)
        trailNode(topfun)
        topfun.args[topfun.value] = cn
      end

    elseif cn.tag == 2  # choice node? create choice point and follow left:
      currentchoicenum += 1
      #println("CREATE CHOICE POINT " * string(currentchoicenum))
      # create choice point:
      savedstack = copy(stack)
      push!(savedstack, cn)
      newcp = ChoicePoint(currentchoicenum, savedstack, [])
      push!(backtrackstack, newcp)
      numchoices += 1
      setForwardNode(cn,cn.args[1])
      cn = cn.args[1]
      if !isempty(stack)
        topfun = last(stack)
        trailNode(topfun)
        topfun.args[topfun.value] = cn
      end

    elseif cn.tag == 3  # failure node?
      cp = pop!(backtrackstack)
      if length(cp.stack) == 0 # last choice point?
        stack = []
      else
        #println("BACKTRACK: RESTORE CHOICE POINT " * string(cp.choicenum))
        currentchoicenum = last(backtrackstack).choicenum
        # restore trailed changes of choice point:
        while length(cp.trail) > 0
          tn = pop!(cp.trail)
          tn1 = tn[1]
          tn2 = tn[2]
          tn1.tag    = tn2.tag
          tn1.value  = tn2.value
          tn1.otsk   = tn2.otsk
          tn1.fcode  = tn2.fcode
          tn1.args   = tn2.args
          tn1.symbol = tn2.symbol
          #print("UNTRAIL NODE: "); printNode(tn1); println()
          #traceEval(isempty(stack) ? cn : stack[1])
        end
        stack = cp.stack
        cn = pop!(stack)
        #println("GO RIGHT IN CHOICE POINT " * string(cp.choicenum))
        #traceEval(isempty(stack) ? cn : stack[1])
        setForwardNode(cn, cn.args[2])
        cn = cn.args[1]
        if !isempty(stack)
          topfun = last(stack)
          trailNode(topfun)
          topfun.args[topfun.value] = cn
        end
      end

    # final alternative: cn.tag == 0, i.e., operation
    else
      if cn.value == 0 # no argument demanded?
        cn.fcode(cn) # apply rewrite rule
      else
        push!(stack,cn)   # push current function on stack
        cn = cn.args[cn.value] # proceed with demanded argument
      end
    end
  end
  if !quiet
    println("NUMBER OF TRAILED ELEMENTS: " * string(trailcounter))
    println("NUMBER OF CREATED CHOICE POINTS: " * string(numchoices))
  end
end

##############################################################################
# Auxiliary operations for dealing with choice nodes

# Create a choice node with a fresh choice identifier
function makeChoice(otsk, l, r) :: Node
  return Node(2, 0, otsk, nothing, [l,r], "?")
end

# Replace argument node by a choice node
function setChoice(root :: Node, l :: Node, r :: Node)
  trailNode(root)
  root.symbol = "?"
  root.tag    = 2
  root.value  = 0
  empty!(root.args); push!(root.args,l); push!(root.args,r)
end

##############################################################################
# Auxiliary operations

# Node representing a void value. Used to initialize Node variables.
function voidNode(otsk)
  return Node(10, 0, otsk, nothing, [], "()")
end

# Create a free variable node
function makeFree(otsk) :: Node
  return makeNode(4, 0, nothing, [], "free")
end

# Replace first argument node by content of second argument node.
# Used to instantiate free variable nodes
function makeNode(ntag, nvalue, nfcode, nargs, nsymbol)
  global currentchoicenum
  return Node(ntag, nvalue, currentchoicenum, nfcode, nargs, nsymbol)
end

# Set the first argument node (which is the root of a redex)
# to the right-hand side provided as the remaining arguments.
function setRHS(root :: Node, ntag, nvalue, nfcode, nargs, nsymbol)
  global currentchoicenum
  trailNode(root)
  root.otsk    = currentchoicenum
  root.tag     = ntag
  root.value   = nvalue
  root.fcode   = nfcode
  root.args    = nargs
  root.symbol  = nsymbol
end

# Set the first argument node (which is the root of a redex)
# to a constant value provided as the remaining arguments.
function setConst(root :: Node, ntag, nvalue)
  global currentchoicenum
  trailNode(root)
  root.otsk    = currentchoicenum
  root.tag     = ntag
  root.value   = nvalue
end

# Set the first argument node (which is the root of a redex)
# to a float value provided as the remaining arguments.
function setFloatConst(root :: Node, fvalue)
  global currentchoicenum
  trailNode(root)
  root.otsk    = currentchoicenum
  root.tag     = 6
  root.symbol  = fvalue
end

# Set the first argument node to a forward node
function setForwardNode(root :: Node, fnode :: Node)
  trailNode(root)
  root.tag  = 1
  empty!(root.args); push!(root.args, fnode)
end

# Bind the first argument node (a free variable) to another node.
# Implemented by setting a forward node.
function bindVarNode(varnode :: Node, n :: Node)
  setForwardNode(varnode, n)
end

##############################################################################
# ADDITIONAL STUFF TO SUPPORT THE REPL:
##############################################################################
# Normal form computation.

# Wrap a node with the `normalForm` operation.
function makeNormalForm(root :: Node)
  return makeNode(0, 1, rew_nf, [root], "normalForm")
end

# Operation to compute the normal form of its argument.
# Implemented by wrapping all arguments with `normalForm` and
# evaluating successively all arguments (by operation `rew_nfArgs`).
function rew_nf(root :: Node)
  trailNode(root)
  x = root.args[1]
  if x.tag == 4 # free variable: create forward node
    root.tag  = 1
    empty!(root.args); push!(root.args,x)
  elseif x.tag < 10 || length(x.args) == 0 # no argument evaluation necessary
    root.tag    = x.tag
    root.symbol = x.symbol
    root.value  = x.value
    root.otsk   = x.otsk
    root.fcode  = nothing
    root.args   = x.args
  else
    root.value = 1
    root.fcode = rew_nfArgs
    root.args  = vcat(map(makeNormalForm,x.args), [x]) # remember constructor
  end
end

function rew_nfArgs(root :: Node)
  trailNode(root)
  darg = root.value + 1
  if darg == length(root.args) # all constructor arguments evaluated
    cterm       = pop!(root.args)
    root.tag    = cterm.tag
    root.value  = cterm.value
    root.otsk   = cterm.otsk
    root.fcode  = nothing
    root.symbol = cterm.symbol
  else
    root.value = darg
  end
end

##############################################################################
# Auxiliaries for printing results and tracing

# Print a graph node as a term.
function printNode(n :: Node)
  if n.tag == 1 # forward node
    print("^") ; printNode(n.args[1])
  elseif n.tag == 2 # choice node
    print("(")
    printNode(n.args[1])
    print(" ?" * string(n.value) * " ")
    printNode(n.args[2])
    print(")")
  elseif n.tag == 3 # fail node
    print("!")
  elseif n.tag == 4 # free node
    #print("FREE")
    #print("_x" * string(objectid(n)))
    print("_x" * string(objectid(n))[end-3:end])
    #printNodeSymbol(n)
  elseif n.tag == 5 # integer node
    print(string(n.value))
  elseif n.tag == 7 # character node
    if n.value > 31
      print("'" * Char(n.value) * "'")
    else
      print("(char " * string(n.value) * ")")
    end
  elseif length(n.args) == 0
    printNodeSymbol(n)
  elseif length(n.args) == 2 && !isletter(n.symbol[1]) # infix operator
    print("(")
    printNode(n.args[1])
    print(" ") ; printNodeSymbol(n) ; print(" ")
    printNode(n.args[2])
    print(")")
  else
    print("(")
    printNodeSymbol(n)
    for i = 1 : length(n.args)
      print(" ") ; printNode(n.args[i])
    end
    print(")")
  end
end

# Print a node symbol
function printNodeSymbol(n :: Node)
  #print(n.symbol * (n.tag == 0 ? "[" * string(objectid(n)) * "]" : ""))
  #print(n.symbol * (n.tag == 0 ? "[" * string(n.otsk) * "]" : ""))
  print(n.symbol)
end

function traceEval(cn :: Node)
  print("ROOT: ")
  printNode(cn)
  println()
end

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

end

max_stars_repo_name: cau-placc/julia-curry

max_issues_repo_name: cau-placc/julia-curry

max_forks_repo_name: cau-placc/julia-curry


text:
abstract type AbstractCell{dim,V,F} end

# Vertices
struct Vertex{dim,T}
    x::Tensors.Vec{dim, T}
end

Base.eltype(::Vertex{dim,T}) where {dim,T} = T

@inline get_coordinates(vertex::Vertex) = vertex.x

#--------------- Cells
struct Cell{dim, N, M, L}
    vertices::NTuple{N,Int}
end

#Common cell types
const LineCell = Cell{1,2,1,1}
@inline get_cell_name(::LineCell) = "Line"

const TriangleCell = Cell{2,3,3,1}
@inline get_cell_name(::TriangleCell) = "Triangle"
@inline reference_edge_vertices(::TriangleCell) = ((2,3),(3,1),(1,2))

const RectangleCell = Cell{2,4,4,1}
@inline get_cell_name(::RectangleCell) = "Rectangle"
@inline reference_edge_vertices(::RectangleCell) = ((1,2),(2,3),(3,4),(4,1))

@inline reference_edge_vertices(::Cell{2,N,N,1})  where {N} = Tuple((i,mod1(i+1,N)) for i in 1:N)

const HexagonCell = Cell{2,6,6,1}
@inline get_cell_name(::HexagonCell) = "Hexagon"

const TetrahedronCell = Cell{3,4,6,4}
@inline get_cell_name(::TetrahedronCell) = "Tetrahedron"
reference_facet_vertices(::Type{TetrahedronCell}) = reference_face_nodes(Tetrahedron)

const HexahedronCell = Cell{3,8,12,6}
@inline get_cell_name(::HexahedronCell) = "Hexahedron"
reference_facet_vertices(::Type{HexahedronCell}) = reference_face_nodes(Hexahedron)

# API
@inline getnvertices(cell::Cell{dim,N}) where {dim,N} = N
@inline getnedges(cell::Cell{dim,N,M}) where {dim,N,M} = M
@inline getnfaces(cell::Cell{dim,N,M,P}) where {dim,N,M,P} = P
gettopology(cell::Cell{1,N,M,P}) where {N,M,P} = Dict(0=>N,1=>M)
gettopology(cell::Cell{2,N,M,P}) where {N,M,P} = Dict(0=>N,1=>M,2=>P)
gettopology(cell::Cell{3,N,M,P}) where {N,M,P} = Dict(0=>N,1=>M,2=>P,3=>1)

# ----------------- Mesh
struct PolytopalMesh{dim,T,C} <: AbstractPolytopalMesh{dim,T}
    cells::Vector{C}
    vertices::Vector{Vertex{dim,T}}
    # Sets
    cellsets::Dict{String,Set{Int}}
    facesets::Dict{String,Set{FaceIndex}}
    edgesets::Dict{String,Set{EdgeIndex}}
    vertexsets::Dict{String,Set{Int}}
end

function PolytopalMesh(cells,
              vertices;
              cellsets::Dict{String,Set{Int}}=Dict{String,Set{Int}}(),
              facesets::Dict{String,Set{FaceIndex}}=Dict{String,Set{FaceIndex}}(),
              edgesets::Dict{String,Set{EdgeIndex}}=Dict{String,Set{EdgeIndex}}(),
              vertexsets::Dict{String,Set{Int}}=Dict{String,Set{Int}}()) where {dim}
    return PolytopalMesh(cells, vertices, cellsets, facesets, edgesets, vertexsets)
end

# Generic Interface
getfacet(mesh::PolytopalMesh, facet) = facet
getdim(mesh::PolytopalMesh{dim}) where {dim} = dim
getncellvertices(mesh::PolytopalMesh, cell_idx::Int) = getnvertices(mesh.cells[cell_idx])
gettopology(mesh::PolytopalMesh, cell::Cell) = gettopology(cell)
getnvertices(mesh::PolytopalMesh, cell::Cell) = getnvertices(cell)
getnedges(mesh::PolytopalMesh, cell::Cell) = getnedges(cell)
getcell(mesh::PolytopalMesh, idx::Int) = mesh.cells[idx]
getCellType(mesh::PolytopalMesh, idx::Int) = typeof(getcell(mesh,idx))

function getcellsubentities(mesh::PolytopalMesh{2},cellidx::Int,element::Int)
  if element == 0
      return mesh.cells[cellidx].vertices
  elseif element == 1
      return Tuple(EdgeIndex(cellidx,i) for i in 1:getnedges(mesh.cells[cellidx]))
  else
      throw("Topology element of order $element not available for cell type")
  end
end

function entityeltype(mesh::PolytopalMesh{2},dim)
    if dim == 0
        return Int
    elseif dim == 1
        return EdgeIndex
    else
        error("mesh $mesh has not entity of dim $dim")
    end
end

reference_edge_vertices(mesh::PolytopalMesh, cell::Cell) = reference_edge_vertices(cell)

function getverticescoords(mesh::PolytopalMesh{dim,T}, cell::Cell) where {dim,T}
    N = getnvertices(cell)
    coords = Vector{Tensors.Vec{dim,T}}(undef, N)
    for (i,j) in enumerate(cell.vertices)
        coords[i] = mesh.vertices[j].x
    end
    return coords
end

# API

@inline getncells(mesh::PolytopalMesh) = length(mesh.cells)
@inline getnvertices(mesh::PolytopalMesh) = length(mesh.vertices)
@inline getverticesidx(mesh::PolytopalMesh, cell_idx) = mesh.cells[cell_idx].vertices
@inline getvertexset(mesh::PolytopalMesh, set::String) = mesh.vertexsets[set]
@inline getedgeset(mesh::PolytopalMesh, set::String) = mesh.edgesets[set]
getcells(mesh::PolytopalMesh) = mesh.cells

"""
function getcoords(mesh, vertex_idx::Int)
Return a Tensor.Vec with the coordinates of vertex with index `vertex_idx`
"""
getvertexcoords(mesh::PolytopalMesh, vertex_idx::Int) = mesh.vertices[vertex_idx].x

"""
    getverticescoords(mesh::PolytopalMesh, cell_idx)
Return a vector with the coordinates of the vertices of cell number `cell`.
"""
function getverticescoords(mesh::PolytopalMesh{dim,T}, cell_idx::Int) where {dim,T}
    N = getnvertices(mesh.cells[cell_idx])
    coords = Vector{Tensors.Vec{dim,T}}(undef, N)
    for (i,j) in enumerate(mesh.cells[cell_idx].vertices)
        coords[i] = mesh.vertices[j].x
    end
    return coords
end

function getvertexcoords(mesh::PolytopalMesh{dim,T}, cell::Cell, vidx::Int) where {dim,T}
    return mesh.vertices[cell.vertices[vidx]].x
end

function getverticescoords(mesh::PolytopalMesh{dim,T}, edge_idx::EdgeIndex) where {dim,T}
    cell = getcells(mesh)[edge_idx.cellidx]
    ref_edge = reference_edge_vertices(mesh,cell)[edge_idx.idx]
    coords = Vector{Tensors.Vec{dim,T}}(undef, 2)
    for i in 1:2
        coords[i] = mesh.vertices[cell.vertices[ref_edge[i]]].x
    end
    return coords
end

function getverticesindices(mesh::PolytopalMesh{dim,T}, edge_idx::EdgeIndex) where {dim,T}
    cell = getcells(mesh)[edge_idx.cellidx]
    ref_edge = reference_edge_vertices(mesh,cell)[edge_idx.idx]
    return [cell.vertices[ref_edge[i]] for i in 1:2]
end

function get_cell_connectivity_list(mesh::PolytopalMesh{dim,T,C}) where {dim,T,C}
    cells_m = Vector()
    for k = 1:getncells(mesh)
        push!(cells_m,mesh.cells[k].vertices)
    end
    cells_m
end

getverticesindices(mesh::PolytopalMesh,cell::Cell) = cell.vertices

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
"""

Structure that contains a two-dimensional point, subtype of FieldVector from StaticArrays.

"""
struct Point2D <: FieldVector{2,Float64}
  x::Float64
  y::Float64
end
norm2(p::Point2D) = (p.x^2 + p.y^2)
norm(p::Point2D) = sqrt(norm2(p))

"""

Structure that contains a three-dimensional point, subtype of FieldVector from StaticArrays.

"""
struct Point3D <: FieldVector{3,Float64}
  x::Float64
  y::Float64
  z::Float64
end
norm2(p::Point3D) = (p.x^2 + p.y^2 + p.z^2)
norm(p::Point3D) = sqrt(norm2(p))

"""

```
dim(T::DataType)
```

Returns the dimension associated to a `Point2D` (`2`) or `Point3D` (`3`).

"""
function dim(T::DataType)
  T == Point2D && return 2
  T == Point3D && return 3
end


text:
using LinearAlgebra, Random, DataFrames, CSV, Printf, CPLEX, Gurobi;
using Distributed
addprocs(6)
@everywhere using JuMP, BARON
include("DecompositionAlgorithm_AdaptiveSampling.jl")
include("DecompositionAlgorithm_RandomSampling.jl")
include("Utilities.jl")

@everywhere root=("//panfs//roc//groups//10//qizh//rishabh//ProdPlanning1//Adaptive//")
@everywhere include("ASP.jl")

function ModelParameters(II, JJ, TT, Mu, D)
    Pmax = 400*ones(Int16, 38)
    Q0 = zeros(Int16, 28)
    Q0[10:26] = [10 5 0 5 5 7 10 2 3 5 10 7 5 5 3 3 0]
    Qmin = zeros(Int16, 28)
    Qmax = 200*ones(Int16, 28)
    Wmax = zeros(Int16, 28, TT)
    Wmax[1:10, 1] = [200 125 200 100 50 250 600 40 1000 20]
    for t = 1:TT
        Wmax[1:10, t] = 400*ones(Int16, 10)
    end
    IhatJ = zeros(Int16, 28, 38)
    IbarJ = zeros(Int16, 28, 38)

    for j = 1:JJ
        for i = 1:II
            if Mumax[i,j] <= -0.001
                IhatJ[j,i] = 1
            elseif Mumax[i,j] >= 0.001
                IbarJ[j,i] = 1
            end
        end
    end

    A = zeros(Float64, (4*JJ+2*II)*TT, (II+JJ)*TT)
    b = zeros(Float64, (4*JJ+2*II)*TT)

    for i = 1:II
        for j = 1:JJ
            Mu[i,j] = abs(Mu[i,j])
        end
    end

    # Constraint Matrix
    for t = 1:TT
        for j = 1:JJ
            b[JJ*(t-1)+j] = -(Qmin[j] - Q0[j] + sum(D[j, tprime] for tprime = 1:t))
            b[JJ*TT + JJ*(t-1)+j] = Qmax[j] - Q0[j] + sum(D[j, tprime] for tprime = 1:t)
            b[2*JJ*TT+JJ*(t-1)+j] = Wmax[j,t]
            b[3*JJ*TT+JJ*(t-1)+j] = 0
            for tprime = 1:t
                for i = 1:II
                    if IhatJ[j,i] == 1
                        A[JJ*(t-1)+j, II*(tprime-1)+i] = -Mu[i,j]
                        A[JJ*TT+JJ*(t-1)+j, II*(tprime-1)+i] = Mu[i,j]
                    elseif IbarJ[j,i] == 1
                        A[JJ*(t-1)+j, II*(tprime-1)+i] = Mu[i,j]
                        A[JJ*TT+JJ*(t-1)+j, II*(tprime-1)+i] = -Mu[i,j]
                    end
                    A[JJ*(t-1)+j, II*TT+JJ*(tprime-1)+j] = -1
                    A[JJ*TT+JJ*(t-1)+j, II*TT+JJ*(tprime-1)+j] = 1
                    A[4*JJ*TT+II*(t-1)+i, II*(t-1)+i] = 1
                    A[4*JJ*TT+II*TT+II*(t-1)+i, II*(t-1)+i] = -1
                    b[4*JJ*TT+II*(t-1)+i] = Pmax[i]
                    b[4*JJ*TT+II*TT+II*(t-1)+i] = 0
                end
            end
            A[2*JJ*TT+JJ*(t-1)+j, II*TT+JJ*(t-1)+j] = 1
            A[3*JJ*TT+JJ*(t-1)+j, II*TT+JJ*(t-1)+j] = -1
        end
    end
    return A, b
end

global rng = MersenneTwister(1234)

II = 38 # No of Processes
JJ = 28 # No of Chemicals involved
H = 1 # No of time periods
TT = H # This code uses TT isntead of H as the total number of time periods
S = 100 # number of samples for heuristic adative sampling
dim = (II+JJ)*TT
global n = dim

global Nominal_Demand = zeros(Int16, 28)
Nominal_Demand[10:25] = [40 150 30 75 60 30 30 50 150 70 30 75 100 75 250 50]
global Mumax = zeros(Float64, 38, 28)
Mumax[1,1] = 0.58; Mumax[1,6] = 0.63; Mumax[1,11] = -1;
Mumax[2,6] = 0.64; Mumax[2,12] = -1;
Mumax[3, 1] = -.055; Mumax[3,2] = 1.25; Mumax[3,11] = -1;
Mumax[4,2] = 0.40; Mumax[4,3] = 0.69; Mumax[4, 14] = -1;
Mumax[5, 13] = -1; Mumax[5,14] = 2.3; Mumax[5,17] = -1.7;
Mumax[6, 2] = 0.74; Mumax[6,15] = -1;
Mumax[7, 13] = -1; Mumax[7,15] = 1.1;
Mumax[8, 3] = 1; Mumax[8,16] = -1;
Mumax[9,16] = 1.26; Mumax[9,17] = -1;
Mumax[10,13] = 1.57; Mumax[10,27] = -1;
Mumax[11, 3] = 1.01; Mumax[11,17] = -1;
Mumax[12,3] = 0.76; Mumax[12, 4] = 0.28; Mumax[12, 8] = -1;
Mumax[13, 8] = 1.14; Mumax[13, 18] = -1;
Mumax[14,2] = 0.78; Mumax[14,13] = -1;
Mumax[15, 12] = -1; Mumax[15, 19] = 1.34;
Mumax[16, 4] = 0.60; Mumax[16,19] = -1;
Mumax[17, 4] = 0.67; Mumax[17, 12] = -1;
Mumax[18, 12] = 1.1; Mumax[18,20] = -1;
Mumax[19, 19] = 0.98; Mumax[19, 20] = -1;
Mumax[20, 4] = 0.35; Mumax[20,20] = 0.71; Mumax[20,21] = -1;
Mumax[21, 6] = 0.32; Mumax[21,20] = 0.72; Mumax[21,21] = -1;
Mumax[22, 4] = 0.88; Mumax[22, 5] = -1; Mumax[22, 24] = -0.03;
Mumax[23, 1] = 0.56; Mumax[23, 5] = 0.92; Mumax[23, 11] = -1;
Mumax[24, 4] = 0.39; Mumax[24, 28] = -1;
Mumax[25, 5] = -1; Mumax[25, 28] = -1;
Mumax[26, 4] = 0.3; Mumax[26, 23] = -1;
Mumax[27, 20] = 0.65; Mumax[27, 22]= -1; Mumax[27, 27] = 0.46;
Mumax[28, 7] = 0.56; Mumax[28, 10] = 0.56; Mumax[28, 20] = -1;
Mumax[29, 12] = 1.2; Mumax[29, 22] = -1;
Mumax[30, 4] = 1.17; Mumax[30, 25] = -1;
Mumax[31, 5] = 0.75; Mumax[31, 24] = -1;
Mumax[32, 4] = 0.53; Mumax[32, 24] = -1;
Mumax[33, 12] = 0.6; Mumax[33, 20] = 0.82; Mumax[33, 21] = -1;
Mumax[34, 10] = 0.42; Mumax[34, 25] = -1;
Mumax[35, 9] = 0.5; Mumax[35, 10] = -1;
Mumax[36, 7] = 0.53; Mumax[36, 24] = -1; Mumax[36, 25] = 0.57;
Mumax[37, 24] = -1; Mumax[37, 28] = 1.44;
Mumax[38, 2] = -0.38; Mumax[38, 3] = -0.22; Mumax[38, 4] = -1; Mumax[38, 9] = 3.08; Mumax[38, 26] = -1.81;
Mu = deepcopy(Mumax)
D = zeros(Float64, JJ, TT)
for t = 1:TT
    D[:, t] = Nominal_Demand[1:JJ]
end

A, b = ModelParameters(II, JJ, TT, Mu, D)

obj = rand(rng,1:1000,dim)
obj = -1*obj/ norm(obj, 1)
noisy_p = rand(rng,1.0:1000.0, 10, dim)
for i = 1:10
    noisy_p[i, :] = -1*noisy_p[i,:]/norm(noisy_p[i,:], 1)
end

global c = length(A[:,1])
TS = 100
test_x  = zeros(Float16, TS, n)
b_test = zeros(Float64, TS, c)
A_test = zeros(Float64, c, n, TS)

for test = 1:TS
    mutest = zeros(Float64, II, JJ)
    Dtest = zeros(Float64, JJ, TT)
    for i = 1:II
        for j = 1:JJ
            mutest[i,j] = (rand(rng, 50:100)/100)*Mumax[i,j]
            for t = 1:TT
                Dtest[j, t] = (rand(rng, 90:110)/100)*Nominal_Demand[j]
            end
        end
    end
    A_test[:,:,test], b_test[test, :] = ModelParameters(II,JJ,TT,mutest,Dtest)
    test_x[test, :] = SolveModel(A_test[:, :, test], obj, b_test[test, :], n)
end

instance = 1
adaptive_sampling(A,b,obj,noisy_p,A_test,b_test,test_x,instance,S)
random_sampling(A,b,obj,noisy_p,A_test,b_test,test_x,instance)

max_stars_repo_name: ddolab/InvLP

max_issues_repo_name: ddolab/InvLP

max_forks_repo_name: ddolab/InvLP


text:
function parseMask(line)
    # result = "XXXXXXXXXXXXXXXXXXXXXXXXXXXX" * String(split(line,"mask = ")[2])
    result = "0000000000000000000000000000" * String(split(line,"mask = ")[2])
    return result
end

struct Instruction
    pos::Int
    value::Int
end

function parseInstruction(line)
    insideSquareBrackets = r"\[([^()]*)\]"
    pos = match(insideSquareBrackets, line).captures[1]
    value = split(line, " = ")[2]
    result = Instruction(parse(Int, pos), parse(Int, value))
    return result
end

function applyCondition(x::Tuple{Char,Char})
    m = x[1]
    c = x[2]
    if m == 'X'
        return 'X'
    elseif m == '0'
        return c
    elseif m == '1'
        return '1'
    else
        throw("wrong character for mask")
    end
end

function applyMask(mask::String, pos::Int)
    println(mask)
    posBstr = bitstring(pos)
    println(posBstr)
    c = zip(mask, posBstr)
    resultStr = String(map(x -> applyCondition(x), c))
    #result = parse(Int, resultStr, base=2)
    return resultStr
end

function applyPermutation(floatingPos::String, permutation)
    result = collect(Char, floatingPos)
    for tuple in permutation
        index = tuple[1]
        value = tuple[2]
        result[index] = value
    end
    return String(result)
end

function execute(mask, instruction, memory)
    floatingPos = applyMask(mask, instruction.pos)
    xs = findall(x -> x == 'X', floatingPos)
    permutations = (2^(length(xs)))
    println(permutations)
    for i in 0:permutations-1
        bits = string(i, base = 2, pad = length(xs))
        permutation = zip(xs, bits)
        newPosStr = applyPermutation(floatingPos, permutation)
        # println(newPosStr)
        newPos = parse(Int, newPosStr, base = 2)
        memory[newPos] = instruction.value
    end
end

function runProgram(inputFile)
    memory = Dict()
    mask = ""
    for line in readlines(inputFile)
        if startswith(line,"mask")
            mask = parseMask(line)
        else
            instruction = parseInstruction(line)
            execute(mask, instruction, memory)
        end
    end
    println(reduce(+,values(memory)))
end

runProgram("data.txt")
# memory = zeros(Int, 100000)
# mask = parseMask("mask = X00XX111000XX011110110X1X0101000101X")
# instruction = Instruction(0,0)
# execute(mask, instruction, memory)


text:
groupStart() = @apicall(:ncclGroupStart, ())
groupEnd() = @apicall(:ncclGroupEnd, ())
function group(f)
    groupStart()
    f()
    groupEnd()
end


text:
@testset "Library" begin
    @testset "Camel Case Splitting" begin
        @test split(camel_split_l_c1) == ["my", "Camel", "2"]
        @test split(camel_split_l_c2) == ["my", "Camel", "2", "KHR", "Ext"]
        @test split(camel_split_l_c3) == ["my", "Camel", "2", "Ext", "4"]
        @test split(camel_split_u_c1) == ["My", "Camel", "2"]
        @test split(camel_split_u_c2) == ["My", "Camel", "2", "KHR", "Ext"]
        @test split(camel_split_u_c3) == ["My", "Camel", "2", "Ext", "4"]
    end
    @testset "Conversions" begin
        @test Base.convert(CamelCaseLower, const_s1) == const_c1
        @test Base.convert(CamelCaseLower, const_s2) == const_c1
        @test Base.convert(CamelCaseUpper, const_s1) == const_c2
        @test Base.convert(CamelCaseUpper, const_s2) == const_c2
        @test Base.convert(SnakeCaseLower, const_c1) == const_s1
        @test Base.convert(SnakeCaseLower, const_c2) == const_s1
        @test Base.convert(SnakeCaseUpper, const_c1) == const_s2
        @test Base.convert(SnakeCaseUpper, const_c2) == const_s2
    end
    @testset "Remove parts" begin
        @test remove_parts(longstr, [1, 3]).value == "snake_with_many_words"
        @test remove_parts(longstr, 1).value == "snake_case_with_many_words"
    end
end


text:
using ImageView, Images   # TestImages  just a reminder poop stinks!
using Gtk.ShortNames

grid, frames, canvases = canvasgrid((1,1))  # 1 row, 2 columns

#
#img = load("/home/travis/Documents/Mael Engine/sprites/cloud.png")
cloud = load("/home/travis/Documents/Mael Engine/sprites/cloud.png")
cloud2 = load("/home/travis/Documents/Mael Engine/sprites/cloud2.png")
LeftPlayer = load("/home/travis/Documents/Mael Engine/sprites/leftplayer.png")
MiddlePlayer = load("/home/travis/Documents/Mael Engine/sprites/middleplayer.png")
img = fill(RGBA(0.65625,0.796875,1,1), 500, 800)

imshow(canvases[1,1], img )
#win = @Window(canvas, "Grinding Gears", 500, 800)

win = Window(grid, "Grinding Gears", 800, 500)
showall(win)
# ==============================================================================
# This draws a line to a picture!
# ==============================================================================
function line(image, x1, y1, x2, y2, color)
        dx, dy = x2-x1,  y2-y1

        for x in x1:x2
              y = y1 + dy * (x - x1) / dx
              image[ round(Int32, y), round(Int32, x) ] = color
        end
end
# ==============================================================================
# This draws a rectrangle to a picture!
# ==============================================================================
function rect(image, left, top, width, height, color)
    for y in 1:height, x in 1:width
        pixel(image, top+y, left+x, color )
    end
end
# ==============================================================================
function clearScreen(image, color)
    height, width = size(image)
    for y in 1:height, x in 1:width
        image[ y, x ] = color
    end
end
# ==============================================================================
# This draws a rectrangle to a picture! 1="centered" 2= topCenter
# ==============================================================================
function drawSprite(image, sprite, x, y, position=1)
    # width, height = size(image)
    SpHeight, SpWidth = size(sprite)
    if position==1
        cx,cy = round(Int32, SpWidth/2), round(Int32, SpHeight/2)
    elseif  position==2
        cx,cy = round(Int32, SpWidth/2), SpHeight
    elseif position==3
        cx,cy = round(Int32, SpWidth/2), 0
    end
    left, top = x-cy, y-cx

    for y in 1:SpWidth,  x in 1:SpHeight
        pixel(image, left+x, top+y, sprite[ x, y ] )
    end
end
# ==============================================================================
function drawSprite(image, sprite) # sprite.image, sprite.x, sprite.y, sprite.position=1
    SpHeight, SpWidth = size(sprite.image)
    if sprite.position==1
        cx,cy = round(Int32, SpWidth/2), round(Int32, SpHeight/2)
    elseif  sprite.position==2
        cx,cy = round(Int32, SpWidth/2), SpHeight
    elseif sprite.position==3
        cx,cy = round(Int32, SpWidth/2), 0
    end
    left, top =  sprite.y-cy, sprite.x-cx

    for y in 1:SpWidth,  x in 1:SpHeight
        pixel(image, left+x, top+y, sprite.image[ x, y ] )
    end
end
# ==============================================================================
# This draws a pixel with Alpha chanel
# ==============================================================================
function pixel(image, X, Y, color )
    a = alpha(color)
    width, height = size(image)
    if X>0 && X < width && Y>0 && Y < height
        image[ X,Y ] = image[ X,Y ]*(1-a) + (color*a)
    end
end
# ==============================================================================
# This calculates the positions of all sprites draws them.
# ==============================================================================
function Paint(spritList)
     height, width  = size(img)
    clearScreen(img, RGBA(0.65625,0.796875,1,1))
    for sprit in spritList
        drawSprite(img, sprit)
        if sprit.direction == 0
        elseif sprit.direction == 1
            if sprit.x < width
                sprit.x += 4
            else
                sprit.direction =2
            end
        elseif sprit.direction == 2
            if sprit.x > 0
                sprit.x -= 4
            else
                sprit.direction = 1
            end
        elseif sprit.direction == 3
            if sprit.y < height
                sprit.y += 4
            else
                sprit.direction =4
            end
        elseif sprit.direction == 4
            if sprit.y > 0
                sprit.y -= 4
            else
                sprit.direction = 3
            end
        end
    end
    imshow(canvases[1,1], img)
end
# ==============================================================================
# Data structure to hold sprites
# ==============================================================================
mutable struct Sprite
    image           # Image or image class. for example it could be mario but mario could change depending on what he is doing
    x::Int          # Coordinates
    y::Int          # Coordinates
    direction::Int  # direction of movement.
    movement::Int   # Type of movement: walk, jump, float, fall...
    position::Int   # the positioning: centered, above, below...
end
# ==============================================================================

spritelist = []

push!(spritelist, Sprite(cloud, 100, 100, 1, 3, 2) )
push!(spritelist, Sprite(cloud, 100, 600, 1, 3, 2) )
push!(spritelist, Sprite(cloud2, 100, 400, 2, 3, 2) )
push!(spritelist, Sprite(cloud, 150, 200, 1, 3, 2) )

push!(spritelist, Sprite(LeftPlayer, 100, 150, 1, 3, 2) )
push!(spritelist, Sprite(LeftPlayer, 200, 150, 1, 3, 2) )
# mario = Sprite(LeftPlayer, 100, 150, 1, 3, 2)
# push!(spritelist, mario)

# key events in GTK -> gdk.jl
const Left = 0xff51
const Up = 0xff52
const Right = 0xff53
const Down = 0xff54

signal_connect(win, "key-press-event") do widget, event
    if event.keyval == Left
        spritelist[5].direction = 2
    elseif event.keyval == Up
        spritelist[5].direction = 4
    elseif event.keyval == Right
        spritelist[5].direction = 1
    elseif event.keyval == Down
        spritelist[5].direction = 3
    end
    if event.keyval == 97
        spritelist[6].direction = 2
    elseif event.keyval == 119
        spritelist[6].direction = 4
    elseif event.keyval == 100
        spritelist[6].direction = 1
    elseif event.keyval == 115
        spritelist[6].direction = 3
    end
    print(event.keyval)
end
signal_connect(win, "key-release-event") do widget, event
    #c = Char(event.keyval)
 # println("You released key $(c)") #.keyval
end

u1 = Timer(u1 -> Paint(spritelist), 0.0003, 0.001)

max_stars_repo_name: itsKriptunium/MaelEngine

max_issues_repo_name: itsKriptunium/MaelEngine

max_forks_repo_name: itsKriptunium/MaelEngine