text:
""" @ch_struct

Macro for defining CH read/writable structs with the support for
restricting fields by server revision.
"""
macro ch_struct(expr)
    expr isa Expr && expr.head === :struct || error("Invalid usage of @ch_struct")

    struct_name = expr.args[2]
    cleared_block = Expr(:block)
    args = parse_struct_args!(expr.args[3], cleared_block)
    args_names = []

    expr.args[3] = cleared_block
    reader = _make_reader(struct_name, args)
    writer = _make_writer(struct_name, args)
    return Expr(:block, expr, esc(reader), esc(writer))
end

function _make_reader(struct_name, args)

    :(function chread(sock::ClickHouseSock, ::Type{$struct_name})
        $(_field_reader.(args)...)
        return $struct_name($(
                _arg_name.(
                    getindex.(args, 1)
                )...
            )
            )
    end)
end
function _make_writer(struct_name, args)

    :(function chwrite(sock::ClickHouseSock, packet::$struct_name)
        $(_field_writer.(args)...)
    end)

end

function _field_reader(arg_tuple)
    (arg, check) = arg_tuple
    if isnothing(check)
        :($(_arg_name(arg)) = chread(sock, $(_arg_type(arg))))
    else
        :(
            $(_arg_name(arg)) = $(check)(sock.server_rev) ?
                chread(sock, $(_arg_type(arg))) :
                $(_arg_def(arg))
        )
    end
end
function _field_writer(arg_tuple)
    (arg, check) = arg_tuple
    if isnothing(check)
        :(chwrite(sock, packet.$(_arg_name(arg))))
    else
        :(
            if $(check)(sock.server_rev)
                chwrite(sock, packet.$(_arg_name(arg)))
            end
        )
    end
end

_arg_name(arg) = arg.head == :(::) ?
        arg.args[1] :
        arg.args[1].args[1]

_arg_type(arg) = arg.head == :(::) ?
        arg.args[2] :
        arg.args[1].args[2]

_arg_def(arg) = arg.head == :(::) ?
        nothing :
        arg.args[2]

function parse_struct_args!(expr, cleared_expr)
    result = []
    last_rev::Int = 0

    for field in expr.args
        !isa(field, Expr) && continue
        arg::Union{Expr, Nothing} = nothing
        check::Union{Symbol, Nothing} = nothing

        if field.head == :(::)
           arg = field
           push!(cleared_expr.args, arg)
        elseif field.head == :macrocall
            check = Symbol(string(field.args[1])[2:end])
            arg = field.args[3]
            arg.head != :(=) &&
                error("Version restricted fields must have default value")
            arg.head = :kw
            push!(cleared_expr.args, arg.args[1]) #remove macrocall and default value
        else
            error("Invalid usage of @ch_struct")
        end
        push!(result, (
            arg = arg,
            check = check
        ))
    end
    return result
end

macro reg_packet(code, packet_struct)
    result = Expr(:block)
    func = :(packet_struct(::Val{$code}) = $packet_struct)
    push!(result.args, esc(func))
    func = :(packet_code(::Type{$packet_struct}) = $code)
    push!(result.args, esc(func))
    return result
end


text:
"""
    canibuy(vehicle, price, monthly_budget)

Return a message describing if you can afford the monthly payments on a given vehicle.
"""
function canibuy(vehicle, price, monthly_budget)
    if price / 60 < monthly_budget
        "Yes! I'm getting a $vehicle."
    elseif price / 60 > monthly_budget + 10
        "Damn! No $vehicle for me."
    else
        "I'll have to be frugal if I want a $vehicle."
    end
end

"""
    licence(vehicle, kind)

Return a message that describes if `vehicle` requires a licence to operate based on its `kind`.
"""
function licence(vehicle, kind)
    s = kind == "car" ? "a" : "no"
    "The $vehicle requires $s licence to operate."
end

"""
    fee(msrp, age, kind)

Return the registration fee for a given vehicle.
"""
function fee(msrp, age, kind)
    kind == "bike" && return 0
    age >= 10 && return 25

    value = max(msrp, 25000)
    percent = 1 - age / 10

    return value * percent / 100
end

ext: jl


text:
"""
    constrain_objective_value(tolerance)

Limit the objective value to `tolerance`-times the current objective value, as
with [`objective_bounds`](@ref).
"""
constrain_objective_value(tolerance) =
    (_, opt_model) -> begin
        lambda_min, lambda_max = objective_bounds(tolerance)(objective_value(opt_model))
        old_objective = objective_function(opt_model)
        @constraint(opt_model, lambda_min <= sum(old_objective) <= lambda_max)
    end

"""
    change_constraint(id::String; lb=nothing, ub=nothing)

Change the lower and upper bounds (`lb` and `ub` respectively) of reaction `id` if supplied.
"""
change_constraint(id::String; lb = nothing, ub = nothing) =
    (model, opt_model) -> begin
        ind = first(indexin([id], reactions(model)))
        isnothing(ind) && throw(DomainError(id, "No matching reaction was found."))
        set_optmodel_bound!(ind, opt_model, lb = lb, ub = ub)
    end

"""
    change_objective(new_objective::Union{String,Vector{String}}; weights=[], sense=MAX_SENSE)

Modification that changes the objective function used in a constraint based
analysis function.  `new_objective` can be a single reaction identifier, or an
array of reactions identifiers.

Optionally, the objective can be weighted by a vector of `weights`, and a
optimization `sense` can be set to either `MAX_SENSE` or `MIN_SENSE`.
"""
change_objective(
    new_objective::Union{String,Vector{String}};
    weights = [],
    sense = MAX_SENSE,
) =
    (model, opt_model) -> begin

        # Construct objective_indices array
        if typeof(new_objective) == String
            objective_indices = indexin([new_objective], reactions(model))
        else
            objective_indices =
                [first(indexin([rxnid], reactions(model))) for rxnid in new_objective]
        end

        any(isnothing.(objective_indices)) && throw(
            DomainError(new_objective, "No matching reaction found for one or more ids."),
        )

        # Initialize weights
        opt_weights = spzeros(n_reactions(model))

        isempty(weights) && (weights = ones(length(objective_indices))) # equal weights

        for (j, i) in enumerate(objective_indices)
            opt_weights[i] = weights[j]
        end

        v = opt_model[:x]
        @objective(opt_model, sense, sum(opt_weights[i] * v[i] for i in objective_indices))
    end


text:
struct BoxmapTile
    x::Float64
    y::Float64
    w::Float64
    h::Float64
end

"""
    buildrow(A, x, y, w, h)

Make a row of tiles from A.
"""
function buildrow(A, x, y, w, h)
    covered = sum(A)
    width = covered/h
    tiles = BoxmapTile[]
    for value in A
        push!(tiles, BoxmapTile(x, y, width, value/width))
        y += value/width
    end
    return tiles
end

"""
    buildcolumn(A, x, y, w, h)

Make a column of tiles from A.
"""
function buildcolumn(A, x, y, w, h)
    covered = sum(A)
    height = covered/w
    tiles = BoxmapTile[]
    for value in A
        push!(tiles, BoxmapTile(x, y, value/height, height))
        x += value/height
    end
    return tiles
end

"""
    layout(A, x, y, w, h)

From A, make a row of tiles (if wider than tall) or a column of tiles (if taller than wide).
"""
function layout(A, x, y, w, h)
    if w >= h
        return buildrow(A, x, y, w, h)
    else
        return buildcolumn(A, x, y, w, h)
    end
end

function rowlayoutleftover(A, x, y, w, h)
    covered = sum(A)
    width = covered/h
    layoutleftoverx = x + width
    layoutleftoverw = w - width
    return layoutleftoverx, y, layoutleftoverw, h
end

function columnlayoutleftover(A, x, y, w, h)
    covered = sum(A)
    height = covered/w
    layoutleftovery = y + height
    layoutleftoverh = h - height
    return x, layoutleftovery, w, layoutleftoverh
end

function layoutleftover(A, x, y, w, h)
    if w >= h
        return rowlayoutleftover(A, x, y, w, h)
    else
        return columnlayoutleftover(A, x, y, w, h)
    end
end

"""
    highestaspectratio()

Find the highest aspect ratio of a list of rectangles, given the length
of the side along which they are to be laid out.
"""
function highestaspectratio(A, x, y, w, h)
    return maximum([max(tile.w/tile.h, tile.h/tile.w) for tile in layout(A, x, y, w, h)])
end

function buildboxmap(A, x, y, w, h)
    length(A) == 0 && return []
    length(A) == 1 && return layout(A, x, y, w, h)
    i = 1
    while i < length(A) &&
        (highestaspectratio(A[1:i], x, y, w, h) >= highestaspectratio(A[1:i+1], x, y, w, h))
        i += 1
    end
    layoutleftoverx, layoutleftovery, layoutleftoverw, layoutleftoverh = layoutleftover(A[1:i], x, y, w, h)
    return append!(layout(A[1:i], x, y, w, h), buildboxmap(A[i+1:end], layoutleftoverx, layoutleftovery, layoutleftoverw, layoutleftoverh))
end

"""
    boxmap(A::Array, pt, w, h)

Build a box map of the values in `A` with one corner at `pt` and width `w` and
height `h`. There are `length(A)` boxes. The areas of the boxes are proportional
to the original values, scaled as necessary.

The return value is an array of BoxmapTiles. For example:

```
[BoxmapTile(0.0, 0.0, 10.0, 20.0)
 BoxmapTile(10.0, 0.0, 10.0, 13.3333)
 BoxmapTile(10.0, 13.3333, 10.0, 6.66667)]
```

with each tile containing `(x, y, w, h)`. `box()` and `BoundingBox()` can work with
BoxmapTiles as well.

# Example

```
using Luxor
@svg begin
    fontsize(16)
    fontface("HelveticaBold")
    pt = Point(-200, -200)
    a = rand(10:200, 15)
    tiles = boxmap(a, Point(-200, -200), 400, 400)
    for (n, t) in enumerate(tiles)
        randomhue()
        bb = BoundingBox(t)
        box(bb - 2, :stroke)
        box(bb - 5, :fill)
        sethue("white")
        text(string(n), midpoint(bb[1], bb[2]), halign=:center)
    end
end 400 400 "boxmap.svg"
```
"""
function boxmap(A::Array, pt::Point, w, h)
    # algorithm basically from Bruls, Huizing, van Wijk, "Squarified Treemaps"
    # without the silly name
    !all(n -> n > 0, A) && error("all values must be positive")
    sort!(A, rev=true)
    totalvalue = sum(A)
    totalarea = w * h
    normalizedA = A .* (totalarea/totalvalue)
    return buildboxmap(normalizedA, pt.x, pt.y, w, h)
end

"""
    box(tile::BoxmapTile, action::Symbol=:none; vertices=false)

Use a Boxmaptile to make or draw a rectangular box. Use `vertices=true` to obtain
the coordinates.

Create boxmaps using `boxmap()`.
"""
function box(tile::BoxmapTile, action::Symbol=:none; vertices=false)
    if vertices
        return [Point(tile.x,           tile.y + tile.h),
                Point(tile.x,           tile.y),
                Point(tile.x + tile.w,  tile.y),
                Point(tile.x + tile.w,  tile.y + tile.h)]
    end
    rect(tile.x, tile.y, tile.w, tile.h, action)
end

"""
    BoundingBox(tile::BoxmapTile)

Return a BoundingBox of a BoxmapTile (as created with `boxmap()`).
"""
function BoundingBox(tile::BoxmapTile)
    lcorner = Point(tile.x, tile.y)
    ocorner = Point(lcorner.x + tile.w, lcorner.y + tile.h)
    return BoundingBox(lcorner, ocorner)
end


text:
using Statistics, Biomechanics.SpatiotemporalMetrics, Random
using Biomechanics: interleave

Random.seed!(0xBEEF)

stride = 1.3
lstep = 0.6*stride
rstep = stride - lstep

rswing_σ = randn(100) .* .001;
lswing_σ = randn(100) .* .001;
rds_σ = randn(100) .* .001;
lds_σ = randn(100) .* .001;

@testset "Walking" begin
    lstance = 0.70*stride
    rstance = 0.75*stride
    lswing = stride - lstance

    rswing = stride - rstance
    lswings = lswing .+ lswing_σ;
    rswings = rswing .+ rswing_σ;
    ldss = lstep .- lswings .+ lds_σ;
    rdss = rstep .- rswings .+ rds_σ;

    lsteps = lswings + ldss;
    rsteps = rswings + rdss;

    strides = lsteps + rsteps;

    ge = cumsum([0; interleave(ldss, lswings, rdss, rswings)]);

    rfs = ge[1:4:end];
    lfo = ge[2:4:end];
    lfs = ge[3:4:end];
    rfo = ge[4:4:end];

    @test stridetimes(rfs) ≈ strides

    @test collect(steptimes(;lfs, rfs)) ≈ [lsteps, rsteps]

    @test collect(swingstance(rfs, rfo; normalize=false)) ≈ [rswings, strides - rswings]

    @test swing(rfs, rfo; normalize=false) ≈ rswings
    @test stance(rfs, rfo; normalize=false) ≈ strides - rswings

    @test hcat(swingstance(rfs, rfo)...) ≈ [rswings strides - rswings] ./ strides

    @test swing(rfs, rfo) ≈ rswings ./ strides
    @test stance(rfs, rfo) ≈ (strides - rswings) ./ strides

    @test singlesupport(;lfs, lfo, rfs, rfo, normalize=false) ≈ (lswings + rswings)[1:99]
    @test doublesupport(;lfs, lfo, rfs, rfo, normalize=false) ≈ (strides - lswings - rswings)[1:100]

    @test singlesupport(;lfs, lfo, rfs, rfo, normalize=true) ≈
        (lswings + rswings)[1:99] ./ strides[1:99]
    @test doublesupport(;lfs, lfo, rfs, rfo, normalize=true) ≈
        (strides - lswings - rswings)[1:100] ./ strides[1:100]
end

@testset "Running" begin
    lstance = 0.30*stride
    rstance = 0.35*stride
    lswing = stride - lstance
    rswing = stride - rstance
    lstances = lstance .+ lswing_σ;
    rstances = rstance .+ rswing_σ;
    lfloats = lstep .- lstances .+ lds_σ;
    rfloats = rstep .- rstances .+ rds_σ;

    lsteps = rstances + lfloats;
    rsteps = lstances + rfloats;

    strides = lsteps + rsteps;

    ge = cumsum([0; interleave(rstances, lfloats, lstances, rfloats)]);

    rfs = ge[1:4:end];
    rfo = ge[2:4:end];
    lfs = ge[3:4:end];
    lfo = ge[4:4:end];

    @test stridetimes(rfs) ≈ strides

    @test collect(steptimes(;lfs, rfs)) ≈ [lsteps, rsteps]

    @test collect(swingstance(rfs, rfo; normalize=false)) ≈ [strides - rstances, rstances]

    @test swing(rfs, rfo; normalize=false) ≈ strides - rstances
    @test stance(rfs, rfo; normalize=false) ≈ rstances

    @test hcat(swingstance(rfs, rfo)...) ≈ [strides - rstances rstances] ./ strides

    @test swing(rfs, rfo) ≈ (strides - rstances) ./ strides
    @test stance(rfs, rfo) ≈ rstances ./ strides

    @test singlesupport(;lfs, lfo, rfs, rfo, normalize=false) ≈ (lstances + rstances)[1:99]
    @test (doublesupport(;lfs, lfo, rfs, rfo, normalize=false) .== 0) |> prod

    @test singlesupport(;lfs, lfo, rfs, rfo, normalize=true) ≈ (lstances + rstances)[1:99] ./ strides[1:99]
    @test (doublesupport(;lfs, lfo, rfs, rfo, normalize=true) .== 0) |> prod
end


text:
include("SLPDE.jl")
include("etdrk4.jl")
include("evaluate_lambda.jl")
include("sphcesaro.jl")

using PlotlyJS

n = 100

θ = (0.5:n-0.5)/n
φ = (0:2n-2)*2/(2n-1)
x = [cospi(φ)*sinpi(θ) for θ in θ, φ in φ]
y = [sinpi(φ)*sinpi(θ) for θ in θ, φ in φ]
z = [cospi(θ) for θ in θ, φ in φ]

RHSF = -exp.(-30.*((x-1/4).^2+(y-sqrt(11)/4).^2+(z-1/2).^2))-exp.(-50.*z.^2)

RHSV = zero(RHSF)
A_mul_B!(RHSV, FastTransforms.plan_analysis(RHSF), RHSF)
RHSU = fourier2sph(RHSV)

θ = [0;(0.5:n-0.5)/n;1]
φ = [(0:2n-2)*2/(2n-1);2]
x = [cospi(φ)*sinpi(θ) for θ in θ, φ in φ]
y = [sinpi(φ)*sinpi(θ) for θ in θ, φ in φ]
z = [cospi(θ) for θ in θ, φ in φ]

RHSF = [mean(RHSF[1,:])*ones(1, size(RHSF, 2)); RHSF; mean(RHSF[end,:])*ones(1, size(RHSF, 2))]
RHSF = [RHSF RHSF[:,1]]

s = surface(x=x, y=y, z=z, colorscale = "Viridis", surfacecolor = RHSF, cmin = minimum(RHSF), cmax = maximum(RHSF), showscale = false)
ax = attr(visible = false)
cam = attr(up = attr(x=0,y=0,z=1), center = attr(x=0,y=0,z=0), eye = attr(x=0.75,y=0.75,z=0.75))
layout = Layout(width = 500, height = 500, autosize = false, margin = attr(l = 0, r = 0, b = 0, t = 0), scene = attr(xaxis = ax, yaxis = ax, zaxis = ax, camera = cam))
p = plot(s, layout)
savefig(p,"/Users/Mikael/PRHS.pdf";js=:remote)

L = create_linear_operator(LaplaceBeltrami(1.0), size(RHSU, 1), size(RHSU, 2)÷2 + 1)
L[1] = 1

LHSU = FastTransforms.sph_zero_spurious_modes!(RHSU./L)
LHSV = sph2fourier(LHSU)
LHSF = zero(LHSV)
A_mul_B!(LHSF, FastTransforms.plan_synthesis(LHSV), LHSV)

LHSF = [mean(LHSF[1,:])*ones(1, size(LHSF, 2)); LHSF; mean(LHSF[end,:])*ones(1, size(LHSF, 2))]
LHSF = [LHSF LHSF[:,1]]

s = surface(x=x, y=y, z=z, colorscale = "Viridis", surfacecolor = LHSF, cmin = minimum(LHSF), cmax = maximum(LHSF), showscale = false)
ax = attr(visible = false)
cam = attr(up = attr(x=0,y=0,z=1), center = attr(x=0,y=0,z=0), eye = attr(x=0.75,y=0.75,z=0.75))
layout = Layout(width = 500, height = 500, autosize = false, margin = attr(l = 0, r = 0, b = 0, t = 0), scene = attr(xaxis = ax, yaxis = ax, zaxis = ax, camera = cam))
p = plot(s, layout)
savefig(p,"/Users/Mikael/LPLHS.pdf";js=:remote)

NL = create_linear_operator(NonlocalLaplaceBeltrami(1.0, 0.0, 1.5), size(RHSU, 1), size(RHSU, 2)÷2 + 1)
NL[1] = 1

LHSU = FastTransforms.sph_zero_spurious_modes!(RHSU./NL)
LHSV = sph2fourier(LHSU)
LHSF = zero(LHSV)
A_mul_B!(LHSF, FastTransforms.plan_synthesis(LHSV), LHSV)

LHSF = [mean(LHSF[1,:])*ones(1, size(LHSF, 2)); LHSF; mean(LHSF[end,:])*ones(1, size(LHSF, 2))]
LHSF = [LHSF LHSF[:,1]]

s = surface(x=x, y=y, z=z, colorscale = "Viridis", surfacecolor = LHSF, cmin = minimum(LHSF), cmax = maximum(LHSF), showscale = false)
ax = attr(visible = false)
cam = attr(up = attr(x=0,y=0,z=1), center = attr(x=0,y=0,z=0), eye = attr(x=0.75,y=0.75,z=0.75))
layout = Layout(width = 500, height = 500, autosize = false, margin = attr(l = 0, r = 0, b = 0, t = 0), scene = attr(xaxis = ax, yaxis = ax, zaxis = ax, camera = cam))
p = plot(s, layout)
savefig(p,"/Users/Mikael/NLPLHS.pdf";js=:remote)


text:
function op_Vpot( pw::PWGrid, Vpot, psi::Array{ComplexF64,2}, ik::Int )
  #
  Ns = pw.Ns
  Ω  = pw.Ω
  Npoints = prod(Ns)
  Ncols   = size(psi)[2]

  ctmp = zeros(ComplexF64, Npoints, Ncols)
  idx = pw.gkvec.idx_gk[ik]
  for ic = 1:Ncols
    ctmp[idx,ic] = psi[:,ic]
  end

  # get values of psi in real space grid via forward transform
  ctmp = G_to_R( Ns, ctmp )

  cVpsi = R_to_G( Ns, Diagprod(Vpot, ctmp) )
  return cVpsi[idx,:]
end

# B is usually consists of more than one-column
function Diagprod( a,B )
  Ncol    = size(B)[2]
  Npoints = size(B)[1]
  out = zeros( ComplexF64, size(B) )
  for ic = 1:Ncol
    for ip = 1:Npoints
      out[ip,ic] = a[ip]*B[ip,ic]
    end
  end
  return out
end


text:
ThreadStats() = ThreadStats(0, 0, 0, 0)

Thread() = Thread(TimeManagement(),
    Board(),
    [MoveStack(MAX_PLY + 1) for i in 0:MAX_PLY+1],
    ThreadStats(),
    [MoveOrder() for i in 0:MAX_PLY+1],
    MoveStack(MAX_PLY + 1),
    PieceStack(MAX_PLY + 1),
    zeros(Int, MAX_PLY + 1),
    [MoveStack(MAX_QUIET_TRACK) for i in 0:MAX_PLY + 1],
    ButterflyHistTable(),
    CounterHistTable(),
    CounterHistTable(),
    MoveStack(MAX_PLY + 1),
    MoveStack(MAX_PLY + 1),
    CounterTable(),
    PawnKingTable(),
    false)

const ThreadPool = Vector{Thread}
ThreadPool(n::Int) = [Thread() for i in 1:n]

# we can't increase the number of threads in julia while the process is running...
# so we set an atexit condition to reload the engine & process from scratch
function createthreadpool(board::Board, num_threads::Int)
    if Threads.nthreads() !== num_threads
        fen = exportfen(board)
        atexit(() -> withenv(() -> run(`julia AstellarnEngine.jl $fen "info string set Threads to $num_threads"`), "JULIA_NUM_THREADS" => num_threads))
        exit()
    end
end

function setthreadpoolboard!(threads::ThreadPool, board::Board)
    for thread in threads
        copy!(thread.board, board)
    end
end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
@testset "gradient" begin

N = [4,5]
M = prod(N)
r = reshape(collect(1:2M), M, 2)'[:]  # index mapping from block matrix to narrowly banded matrix
Z = spzeros(M,M)

G = rand(ComplexF64, N..., 2)  # output array of vector field; make sure not to be confused with gradient operator
f = zeros(ComplexF64, N...)  # input array of scalar
g = zeros(ComplexF64, 2M)  # column vector representation of G

@testset "create_grad and apply_grad!" begin
    for ci = CartesianIndices((false:true,false:true))
        # Construct Grad for a uniform grid and Bloch boundaries.
        isfwd = Vector{Bool}([ci.I...])
        Grad = create_grad(isfwd, N, order_cmpfirst=false)

        # Test the overall coefficients.
        @test size(Grad) == (2M,M)
        @test all(any(Grad.≠0, dims=1))  # no zero columns
        @test all(any(Grad.≠0, dims=2))  # no zero rows
        @test all(sum(Grad, dims=1) .== 0)  # all column sums are zero, because each input field to Grad is used twice in each Cartesian direction, once multiplied with +1 and once with -1
        @test all(sum(Grad, dims=2) .== 0)  # all row sums are zero, because Grad * ones(M) = 0
        @test all(sum(abs.(Grad), dims=1) .== 4)  # each column of Grad has six nonzero entries, which are ±1's
        @test all(sum(abs.(Grad), dims=2) .== 2)  # each row of Grad has two nonzero entries, which are ±1's

        ∂x = (nw = 1; create_∂(nw, isfwd[nw], N))
        ∂y = (nw = 2; create_∂(nw, isfwd[nw], N))
        @test Grad == [∂x; ∂y]

        # Construct Grad for a nonuniform grid and general boundaries.
        ∆l⁻¹ = rand.(tuple(N...))  # isfwd = true (false) uses ∆l⁻¹ at dual (primal) locations
        isbloch = [true, false]
        e⁻ⁱᵏᴸ = rand(ComplexF64, 2)

        Grad = create_grad(isfwd, ∆l⁻¹, isbloch, e⁻ⁱᵏᴸ, order_cmpfirst=false)

        # Test Grad.
        ∂x = (nw = 1; create_∂(nw, isfwd[nw], N, ∆l⁻¹[nw], isbloch[nw], e⁻ⁱᵏᴸ[nw]))
        ∂y = (nw = 2; create_∂(nw, isfwd[nw], N, ∆l⁻¹[nw], isbloch[nw], e⁻ⁱᵏᴸ[nw]))
        @test Grad == [∂x; ∂y]

        # Test Cartesian-component-first ordering.
        Grad_cmpfirst = create_grad(isfwd, ∆l⁻¹, isbloch, e⁻ⁱᵏᴸ, order_cmpfirst=true)
        @test Grad_cmpfirst == Grad[r,:]

        # Test apply_grad!.
        fvec = f[:]
        mul!(g, Grad, fvec)
        apply_grad!(G, f, Val(:(=)), isfwd, ∆l⁻¹, isbloch, e⁻ⁱᵏᴸ)  # no keyword argument order_cmpfirst in apply_XXX!()
        @test G[:] ≈ g

        mul!(g, Grad_cmpfirst, fvec)
        apply_grad!(G, f, Val(:(=)), isfwd, ∆l⁻¹, isbloch, e⁻ⁱᵏᴸ)  # no keyword argument order_cmpfirst in apply_XXX!()
        @test G[:][r] ≈ g
    end
end  # @testset "create_grad and apply_grad!"

@testset "curl of gradient" begin
    # Construct Curl and Grad for a uniform grid and periodic boundaries.
    N = [3,4,5]
    M = prod(N)
    isfwd = [true, true, true]  # curl(U) and gradient to generate U are differentiated forward

    Curl = create_curl(isfwd, N, order_cmpfirst=false)
    Grad = create_grad(isfwd, N, order_cmpfirst=false)

    # Construct Curl * Grad
    A = Curl * Grad

    # Test Divergence of curl.
    @test size(A) == (3M,M)
    @test all(A.==0)
end  # @testset "curl of gradent"

end  # @testset "gradient"


text:
struct TimeSeriesData <: InfrastructureSystemsType
    data::TimeSeries.TimeArray
    internal::InfrastructureSystemsInternal

    function TimeSeriesData(data, internal)
        if length(data) < 2
            throw(ArgumentError("time array length must be at least 2"))
        end

        return new(data, internal)
    end
end

function TimeSeriesData(data)
    return TimeSeriesData(data, InfrastructureSystemsInternal())
end

function Base.summary(data::TimeSeriesData)
    return "TimeSeriesData"
end

function Base.show(io::IO, ::MIME"text/plain", ts::TimeSeriesData)
    println(io, "UUID=$(get_uuid(ts))")
    println(io, "data=$(ts.data)")
end

Base.length(ts::TimeSeriesData) = length(ts.data)
get_initial_time(ts::TimeSeriesData) = TimeSeries.timestamp(ts.data)[1]
get_horizon(ts::TimeSeriesData) = length(ts.data)
get_resolution(ts::TimeSeriesData) = TimeSeries.timestamp(ts.data)[2] -
                                     TimeSeries.timestamp(ts.data)[1]

#function get_time_series(data::TimeSeriesData)
#    if !isnothing(data)
#        return data[Symbol(data.component_name)]
#    end
#
#    ta = read_time_series(data.file_path)
#    return ta[Symbol(data.component_name)]
#end


text:
abstract type AbstractPotentialSimulationSetup{T, N} end

# """
#     PotentialSimulationSetup{T, N, S} <: AbstractPotentialSimulationSetup{T, N}
#
# Collection struct. It holds the grid, the potential, the point types, the charge density and the dielectric distribution.
# """
struct PotentialSimulationSetup{T, N, S} <: AbstractPotentialSimulationSetup{T, N}
    grid::Grid{T, N, S}
    potential::Array{T, N}
    point_types::Array{PointType, N}
    q_eff_imp::Array{T, N}
    q_eff_fix::Array{T, N}
    ϵ_r::Array{T, N}
end

struct PotentialSimulationSetupRB{
            T, S, N1, DATN1<:AbstractArray{T,N1},
            N2, AT,
            DATN2<:AbstractArray{T,N2},
            DATPT<:AbstractArray{PointType,N2},
            DATGW<:AbstractArray{T,2},
            DATSOR<:AbstractArray{T,1}
        } <: AbstractPotentialSimulationSetup{T, N1}
    grid::Grid{T, N1, S, AT}
    potential::DATN2 # Array{T, N2} or e.g. CuArray{T, N2}
    point_types::DATPT # Array{PointType, N2} or ...
    volume_weights::DATN2 # Array{T, N2}
    q_eff_imp::DATN2 # Array{T, N2}
    q_eff_fix::DATN2 # Array{T, N2}
    ϵ_r::DATN1 # Array{T, N1}
    geom_weights::NTuple{3, DATGW} # NTuple{3, Array{T, 2}}
    sor_const::DATSOR # Array{T, 1}
    bias_voltage::T
    maximum_applied_potential::T
    minimum_applied_potential::T
    depletion_handling_potential_limit::T
    grid_boundary_factors::NTuple{3, NTuple{2, T}}
end

function Adapt.adapt_structure(to, pssrb::PotentialSimulationSetupRB{T, S, 3}) where {T, S}
    PotentialSimulationSetupRB(
        pssrb.grid,
        adapt(to, pssrb.potential),
        adapt(to, pssrb.point_types),
        adapt(to, pssrb.volume_weights),
        adapt(to, pssrb.q_eff_imp),
        adapt(to, pssrb.q_eff_fix),
        adapt(to, pssrb.ϵ_r),
        adapt(to, pssrb.geom_weights),
        adapt(to, pssrb.sor_const),
        pssrb.bias_voltage,
        pssrb.maximum_applied_potential,
        pssrb.minimum_applied_potential,
        pssrb.depletion_handling_potential_limit,
        pssrb.grid_boundary_factors
    )
end

include("BoundaryConditions/BoundaryConditions.jl")
include("PotentialSimulationSetupRBCylindrical.jl")
include("PotentialSimulationSetupRBCartesian3D.jl")


text:
using Test
using DataWrangler

@testset "boxcox" begin

    x = rand(100)
    bc = boxcox(x)
    @test iboxcox(bc[:x],bc[:λ]) ≈ x

    x = rand(100) .- 0.5
    bc = boxcox(x)
    @test iboxcox(bc[:x],bc[:λ]) ≈ x

    x = rand(100)
    bc = boxcox(x)
    @test iboxcox(bc[:x],bc[:λ]) ≈ x

    x =  Int64.(round.(100*rand(100)))
    bc = boxcox(x)
    @test iboxcox(bc[:x],bc[:λ]) ≈ x

    x =  Int64.(round.(100*rand(100) .- 50))
    bc = boxcox(x)
    @test iboxcox(bc[:x],bc[:λ]) ≈ x

end


text:
#
# types.jl --
#
# Definitions fo types for IPC module of Julia.
#
#------------------------------------------------------------------------------
#
# This file is part of IPC.jl released under the MIT "expat" license.
# Copyright (C) 2016-2018, Éric Thiébaut (https://github.com/emmt/IPC.jl).
#

"""
`TimeoutError` is used to throw a timeout exception.
"""
struct TimeoutError <: Exception; end

struct TimeVal
    sec::_typeof_timeval_sec    # seconds
    usec::_typeof_timeval_usec  # microseconds
end

struct TimeSpec
    sec::_typeof_timespec_sec    # seconds
    nsec::_typeof_timespec_nsec  # nanoseconds
end

struct ProcessId
    value::_typeof_pid_t
end

struct UserId
    value::_typeof_uid_t
end

# Counterpart of C `sigset_t` structure.  Must be mutable to avoid wrapping it
# in a Ref when passed by address.
mutable struct SigSet
    bits::_typeof_sigset
    SigSet() = fill!(new(), false)
end

# Counterpart of C `struct sigaction` structure.  Must be mutable to avoid
# wrapping it in a Ref when passed by address.
mutable struct SigAction
    handler::Ptr{Cvoid}
    mask::SigSet
    flags::_typeof_sigaction_flags
end

# Counterpart of C `siginfo_t` structure; although only a pointer of this is
# ever used by signal handlers.  Must be mutable to avoid wrapping it in a Ref
# when passed by address.
mutable struct SigInfo
    bits::_typeof_siginfo
    function SigInfo()
        obj = new()
        ccall(:memset, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t),
              pointer_from_objref(obj), 0, sizeof(obj))
        return obj
    end
end

struct Key
    value::_typeof_key_t
    Key(value::Integer) = new(value)
end

struct ShmId
    value::Cint
    ShmId(value::Integer) = new(value)
end

mutable struct ShmInfo
    atime::UInt64 # last attach time
    dtime::UInt64 # last detach time
    ctime::UInt64 # last change time
    segsz::UInt64 # size of the public area
    id::Int32     # shared memory identifier
    cpid::Int32   # process ID of creator
    lpid::Int32   # process ID of last operator
    nattch::Int32 # no. of current attaches
    mode::UInt32  # lower 9 bits of access modes
    uid::UInt32   # effective user ID of owner
    gid::UInt32   # effective group ID of owner
    cuid::UInt32  # effective user ID of creator
    cgid::UInt32  # effective group ID of creator
    ShmInfo() = new(0,0,0,0,0,0,0,0,0,0,0,0,0)
end

mutable struct SemInfo
    otime::UInt64 # last semop time
    ctime::UInt64 # last change time
    nsems::Int32  # number of semaphores in set
    id::Int32     # semaphore set identifier
    uid::UInt32   # effective user ID of owner
    gid::UInt32   # effective group ID of owner
    cuid::UInt32  # effective user ID of creator
    cgid::UInt32  # effective group ID of creator
    mode::UInt32  # lower 9 bits of access modes
    SemInfo() = new(0,0,0,0,0,0,0,0,0)
end

# Must be mutable to allow for finalizing.
mutable struct FileDescriptor
    fd::Cint
end

abstract type MemoryBlock end

mutable struct DynamicMemory <: MemoryBlock
    ptr::Ptr{Cvoid}
    len::Int
    function DynamicMemory(len::Integer)
        @assert len ≥ 1
        ptr = Libc.malloc(len)
        ptr != C_NULL || throw(OutOfMemoryError())
        return finalizer(_destroy, new(ptr, len))
    end
end

mutable struct SharedMemory{T<:Union{String,ShmId}} <: MemoryBlock
    ptr::Ptr{Cvoid} # mapped address of shared memory segment
    len::Int        # size of shared memory segment (in bytes)
    volatile::Bool  # true if shared memory is volatile (only for the creator)
    id::T           # identifier of shared memory segment
    function SharedMemory{T}(ptr::Ptr{Cvoid},
                             len::Integer,
                             volatile::Bool,
                             id::T) where {T<:Union{String,ShmId}}
        return finalizer(_destroy, new(ptr, len, volatile, id))
    end
end

struct WrappedArray{T,N,M} <: DenseArray{T,N}
    # All members shall be considered as private.
    arr::Array{T,N}  # wrapped Julia array
    mem::M           # object providing the memory
    function WrappedArray{T,N,M}(ptr::Ptr{T},
                                 dims::NTuple{N,<:Integer},
                                 mem::M) where {T,N,M}
        arr = unsafe_wrap(Array, ptr, dims)
        return new{T,N,M}(arr, mem)
    end
end

const WrappedVector{T,M} = WrappedArray{T,1,M}
const WrappedMatrix{T,M} = WrappedArray{T,2,M}

const ShmArray{T,N,M<:SharedMemory} = WrappedArray{T,N,M}
const ShmVector{T,M} = ShmArray{T,1,M}
const ShmMatrix{T,M} = ShmArray{T,2,M}

# Header for saving a minimal description of a wrapped array.  The layout is:
#
#   Name   Size
#   --------------
#   magic  4 bytes
#   etype  2 bytes
#   ndims  2 bytes
#   offset 8 bytes
#
# This header is supposed to be directly followed by the array dimensions
# stored as 8-byte signed integers.
#
struct WrappedArrayHeader
    magic::UInt32 # magic number to check correctness
    etype::UInt16 # identifier of element type
    ndims::UInt16 # number of dimensions
    offset::Int64 # total size of header
end
@assert rem(sizeof(WrappedArrayHeader), sizeof(Int64)) == 0

mutable struct Semaphore{T}
    ptr::Ptr{Cvoid}
    lnk::T
    # Provide inner constructor to force fully qualified calls.
    Semaphore{T}(ptr::Ptr, lnk) where {T} = new{T}(ptr, lnk)
end


text:
"""
    @test_reference filename expr [by] [kw...]

Tests that the expression `expr` with reference `filename` using
equality test strategy `by`.

The pipeline of `test_reference` is:

1. preprocess `expr`
2. read and preprocess `filename` via FileIO.jl
3. compare the results using `by`
4. if test fails in an interactive session (e.g, `include(test/runtests.jl)`), an interactive dialog will be trigered.

Arguments:

* `filename::String`: _relative_ path to the file that contains the macro invocation.
* `expr`: the actual content used to compare.
* `by`: the equality test function. By default it is `isequal` if not explicitly stated.
* `format`: Force reading the file using a specific format

# Types

The file-extension of `filename`, as well as the type of the
result of evaluating `expr`, determine how contents are processed and
compared. The contents is treated as:

* Images when `expr` is an image type, i.e., `AbstractArray{<:Colorant}`;
* SHA256 when `filename` endswith `*.sha256`;
* Any file-type which FileIO.jl handles and with the proper backend loaded;
* Text as a fallback.

## Any FileIO.jl handled filetype

Any file-types which can be read by [FileIO.jl](https://github.com/JuliaIO/FileIO.jl) can be used.
Note that most Julia values can be stored using packages such as, e.g.,
[BSON.jl](https://github.com/JuliaIO/BSON.jl) or [JLD](https://github.com/JuliaIO/JLD.jl) and can thus be used in reference-tests.

## Images

Images are compared _approximately_ using a different `by` to ignore most encoding
and decoding errors. The default function is generated from [`psnr_equality`](@ref).

The file can be either common image files (e.g., `*.png`), which are handled by
`FileIO`, or text-coded `*.txt` files, which is handled by `ImageInTerminal`.
Text-coded image has less storage requirements and also allows to view the
reference file in a simple terminal using `cat`.

## SHA256

The hash of the `expr` and content of `filename` are compared.

!!! tip

    This is useful for a convenient low-storage way of making
    sure that the return value doesn't change for selected test
    cases.

## Fallback

Simply test the equality of `expr` and the contents of `filename` without any preprocessing.
Note that reading `filename` will return a `String`.

# Examples

```julia
# compare text-file against string representation of a value
@test_reference "stringtest1.txt" collect(1:20)

# test number with absolute tolerance 10
@test_reference "references/string3.txt" 1338 by=(ref, x)->isapprox(ref, x; atol=10)

# store a floating point array ar in BSON-file and compare it back.
# Note that a Dict is needed and the custom comparison function.
using BSON
comp(d1, d2) = keys(d1)==keys(d2) &&
    all([ v1≈v2 for (v1,v2) in zip(values(d1), values(d2))])
@test_reference "reftest-files/X.bson" Dict(:ar=>[1, pi, 4.5]) by=comp

# store as string using ImageInTerminal with encoding size (5,10)
using TestImages
@test_reference "camera.txt" testimage("cameraman") size=(5,10)

# using folders in the relative path is allowed
@test_reference "references/camera.png" testimage("cameraman")

# Images can also be stored as hash. Note however that this
# can only check for equality (no tolerance possible)
@test_reference "references/camera.sha256" testimage("cameraman")

# test images with custom Peak Signal-to-Noise Ratio (psnr) threshold
@test_reference "references/camera.png" testimage("cameraman") by=psnr_equality(20)
```
"""
macro test_reference(reference, actual, kws...)
    dir = Base.source_dir()
    expr = :(test_reference(abspath(joinpath($dir, $(esc(reference)))), $(esc(actual))))
    for kw in kws
        (kw isa Expr && kw.head == :(=)) || error("invalid signature for @test_reference")
        k, v = kw.args
        push!(expr.args, Expr(:kw, k, esc(v)))
    end
    expr
end

function test_reference(
    filename::AbstractString, raw_actual;
    by = nothing, render = nothing, format = nothing, kw...)
    format isa AbstractString && (format = FileIO.DataFormat{Symbol(format)})
    reference_file = format === nothing ? query_extended(filename) : File{format}(filename)
    test_reference(reference_file, raw_actual, by, render; kw...)
end

function test_reference(
    reference_file::File{F},
    raw_actual::T,
    equiv=nothing,
    rendermode=nothing;
    kw...) where {F <: DataFormat, T}

    reference_path = reference_file.filename
    reference_dir, reference_filename = splitdir(reference_path)

    actual = _convert(F, raw_actual; kw...)

    # infer the default rendermode here
    # since `nothing` is always passed to this method from
    # test_reference(filename::AbstractString, raw_actual; kw...)
    if rendermode === nothing
        rendermode = default_rendermode(F, actual)
    end

    if !isfile(reference_path)  # when reference file doesn't exists
        mkpath(reference_dir)
        savefile(reference_file, actual)
        @info(
            "Reference file for \"$reference_filename\" did not exist. It has been created",
            new_reference=reference_path,
        )

        # TODO: move encoding out from render
        render(rendermode, actual)

        @info("Please run the tests again for any changes to take effect")
        return nothing # skip current test case
    end

    # file exists
    reference = loadfile(typeof(actual), reference_file)

    if equiv === nothing
        # generally, `reference` and `actual` are of the same type after preprocessing
        equiv = default_equality(reference, actual)
    end

    if equiv(reference, actual)
        @test true # to increase test counter if reached
    else  # When test fails
        # Saving actual file so user can look at it
        actual_path = joinpath(mismatch_staging_dir(), reference_filename)
        actual_file = typeof(reference_file)(actual_path)
        savefile(actual_file, actual)

        # Report to user.
        @info(
            "Reference Test for \"$reference_filename\" failed.",
            reference=reference_path,
            actual=actual_path,
        )
        render(rendermode, reference, actual)

        if !isinteractive() && !force_update()
            error("""
            To update the reference images either run the tests interactively with 'include(\"test/runtests.jl\")',
            or to force-update all failing reference images set the environment variable `JULIA_REFERENCETESTS_UPDATE`
            to "true" and re-run the tests via Pkg.
            """)
        end

        if force_update() || input_bool("Replace reference with actual result?")
            mv(actual_path, reference_path; force=true)  # overwrite old file it
            @info("Please run the tests again for any changes to take effect")
        else
            @test false
        end
    end
end

force_update() = tryparse(Bool, get(ENV, "JULIA_REFERENCETESTS_UPDATE", "false")) === true

"""
    mismatch_staging_dir()

The directory where we store files that don't match so user can look at them.
If the enviroment variable `REFERENCE_TESTS_STAGING_PATH` is set then we will use that directory.
If not a temporary directory will be created. Note that this temporary directory will not be
deleted when julia exists. Since in that case the files would may be deleted before you can
look at them. You should use your operating systems standard mechanisms to clean up excess
temporary directories.
"""
function mismatch_staging_dir()
    return mkpath(expanduser(get(ENV,
        "REFERENCE_TESTS_STAGING_PATH",
        VERSION >= v"1.3" ? mktempdir(; cleanup=false) : mktempdir()
    )))
end


text:
VALIDDISTRIBUTIONS = [i for i in names(Distributions) if uppercasefirst(string(i)) == string(i)]
VALIDDISTRIBUTIONSCUDA = [:Normal,:Uniform]

#Random distribution transformations for cuda capabilities
NormalCUDA(x,μ,σ) = σ*CUDA.sqrt(2.)*CUDA.erfinv(2*(x-.5))+μ
UniformCUDA(x,l0,l1) = (l1-l0)*x+l0

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
"""
    flux_em!(femL, femR, ELL, BLL, EL, BL, ER, BR, ERR, BRR, ϕL, ϕR, ψL, ψR, dxL, dxR, Ap, An, D, sol, χ, ν, dt)

Wave propagation method for Maxwell's equations

* @args: {E, B, ϕ, ψ} in left-left, left, right, and right-right cells
* @args: eigenmatrix (A -> A+ & A-), eigenvalue (D)
* @args: full size of left & right cells
* @args: speed of light (sol)
* @args: auxiliary parameters (χₑ, νᵦ)

"""
function flux_em!(
    femL::T1,
    femR::T1,
    ELL::T2,
    BLL::T2,
    EL::T2,
    BL::T2,
    ER::T2,
    BR::T2,
    ERR::T2,
    BRR::T2,
    ϕL,
    ϕR,
    ψL,
    ψR,
    dxL, # NO need to multiply 0.5
    dxR, # NO need to multiply 0.5
    Ap::T3,
    An::T3,
    D::T4,
    sol,
    χ,
    ν,
    dt,
) where {
    T1<:AbstractArray{<:AbstractFloat,1},
    T2<:AbstractArray{<:Real,1},
    T3<:AbstractArray{<:Real,2},
    T4<:AbstractArray{<:Real,1},
}

    @assert length(femL) == length(femR) == 8

    slop = zeros(8, 8)
    slop[3, 1] = -0.5 * sol^2 * (BR[2] - BL[2]) + 0.5 * sol * (ER[3] - EL[3])
    slop[5, 1] = 0.5 * sol * (BR[2] - BL[2]) - 0.5 * (ER[3] - EL[3])
    slop[2, 2] = 0.5 * sol^2 * (BR[3] - BL[3]) + 0.5 * sol * (ER[2] - EL[2])
    slop[6, 2] = 0.5 * sol * (BR[3] - BL[3]) + 0.5 * (ER[2] - EL[2])
    slop[1, 3] = 0.5 * sol * χ * (ER[1] - EL[1])
    slop[7, 3] = 0.5 * χ * (ER[1] - EL[1])
    slop[4, 4] = 0.0
    slop[8, 4] = 0.0
    slop[3, 5] = -0.5 * sol^2 * (BR[2] - BL[2]) - 0.5 * sol * (ER[3] - EL[3])
    slop[5, 5] = -0.5 * sol * (BR[2] - BL[2]) - 0.5 * (ER[3] - EL[3])
    slop[2, 6] = 0.5 * sol^2 * (BR[3] - BL[3]) - 0.5 * sol * (ER[2] - EL[2])
    slop[6, 6] = -0.5 * sol * (BR[3] - BL[3]) + 0.5 * (ER[2] - EL[2])
    slop[1, 7] = -0.5 * sol * χ * (ER[1] - EL[1])
    slop[7, 7] = 0.5 * χ * (ER[1] - EL[1])
    slop[4, 8] = 0.0
    slop[8, 8] = 0.0

    limiter = zeros(8, 8)
    limiter[3, 1] = -0.5 * sol^2 * (BL[2] - BLL[2]) + 0.5 * sol * (EL[3] - ELL[3])
    limiter[5, 1] = 0.5 * sol * (BL[2] - BLL[2]) - 0.5 * (EL[3] - ELL[3])
    limiter[2, 2] = 0.5 * sol^2 * (BL[3] - BLL[3]) + 0.5 * sol * (EL[2] - ELL[2])
    limiter[6, 2] = 0.5 * sol * (BL[3] - BLL[3]) + 0.5 * (EL[2] - ELL[2])
    limiter[1, 3] = 0.5 * sol * χ * (EL[1] - ELL[1])
    limiter[7, 3] = 0.5 * χ * (EL[1] - ELL[1])
    limiter[4, 4] = 0.0
    limiter[8, 4] = 0.0
    limiter[3, 5] = -0.5 * sol^2 * (BRR[2] - BR[2]) - 0.5 * sol * (ERR[3] - ER[3])
    limiter[5, 5] = -0.5 * sol * (BRR[2] - BR[2]) - 0.5 * (ERR[3] - ER[3])
    limiter[2, 6] = 0.5 * sol^2 * (BRR[3] - BR[3]) - 0.5 * sol * (ERR[2] - ER[2])
    limiter[6, 6] = -0.5 * sol * (BRR[3] - BR[3]) + 0.5 * (ERR[2] - ER[2])
    limiter[1, 7] = -0.5 * sol * χ * (ERR[1] - ER[1])
    limiter[7, 7] = 0.5 * χ * (ERR[1] - ER[1])
    limiter[4, 8] = 0.0
    limiter[8, 8] = 0.0

    for i = 1:8
        limiter_theta = sum(slop[:, i] .* limiter[:, i]) / (sum(slop[:, i] .^ 2) + 1.e-7)
        slop[:, i] .*=
            max(0.0, min(min((1.0 + limiter_theta) / 2.0, 2.0), 2.0 * limiter_theta))
    end

    for i = 1:8
        femL[i] =
            sum(An[i, 1:3] .* (ER .- EL)) +
            sum(An[i, 4:6] .* (BR .- BL)) +
            An[i, 7] * (ϕR - ϕL) +
            An[i, 8] * (ψR - ψL) +
            0.5 * sum(
                fortsign.(1.0, D) .* (1.0 .- dt ./ (0.5 * (dxL + dxR)) .* abs.(D)) .*
                slop[i, :],
            )
        femR[i] =
            sum(Ap[i, 1:3] .* (ER .- EL)) +
            sum(Ap[i, 4:6] .* (BR .- BL)) +
            Ap[i, 7] * (ϕR - ϕL) +
            Ap[i, 8] * (ψR - ψL) -
            0.5 * sum(
                fortsign.(1.0, D) .* (1.0 .- dt ./ (0.5 * (dxL + dxR)) .* abs.(D)) .*
                slop[i, :],
            )
    end

    return nothing

end

"""
    flux_emx!(femL, femR, femLU, femLD, femRU, femRD,
    ELL, BLL, EL, BL, ER, BR, ERR, BRR, ϕL, ϕR, ψL, ψR,
    dxL, dxR, A1p, A1n, A2p, A2n, D, sol, χ, ν, dt)

Wave propagation method for 2D Maxwell's equations

* @args: {E, B, ϕ, ψ} in left-left, left, right, and right-right cells
* @args: eigenmatrix (A -> A+ & A-), eigenvalue (D)
* @args: full size of left & right cells
* @args: speed of light (sol)
* @args: auxiliary parameters (χₑ, νᵦ)

"""
function flux_emx!(
    femL::X,
    femR::X,
    femLU::X,
    femLD::X,
    femRU::X,
    femRD::X,
    ELL::Y,
    BLL::Y,
    EL::Y,
    BL::Y,
    ER::Y,
    BR::Y,
    ERR::Y,
    BRR::Y,
    ϕL,
    ϕR,
    ψL,
    ψR,
    dxL, # NO need to multiply 0.5
    dxR, # NO need to multiply 0.5
    A1p::A,
    A1n::A,
    A2p::A,
    A2n::A,
    D::B,
    sol,
    χ,
    ν,
    dt,
) where {
    X<:AbstractArray{<:AbstractFloat,1},
    Y<:AbstractArray{<:Real,1},
    A<:AbstractArray{<:Real,2},
    B<:AbstractArray{<:Real,1},
}

    slop = zeros(8, 8)
    slop[3, 1] = -0.5 * sol^2 * (BR[2] - BL[2]) + 0.5 * sol * (ER[3] - EL[3])
    slop[5, 1] = 0.5 * sol * (BR[2] - BL[2]) - 0.5 * (ER[3] - EL[3])
    slop[2, 2] = 0.5 * sol^2 * (BR[3] - BL[3]) + 0.5 * sol * (ER[2] - EL[2])
    slop[6, 2] = 0.5 * sol * (BR[3] - BL[3]) + 0.5 * (ER[2] - EL[2])
    slop[1, 3] = 0.5 * sol * χ * (ER[1] - EL[1])
    slop[7, 3] = 0.5 * χ * (ER[1] - EL[1])
    slop[4, 4] = 0.0
    slop[8, 4] = 0.0
    slop[3, 5] = -0.5 * sol^2 * (BR[2] - BL[2]) - 0.5 * sol * (ER[3] - EL[3])
    slop[5, 5] = -0.5 * sol * (BR[2] - BL[2]) - 0.5 * (ER[3] - EL[3])
    slop[2, 6] = 0.5 * sol^2 * (BR[3] - BL[3]) - 0.5 * sol * (ER[2] - EL[2])
    slop[6, 6] = -0.5 * sol * (BR[3] - BL[3]) + 0.5 * (ER[2] - EL[2])
    slop[1, 7] = -0.5 * sol * χ * (ER[1] - EL[1])
    slop[7, 7] = 0.5 * χ * (ER[1] - EL[1])
    slop[4, 8] = 0.0
    slop[8, 8] = 0.0

    limiter = zeros(8, 8)
    limiter[3, 1] = -0.5 * sol^2 * (BL[2] - BLL[2]) + 0.5 * sol * (EL[3] - ELL[3])
    limiter[5, 1] = 0.5 * sol * (BL[2] - BLL[2]) - 0.5 * (EL[3] - ELL[3])
    limiter[2, 2] = 0.5 * sol^2 * (BL[3] - BLL[3]) + 0.5 * sol * (EL[2] - ELL[2])
    limiter[6, 2] = 0.5 * sol * (BL[3] - BLL[3]) + 0.5 * (EL[2] - ELL[2])
    limiter[1, 3] = 0.5 * sol * χ * (EL[1] - ELL[1])
    limiter[7, 3] = 0.5 * χ * (EL[1] - ELL[1])
    limiter[4, 4] = 0.0
    limiter[8, 4] = 0.0
    limiter[3, 5] = -0.5 * sol^2 * (BRR[2] - BR[2]) - 0.5 * sol * (ERR[3] - ER[3])
    limiter[5, 5] = -0.5 * sol * (BRR[2] - BR[2]) - 0.5 * (ERR[3] - ER[3])
    limiter[2, 6] = 0.5 * sol^2 * (BRR[3] - BR[3]) - 0.5 * sol * (ERR[2] - ER[2])
    limiter[6, 6] = -0.5 * sol * (BRR[3] - BR[3]) + 0.5 * (ERR[2] - ER[2])
    limiter[1, 7] = -0.5 * sol * χ * (ERR[1] - ER[1])
    limiter[7, 7] = 0.5 * χ * (ERR[1] - ER[1])
    limiter[4, 8] = 0.0
    limiter[8, 8] = 0.0

    for i = 1:8
        limiter_theta = sum(slop[:, i] .* limiter[:, i]) / (sum(slop[:, i] .^ 2) + 1.e-7)
        slop[:, i] .*=
            max(0.0, min(min((1.0 + limiter_theta) / 2.0, 2.0), 2.0 * limiter_theta))
    end

    for i = 1:8
        femL[i] =
            sum(A1n[i, 1:3] .* (ER .- EL)) +
            sum(A1n[i, 4:6] .* (BR .- BL)) +
            A1n[i, 7] * (ϕR - ϕL) +
            A1n[i, 8] * (ψR - ψL) +
            0.5 * sum(
                fortsign.(1.0, D) .* (1.0 .- dt ./ (0.5 * (dxL + dxR)) .* abs.(D)) .*
                slop[i, :],
            )
        femR[i] =
            sum(A1p[i, 1:3] .* (ER .- EL)) +
            sum(A1p[i, 4:6] .* (BR .- BL)) +
            A1p[i, 7] * (ϕR - ϕL) +
            A1p[i, 8] * (ψR - ψL) -
            0.5 * sum(
                fortsign.(1.0, D) .* (1.0 .- dt ./ (0.5 * (dxL + dxR)) .* abs.(D)) .*
                slop[i, :],
            )
    end

    # high order correction
    for i = 1:8
        femL[i] +=
            0.5 *
            sum(fortsign.(1.0, D) .* (1.0 .- dt / (dxL + dxR) * abs.(D)) .* slop[i, :])
        femR[i] -=
            0.5 *
            sum(fortsign.(1.0, D) .* (1.0 .- dt / (dxL + dxR) * abs.(D)) .* slop[i, :])
    end

    # transverse correction
    for i = 1:8
        femRU[i] = sum(A2p[i, :] .* femR)
        femRD[i] = sum(A2n[i, :] .* femR)
        femLU[i] = sum(A2p[i, :] .* femL)
        femLD[i] = sum(A2n[i, :] .* femL)
    end
    femRU .*= -0.5 * dt / (dxL + dxR)
    femRD .*= -0.5 * dt / (dxL + dxR)
    femLU .*= -0.5 * dt / (dxL + dxR)
    femLD .*= -0.5 * dt / (dxL + dxR)

    return nothing

end

"""
    flux_emx!(femL, femR, femLU, femLD, femRU, femRD,
    ELL, BLL, EL, BL, ER, BR, ERR, BRR, ϕL, ϕR, ψL, ψR,
    dxL, dxR, A1p, A1n, A2p, A2n, D, sol, χ, ν, dt)

Wave propagation method for 2D Maxwell's equations

* @args: {E, B, ϕ, ψ} in left-left, left, right, and right-right cells
* @args: eigenmatrix (A -> A+ & A-), eigenvalue (D)
* @args: full size of left & right cells
* @args: speed of light (sol)
* @args: auxiliary parameters (χₑ, νᵦ)

"""
function flux_emy!(
    femL::AbstractArray{<:Real,1},
    femR::AbstractArray{<:Real,1},
    femRU::AbstractArray{<:Real,1},
    femRD::AbstractArray{<:Real,1},
    femLU::AbstractArray{<:Real,1},
    femLD::AbstractArray{<:Real,1},
    ELL::AbstractArray{<:Real,1},
    BLL::AbstractArray{<:Real,1},
    EL::AbstractArray{<:Real,1},
    BL::AbstractArray{<:Real,1},
    ER::AbstractArray{<:Real,1},
    BR::AbstractArray{<:Real,1},
    ERR::AbstractArray{<:Real,1},
    BRR::AbstractArray{<:Real,1},
    ϕL::Real,
    ϕR::Real,
    ψL::Real,
    ψR::Real,
    dxL::Real,
    dxR::Real,
    A1p::AbstractArray{<:Real,2},
    A1n::AbstractArray{<:Real,2},
    A2p::AbstractArray{<:Real,2},
    A2n::AbstractArray{<:Real,2},
    D::AbstractArray{<:Real,1},
    sol::Real,
    χ::Real,
    ν::Real,
    dt::Real,
)

    slop = zeros(8, 8)
    slop[3, 1] = 0.5 * sol^2 * (BR[1] - BL[1]) + 0.5 * sol * (ER[3] - EL[3])
    slop[4, 1] = 0.5 * sol * (BR[1] - BL[1]) + 0.5 * (ER[3] - EL[3])
    slop[1, 2] = -0.5 * sol^2 * (BR[3] - BL[3]) + 0.5 * sol * (ER[1] - EL[1])
    slop[6, 2] = 0.5 * sol * (BR[3] - BL[3]) - 0.5 * (ER[1] - EL[1])
    slop[2, 3] = 0.5 * sol * χ * (ER[2] - EL[2])
    slop[7, 3] = 0.5 * χ * (ER[2] - EL[2])
    slop[5, 4] = 0.0
    slop[8, 4] = 0.0
    slop[3, 5] = 0.5 * sol^2 * (BR[1] - BL[1]) - 0.5 * sol * (ER[3] - EL[3])
    slop[4, 5] = -0.5 * sol * (BR[1] - BL[1]) + 0.5 * (ER[3] - EL[3])
    slop[1, 6] = -0.5 * sol^2 * (BR[3] - BL[3]) - 0.5 * sol * (ER[1] - EL[1])
    slop[6, 6] = -0.5 * sol * (BR[3] - BL[3]) - 0.5 * (ER[1] - EL[1])
    slop[2, 7] = -0.5 * sol * χ * (ER[2] - EL[2])
    slop[7, 7] = 0.5 * χ * (ER[2] - EL[2])
    slop[5, 8] = 0.0
    slop[8, 8] = 0.0

    limiter = zeros(8, 8)
    limiter[3, 1] = 0.5 * sol^2 * (BL[1] - BLL[1]) + 0.5 * sol * (EL[3] - ELL[3])
    limiter[4, 1] = 0.5 * sol * (BL[1] - BLL[1]) + 0.5 * (EL[3] - ELL[3])
    limiter[1, 2] = -0.5 * sol^2 * (BL[3] - BLL[3]) + 0.5 * sol * (EL[1] - ELL[1])
    limiter[6, 2] = 0.5 * sol * (BL[3] - BLL[3]) - 0.5 * (EL[1] - ELL[1])
    limiter[2, 3] = 0.5 * sol * χ * (EL[2] - ELL[2])
    limiter[7, 3] = 0.5 * χ * (EL[2] - ELL[2])
    limiter[5, 4] = 0.0
    limiter[8, 4] = 0.0
    limiter[3, 5] = 0.5 * sol^2 * (BRR[1] - BR[1]) - 0.5 * sol * (ERR[3] - ER[3])
    limiter[4, 5] = -0.5 * sol * (BRR[1] - BR[1]) + 0.5 * (ERR[3] - ER[3])
    limiter[1, 6] = -0.5 * sol^2 * (BRR[3] - BR[3]) - 0.5 * sol * (ERR[1] - ER[1])
    limiter[6, 6] = -0.5 * sol * (BRR[3] - BR[3]) - 0.5 * (ERR[1] - ER[1])
    limiter[2, 7] = -0.5 * sol * χ * (ERR[2] - ER[2])
    limiter[7, 7] = 0.5 * χ * (ERR[2] - ER[2])
    limiter[5, 8] = 0.0
    limiter[8, 8] = 0.0

    for i = 1:8
        limiter_theta = sum(slop[:, i] .* limiter[:, i]) / (sum(slop[:, i] .^ 2) + 1.e-7)
        slop[:, i] .*=
            max(0.0, min(min((1.0 + limiter_theta) / 2.0, 2.0), 2.0 * limiter_theta))
    end

    for i = 1:8
        femL[i] =
            sum(A2n[i, 1:3] .* (ER .- EL)) +
            sum(A2n[i, 4:6] .* (BR .- BL)) +
            A2n[i, 7] * (ϕR - ϕL) +
            A2n[i, 8] * (ψR - ψL) +
            0.5 * sum(
                fortsign.(1.0, D) .* (1.0 .- dt ./ (0.5 * (dxL + dxR)) .* abs.(D)) .*
                slop[i, :],
            )
        femR[i] =
            sum(A2p[i, 1:3] .* (ER .- EL)) +
            sum(A2p[i, 4:6] .* (BR .- BL)) +
            A2p[i, 7] * (ϕR - ϕL) +
            A2p[i, 8] * (ψR - ψL) -
            0.5 * sum(
                fortsign.(1.0, D) .* (1.0 .- dt ./ (0.5 * (dxL + dxR)) .* abs.(D)) .*
                slop[i, :],
            )
    end

    # high order correction
    for i = 1:8
        femL[i] +=
            0.5 *
            sum(fortsign.(1.0, D) .* (1.0 .- dt / (dxL + dxR) .* abs.(D)) .* slop[i, :])
        femR[i] -=
            0.5 *
            sum(fortsign.(1.0, D) .* (1.0 .- dt / (dxL + dxR) .* abs.(D)) .* slop[i, :])
    end

    # transverse correction
    for i = 1:8
        femRU[i] = sum(A1p[i, :] .* femR)
        femRD[i] = sum(A1n[i, :] .* femR)
        femLU[i] = sum(A1p[i, :] .* femL)
        femLD[i] = sum(A1n[i, :] .* femL)
    end
    femRU .*= -0.5 * dt / (dxL + dxR)
    femRD .*= -0.5 * dt / (dxL + dxR)
    femLU .*= -0.5 * dt / (dxL + dxR)
    femLD .*= -0.5 * dt / (dxL + dxR)

    return nothing

end


text:
#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2022-Jan-27: define the function to grow a netcdf file
#     2022-Jan-28: add documentation
#
#######################################################################################################################################################################################################
"""
This function grows variables within a netcdf file, and the supported methods are

$(METHODLIST)

"""
function grow_nc! end

#######################################################################################################################################################################################################
#
# Changes to the method
# General
#     2022-Jan-28: add the function to grow semi-automatically with respect to data dimensions
#     2022-Jan-28: add documentation
#
#######################################################################################################################################################################################################
"""

    grow_nc!(ds::Dataset, var_name::String, in_data::Union{AbstractFloat,Array,Int,String}, pending::Bool)

Grow the netcdf dataset, given
- `ds` A `NCDatasets.Dataset` type dataset
- `var_name` New variable name to write to
- `in_data` New data to grow, can be integer, float, and string with N dimens
- `pending` If true, the new data is appened to the end (growth); if false, the data will replace the ones from the bottom (when dimension has already growed)

Note that if there are more variables to grow at the same time, set `pending` to `true` only for the first time you call this function, and set `pending` to `false` for the rest variables.

---
# Examples
```julia
create_nc!("test.nc", String["lon", "lat", "ind"], [36, 18, 0]);
dset = Dataset("test.nc", "a");
append_nc!(dset, "lat", collect(1:18), Dict("longname" => "latitude"), ["lat"]);
append_nc!(dset, "lon", collect(1:36), Dict("longname" => "longitude"), ["lon"]; compress=4);
append_nc!(dset, "ind", collect(1:5), Dict("longname" => "index"), ["ind"]);
append_nc!(dset, "d2d", rand(36,18), Dict("longname" => "a 2d dataset"), ["lon", "lat"]);
append_nc!(dset, "d3d", rand(36,18,5), Dict("longname" => "a 3d dataset"), ["lon", "lat", "ind"]);
grow_nc!(dset, "ind", 6, true);
grow_nc!(dset, "d3d", rand(36,18), false);
grow_nc!(dset, "d3d", rand(36,18), true);
grow_nc!(dset, "ind", 7, false);
close(dset);
```
"""
grow_nc!(ds::Dataset, var_name::String, in_data::Union{AbstractFloat,Array,Int,String}, pending::Bool) = (
    # make sure the data to grow has -1 or the same dimensions as the target, e.g., a 3D dataset can grow with 2D or 3D input
    _dim_ds = length(size(ds[var_name]));
    _dim_in = length(size(in_data));
    @assert _dim_in in [_dim_ds, _dim_ds - 1] "Data to grow must have same or -1 dimensions compared to data in the netcdf file!";
    @assert _dim_ds <= 3 "This function only supports 1D to 3D datasets!";

    # calculate how many layers to add
    if _dim_in < _dim_ds
        _n = 1;
    else
        _n = (typeof(in_data) <: Array ? size(in_data)[end] : 1);
    end;

    # if the data need to pend to the end (grow in unlimited dimension)
    if pending
        if _dim_ds == 1
            ds[var_name][end+1:end+_n] = in_data;
        elseif _dim_ds == 2
            ds[var_name][:,end+1:end+_n] = in_data;
        elseif _dim_ds == 3
            ds[var_name][:,:,end+1:end+_n] = in_data;
        end;
    end;

    # if the unlimited dimension has grown already
    if !pending
        if _dim_ds == 1
            ds[var_name][end+1-_n:end] = in_data;
        elseif _dim_ds == 2
            ds[var_name][:,end+1-_n:end] = in_data;
        elseif _dim_ds == 3
            ds[var_name][:,:,end+1-_n:end] = in_data;
        end;
    end;

    return nothing
);

#######################################################################################################################################################################################################
#
# Changes to the method
# General
#     2022-Jan-28: add the function to grow semi-automatically with respect to data dimensions
#     2022-Jan-28: add documentation
#
#######################################################################################################################################################################################################
"""

    grow_nc!(file::String, var_name::String, in_data::Union{AbstractFloat,Array,Int,String}, pending::Bool)

Grow the netcdf dataset, given
- `file` Path of the netcdf dataset
- `var_name` New variable name to write to
- `in_data` New data to grow, can be integer, float, and string with N dimens
- `pending` If true, the new data is appened to the end (growth); if false, the data will replace the ones from the bottom (when dimension has already growed)

Note that if there are more variables to grow at the same time, set `pending` to `true` only for the first time you call this function, and set `pending` to `false` for the rest variables.

---
# Examples
```julia
create_nc!("test.nc", String["lon", "lat", "ind"], [36, 18, 0]);
dset = Dataset("test.nc", "a");
append_nc!(dset, "lat", collect(1:18), Dict("longname" => "latitude"), ["lat"]);
append_nc!(dset, "lon", collect(1:36), Dict("longname" => "longitude"), ["lon"]; compress=4);
append_nc!(dset, "ind", collect(1:5), Dict("longname" => "index"), ["ind"]);
append_nc!(dset, "d2d", rand(36,18), Dict("longname" => "a 2d dataset"), ["lon", "lat"]);
append_nc!(dset, "d3d", rand(36,18,5), Dict("longname" => "a 3d dataset"), ["lon", "lat", "ind"]);
close(dset);
grow_nc!("test.nc", "ind", 6, true);
grow_nc!("test.nc", "d3d", rand(36,18), false);
grow_nc!("test.nc", "d3d", rand(36,18), true);
grow_nc!("test.nc", "ind", 7, false);
```
"""
grow_nc!(file::String, var_name::String, in_data::Union{AbstractFloat,Array,Int,String}, pending::Bool) = (
    _dset = Dataset(file, "a");
    grow_nc!(_dset, var_name, in_data, pending);
    close(_dset);

    return nothing
);


text:
function solve_JuMPmodel!(JuMPmodel::JuMP.Model, save_path::String; kwargs...)

    timed_log = Dict{Symbol, Any}()

    if JuMPmodel.moi_backend.state == MOIU.NO_OPTIMIZER

        if !(:optimizer in keys(kwargs))
            error("No Optimizer has been defined, can't solve the operational problem")
        end

        _, timed_log[:timed_solve_time],
        timed_log[:solve_bytes_alloc],
        timed_log[:sec_in_gc] = @timed JuMP.optimize!(JuMPmodel,
                                                        kwargs[:optimizer])

    else

        _, timed_log[:timed_solve_time],
        timed_log[:solve_bytes_alloc],
        timed_log[:sec_in_gc] = @timed JuMP.optimize!(JuMPmodel)

    end

    model_status = JuMP.primal_status(JuMPmodel)
    if model_status != MOI.FEASIBLE_POINT::MOI.ResultStatusCode
        error("Status is $(model_status)")
    end

    #creating the results to print to memory
    vars_result = get_model_result(JuMPmodel)
    optimizer_log = get_optimizer_log(JuMPmodel)
    obj_value = Dict(:OBJECTIVE_FUNCTION => JuMP.objective_value(JuMPmodel))
    merge!(optimizer_log, timed_log)

    #results to be printed to memory
    results = Dict(:vars => vars_result,
                   :obj_value => obj_value,
                   :optimizer_log => optimizer_log)

    write_model_results(results, save_path)

    return

end

function make_initial_conditions_from_data(sys::PowerSystems.System)
    inital_conditions = Dict{String, Any}()
    for gen in PowerSystems.get_components(PowerSystems.ThermalStandard, sys)
        ini_g = Dict{Symbol, Float64}()
        name = PowerSystems.get_name(gen)
        ini_g[:power_output_t0] =  PowerSystems.get_activepower(gen)
        activepowerlimits = PowerSystems.get_tech(gen) |> PowerSystems.get_activepowerlimits
        ini_g[:min_power] = activepowerlimits.min
        ini_g[:unit_on_t0] = 1.0*(ini_g[:power_output_t0] > 0)
        ini_g[:time_down_t0] =  999.0*(1.0 - (ini_g[:power_output_t0] > 0))
        ini_g[:time_up_t0] =  999.0*(ini_g[:power_output_t0] > 0)
        inital_conditions[name] = ini_g
    end
    return inital_conditions
end

function update_commitment_status!(ed_model, uc_model, t)

    status_solutions = JuMP.value.(uc_model.obj_dict[:ug])[:,t]
    for (ix ,v) in enumerate(status_solutions)
        name = axes(status_solutions)[1][ix]
        JuMP.fix.(ed_model.obj_dict[:ug][name, :], v; force = true)
    end
    return
end

function update_initial_conditions!(ic, uc_m, ed_m, t)

    pg = ed_m[:pg]
    ug = uc_m[:ug]

    for (k, v) in ic
        #@assert ug[k, t] == 1.0*(JuMP.value(pg[k, t]) > eps())
        v[:power_output_t0] = JuMP.value(pg[k, end]) + JuMP.value(ug[k, t])*v[:min_power]
        # If the unit was on, add to the count
        if v[:unit_on_t0] > 0
            if v[:unit_on_t0] == JuMP.value(ug[k, t])
                v[:time_up_t0] += 1.0
            elseif v[:unit_on_t0] != JuMP.value(ug[k, t])
                v[:time_up_t0] = 0.0
                v[:time_down_t0] = 1.0
            end
            v[:unit_on_t0] = JuMP.value(ug[k, t])
            continue
        end

        # If the unit was off, add to the count
        if v[:unit_on_t0] < 1
            if v[:unit_on_t0] == JuMP.value(ug[k, t])
                v[:time_down_t0] += 1.0
            elseif v[:unit_on_t0] != JuMP.value(ug[k, t])
                v[:time_up_t0] = 1.0
                v[:time_down_t0] = 0.0
            end
            v[:unit_on_t0] = JuMP.value(ug[k, t])
            continue
        end
    end

end


text: ~n=1 in[i^2+j^2==n for i=0:n,j=0:n]

max_stars_repo_name: pgrepds/Codegolf

max_issues_repo_name: pgrepds/Codegolf

max_forks_repo_name: pgrepds/Codegolf


text:
# Simple tests for the functions used in the simulation
cd("D:\\GIT\\correlatedtraits\\inst")
@everywhere include("simulation.jl")
using Base.Test

# set some parameters for testing
n = 20               # number of individuals
M = [1.60, 2.70]     # phenotype means
V = [0.40, 0.30]     # total phenotypic variances
b = 4.05             # Beverton-Holt parameter
ρ = [0.9, 0.0]       # additive and environmental covariances
H = [0.8, 0.8]       # trait heritabilities
p = 3                # population number
r = 5                # replicate number
ngens = 20           # number of generations to simulate

srand(42)
# ------------------------------------------------------------------------------
# from utilities.jl ------------------------------------------------------------
# ------------------------------------------------------------------------------

# rep_each functions ###########################################################
A = Vector{Int64}(7)

rep_each!(A, [1, 2, 3, 4], [1, 2, 0, 4]); @test A == [1, 2, 2, 4, 4, 4, 4]
@test rep_each([1.1, 2.2, 3.3], [1, 0, 3])        == [1.1, 3.3, 3.3, 3.3]
@test rep_eachi([1, 2, 3], [1, 0, 3])             == [1, 3, 3, 3]

# vector math functions ########################################################
A = rand(100); B = rand(100)
C = Vector{Float64}(100); C[:] = A

@test vector_add(A, B)  == A .+ B
vector_add!(A,  B); @test A == C .+ B
vector_add!(A, -B); @test A == C
vector_scalar_add!(A,  0.5); @test A == C .+ 0.5
vector_scalar_add!(A, -0.5); @test A == C
@test vector_mult(A, B) == A .* B
vector_scalar_mult!(A, 0.5); @test A == C .* 0.5

μ1  = rand(50) .+ 5.0
μ2  = rand(50) .+ 3.0
ix  = vcat(collect(1:50), collect(1:50))
D   = zeros(Float64, 50)

sum_prod_deviates!(D, μ1, μ2, A, B, ix);
@test D == sum(reshape((A .- vcat(μ1, μ1)) .* (B .- vcat(μ2, μ2)), 50, 2), 2)[:]

# functions that refer to reference vectors to make calculations ###############
A = zeros(Int64, 3); B = zeros(Float64, 3)

running_count!(A, [1, 1, 1, 2, 2, 3]);               @test A == [3, 2, 1]
running_sum!(B, [1.1, 2.2, 1.1, 3.3], [1, 2, 1, 3]); @test B == [2.2, 2.2, 3.3]

map_vals!(B, [3.0, 2.0, 1.0], [3, 2, 1]);  @test B == [1.0, 2.0, 3.0]
@test map_vals([3.0, 2.0, 1.0], [3, 2, 1])         == [1.0, 2.0, 3.0]
map_valsi!(A, [3, 2, 1], [3, 2, 1]);       @test A == [1, 2, 3]
@test map_valsi([3, 2, 1], [3, 2, 1])              == [1, 2, 3]

ix = Vector{Int64}(5); dict = Dict(1=>4, -1=>2, 2=>3, 0=>1)
# x only includes entries in dictionary
map_dict!(ix, [-1, -1, 0, 1, 2], dict); @test ix == [2, 2, 1, 4, 3]
# x includes some entries not in dictionary
map_dict!(ix, [-1, -3, 0, 1, 2], dict); @test ix == [2, 0, 1, 4, 3]

# functions that operate on arrays #############################################
A = cor2cov([1.0 1.0; 1.0 1.0], [1, 2])
pos_def_check!(A); @test A == [1.0 1.996; 1.996 4.0] # A should change
pos_def_check!(A); @test A == [1.0 1.996; 1.996 4.0] # A should not change

A = [3 1; 1 2; 2 4; 2 3; 4 5]
sortrowsi!(A); @test A == [1 2; 2 3; 2 4; 3 1; 4 5]

B = Array{Int64}(4, 2)
group_index!(B, A, collect(1:4), 4, 5); @test B == [1 1; 2 3; 4 4; 5 5]

# Et cetera #################################################################### ?r2
@test clocktime(4242.676) == "01:10:42:676"

@test r1(pi)  ==  3.1
@test r1(-pi) == -3.1
@test r1(0)   ==  0

# test random seed generator
seed_test = Vector{Int64}(21 * 21 * 11 * 11 * 10 * 20)
c = 0
for i in vcat(-0.99, -0.9:0.1:0.9, 0.99)
  for j in vcat(-0.99, -0.9:0.1:0.9, 0.99)
    for k in vcat(0.01, 0.1:0.1:0.9, 0.99)
      for l in vcat(0.01, 0.1:0.1:0.9, 0.99)
        for m in 1:10
          for n in 1:20
            c += 1
            seed_test[c] = plant_seed([i j], [k l], m, n)
          end
        end
      end
    end
  end
end
@test length(seed_test) == length(unique(seed_test))

# ------------------------------------------------------------------------------
# from simulation.jl -----------------------------------------------------------
# ------------------------------------------------------------------------------

# helper functions #############################################################
@test outname(p, r) == "correlated_traits_3_5.csv"

# growth functions #############################################################
@test growth(0,  exp.(M[2]), b) == 0
@test growth(3,  exp.(M[2]), b) ≈ 1.1315385342
@test growth(10, exp.(M[2]), b) ≈ 0.3585477524
@test growth(50, exp.(M[2]), b) ≈ 0.0731190748

# genetic variance functions ###################################################

# test the popmatrix function --------------------------------------------------
pmx = init_inds(10^7, M, V, ρ, H)

# check basic qualities of the data frame
@test nrow(pmx) == 10^7
@test ncol(pmx) == 5
@test unique(pmx[:X]) == [0]

# check that mean phenotypes are correct
@test r1(mean(pmx[:PD])) ≈ r1(M[1])
@test r1(mean(pmx[:Pr])) ≈ r1(M[2])

# check that variances are correct
@test r1(var(pmx[:PD])) ≈ V[1]                            # total
@test r1(var(pmx[:Pr])) ≈ V[2]
@test r1(var(pmx[:D])) ≈ r1(V[1] * H[1])                  # additive
@test r1(var(pmx[:r])) ≈ r1(V[2] * H[2])
@test r1(var(pmx[:PD] - pmx[:D])) ≈ r1(V[1] * (1 - H[1])) # environmental
@test r1(var(pmx[:Pr] - pmx[:r])) ≈ r1(V[2] * (1 - H[2]))

# check that correlations are correct
@test r1(cor(pmx[:D], pmx[:r])) ≈ ρ[1]
@test r1(cor(pmx[:PD] - pmx[:D], pmx[:Pr] - pmx[:r])) ≈ ρ[2]

# test the covariance matrix functions -----------------------------------------
@test    VE(V, 0, [1, 1]) == zeros(Float64, 2, 2)
@test    VE(V, 1, [1, 1]) == zeros(Float64, 2, 2)
@test    VE(V, 0, [0, 0])  ≈ [V[1]   0; 0   V[2]]
@test r1(VE(V, 1, [0, 0])) ≈ [V[1] 0.3; 0.3 V[2]]

@test    VA(V, 0, [0, 0]) == zeros(Float64, 2, 2)
@test    VA(V, 1, [0, 0]) == zeros(Float64, 2, 2)
@test    VA(V, 0, [1, 1])  ≈  [V[1]   0; 0   V[2]]
@test r1(VA(V, 1, [1, 1])) == [V[1] 0.3; 0.3 V[2]]

# accordian_bins ###############################################################
@test length(accordian_bins(collect(-2000:2000))) == 2301
@test accordian_bins(collect(-2000:2000))[1:100]     == collect(-2000:-1901)
@test accordian_bins(collect(-2000:2000))[2202:2301] == collect(1901:2000)
@test accordian_bins(collect(-2000:2000))[1101:1201] == collect(-50:50)
@test accordian_bins(collect(-1:1)) == collect(-1150:1150)

# metrics ######################################################################

# occupied patches, genotype and phenotype metrics -----------------------------
pmx = init_inds(10^6, M, V, ρ, H)
A = metrics(pmx, M, false, "occupied")
@test nrow(A)   == 1
@test A[1, :X]  == 0
@test A[1, :Nx] == 10^6
@test r1(A[1, :Mean_D]) == r1(M[1])
@test r1(A[1, :Mean_r]) == r1(M[2])
@test r1(A[1, :V_D])    == r1(V[1] * H[1])
@test r1(A[1, :V_r])    == r1(V[2] * H[2])
@test r1(A[1, :C_Dr])   == r1(cov(pmx[:D], pmx[:r]))

A = metrics(pmx, M, true, "occupied")
@test nrow(A)   == 1
@test A[1, :X]  == 0
@test A[1, :Nx] == 10^6
@test r1(A[1, :Mean_D]) == r1(M[1])
@test r1(A[1, :Mean_r]) == r1(M[2])
@test r1(A[1, :V_D])    == r1(V[1] * (1 - H[1]))
@test r1(A[1, :V_r])    == r1(V[2] * (1 - H[2]))
@test r1(A[1, :C_Dr])   == r1(cov(pmx[:PD] - pmx[:D], pmx[:Pr] - pmx[:r]))

# accordian patches, genotype and phenotype metrics ----------------------------
A = metrics(pmx, M, false, "accordian")
@test nrow(A)   == 2301
@test A[:X]  == collect(-1150:1150)
@test A[1151, :Nx] == 10^6
@test A[vcat(collect(1:1150), collect(1152:2301)), :Nx] == zeros(Int64, 2300)
@test r1(A[1151, :Mean_D]) == r1(M[1])
@test r1(A[1151, :Mean_r]) == r1(M[2])
@test r1(A[1151, :V_D])    == r1(V[1] * H[1])
@test r1(A[1151, :V_r])    == r1(V[2] * H[2])
@test r1(A[1151, :C_Dr])   == r1(cov(pmx[:D], pmx[:r]))

A = metrics(pmx, M, true, "accordian")
@test nrow(A)   == 2301
@test A[:X]  == collect(-1150:1150)
@test A[1151, :Nx] == 10^6
@test A[vcat(collect(1:1150), collect(1152:2301)), :Nx] == zeros(Int64, 2300)
@test r1(A[1151, :Mean_D]) == r1(M[1])
@test r1(A[1151, :Mean_r]) == r1(M[2])
@test r1(A[1151, :V_D])    == r1(V[1] * (1 - H[1]))
@test r1(A[1151, :V_r])    == r1(V[2] * (1 - H[2]))
@test r1(A[1151, :C_Dr])   == r1(cov(pmx[:PD] - pmx[:D], pmx[:Pr] - pmx[:r]))

# edge patches, genotype and phenotype metrics ---------------------------------
A = metrics(pmx, M, false, "edges")
@test nrow(A)   == 3
@test A[:X]  == collect(-1:1)
@test A[2, :Nx] == 10^6
@test A[[1, 3], :Nx] == zeros(Int64, 2)
@test r1(A[2, :Mean_D]) == r1(M[1])
@test r1(A[2, :Mean_r]) == r1(M[2])
@test r1(A[2, :V_D])    == r1(V[1] * H[1])
@test r1(A[2, :V_r])    == r1(V[2] * H[2])
@test r1(A[2, :C_Dr])   == r1(cov(pmx[:D], pmx[:r]))

A = metrics(pmx, M, true, "edges")
@test nrow(A)   == 3
@test A[:X]  == collect(-1:1)
@test A[2, :Nx] == 10^6
@test A[[1, 3], :Nx] == zeros(Int64, 2)
@test r1(A[2, :Mean_D]) == r1(M[1])
@test r1(A[2, :Mean_r]) == r1(M[2])
@test r1(A[2, :V_D])    == r1(V[1] * (1 - H[1]))
@test r1(A[2, :V_r])    == r1(V[2] * (1 - H[2]))
@test r1(A[2, :C_Dr])   == r1(cov(pmx[:PD] - pmx[:D], pmx[:Pr] - pmx[:r]))

### LIFE HISTORY FUNCTIONS #####################################################

# mate pairing -----------------------------------------------------------------
pai = [1 1; 2 3; 4 4; 5 5]
sx  = [1 2; 2 3; 2 4; 3 1; 4 5]
m = Vector{Int64}(5)
find_mates!(m, pai, sx, 4); @test m == [-99, -99, 4, 3, -99]

pmx = init_inds(5, M, V, ρ, H)
pmx[:X] = [3, 1, 2, 2, 4]
@test mate(pmx) == [-99, -99, 4, 3, -99]

# reproduction -----------------------------------------------------------------
Nx = ones(Int64, 5)
r  = zeros(Float64, 5)
@test reproduce(Nx, r, m, b)[[1, 2, 5]] == [0, 0, 0]

Nx = rep_eachi([2],  [10^6])
r  = rep_each([M[2]], [10^6])
m  = ones(Int64, 10^6)
@test r1(mean(reproduce(Nx, r, m, b))) ≈ 1.6
@test r1( var(reproduce(Nx, r, m, b))) ≈ 1.6

Nx = rep_eachi( [3],  [10^6])
r  = rep_each([M[2]], [10^6])
m  = ones(Int64, 10^6)
@test r1(mean(reproduce(Nx, r, m, b))) ≈ 1.1
@test r1( var(reproduce(Nx, r, m, b))) ≈ 1.1

# dispersal --------------------------------------------------------------------
# this requires fairly large sample sizes for the variance to work out correctly
x = zeros(Int64, 10^7)
λ = rep_each([exp.(M[1])], [10^7])
@test r1(mean(    disperse(x, λ)))  == 0
@test r1(mean(abs(disperse(x, λ)))) == r1(exp.(M[1]))
@test r1( var(abs(disperse(x, λ)))) == r1(exp.(M[1]))

@test abs(minimum(disperse(x, λ))) < 75
@test abs(maximum(disperse(x, λ))) < 75

x = floor(Int64, vcat(1:(10^6)/2, 1:(10^6)/2))
μ = rep_each([exp.(M[1])], [10^6])
@test r1(mean(abs(disperse(x, λ)  .- x))) == r1(exp.(M[1]))

# combined reproduction and dispersal ------------------------------------------

# reduce sources of variance for testing
V1 = [0.0001, 0.0001]
H1 = [0.01,   0.01]
ρ1 = [0.0,    0.0]
pmx = init_inds(10^6, M, V1, ρ1, H1)
# ensure 2 individuals/patch so matings are successful
pmx[:X] = vcat(1:(nrow(pmx)/2), 1:(nrow(pmx)/2))
pmx1 = repro_disp(pmx, M, V1, ρ1, H1, b, true)

# test traits
@test r1(mean(pmx1[:PD])) == r1(mean(pmx[:PD]))
@test r1(mean(pmx1[:Pr])) == r1(mean(pmx[:Pr]))
@test r1(mean(pmx1[:D]))  == r1(mean(pmx[:D]))
@test r1(mean(pmx1[:r]))  == r1(mean(pmx[:r]))

# test reproduction
@test r1(nrow(pmx1) / nrow(pmx)) ≈ 1.6

# test dispersal
@test r1(mean(abs(pmx1[:X] - pmx1[:dX]))) ≈ r1(exp.(M[1]))
@test r1( var(abs(pmx1[:X] - pmx1[:dX]))) ≈ r1(exp.(M[1]))
@test abs(maximum(pmx1[:X]) -  maximum(pmx[:X])) < 75
@test abs(minimum(pmx1[:X]) -  minimum(pmx[:X])) < 75

# more detailed test of traits using the usual parameters
pmx = init_inds(40*10000, M, V, ρ, H)
pmx[:X] = rep_eachi(vcat(1:10000), rep_eachi([40], [10000]))
pmx1 = repro_disp(pmx, M, V, ρ, H, b)

# test (co)variances
CA = cov(hcat(convert(Vector{Float64}, pmx1[:D]),
              convert(Vector{Float64}, pmx1[:r])))
ΣA = VA(V, ρ[1], H)
@test r1(CA[1]) == r1(ΣA[1])
@test r1(CA[2]) == r1(ΣA[2])
@test r1(CA[4]) == r1(ΣA[4])

CE = cov(hcat(convert(Vector{Float64}, pmx1[:PD] - pmx1[:D]),
              convert(Vector{Float64}, pmx1[:Pr] - pmx1[:r])))
ΣE = VE(V, ρ[2], H)
@test r1(CE[1]) == r1(ΣE[1])
@test r1(CE[2]) == r1(ΣE[2])
@test r1(CE[4]) == r1(ΣE[4])

@test r1(ΣA  ./ (ΣA .+ ΣE)) == r1(CA  ./ (CA .+ CE))

### SIMULATION FUNCTIONS #######################################################

# generate output --------------------------------------------------------------
i = p = r = 2
pmx = init_inds(20, M, V, ρ, H)
pars = [p, r, M[1], M[2], H[1], H[2], ρ[1], ρ[2], i]

A  = generate_output(pmx, p, r, M, H, ρ, i, "full")
@test size(A)    == (12, 2315)
@test A[1, 2:11] == vcat(pars, nrow(pmx))
@test A[2, 1165] == nrow(pmx)
@test A[2, 13]   == 20

A  = generate_output(pmx, p, r, M, H, ρ, i, "zero")
@test size(A)       == (12, 2315)
@test A[1, 2:11]    == vcat(pars, 0)
@test A[2, 12:2315] == rep_each([0.0], [2304])
@test A[2, 13]      == 0

A  = generate_output(pmx, p, r, M, H, ρ, i, "timeout")
@test size(A)       == (12, 2315)
@test A[1, 2:11]    == vcat(pars, -i)
@test A[2, 15:2315] == rep_each([0.0], [2301])
@test A[2, 13]      == -999

#runsim ------------------------------------------------------------------------
p = r = 2
pmx = init_inds(20, M, V, ρ, H)
pars = [p, r, M[1], M[2], H[1], H[2], ρ[1], ρ[2]]

A = runsim(pmx, 5, M, V, ρ, H, b, p, r, time())
@test size(A)          == (60, 2315)
@test A[1:12:60, 2:9]  == repmat(pars', 5, 1)
@test A[1:12:60, 10]   == vcat(1.0:5.0)

A = runsim(pmx, 5, M, V, ρ, H, b, p, r, time() - 3600*40)
@test size(A)          == (60, 2315)
@test A[1:12:60, 2:9]  == repmat(pars', 5, 1)
@test A[1:12:60, 10]   == vcat(1.0:5.0)
@test A[1:12:60, 11]   == -vcat(1.0:5.0)
@test A[2, 15:2315]    == rep_each([0.0], [2301])
@test A[2, 12:14]      == [-999, -999, -999]

pmx = init_inds(0, M, V, ρ, H)
A = runsim(pmx, 5, M, V, ρ, H, b, p, r, time())
@test size(A)          == (60, 2315)
@test A[1:12:60, 2:9]  == repmat(pars', 5, 1)
@test A[1:12:60, 10]   == vcat(1.0:5.0)
@test A[1:12:60, 11]   == rep_each([0.0], [5])
@test A[2, 15:2315]    == rep_each([0.0], [2301])
@test A[2, 12:14]      == [0, 0, 0]

max_stars_repo_name: bochocki/correlatedtraits

max_issues_repo_name: bochocki/correlatedtraits

max_forks_repo_name: bochocki/correlatedtraits


text:
using CSV

"""
    save_solution(gm::GlobalModel; dir::String = SOL_DIR)

Saves the optimal solution of GlobalModel as a CSV.
"""
function save_solution(gm::GlobalModel; name::String = gm.name, dir::String = SOL_DIR)
    open(dir * name * ".txt", "w") do file
        write(file, solution(gm))
    end
end


text:
using Pkg
Pkg.activate(".")

using PowerDynamics
using OrdinaryDiffEq
using DifferentialEquations: DynamicSS

## set per unit MV
const base_power = 1E6 # 1MW
const base_voltage = 20E3 # 20kV
const base_current = base_power / base_voltage # 50A
const base_admittance = base_power / base_voltage^2 # 0.0025Ω^-1
const ω = 2 * π * 50.0 # 314.1593rad/s
# per unit HV
const base_voltage_HV = 110E3 # 110kV

## some additional functions that might develop into new PD.jl features
include("PDpatches.jl")
# control layer (needs ControlledPowerGrid from PDpatches)
include("components/control.jl")
# custom node type
include("components/DGUnit.jl")
include("components/OLTC.jl")
# TODO PD systemsize issue. base State type on GraphData??
include("components/RLLine.jl")
# load actual data
include("../example_cases/CIGRE/cigre_static.jl")

## define test cases

function testbed_chain(
    device,
    S,
    U = complex(1.0, 0.0),
    Y = complex(270.0, -380.0),
    P_ref = t -> -2,
    Q_ref = t -> 0.0,
)
    busses = [SlackAlgebraic(; U = U), device, PQAlgebraic(; S = S)]
    lines = [
        StaticLine(; from = 1, to = 2, Y = Y),
        StaticLine(; from = 2, to = 3, Y = Y),
    ]
    pg = PowerGrid(busses, lines)
    pg, ADN(pg, typeof(device), P_ref, Q_ref)
end

function testbed_CIGRE(
    DG_unit_type;
    idx=1:3,
    P_ref = t -> -2,
    Q_ref = t -> 0.0,
)
    busses_static, lines, T, elist, Zs, Yshs = CIGRE_static()
    busses = copy(busses_static)
    DG_locs = 2:12
    for i in DG_locs
        S_bkgrnd = zero(im)
        try
            S_bkgrnd = busses_static[i].S
        catch
            S_bkgrnd = complex(busses_static[i].P, busses_static[i].Q)
        end
        busses[i] = eval(Symbol(DG_unit_type))(;
            K_pll = 0.1base_voltage, #Hz/pu
            K_PT1 = 1.0, # unit: [y]/[u] = [P]/[P] = 1
            T_PT1 = 0.1, # unit: s
            S_pq = S_bkgrnd, # pu
            Y_n = 0.0,
        )
    end
    pg = PowerGrid(busses[idx], lines[1:length(idx)-1])
    pg, ADN(pg, DG_unit_type, P_ref, Q_ref)
end

function testbed(
    device1,
    device2,
    U = complex(1.0, 0.0),
    Y = complex(262, -375);
    P_ref = t -> -2,
    Q_ref = t -> 0.0,
)
    T = OLTC(
        from = 1,
        to = 2,
        Uos = 110E3 / base_voltage_HV, # Bemessungsspannung Oberspannungsseite in kV
        Uus = 20E3 / base_voltage,# Bemessungsspannung Unterspannungsseite in kV
        k = 0, # Kennzahl der Schaltgruppe, Bsp: YD5-> 5
        ssp = 10.0, # Stufenschalterposition
        stufe_us = 0.625, # Stufung pro Stufenschalterstellung in %
        Sr = 25E6 / base_power, #Bemessungsscheinleistung in MVA
        uk = 12.0, # Kurzschlussspannung in %
        Pvk = 25E3 / base_power, # Kupferverluste in MW
        Pvl = 0.0, # Eisenverluste in kW
        iLeer = 0.0, # Leerlaufstrom in %
    )
    L = PiModelLine(
        from = 2,
        to = 3,
        y = Y,
        y_shunt_km = complex(0., 0.01),
        y_shunt_mk = complex(0., 0.01),
    )
    busses = [SlackAlgebraic(; U = U), device1, device2]
    lines = [T, L]
    pg = PowerGrid(busses, lines)
    pg, ADN(pg, typeof(device), P_ref, Q_ref)
end

function testbed_dynamicline(
    device,
    U = complex(1.0, 0.0),
    R = 0.5,
    L = 0.7;
    P_ref = t -> -2,
    Q_ref = t -> 0.0,
)
    busses = [SlackAlgebraic(; U = U), device, device]
    lines = [
        RLLine(; from = 1, to = 2, R = R, L = L, ω0 = ω),
        RLLine(; from = 1, to = 2, R = R, L = L, ω0 = ω),
    ]
    pg = PowerGrid(busses, lines)
    pg, ADN(pg, typeof(device), P_ref, Q_ref)
end

function testbed_line(
    device,
    U = complex(1.0, 0.0),
    Y = complex(270.0, -380.0);
    P_ref = t -> -2,
    Q_ref = t -> 0.0,
)
    busses = [SlackAlgebraic(; U = U), device, device]
    lines = [
        StaticLine(; from = 1, to = 2, Y = Y),
        StaticLine(; from = 2, to = 3, Y = Y),
    ]
    pg = PowerGrid(busses, lines)
    pg, ADN(pg, typeof(device), P_ref, Q_ref)
end

# all values in p.u.
device = DGUnitPLLPQTracking(;
    K_pll = 0.1, #Hz/pu
    K_PT1 = 1.0, # unit: [y]/[u] = [P]/[P] = 1
    T_PT1 = 10.0, # unit: s
    S_pq = complex(-1.0, -0.25), # pu
    Y_n = 0.0,
)

## setup system

# define events
t_step = 10.0 # step in the power flow reference
t_fault = 1.0 # onset of node short circuit
t_duration = 0.1 # duration of node short circuit

P_ref(t) = t > t_step ? -1.0 : 0
Q_ref(t) = t > t_step ? -0.25 : 0

# edge data
# ode = rhs(pg)
# ode.f.graph_data.e

## find operating point

pg, cpg = testbed_CIGRE(DGUnitPLLPQTracking; idx=1:12, P_ref = t -> P_ref(0.0), Q_ref = t -> Q_ref(0.0))
ode = rhs(cpg)
state(vec) = State(pg, vec)
solution(vec) = PowerGridSolution(vec, pg)

op_guess = initial_guess(cpg)

#initial_error = zeros(2)
op_prob = ODEProblem(ode, op_guess, (0.0, t_step))#, initial_error)
_op = solve(
    SteadyStateProblem(op_prob),
    DynamicSS(Rodas5(); abstol = 1e-8, reltol = 1e-6, tspan=Inf)
)
op = _op |> state

## check dynamics

pg, cpg = testbed_CIGRE(DGUnitPLLPQTracking; idx=1:12, P_ref = P_ref, Q_ref = Q_ref)
ode = rhs(cpg)

# TODO: debug perturb
# TODO: why is there no change???
perturb = Perturbation(3, :θ, Dec(0.0))
u0 = perturb(op)

prob = ODEProblem(ode, u0.vec, (0.0, 30.0))
_sol = solve(
    prob,
    Rodas4(),
    d_discontinuities = [t_step, t_fault, t_fault + t_duration],
)
sol = _sol |> solution

## node short circuit

nsc = NodeShortCircuit(3, 10base_admittance, (t_fault, t_fault+t_duration)) # 10Ω

_sol = simulate_nsc(nsc, cpg, _op, (0., 30.))
sol = _sol |> solution

## plots

using Plots
plotly()

inc = (last(sol.dqsol.prob.tspan) - first(sol.dqsol.prob.tspan)) / 100.
t_array = range(first(sol.dqsol.prob.tspan); stop=last(sol.dqsol.prob.tspan), step=inc)

# TODO: bugs in solution object, e.g. (sol, 2, :v), (sol, 2, :φ)
# (sol, 2, :iabs) triggers call to rhs and hangs
# --> compare with indexing sol

plot(sol, :, :v, c=:gray, label=false)
plot!(sol, 2:3, :v)
hline!(op[2:3, :v], label = "op")
vline!([t_step,], label="step", c=:black)
vline!([t_fault, t_fault + t_duration], label="fault", c=:black)

plot(sol, :, :φ, c=:gray, label=false)
plot!(sol, 2:3, :φ)
plot!(sol, 2:3, :θ, c = "black", label = "PLL")
hline!(op[2:3, :φ], label = "op")
vline!([t_step,], label="step", c=:black)
vline!([t_fault, t_fault + t_duration], label="fault", c=:black)
#ylims!(-0.01, 0.01)

plot(sol, 2:3, :P_g)
hline!(op[2:3, :P_g], label = "op")
plot!(sol, 2:3, :Q_g)
hline!(op[2:3, :Q_g], label = "op")
vline!([t_step,], label="step", c=:black)
vline!([t_fault, t_fault + t_duration], label="fault", c=:black)

plot(sol, 2:3, :P_err)
plot!(sol, 2:3, :Q_err)

Δ = [
    cpg.controller(u, nothing, t) |> first
    for (t, u) in zip(sol.dqsol.t, sol.dqsol.u)
]

plot(sol.dqsol.t, first.(Δ))
plot!(sol.dqsol.t, last.(Δ))
vline!([t_step,], label="step", c=:black)
vline!([t_fault, t_fault + t_duration], label="fault", c=:black)

# currents
node = 3
plot(
    t_array,
    sol.(t_array, node, :iabs),
    label="i_$node"
)
vline!([t_step,], label="step", c=:black)
vline!([t_fault, t_fault + t_duration], label="fault", c=:black)
# current set point
I_r = [
    complex(sol(t, node, :P_g), -sol(t, node, :Q_g)) / sol(t, node, :v)
    for t in t_array
]
plot!(t_array, abs.(I_r), label="I_r_$node")

P_err = sol.(t_array, 2:3, :P_err)
Q_err = sol.(t_array, 2:3, :Q_err)

plot(first.(P_err), first.(Q_err))
plot!(last.(P_err), last.(Q_err))

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
module UpdateWindow

using Gtk, Gtk.ShortNames, Pkg
using ..Ahorn
using Ahorn.PackageHelper

# Modified from stream.jl
function take_stdout(f::Function)
    local oldout = Base.stdout
    local (rd, wr) = redirect_stdout()
    f()
    redirect_stdout(oldout)
    close(wr)
    local res = readavailable(rd)
    close(rd)
    return String(res)
end

function ask_gc_and_exit()
    if ask_dialog("""Ahorn updated!
            But the old version of Ahorn still lingers around.
            We can run gc for you to clean up installed package versions,
            but that is going to affect all your julia environments.
            If you do not have any other Julia packages installed,
            this should be safe to do and save you some disk space.""",
        "Skip", "Run gc", Ahorn.window)
        Pkg.gc()
    end
    info_dialog("""Ahorn will now close.
                Start it again to use the updated version.""",
                Ahorn.window)
    exit()
end

function do_ask_dialog(text::String; no::String="Cancel", yes::String="Update")
    if ask_dialog(text, no, yes, Ahorn.window)
        info_dialog("""Ahorn will now be updated.
            This will freeze the window during the process.
            Do not close Ahorn, you will get notified when the update completes.""", Ahorn.window)

        Pkg.update(["Maple", "Ahorn"])

        ask_gc_and_exit()
    end
end

function updateAhorn(widget::Union{Ahorn.MenuItemsTypes, Nothing}=nothing)
    ask_dialog("Would you like to check for updates?",
            "Cancel", "Check for Updates", Ahorn.window) || return
    try
        localhash = pkghash()
        remotehash = nothing
        if pkghash !== nothing
            remotehash = latesthash()
        else
            localhash = pkghash_tree()
            if pkghash !== nothing
                remotehash = latesthash_tree()
            end
        end

        if localhash === nothing || remotehash === nothing
            updateAhorn_fallback()
            return
        end

        if localhash != remotehash
            do_ask_dialog("""A new version is available!
                    Do you wish to update Ahorn?
                    This will close the program afterwards and you will have to rerun it.""")
        else
            h = pkghash_short()
            do_ask_dialog("""Ahorn seems to be up-to-date$(h !== nothing ? " at hash $h" : "").
                    Do you wish to try updating anyway?
                    This will close the program afterwards and you will have to rerun it.""",
                    yes="Try Updating Anyway")
        end
    catch e
        updateAhorn_fallback()
        return
    end
end

function updateAhorn_fallback()
    ask_dialog("""The update check failed for some reason.
            Do you wish to try updating Ahorn anyway?
            This will download files required for the update if there is one.
            The window might also freeze for a bit.
            Make sure not to close Ahorn during the process!""",
            "Cancel", "Update", Ahorn.window) || return
    try
        info_dialog("""Ahorn will now be updated.
            This will freeze the window during the process.
            Do not close Ahorn, you will get notified when the update completes.""",
            Ahorn.window)
        sim = take_stdout() do
            Pkg.update(["Maple", "Ahorn"])
        end

        println(Base.stdout, sim)

        if occursin(r"(Ahorn|Maple)[^\n]+⇒", sim) || occursin(r"~ (Ahorn|Maple)", sim)
            ask_gc_and_exit()
        else
            h = pkghash_short()
            info_dialog("""Ahorn seems to be up-to-date$(h !== nothing ? " at hash $h" : "").
                    No new update was found.""",
                    Ahorn.window)
        end
    catch e
        println(Base.stderr, "Update check failed")
        for (exc, bt) in Base.catch_stack()
            showerror(Base.stderr, exc, bt)
            println(Base.stderr, "")
        end
        info_dialog("""Something went wrong during the update.
                Please check your error.log file in the Ahorn config directory.""",
                Ahorn.window)
    end
end

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
"""
    DMD(X; alg)
    DMD(X, Y; alg)

Approximates a 'LinearKoopman' with the `AbstractKoopmanAlgorithm` 'alg' from the data matrices `X` or `X` and `Y`, respectively.
If only `X` is given, the data is split into `X[:, 1:end-1]` and `X[:, 2:end]`.

# Example

```julia
alg = DMDPINV()
koopman = DMD(X, alg = alg)

koopman = DMD(X[:, 1:end-1], X[:, 2:end], alg = alg)
```
"""
function DMD(X::AbstractArray; alg::AbstractKoopmanAlgorithm = DMDPINV())
    return DMD(X[:, 1:end-1], X[:, 2:end], alg = alg)
end

function DMD(X::AbstractArray, Y::AbstractArray; alg::AbstractKoopmanAlgorithm = DMDPINV())
    @assert size(X)[2] .== size(Y)[2] "Provide consistent dimensions for data"
    @assert size(Y)[1] .<= size(Y)[2] "Provide consistent dimensions for data"

    # Best Frob norm approximator
    A = alg(X, Y)

    return LinearKoopman(A, zero(eltype(A))*I(size(A,1)), Y*X', X*X', true)
end

"""
    gDMD(X, Y; alg)
    gDMD(t, X ; dt, alg, fdm, itp)

Approximates a 'LinearKoopman' with the `AbstractKoopmanAlgorithm` 'alg' from the data matrices `X` and `Y`.
`X` should contain the state trajectory and `Y` the differential state trajectory.

If no measurements of the differential state are available, `gDMD` can be called with measurement time points `t` as the first argument.
It will then create an interpolation using the interpolation method from `DataInterpolations.jl` defined in `itp`. The trajectory will then be resampled
to equidistant measurements over time corresponding to the mean of `diff(t)` or `dt`, if given.
The differential state measurements will be computed via 'FiniteDifferences.jl', given a `FiniteDifferenceMethod` in `fdm`.

# Example

```julia
koopman = gDMD(X, Y)

fdm = backward_fdm(5,1)
itp = CubicSpline
koopman = gDMD(t, X, fdm = fdm, itp = itp)
```
"""
function gDMD(X::AbstractArray, Y::AbstractArray; alg::AbstractKoopmanAlgorithm = DMDPINV())
    @assert size(X)[2] .== size(Y)[2] "Provide consistent dimensions for data"
    @assert size(Y)[1] .<= size(Y)[2] "Provide consistent dimensions for data"

    # Best Frob norm approximator
    A = alg(X, Y)

    return LinearKoopman(A, zero(eltype(A))*I(size(A,1)), Y*X', X*X', false)
end

function gDMD(t::AbstractVector, X::AbstractArray ; dt::Real = 0.0, alg::DataDrivenDiffEq.AbstractKoopmanAlgorithm = DMDPINV(), fdm::FiniteDifferences.FiniteDifferenceMethod = backward_fdm(5, 1), itp = CubicSpline)
    @assert size(X, 2) == length(t) "Sample size must match."
    @assert test_comp = begin

        if itp ∈ [LinearInterpolation, QuadraticInterpolation] && !isa(fdm, FiniteDifferences.Backward{typeof(fdm.grid), typeof(fdm.coefs)})
            false
        else
            true
        end
    end "LinearInterpolation and QuadraticInterpolation need to a backward finite difference method."

    Δt = iszero(dt) ? mean(diff(t)) : dt
    t̂ = t[1]:Δt:t[end]

    X̂ = zeros(eltype(X), size(X, 1), length(t̂))
    Y = similar(X̂)
    for (i, xi) in enumerate(eachrow(X))
        itp_ = itp(xi, t)
        dx(t) = FiniteDifferences.fdm(fdm, itp_, t)
        X̂[i, :] .= itp_.(t̂)
        Y[i, :] .= dx.(t̂)
    end
    return gDMD(X̂, Y, alg = alg)
end


text:
# We are calculating the energy to check <Psi|H|Psi>
function calculate_energy(wfc, H_k, H_r, dx)
    # Creating momentum and conjugate wavefunctions
    wfc_k = fft(wfc)
    wfc_c = conj(wfc)

    # Finding the momentum and real-space energy terms
    energy_k = wfc_c.*ifft((H_k) .* wfc_k)
    energy_r = wfc_c.* H_r .* wfc

    # Integrating over all space
    energy_final = 0
    for i = 1:length(energy_k)
        energy_final += real(energy_k[i] + energy_r[i])
    end

    return energy_final*dx
end

max_stars_repo_name: atocil/algorithm-archive

max_issues_repo_name: atocil/algorithm-archive

max_forks_repo_name: atocil/algorithm-archive


text:
"""
Computes the tendency to θ_liq_ice due to qt moving between
the working fluid and precipitation
"""
function θ_liq_ice_helper(ts, qt_tendency::FT) where {FT}
    L = TD.latent_heat_vapor(ts)
    Π = TD.exner(ts)
    return -L * qt_tendency / Π / FT(CPP.cp_d(ts.param_set))
end

"""
Computes the tendencies to qt and θ_liq_ice due to precipitation formation
"""
function precipitation_formation(param_set::APS, precipitation_model, qr, qs, area, ρ0, dt, ts)

    qr_tendency = 0.0

    if area > 0.0

        q = TD.PhasePartition(ts)

        if precipitation_model == "clima_1m"
            qr_tendency = min(
                q.liq / dt,
                (
                    CM1.conv_q_liq_to_q_rai(param_set, q.liq) +
                    CM1.accretion(param_set, liq_type, rain_type, q.liq, qr, ρ0)
                ),
            )
        end
        if precipitation_model == "cutoff"
            qsat = TD.q_vap_saturation(ts)
            qr_tendency = min((q.liq + q.ice) / dt, -CM0.remove_precipitation(param_set, q, qsat))
        end
    end

    qt_tendency = -qr_tendency
    qs_tendency = 0.0
    θ_liq_ice_tendency = θ_liq_ice_helper(ts, qt_tendency)

    return PrecipFormation(θ_liq_ice_tendency, qt_tendency, qr_tendency, qs_tendency)
end


text:
using ClimateMachine.Writers

function ivdc_default_diagnostics_init()
end
function ivdc_default_diagnostics_fini()
end
function ivdc_default_diagnostics_collect()
    mpicomm = Diagnostics.Settings.mpicomm
    dg1 = Diagnostics.Settings.dg.modeldata.diagsconf.myivdcdg
    dg2 = Diagnostics.Settings.dg.modeldata.diagsconf.myivdcdg2
    Q = Diagnostics.Settings.dg.modeldata.diagsconf.Q
end

function setup_single_column_default_diagnostics(::IVDCConfigType, interval, out_pref; writer=NetCDFWriter(),interpol=nothing)

	return DiagnosticsGroup(
				"IVDCdefault",
				ivdc_default_diagnostics_init,
				ivdc_default_diagnostics_fini,
				ivdc_default_diagnostics_collect,
				interval,
				out_pref,
				writer,
				interpol,
				)

end

max_stars_repo_name: christophernhill/clima-impldiff-testing

max_issues_repo_name: christophernhill/clima-impldiff-testing

max_forks_repo_name: christophernhill/clima-impldiff-testing


text:
#!/usr/bin/env julia

using ArgParse
using Distributed

function compute_pi(a::Float64, b::Float64, n::Int64)::Float64
    result::Float64 = 0.0
    delta = (b - a)/n
    x::Float64 = a
    i::Int64 = 1
    while x < 1.0 && i <= n
        result += sqrt(1.0 - x^2)
        x += delta
        i += 1
    end
    result *= 2.0*delta
    println("($a, $b) -> $result, $n")
    return result
end

function main()
    arg_parser = ArgParseSettings()
    arg_parser.description = "Compute pi using a quadrature on multiple processes"
    @add_arg_table arg_parser begin
        "--nr_points", "-p"
            arg_type = Int64
            help = "number of points for the quadrature method"
            required = true
        "--nr_workers", "-w"
            arg_type = Int
            help = "number of workers to use"
            default = 0
    end
    options = parse_args(ARGS, arg_parser)
    nr_points = options["nr_points"]
    nr_workers = options["nr_workers"]
    println("running with $nr_workers workers")
    if nr_workers == 0
        @time result = compute_pi(-1.0, 1.0, nr_points)
        println(result)
    else
        addprocs(nr_workers)
        @everywhere function compute_pi(a::Float64, b::Float64, n::Int64)::Float64
            result::Float64 = 0.0
            delta = (b - a)/n
            x::Float64 = a
            i::Int64 = 1
            while x < 1.0 && i <= n
                result += sqrt(1.0 - x^2)
                x += delta
                i += 1
            end
            result *= 2.0*delta
            println("($a, $b) -> $result, $n")
            return result
        end

        delta = 2.0/nr_workers
        result = @distributed (+) for i in 1:nr_workers
            @time compute_pi(-1.0 + (i - 1)*delta, -1.0 + i*delta, div(nr_points, nr_workers))
        end
        println(result)
    end
end

main()


text:
using Laplace
using Test

@time @testset "postwidder" begin include("postwidder_tests.jl") end


text:
# Mean of MCChain

mean(c::MCChain, par::Range=1:size(c.samples, 2)) = Float64[mean(c.samples[:, i]) for i = par]

mean(c::MCChain, par::Real) = mean(c, par:par)