text:
# # One-dimensional regression
#
# ## Setup
#
# Loading the necessary packages and setting seed.

using AbstractGPs
using Distributions
using StatsFuns

using Plots
default(legend=:outertopright, size=(700,400))

using Random
Random.seed!(1234)
#md nothing #hide

# Load toy regression
# [dataset](https://github.com/GPflow/docs/blob/master/doc/source/notebooks/basics/data/regression_1D.csv)
# taken from GPFlow examples.

x = [0.8658165855998895, 0.6661700880180962, 0.8049218148148531, 0.7714303440386239,
    0.14790478354654835, 0.8666105548197428, 0.007044577166530286, 0.026331737288148638,
    0.17188596617099916, 0.8897812990554013, 0.24323574561119998, 0.028590102134105955]
y = [1.5255314337144372, 3.6434202968230003, 3.010885733911661, 3.774442382979625,
    3.3687639483798324, 1.5506452040608503, 3.790447985799683, 3.8689707574953,
    3.4933565751758713, 1.4284538820635841, 3.8715350915692364, 3.7045949061144983]
scatter(x, y; xlabel="x", ylabel="y", legend=false)

# We split the observations into train and test data.

x_train = x[1:8]
y_train = y[1:8]
x_test = x[9:end]
y_test = y[9:end]
#md nothing #hide

# We instantiate a Gaussian process with a Matern kernel. The kernel has
# fixed variance and length scale parameters of default value 1.

f = GP(Matern52Kernel())
#md nothing #hide

# We create a finite dimentional projection at the inputs of the training dataset
# observed under Gaussian noise with standard deviation $\sigma = 0.1$, and compute the
# log-likelihood of the outputs of the training dataset.

fx = f(x_train, 0.1)
logpdf(fx, y_train)

# We compute the posterior Gaussian process given the training data, and calculate the
# log-likelihood of the test dataset.

p_fx = posterior(fx, y_train)
logpdf(p_fx(x_test), y_test)

# We plot the posterior Gaussian process along with the observations.

scatter(
    x_train, y_train;
    xlim=(0,1), xlabel="x", ylabel="y",
    title="posterior (default parameters)", label="Train Data",
)
scatter!(x_test, y_test; label="Test Data")
plot!(p_fx, 0:0.001:1; label=false)

# ## Markov Chain Monte Carlo
#
# Previously we computed the log likelihood of the untuned kernel parameters of the GP.
# We now also perform approximate inference over said kernel parameters using different
# Markov chain Monte Carlo (MCMC) methods. I.e., we approximate the posterior distribution
# of the kernel parameters with samples from a Markov chain.
#
# We define a function which returns the log-likelihood of the data for different variance
# and inverse lengthscale parameters of the Matern kernel. We ensure that these parameters are
# positive with the softplus function
# ```math
# f(x) = \log (1 + \exp x).
# ```

struct GPLoglikelihood{X,Y}
    x::X
    y::Y
end

function (ℓ::GPLoglikelihood)(params)
    kernel = ScaledKernel(
        transform(
            Matern52Kernel(),
            ScaleTransform(softplus(params[1]))
        ),
        softplus(params[2]),
    )
    f = GP(kernel)
    fx = f(ℓ.x, 0.1)
    return logpdf(fx, ℓ.y)
end

const loglik_train = GPLoglikelihood(x_train, y_train)
#md nothing #hide

# We define a Gaussian prior for the joint distribution of the two transformed kernel
# parameters. We assume that both parameters are independent with mean 0 and variance 1.

logprior(params) = logpdf(MvNormal(2, 1), params)
#md nothing #hide

# ### Hamiltonian Monte Carlo
#
# We start with a Hamiltonian Monte Carlo (HMC) sampler. More precisely, we use the
# [No-U-Turn sampler (NUTS)](http://www.jmlr.org/papers/volume15/hoffman14a/hoffman14a.pdf),
# which is provided by the Julia packages
# [AdvancedHMC.jl](https://github.com/TuringLang/AdvancedHMC.jl/) and
# [DynamicHMC.jl](https://github.com/tpapp/DynamicHMC.jl/).
#
# #### AdvancedHMC
#
# We start with performing inference with AdvancedHMC.

using AdvancedHMC
using ForwardDiff

# Set the number of samples to draw and warmup iterations.

n_samples = 2_000
n_adapts = 1_000
#md nothing #hide

# Define a Hamiltonian system of the log joint probability.

logjoint_train(params) = loglik_train(params) + logprior(params)
metric = DiagEuclideanMetric(2)
hamiltonian = Hamiltonian(metric, logjoint_train, ForwardDiff)
#md nothing #hide

# Define a leapfrog solver, with initial step size chosen heuristically.

initial_params = rand(2)
initial_ϵ = find_good_stepsize(hamiltonian, initial_params)
integrator = Leapfrog(initial_ϵ)
#md nothing #hide

# Define an HMC sampler, with the following components:
# - multinomial sampling scheme,
# - generalised No-U-Turn criteria, and
# - windowed adaption for step-size and diagonal mass matrix

proposal = NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)
adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator))
#md nothing #hide

# We draw samples from the posterior distribution of kernel parameters. These samples
# are in the unconstrained space $\mathbb{R}^2$.

samples, _ = sample(
    hamiltonian,
    proposal,
    initial_params,
    n_samples,
    adaptor,
    n_adapts;
    progress=false,
)
#md nothing #hide

# We transform the samples back to the constrained space and compute the mean of both
# parameters:

samples_constrained = [map(softplus, p) for p in samples]
mean_samples = mean(samples_constrained)

# We plot a histogram of the samples for the two parameters.
# The vertical line in each graph indicates the mean of the samples.

histogram(
    reduce(hcat, samples_constrained)';
    xlabel="sample", ylabel="counts", layout=2,
    title=["inverse length scale" "variance"], legend=false,
)
vline!(mean_samples'; linewidth=2)

# We approximate the log-likelihood of the test data using the posterior Gaussian processes
# for kernels with the sampled kernel parameters. We can observe that there is a significant
# improvement over the log-likelihood of the test data with respect to the posterior
# Gaussian process with default kernel parameters of value 1.

struct GPPosterior{X,Y}
    x::X
    y::Y
end

function (g::GPPosterior)(p)
    kernel = ScaledKernel(
        transform(
            Matern52Kernel(),
            ScaleTransform(softplus(p[1]))
        ),
        softplus(p[2]),
    )
    f = GP(kernel)
    return posterior(f(g.x, 0.1), g.y)
end

const gp_posterior = GPPosterior(x_train, y_train)

mean(logpdf(gp_posterior(p)(x_test), y_test) for p in samples)

# We sample a function from the posterior GP for the final 100 samples of kernel
# parameters.

plt = scatter(
    x_train, y_train;
    xlim=(0,1), xlabel="x", ylabel="y",
    title="posterior (AdvancedHMC)", label="Train Data",
)
scatter!(plt, x_test, y_test; label="Test Data")
for p in samples[(end-100):end]
    sampleplot!(plt, gp_posterior(p)(0:0.02:1), 1)
end
plt

# #### DynamicHMC
#
# We repeat the inference with DynamicHMC.

using DynamicHMC
using LogDensityProblems

# We have to implement the LogDensityProblems interface for `GPLogLikelihood`.

## Log joint density
LogDensityProblems.logdensity(ℓ::GPLoglikelihood, params) = ℓ(params) + logprior(params)

## The parameter space is two-dimensional
LogDensityProblems.dimension(::GPLoglikelihood) = 2

## `GPLoglikelihood` does not allow to evaluate derivatives of
## the log-likelihood function
function LogDensityProblems.capabilities(::Type{<:GPLoglikelihood})
    LogDensityProblems.LogDensityOrder{0}()
end

# Now we can draw samples from the posterior distribution of kernel parameters with
# DynamicHMC. Again we use [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl)
# to compute the derivatives of the log joint density with automatic differentiation.

samples = mcmc_with_warmup(
    Random.GLOBAL_RNG,
    ADgradient(:ForwardDiff, loglik_train),
    n_samples;
    reporter = NoProgressReport(),
).chain
#md nothing #hide

# We transform the samples back to the constrained space and compute the mean of both
# parameters:

samples_constrained = [map(softplus, p) for p in samples]
mean_samples = mean(samples_constrained)

# We plot a histogram of the samples for the two parameters.
# The vertical line in each graph indicates the mean of the samples.

histogram(
    reduce(hcat, samples_constrained)';
    xlabel="sample", ylabel="counts", layout=2,
    title=["inverse length scale" "variance"], legend=false,
)
vline!(mean_samples'; linewidth=2)

# Again we can observe that there is a significant improvement over the log-likelihood
# of the test data with respect to the posterior Gaussian process with default kernel
# parameters.

mean(logpdf(gp_posterior(p)(x_test), y_test) for p in samples)

# We sample a function from the posterior GP for the final 100 samples of kernel
# parameters.

plt = scatter(
    x_train, y_train;
    xlim=(0,1), xlabel="x", ylabel="y",
    title="posterior (DynamicHMC)", label="Train Data",
)
scatter!(plt, x_test, y_test; label="Test Data")
for p in samples[(end-100):end]
    sampleplot!(plt, gp_posterior(p)(0:0.02:1), 1)
end
plt

# ### Elliptical slice sampling
#
# Instead of HMC, we use
# [elliptical slice sampling](http://proceedings.mlr.press/v9/murray10a/murray10a.pdf)
# which is provided by the Julia package
# [EllipticalSliceSampling.jl](https://github.com/TuringLang/EllipticalSliceSampling.jl/).

using EllipticalSliceSampling

# We draw 2000 samples from the posterior distribution of kernel parameters.

samples = sample(
    ESSModel(
        MvNormal(2, 1), # Gaussian prior
        loglik_train,
    ),
    ESS(),
    n_samples;
    progress=false,
)
#md nothing #hide

# We transform the samples back to the constrained space and compute the mean of both
# parameters:

samples_constrained = [map(softplus, p) for p in samples]
mean_samples = mean(samples_constrained)

# We plot a histogram of the samples for the two parameters.
# The vertical line in each graph indicates the mean of the samples.

histogram(
    reduce(hcat, samples_constrained)';
    xlabel="sample", ylabel="counts", layout=2,
    title=["inverse length scale" "variance"],
)
vline!(mean_samples'; layout=2, labels="mean")

# Again we can observe that there is a significant improvement over the log-likelihood
# of the test data with respect to the posterior Gaussian process with default kernel
# parameters.

mean(logpdf(gp_posterior(p)(x_test), y_test) for p in samples)

# We sample a function from the posterior GP for the final 100 samples of kernel
# parameters.

plt = scatter(
    x_train, y_train;
    xlim=(0,1), xlabel="x", ylabel="y",
    title="posterior (EllipticalSliceSampling)", label="Train Data",
)
scatter!(plt, x_test, y_test; label="Test Data")
for p in samples[(end-100):end]
    sampleplot!(plt, gp_posterior(p)(0:0.02:1), 1)
end
plt

# ## Variational Inference
#
# Sanity check for the Evidence Lower BOund (ELBO) implemented according to
# M. K. Titsias's _Variational learning of inducing variables in sparse Gaussian processes_.

elbo(fx, y_train, f(rand(5)))

# We use the LBFGS algorithm to maximize the given ELBO. It is provided by the Julia
# package [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl).

using Optim

# We define a function which returns the negative ELBO for different variance and inverse
# lengthscale parameters of the Matern kernel and different pseudo-points. We ensure that
# the kernel parameters are positive with the softplus function
# ```math
# f(x) = \log (1 + \exp x),
# ```
# and that the pseudo-points are in the unit interval $[0,1]$ with the logistic function
# ```math
# f(x) = \frac{1}{1 + \exp{(-x)}}.
# ```

struct NegativeELBO{X,Y}
    x::X
    y::Y
end

function (g::NegativeELBO)(params)
    kernel = ScaledKernel(
        transform(
            Matern52Kernel(),
            ScaleTransform(softplus(params[1]))
        ),
        softplus(params[2]),
    )
    f = GP(kernel)
    fx = f(g.x, 0.1)
    return -elbo(fx, g.y, f(logistic.(params[3:end])))
end
#md nothing #hide

# We randomly initialize the kernel parameters and 5 pseudo points, and minimize the
# negative ELBO with the LBFGS algorithm and obtain the following optimal parameters:

x0 = rand(7)
opt = optimize(NegativeELBO(x_train, y_train), x0, LBFGS())

#-

opt.minimizer

# The optimized value of the inverse lengthscale is

softplus(opt.minimizer[1])

# and of the variance is

softplus(opt.minimizer[2])

# We compute the log-likelihood of the test data for the resulting approximate
# posterior. We can observe that there is a significant improvement over the
# log-likelihood with the default kernel parameters of value 1.

opt_kernel = ScaledKernel(
    transform(
        Matern52Kernel(),
        ScaleTransform(softplus(opt.minimizer[1]))
    ),
    softplus(opt.minimizer[2]),
)
opt_f = GP(opt_kernel)
opt_fx = opt_f(x_train, 0.1)
ap = approx_posterior(VFE(), opt_fx, y_train, opt_f(logistic.(opt.minimizer[3:end])))
logpdf(ap(x_test), y_test)

# We visualize the approximate posterior with optimized parameters.

scatter(
    x_train, y_train;
    xlim=(0,1), xlabel="x", ylabel="y",
    title="posterior (VI with sparse grid)", label="Train Data",
)
scatter!(x_test, y_test; label="Test Data")
plot!(ap, 0:0.001:1; label=false)
vline!(logistic.(opt.minimizer[3:end]); label="Pseudo-points")


text:
"""
    DualPerplex{T <: Real} = Nilplexification{Prplex{T}}

A dual-perplex number is a nilplexification of the perplex numbers.
"""
const DualPerplex{T <: Real} = Nilplexification{Prplex{T}}

function DualPerplex(z::Prplex{T}) where T <: Real
    DualPerplex{T}(z, zero(Prplex{T}))
end

function DualPerplex(a::T, b::T, c::T, d::T) where T <: Real
    DualPerplex{T}(Prplex{T}(a, b), Prplex{T}(c, d))
end

function DualPerplex(a::T, b::T, c::T) where T <: Real
    DualPerplex{T}(Prplex{T}(a, b), Prplex{T}(c, zero(T)))
end

function DualPerplex(a::T, b::T) where T <: Real
    DualPerplex{T}(Prplex{T}(a, b), zero(Prplex{T}))
end

function DualPerplex(a::T) where T <: Real
    DualPerplex{T}(Prplex{T}(a, zero(T)), zero(Prplex{T}))
end

function DualPerplex(a::Real, b::Real)
    DualPerplex(promote(a, b)...)
end

function DualPerplex(a::Real, b::Real, c::Real)
    DualPerplex(promote(a, b, c)...)
end

function DualPerplex(a::Real, b::Real, c::Real, d::Real)
    DualPerplex(promote(a, b, c, d)...)
end

"""
    unreal(z::DualPerplex)

The unreal part of a dual-perplex number is a 3-vector.
"""
function unreal(z::DualPerplex)
    vcat(unreal(z.l), array(z.r))
end

function show(io::IO, z::DualPerplex)
    print(io, "[1: ")
    print(io, z.l.l)
    print(io, ", R: ")
    print(io, z.l.r)
    print(io, ", A: ")
    print(io, z.r.l)
    print(io, ", RA: ")
    print(io, z.r.r)
    print(io, "]")
end

function random(::Type{DualPerplex{T}}) where T <: Real
    DualPerplex{T}(random(Prplex{T}), random(Prplex{T}))
end

function barR(z::DualPerplex{T}) where T <: Real
    DualPerplex{T}(barR(z.l), barR(z.r))
end

function barA(z::DualPerplex{T}) where T <: Real
    DualPerplex{T}(z.l, -(z.r))
end

function quadrance(z::DualPerplex)
    (z * barA(z)).l
end

function norm(z::DualPerplex)
    abs2(quadrance(z))
end

function iszerodivisor(z::DualPerplex)
    iszerodivisor(quadrance(z))
end

function inv(z::DualPerplex)
    if iszerodivisor(z)
        error(ZeroDivisorInverse)
    end

    q = quadrance(z)
    barA(z) * barR(q) / abs2(q)
end

ext: jl


text:
include("Networks.jl")
include("DIP.jl")

# use :reco=DIP for dispatch


text:
module Servers
using StructTypes, Distributed, JSON3, HTTP, REPL.TerminalMenus, Parameters
using ..DFControl
using ..Utils
using ..Jobs

export Server

const SERVER_DIR = DFControl.config_path("servers")

@enum Scheduler Slurm=1 Bash=2

#TODO: there is some disconnect between domain and ssh configuration to be used...
@with_kw mutable struct Server
    name::String           = "nothing"
    username::String       = "nothing"
    domain::String         = "nothing"
    port::Int              = 8080
    scheduler::Scheduler   = Bash
    mountpoint::String     = ""
    julia_exec::String     = "julia"
    default_jobdir::String = homedir()
    local_port::Int = 0
    max_concurrent_jobs::Int = 100
end

function Server(s::String; name="")
    server = load_server(s) #First check if previous server exists
    if server !== nothing
        return server
    end
    #TODO load existing config
    # Create new server
    if occursin("@", s)
        username, domain = split(s, "@")
    else
        @info "Server with name $s not found."
        name = s
        username, domain = "", ""
    end
    @info "Creating new Server configuration..."
    while isempty(name)
        print("Please specify the Server's identifying name:")
        name = readline()
    end
    if load_server(name) !== nothing
        @warn "A server with $name was already configured and will be overwritten."
    end
    try
        @info "Trying to pull existing configuration from $username@$domain..."

        tserver = load_remote_config(username, domain)

        local_port = 0
        try
            run(`nc -vz $domain 22`)
            local_port = 0
        catch
            port_str = ""
            while isempty(port_str)
                print("Local tunnel port:")
                port_str = readline()
            end
            local_port = parse(Int, port_str)
        end

        tserver.name       = name
        tserver.domain     = domain
        tserver.local_port = local_port
        println("Server configured as:")
        println(tserver)
        save(tserver)
        return tserver
    catch
        while isempty(username)
            print("Username:")
            username = readline()
        end
        while isempty(domain)
            print("Domain:")
            domain = readline()
        end
        print("Port (default: 8080):")
        port_str = readline()
        port = isempty(port_str) ? 8080 : parse(Int, port_str)

        print("Julia Exec (default: julia):")
        julia_str = readline()
        julia = isempty(julia_str) ? "julia" : julia_str

        @info "Trying to connect to $username@$domain..."
        dir = remotecall_fetch(homedir,
                               Distributed.addprocs([("$username" * "@" * "$domain", 1)], exename=julia)[1])
        print("Jobs top dir (default: $dir):")
        dir = abspath(readline())

        scheduler_choice = request("Please select scheduler:",
                                   RadioMenu([string.(instances(Scheduler))...]))
        scheduler_choice == -1 && return
        scheduler = Scheduler(scheduler_choice)
        mounted_choice = request("Has the server been mounted?", RadioMenu(["yes", "no"]))
        mounted_choice == -1 && return
        if mounted_choice == 1
            print("Please specify mounting point:")
            mountpoint = readline()
        else
            mountpoint = ""
        end
        server = Server(name, username, domain, port, scheduler, mountpoint, julia_str, dir, 0)
        println("Server configured as:")
        println(server)
        save(server)
        return server
    end
end

StructTypes.StructType(::Type{Server}) = StructTypes.Struct()
islocal(s::Server) = s.domain == "localhost"

Base.joinpath(s::Server, p...) = joinpath(s.default_jobdir, p...)
Base.ispath(s::Server, p...) = islocal(s) ? ispath(p...) : JSON3.read(HTTP.get(s, "/get_ispath/" * joinpath(p...)).body, Bool)

Utils.searchdir(s::Server, dir, str) = joinpath.(dir, filter(x->occursin(str, x), readdir(s, dir)))

parse_server_config(config) = JSON3.read(config, Server)
read_server_config(config_file) = parse_server_config(read(config_file, String))

function load_remote_config(username, domain)
    cmd = "cat ~/.julia/config/DFControl/servers/localhost.json"
    return parse_server_config(server_command(username, domain, cmd).stdout)
end

function known_servers(fuzzy = "")
    if ispath(SERVER_DIR)
        servers = [JSON3.read(read(joinpath(SERVER_DIR, s), String), Server)
                   for s in filter(x -> occursin(fuzzy, x), readdir(SERVER_DIR))]
    else
        servers = Server[]
    end
    return servers
end

function load_server(name::String)
    if occursin("@", name)
        return getfirst(x -> x.name == name, known_servers())
    else
        all = known_servers(name)
        if length(all) > 0
            return all[1]
        else
            return nothing
        end
    end
end

function save(s::Server)
    mkpath(SERVER_DIR)
    if ispath(joinpath(SERVER_DIR, s.name * ".json"))
        @info "Updating previously existing configuration for server $s."
    end
    JSON3.write(joinpath(SERVER_DIR, s.name * ".json"),  s)
end

ssh_string(s::Server) = s.username * "@" * s.domain
http_string(s::Server) = s.local_port != 0 ? "http://localhost:$(s.local_port)" : "http://$(s.domain):$(s.port)"

function HTTP.request(method::String, s::Server, url, args...; kwargs...)
    return HTTP.request(method, string(http_string(s), url), args...; kwargs...)
end

for f in (:get, :put, :post, :head)
    str = uppercase(string(f))
    @eval function HTTP.$(f)(s::Server, url, args...; kwargs...)
        return HTTP.request("$($str)", s, url, args...; kwargs...)
    end
end

function Distributed.addprocs(server::Server, nprocs::Int = 1, args...; kwargs...)
    if islocal(server)
        proc = Distributed.addprocs(nprocs, args...; kwargs...)
    else
        proc = Distributed.addprocs([(ssh_string(server), nprocs)], args...; exename=`$(server.julia_exec)`, dir=server.default_jobdir, tunnel=true, kwargs...)
    end
    return proc
end

function remove_server!(name::String)
    return ispath(joinpath(SERVER_DIR, name * ".json")) &&
           rm(joinpath(SERVER_DIR, name * ".json"))
end

function isalive(s::Server)
    try
        p1 = Pipe()
        p2 = Pipe()
        if s.local_port != 0
            run(pipeline(`nc -vz localhost $(s.local_port)`, stdout = p1, stderr=p2))
        else
            run(pipeline(`nc -vz  $(s.domain) $(s.port)`, stdout = p2, stderr=p1))
        end
        close(p1.in)
        close(p2.in)
        resp = HTTP.get(s, "/server_config")
        return JSON3.read(resp.body, Server).username == s.username
    catch
        return false
    end
end

function start(s::Server)
    @info "Starting:\n$s"
    julia_cmd = """$(s.julia_exec) --startup-file=no -t auto -e "using DFControl; DFControl.Resource.run()" &> ~/.julia/config/DFControl/errors.log"""
    if s.domain != "localhost"
        if s.local_port != 0
            t = getfirst(x->occursin("ssh -f -N", x), split(read(pipeline(`ps aux` , stdout = `grep $(s.local_port)`), String), "\n"))

            if t !== nothing
                run(`kill $(split(t)[2])`)
            end
            run(Cmd(`ssh -f -L $(s.local_port):localhost:$(s.port) $(ssh_string(s)) $julia_cmd`, detach=true))
        else
            run(Cmd(`ssh -f $(ssh_string(s)) $julia_cmd`, detach=true))

        end
    else
        scrpt = "using DFControl; DFControl.Resource.run()"
        e = s.julia_exec
        julia_cmd = `$(e) --startup-file=no -t auto -e $(scrpt) '&''>' '~'/.julia/config/DFControl/errors.log '&'`
        run(Cmd(julia_cmd, detach=true), wait=false)
    end

    #TODO: little hack here
    retries = 0
    while !isalive(s) && retries < 60
        tserver = islocal(s) ? read_server_config(DFC.config_path("servers/localhost.json")) : load_remote_config(s.username, s.domain)
        s.port = tserver.port
        sleep(1)
        retries += 1
    end
    if retries == 60
        error("Something went wrong starting the server.")
    else
        if s.local_port == 0
            @info "Daemon on Server $(s.name) started, listening on port $(s.port)."
        else
            @info "Daemon on Server $(s.name) started, listening on local port $(s.local_port)."
        end
    end
    return
end

function maybe_start_server(s::Server)
    if !isalive(s)
        tserver = islocal(s) ? read_server_config(DFC.config_path("servers/localhost.json")) : load_remote_config(s.username, s.domain)
        s.port = tserver.port
        if !isalive(s)
            start(s)
        else
            save(s)
        end
    end
    return s
end
maybe_start_server(j) = (s = Server(j); return maybe_start_server(s))

kill_server(s::Server) = HTTP.put(s, "/kill_server")

function restart_server(s::Server)
    kill_server(s)
    return start(s)
end

function maybe_create_localhost()
    t = load_server("localhost")
    if t === nothing
        scheduler = Sys.which("sbatch") === nothing ? Bash : Slurm
        if !haskey(ENV, "DFCONTROL_PORT")
            port = 8080
        else
            port = parse(Int, ENV["DFCONTROL_PORT"])
        end
        if !haskey(ENV, "DFCONTROL_JOBDIR")
            dir = homedir()
        else
            dir = ENV["DFCONTROL_PORT"]
        end
        julia_exec = joinpath(Sys.BINDIR, "julia")
        out = Server("localhost", ENV["USER"], "localhost", port, scheduler, "", julia_exec,
                     dir, 0)
        save(out)
        return out
    else
        return t
    end
end

"""
    pull(server::Server, server_file::String, local_file::String)

Pulls `server_file` from the server the `local_file`.
"""
function pull(server::Server, server_file::String, filename::String)
    if islocal(server)
        cp(server_file, filename; force = true)
    else
        run(`scp $(ssh_string(server) * ":" * server_file) $filename`)
    end
end

function pull(j::Job, f, t)
    @assert ispath(j, f) "File $f not found in jobdir."
    pull(Server(j.server), joinpath(j, f), t)
end

"""
    push(local_file::String, server::Server, server_file::String)

Pushes the `local_file` to the `server_file` on the server.
"""
function push(filename::String, server::Server, server_file::String)
    if islocal(server)
        cp(filename, server_file; force = true)
    else
        run(`scp $filename $(ssh_string(server) * ":" * server_file)`)
    end
end

function server_command(username, domain, cmd)
    out = Pipe()
    err = Pipe()
    if domain == "localhost"
        process = run(pipeline(ignorestatus(cmd), stdout=out, stderr=err))
    else
        process = run(pipeline(ignorestatus(`ssh $(username * "@" * domain) source /etc/profile '&''&' $cmd`), stdout=out, stderr=err))
    end
    close(out.in)
    close(err.in)

    stdout = String(read(out))
    stderr = String(read(err))
    return (
      stdout = stdout,
      stderr = stderr,
      code = process.exitcode
    )
end

server_command(s::Server, cmd) = server_command(s.username, s.domain, cmd)

function has_modules(s::Server)
    try
        server_command(s, `module avail`).code == 0
    catch
        false
    end
end

function available_modules(s::Server)
    if has_modules(s)
        return server_command(s, `module avail`)
    else
        return String[]
    end
end

function Base.readdir(s::Server, dir::String)
    maybe_start_server(s)
    resp = HTTP.get(s, "/readdir/" * dir)
    return JSON3.read(resp.body, Vector{String})
end

end


text:
using Oceananigans.Operators: Δx, Δy, Δz

#####
##### Halo filling for value and gradient boundary conditions
#####

@inline linearly_extrapolate(c₀, ∇c, Δ) = c₀ + ∇c * Δ

@inline  left_gradient(bc::GBC, c¹, Δ, i, j, args...) = getbc(bc, i, j, args...)
@inline right_gradient(bc::GBC, cᴺ, Δ, i, j, args...) = getbc(bc, i, j, args...)

@inline  left_gradient(bc::VBC, c¹, Δ, i, j, args...) = ( c¹ - getbc(bc, i, j, args...) ) / (Δ/2)
@inline right_gradient(bc::VBC, cᴺ, Δ, i, j, args...) = ( getbc(bc, i, j, args...) - cᴺ ) / (Δ/2)

function _fill_west_halo!(j, k, grid, c, bc::Union{VBC, GBC}, loc, args...)

           #  ↑ x ↑  interior
           #  -----  interior face
    iᴵ = 1 #    *    interior cell
    iᴮ = 1 #  =====  western boundary
    iᴴ = 0 #    *    halo cell

    LX, LY, LZ = loc
    Δ = Δx(iᴮ, j, k, grid, flip(LX), LY, LZ) # Δ between first interior and first west halo point, defined at cell face.
    @inbounds ∇c = left_gradient(bc, c[iᴵ, j, k], Δ, j, k, grid, args...)
    @inbounds c[iᴴ, j, k] = linearly_extrapolate(c[iᴵ, j, k], ∇c, -Δ) # extrapolate westward in -x direction.
end

function _fill_east_halo!(j, k, grid, c, bc::Union{VBC, GBC}, loc, args...)

                     #  ↑ x ↑
    iᴴ = grid.Nx + 1 #    *   halo cell
    iᴮ = grid.Nx + 1 #  ===== eastern boundary
    iᴵ = grid.Nx     #    *   interior cell
                     #  ----- interior face
                     #    ↓   interior

    LX, LY, LZ = loc
    Δ = Δx(iᴮ, j, k, grid, flip(LX), LY, LZ) # Δ between last interior and first east halo point, defined at cell face.
    @inbounds ∇c = right_gradient(bc, c[iᴵ, j, k], Δ, j, k, grid, args...)
    @inbounds c[iᴴ, j, k] = linearly_extrapolate(c[iᴵ, j, k], ∇c, Δ) # extrapolate eastward in +x direction.
end

function _fill_south_halo!(i, k, grid, c, bc::Union{VBC, GBC}, loc, args...)

           #  ↑ y ↑  interior
           #  -----  interior face
    jᴵ = 1 #    *    interior cell
    jᴮ = 1 #  =====  southern boundary
    jᴴ = 0 #    *    halo cell

    LX, LY, LZ = loc
    Δ = Δy(i, jᴮ, k, grid, LX, flip(LY), LZ) # Δ between first interior and first south halo point, defined at cell face.
    @inbounds ∇c = left_gradient(bc, c[i, jᴵ, k], Δ, i, k, grid, args...)
    @inbounds c[i, jᴴ, k] = linearly_extrapolate(c[i, jᴵ, k], ∇c, -Δ) # extrapolate southward in -y direction.
end

function _fill_north_halo!(i, k, grid, c, bc::Union{VBC, GBC}, loc, args...)

                     #  ↑ y ↑
    jᴴ = grid.Ny + 1 #    *   halo cell
    jᴮ = grid.Ny + 1 #  ===== northern boundary
    jᴵ = grid.Ny     #    *   interior cell
                     #  ----- interior face
                     #    ↓   interior

    LX, LY, LZ = loc
    Δ = Δy(i, jᴮ, k, grid, LX, flip(LY), LZ) # Δ between first interior and first north halo point, defined at cell face.
    @inbounds ∇c = right_gradient(bc, c[i, jᴵ, k], Δ, i, k, grid, args...)
    @inbounds c[i, jᴴ, k] = linearly_extrapolate(c[i, jᴵ, k], ∇c, Δ) # extrapolate northward in +y direction.
end

function _fill_bottom_halo!(i, j, grid, c, bc::Union{VBC, GBC}, loc, args...)

           #  ↑ z ↑  interior
           #  -----  interior face
    kᴵ = 1 #    *    interior cell
    kᴮ = 1 #  =====  bottom boundary
    kᴴ = 0 #    *    halo cell

    LX, LY, LZ = loc
    Δ = Δz(i, j, kᴮ, grid, LX, LY, flip(LZ)) # Δ between first interior and first bottom halo point, defined at cell face.
    @inbounds ∇c = left_gradient(bc, c[i, j, kᴵ], Δ, i, j, grid, args...)
    @inbounds c[i, j, kᴴ] = linearly_extrapolate(c[i, j, kᴵ], ∇c, -Δ) # extrapolate downward in -z direction.
end

function _fill_top_halo!(i, j, grid, c, bc::Union{VBC, GBC}, loc, args...)

                     #  ↑ z ↑
    kᴴ = grid.Nz + 1 #    *    halo cell
    kᴮ = grid.Nz + 1 #  =====  top boundary
    kᴵ = grid.Nz     #    *    interior cell
                     #  -----  interior face

    LX, LY, LZ = loc
    Δ = Δz(i, j, kᴮ, grid, LX, LY, flip(LZ)) # Δ between first interior and first top halo point, defined at cell face.
    @inbounds ∇c = right_gradient(bc, c[i, j, kᴵ], Δ, i, j, grid, args...)
    @inbounds c[i, j, kᴴ] = linearly_extrapolate(c[i, j, kᴵ], ∇c, Δ) # extrapolate upward in +z direction.
end


text:
using Compat, ClassicalLaminateTheory

forcesAndMoments =
 @compat Dict(
    :Nx => 1000,
    :Ny => 300,
    :Nxy => 150,
    :Mx => 2.5,
    :My => 5,
    :Mxy => 10)

materialProperties =@compat Dict(
  :E1 => [20.0e6, 13.0*10^6],
  :E2 => [1.5*10^6, 1.0*10^6],
  :G12 => [1.0*10^6, 0.75*10^6],
  :ν12 => [0.34, 0.34],
  :alpha1 => [0.2e-6, 0.2e-5],
  :alpha2 => [0.2e-4, 0.2e-3],
  :beta1 => [0.2e-6, 0.2e-5],
  :beta2 => [0.2e-4, 0.2e-3],
  :thickness => [0.0075, 0.0055])

laminateProperties =@compat Dict(
  :nplies => 4,
  :orientation => [0, 90, 0, 90],
  :materials => [2, 2, 1, 1],
  :repeats => 2,
  :symmetric => true,
  :symmetricrepeats => 2,
  :deltaTemp => -280.0,
  :deltaMoisture => 2
)

createLaminate!(laminateProperties, materialProperties, forcesAndMoments)

fm = [
-3382.29,
-4082.29,
150.0,
2.5,
5.0,
10.0
]

@assert round(laminateProperties[:fTotal], 1) == round(fm, 1)
println("Tests for CLT4 passed!")


text:
function polyDNRR(x::Array{Float64}, degree::Int64, k::Int64, gibbsIter::Int64, burnIn::Int64, dpr1::Float64, dpr2::Float64, wdist::Float64, gEps1::Float64, gEps2::Float64, thLow::Float64,
  thUp::Float64, zLow::Float64, zUp::Float64, papr::Float64, pbpr::Float64, thin::Int64, samplerSeed::Int64)

  srand(samplerSeed)

  filename = "/Results"
  if k .> 0
    savelocation = string(pwd(), filename, "/seed$samplerSeed")
  else
    savelocation = string(pwd(), filename, "/seed$samplerSeed/Reco")
  end
  mkpath(savelocation)
  # initialization
  ss = length(x)

  degree += 1 # {n-deg => (n+1)-parameters}

  nn = Int((gibbsIter - burnIn) / thin)

  theta = zeros(degree) # polynomial coeeficients vector
  sampledTheta = zeros(nn, degree) # matrix to store sample thetas
  thetaprev = copy(theta)

  N = collect(1:1:ss)
  d = collect(1:1:ss)

  noise = zeros(nn)

  clusters = zeros(nn)

  ps = zeros(nn)
  p = 1.

  x0s = zeros(nn)
  sx0 = 0.5

  if k .> 0 #introduce k-length denoised variables

    deltas = zeros(nn)
    delta = 1e-04

    ys = zeros(nn, k)
    sacost = zeros(nn,2)

    dens = zeros(k)
    for j = 1:k
        dens[j] = x[j]#rand(Uniform(zLow, zUp))
    end
   aprobs = zeros(nn)

  end

  toler = 1e-06

  display("Starting MCMC ...")

  ##########################################################################
  # MCMC
  ##########################################################################

  for iter in 1:gibbsIter

      # Compute weights up to ss
      ##########################################################################

      M = maximum(N)
      w = zeros(M)

      for i in 1:M
        w[i] = p * (1 - p)^(i - 1)
      end

      ## Sample precisions
      ##########################################################################

      Lambdas = zeros(M)

      for j in 1:M
        count = 0.0
        temp = 0.0
        if d[1] .== j
          count = count + 1.0
          temp = (x[1] - polyMap(theta, sx0)) ^ 2
        end
        for i in 2:ss
          if d[i] .== j
            count = count + 1.0
            temp = temp + (x[i] - polyMap(theta, x[i-1])) ^ 2
          end
        end
        shapel = gEps1 + 0.5 * count
        ratel = gEps2 + 0.5 * temp
        Lambdas[j] = rand(Gamma(shapel, 1. / ratel))
      end

      ## Sample the indicator variables d(i) for i=1,...,ss
      ##########################################################################

      nc1 = 0.0
      prob1 = 0.0
      for j in 1:1:N[1]
        nc1 = nc1 + Lambdas[j]^0.5 * exp(-0.5 * Lambdas[j] * (x[1] - polyMap(theta, sx0))^2)
      end
      rd = rand()
      for j in 1:1:N[1]
        prob1 = prob1 + Lambdas[j]^0.5 * exp(-0.5 * Lambdas[j] * (x[1] - polyMap(theta, sx0))^2) / nc1
        if rd .< prob1
          d[1] = j
          break
        end
      end

      for i in 2:1:ss
        nc = 0.0
        prob = 0.0
        for j in 1:1:N[i]
          nc = nc + Lambdas[j]^0.5 * exp(-0.5 * Lambdas[j] * (x[i] - polyMap(theta, x[i-1]))^2)
        end
        rd = rand()
        for j in 1:1:N[i]
          prob = prob + Lambdas[j]^0.5 * exp(-0.5 * Lambdas[j] * (x[i] -polyMap(theta, x[i-1]))^2) / nc
          if rd .< prob
            d[i] = j
            break
          end
        end
      end

      # Sample N_i auxilliary variables
      ##########################################################################

      N = etgeornd(p, d)

      ## Number of clusters
      ##########################################################################

      ucl = length(unique(d))
      # println("clusters: $ucl")
      # @printf("uCl: %s \n", ucl)

      ## Sample geometricx probability
      ##########################################################################

      p = rand(Beta(papr + 2 * ss, pbpr + sum(N) - ss))

        ## Sample the vector with the coefficients
        ##########################################################################

        for j in 1:1:length(theta)
          thetaj = copy(theta)
          thetaj[j] = 0.

          tauj = Lambdas[d[1]] * sx0 ^ (2(j-1))
          meanj = Lambdas[d[1]] * sx0 ^ (j-1) * (x[1] - polyMap(thetaj, sx0))
          for i in 2:1:ss
            tauj += Lambdas[d[i]] * x[i-1] ^ (2(j-1))
            meanj += Lambdas[d[i]] * x[i-1] ^ (j-1) * (x[i] - polyMap(thetaj, x[i-1]))
          end
          meanj /= tauj

          vj = -(2. / tauj) * log.(rand()) + (thetaprev[j] - meanj) ^ 2
          theta[j] = rand(Uniform(max(-sqrt(vj)+meanj, thLow), min(sqrt(vj)+meanj, thUp)))

        end

        thetaprev = copy(theta)

        ## Sample initial condition
        ##########################################################################

        # aux = - 2 / Lambdas[d[1]] * log.(rand()) + (x[1] - polyMap(theta, sx0)) ^ 2
        # polyRight = copy(theta)
        # polyRight[1] = polyRight[1] - (x[1] - sqrt(aux))
        # polyRight = polyRight ./ polyRight[end]
        # rootsRight = Polynomials.roots(Poly(polyRight))

        # polyLeft = copy(theta)
        # polyLeft[1] = polyLeft[1] - (x[1] + sqrt(aux))
        # polyLeft = polyLeft ./ polyLeft[end]
        # rootsLeft = Polynomials.roots(Poly(polyLeft))
        # # println("allRoots = $rootsLeft,$rootsRight")
        # allRoots = [rootsLeft rootsRight]
        # #println("allRoots = $allRoots")
        # sel = abs.(imag(allRoots)) .< toler # treat as real the roots with imagine part < ϵ
        # realRoots = sort(real(allRoots[sel]))
        # intervals = rangeIntersection(realRoots, [zLow zUp])
        # sx0 = unifmixrnd(intervals)

        # sample x0
        aux = -(2.0/Lambdas[d[1]])*log(rand()) + (x[1]-polyMap(theta,sx0))^2
        poly_right = copy(theta)
        poly_right[1] = poly_right[1] - (x[1] - sqrt(aux))
        poly_right = poly_right./poly_right[end]
        roots_right = roots(poly_right)
        poly_left = copy(theta)
        poly_left[1] = poly_left[1] - (x[1] + sqrt(aux))
        poly_left = poly_left./poly_left[end]
        roots_left = roots(poly_left)
        allroots = [roots_left roots_right]
        idx = abs.(imag(allroots)).<1e-08
        realroots = sort(real(allroots[idx]))
        intervals = rangeIntersection(realroots,[-20. 20.])
        sx0 = unifmixrnd(intervals)

        if k.> 0

          pvar = 1e-06

          mAccProbs = zeros(k) # vector of acceptance probabilities

          ## Sample xden[1]
          ##########################################################################

          yprop = rand(Normal(dens[1], sqrt(pvar)))

          accProb = min(1, exp(- 0.5 * delta * ( cfy(yprop, sx0, dens[2], theta) - cfy(dens[1], sx0, dens[2], theta) ) - 0.5 * wdist *
                                              ( (yprop - x[1])^2 - (dens[1] - x[1])^2 ) ) )

          if rand() .< accProb
            dens[1] = yprop
          end

          mAccProbs[1] = accProb

          ## Sample xden[2] - xden[ss-1]
          ##########################################################################

          for j in 2:(k - 1)

            yprop = rand(Normal(dens[j], sqrt(pvar)))

            accProb = min(1, exp(- 0.5 * delta * ( cfy(yprop, dens[j-1], dens[j+1], theta) - cfy(dens[j], dens[j-1], dens[j+1], theta) ) - 0.5 * wdist *
                                              ( (yprop - x[j])^2 - (dens[j] - x[j])^2 ) ) )

            if rand() .< accProb
              dens[j] = yprop
            end

            mAccProbs[j] = accProb

          end

          ## Sample xden[k]
          ##########################################################################

          yprop = rand(Normal(dens[k], sqrt(pvar)))

          accProb = min(1, exp(- 0.5 * delta * ( (yprop - polyMap(theta, dens[k-1]))^2 - (dens[k] - polyMap(theta, dens[k-1]))^2 ) - 0.5 * wdist *
                                              ( (yprop - x[k])^2 - (dens[k] - x[k])^2 ) ) )

          if rand() .< accProb
            dens[k] = yprop
          end

          mAccProbs[k] = accProb

          ## Sample denoising parameter delta
          ##########################################################################

          dtemp = (dens[1] - polyMap(theta, sx0)) ^ 2
          for i in 2:ss
            dtemp += (dens[i] - polyMap(theta, dens[i-1])) ^ 2
          end

          drate = dpr2 + 0.5 * dtemp
          delta = rand(Gamma(dpr1 + 0.5 * ss, 1./drate))

        end

        ## After Burn-In period
        ##########################################################################

      if (iter .> burnIn) & ((iter-burnIn) % thin .== 0)

        ii = Int((iter - burnIn)/thin)

        if k .> 0

          ## Calculate cost at each iteration (E_dyn)
          ##########################################################################

          tt = (dens[1] - polyMap(theta, sx0)) ^ 2

          for i in 2:length(dens)
            tt +=  (dens[i] - polyMap(theta, dens[i-1])) ^ 2
          end

          sacost[ii,1] = sqrt(1/length(dens) *tt)

          ## Calculate cost at each iteration (E0)
          ##########################################################################

          tt = 0.

          for i in 1:length(dens)
            tt +=  (dens[i] - x[i]) ^ 2
          end

         sacost[ii,2] = sqrt(1/length(dens) *tt)
       end
          ## Sample noise predictive
          ##########################################################################

          cW = cumsum(w)
          flag = rand()
          if cW[end] .< flag
            pred = rand(Normal(0, sqrt(1 / rand(Gamma(gEps1, 1 / gEps2))))) # draw from the prior
          else
            for j in 1:1:length(cW)
              if flag .< cW[j]
                 pred = rand(Normal(0, sqrt(1 / Lambdas[j])))
                 break
              end
            end
          end

          # Store values
          sampledTheta[ii, :] = thetaprev
          x0s[ii, :] = sx0
          noise[ii] = pred
          clusters[ii] = ucl
          ps[ii] = p
          if k .> 0
            aprobs[ii] = mean(mAccProbs)
            deltas[ii] = delta
            # distws[ii] = wdist
            ys[ii, :] = dens
          end

        end

        if iter % 10000 .== 0
           println("MCMC Iterations: $iter")
        end

  end

  display("MCMC finished !")

  ## Write values in .txt files - specific path
  ##########################################################################

  writedlm(string(savelocation,"/data.txt"), x, '\n')
  writedlm(string(savelocation,"/thetas.txt"), sampledTheta)
  writedlm(string(savelocation,"/x0s.txt"), x0s, '\n')
  writedlm(string(savelocation,"/noise.txt"), noise, '\n')
  writedlm(string(savelocation,"/ucl.txt"), clusters, '\n')
  writedlm(string(savelocation,"/ps.txt"), ps, '\n')

  if k.> 0
    writedlm(string(savelocation,"/probs.txt"), aprobs, '\n')
    writedlm(string(savelocation,"/ys.txt"), ys)
    writedlm(string(savelocation,"/deltas.txt"), deltas, '\n')
    writedlm(string(savelocation,"/cost.txt"), sacost)
  end

  return 0

  end

max_stars_repo_name: KosKal/DNRR

max_issues_repo_name: KosKal/DNRR

max_forks_repo_name: KosKal/DNRR


text:
#####
##### StepRange
#####

# Show we could interop with `StepRange`, which is another contender for a TimeSpans alternative
function Base.StepRange(span::AlignedSpans.AlignedSpan)
    # the rounding here is not ideal
    t = Nanosecond(round(Int, TimeSpans.nanoseconds_per_sample(span.sample_rate)))
    return (span.first_index * t):t:(span.last_index * t)
end

function AlignedSpan(r::StepRange{T,S}) where {T<:Period,S<:Period}
    sample_rate = TimeSpans.NS_IN_SEC / Dates.value(convert(Nanosecond, step(r)))
    i = first(r) / step(r)
    j = last(r) / step(r)
    return AlignedSpan(sample_rate, Int(i), Int(j))
end

@testset "StepRange roundtripping" begin
    for sample_rate in [1.0, 0.5, 100.0, 128.33]
        for (i, j) in [1 => 10, 5 => 20, 3 => 6, 78 => 79]
            span = AlignedSpan(sample_rate, i, j)
            r = StepRange(span)
            span2 = AlignedSpan(r)
            @test span.first_index == span2.first_index
            @test span.last_index == span2.last_index
            @test span.sample_rate ≈ span2.sample_rate rtol = 1e-7 # annoyingly imprecise...
        end
    end
end

#####
##### TimeSpans
#####

@testset "TimeSpans roundtripping" begin
    for sample_rate in [1.0, 0.5, 100.0, 128.33]
        # AlignedSpan -> TimeSpan -> AlignedSpan
        for (i, j) in [1 => 10, 5 => 20, 3 => 6, 78 => 79]
            for mode in (RoundSpanDown, RoundInward, ConstantSamplesRoundingMode(RoundDown))
                as = AlignedSpan(sample_rate, i, j)
                ts = TimeSpan(as)

                # We don't perfectly roundtrip with sample rate 128.33 for two of the rounding modes.
                # However, if we widen the timespan by 1 nanosecond, we can roundtrip in that case.
                # So let's test that, at least.
                if sample_rate == 128.33
                    if mode == RoundInward
                        ts = TimeSpan(max(Nanosecond(0), start(as) - Nanosecond(1)),
                                      stop(as))
                    elseif mode == ConstantSamplesRoundingMode(RoundDown)
                        ts = TimeSpan(start(as), stop(as) + Nanosecond(1))
                    end
                end

                rt = AlignedSpan(sample_rate, ts, mode)
                @test as == rt
            end
        end
    end
end

@test TimeSpans.istimespan(AlignedSpan(1, 1, 1))

#####
##### Onda
#####

@testset "Samples indexing" begin
    for sample_rate in [1.0, 0.5, 100.0, 128.33]
        samples = make_test_samples(sample_rate)
        for (i, j) in [1 => 10, 5 => 20, 3 => 6, 78 => 79]
            @test samples[:, AlignedSpan(sample_rate, i, j)] == samples[:, i:j]
            @test_throws ArgumentError samples[:, AlignedSpan(sample_rate + 1, i, j)]
        end
    end
end

#####
##### JSON/Arrow
#####

spans = [AlignedSpan(1.0, 5, 10), AlignedSpan(111.345, 500, 10000)]
@testset "JSON3 roundtripping" begin
    @test JSON3.read(JSON3.write(spans), Vector{AlignedSpan}) == spans
end

@testset "Arrow roundtripping" begin
    @test Arrow.Table(Arrow.tobuffer((; spans))).spans == spans
end


text:
module PhysicalFESpacesTests

using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Arrays
using Gridap.Fields
using Gridap.FESpaces
using Gridap.Polynomials
using Gridap.CellData
using Gridap.TensorValues
using Test

# domain = (0,1)
# partition = (3,)
domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)
order = 1

trian = get_triangulation(model)
quad = CellQuadrature(trian,order)

u(x) = x[1]
T = Float64

u(x) = VectorValue(x[1],x[2])
T = VectorValue{2,Float64}

# Juno.@run TestFESpace( reffe = :Lagrangian, conformity = :H1, valuetype=T, order = order, model = model, dirichlet_tags = [1,6], dof_space = :physical)
# Juno.@enter TestFESpace( reffe = :Lagrangian, conformity = :H1, valuetype=T, order = order, model = model, dirichlet_tags = [1,6], dof_space = :physical)
# Juno.@enter TestFESpace( reffe = :Lagrangian, conformity = :H1, valuetype=T, order = order, model = model, dirichlet_tags = [1,6])

Vp = TestFESpace(
  reffe = :Lagrangian,
  conformity = :H1,
  valuetype=T,
  order = order,
  model = model,
  dirichlet_tags = [1,6],
  dof_space = :physical)

V = TestFESpace(
  reffe = :Lagrangian,
  conformity = :H1,
  valuetype=VectorValue{2,Float64},
  order = order,
  model = model,
  dirichlet_tags = [1,6])

  V.cell_basis.ref_trait

# Vp = TestFESpace(
#   reffe = :RaviartThomas,
#   conformity = :Hdiv,
#   order = order,
#   model = model,
#   dirichlet_tags = [1,6],
#   dof_space = :physical)
#
# V = TestFESpace(
#   reffe = :RaviartThomas,
#   conformity = :Hdiv,
#   order = order,
#   model = model,
#   dirichlet_tags = [1,6])
# test_single_field_fe_space(V)

Up = TrialFESpace(Vp,u)
U = TrialFESpace(V,u)
U.cell_basis.ref_trait
U.space.cell_basis.ref_trait

# FEM = Gridap.FESpaces
# Juno.@enter TrialFESpace(Vp,u)
# dirichlet_values = FEM.compute_dirichlet_values_for_tags(V,u)

q = get_coordinates(quad)
r = evaluate(U.cell_basis,q)
rp = evaluate(Up.cell_basis,q)

@test all(collect(r) .≈ collect(rp))

@test RefStyle(Up.space.cell_dof_basis) == Val{false}()
@test RefStyle(U.space.cell_dof_basis) == Val{true}()

@test RefStyle(Up.space.cell_basis) == Val{false}()
@test RefStyle(U.space.cell_basis) == Val{true}()

@test RefStyle(Vp.cell_basis) == Val{false}()
@test RefStyle(V.cell_basis) == Val{true}()

uhf = interpolate(u,Up)
uhd = interpolate_dirichlet(u,Up)
uhdf = interpolate_everywhere(u,Up)
uhf.dirichlet_values
uhdf.dirichlet_values
uhd.dirichlet_values

uhf_r = interpolate(u,U)
uhd_r= interpolate_dirichlet(u,U)
uhdf_r = interpolate_everywhere(u,U)
uhd_r.dirichlet_values

@test uhf.free_values == uhf_r.free_values
@test uhf.dirichlet_values == uhf_r.dirichlet_values
@test uhd.dirichlet_values == uhd_r.dirichlet_values
@test uhdf.free_values == uhdf_r.free_values
@test uhdf.dirichlet_values == uhdf_r.dirichlet_values

@test uhf.free_values == uhdf.free_values
@test uhf.dirichlet_values == uhd.dirichlet_values
@test uhdf.free_values == uhf.free_values
@test uhdf.dirichlet_values == uhd.dirichlet_values
@test uhdf.dirichlet_values == uhdf_r.dirichlet_values

e = u - uhf_r

el2 = sqrt(sum(integrate(inner(e,e),trian,quad)))

@test el2 < 1.0e-10

e = u - uhf

el2 = sqrt(sum(integrate(inner(e,e),trian,quad)))

@test el2 < 1.0e-10

end #module


text:
function _add_var_container!(canonical_model::CanonicalModel, var_name::Symbol, ax1, ax2)
    canonical_model.variables[var_name] = _container_spec(canonical_model.JuMPmodel, ax1, ax2)
    return
end

function _add_cons_container!(canonical_model::CanonicalModel, cons_name::Symbol, ax1, ax2)
    canonical_model.constraints[cons_name] = JuMPConstraintArray(undef, ax1, ax2)
    return
end

function _add_param_container!(canonical_model::CanonicalModel, param_reference::UpdateRef, axs...)
    canonical_model.parameters[param_reference] = JuMP.Containers.DenseAxisArray{PJ.ParameterRef}(undef, axs...)
    return
end

function _middle_rename(original::Symbol, split_char::String, addition::String)

    parts = split(String(original),split_char)

    return Symbol(parts[1],"_",addition,"_",parts[2])

end

function _remove_underscore(original::Symbol)

    if !occursin("_", String(original))
        return original
    end

    parts = split(String(original),"_")

    return parts[1]
end

"Replaces the string in `char` with the string`replacement`"
function replace_chars(s::String, char::String, replacement::String)
    return replace(s, Regex("[$char]") => replacement)
end

"Removes the string `char` from the original string"
function remove_chars(s::String, char::String)
    return replace_chars(s::String, char::String, "")
end


text:
function JFVM_test()
# JFVM test script
# This script is supposed to call all the functions of the JFVM package
## Part I: creating an array of different mesh types:
# domain size
Lx= 1.0
Ly= 2*pi
Lz= 2.0
Nx=5
Ny=7
Nz=9
X=[0.01, 0.1, 0.3, 0.5, 0.55, 1.0]
Y= [0.0, 0.1, 1.0, 1.5, 2.9, 3.0, pi, 2*pi]
Z= [0.0, 0.01, 0.1, 0.5, 0.7, 0.95, 1.0, 1.25, 1.39, 2.0]
N_mesh=7
# create nonuniform mesh
mesh_nonuniform= cell(N_mesh)
mesh_nonuniform[1]=createMesh1D(X)
mesh_nonuniform[2]=createMesh2D(X, Y)
mesh_nonuniform[3]=createMesh3D(X, Y, Z)
mesh_nonuniform[4]=createMeshCylindrical1D(X)
mesh_nonuniform[5]=createMeshCylindrical2D(X, Y)
mesh_nonuniform[6]=createMeshCylindrical3D(X, Y, Z)
mesh_nonuniform[7]=createMeshRadial2D(X, Y)
println("Non-uniform mesh created successfully!")
## Part II: create cell and face variables
c_val= 1.0
# nonuniform
c_n=cell(N_mesh)
c_r=cell(N_mesh)
D_n=cell(N_mesh)
for i=1:N_mesh
    c_n[i]= createCellVariable(mesh_nonuniform[i], c_val)
end
for i=1:N_mesh
    D_n[i]= createCellVariable(mesh_nonuniform[i], c_val)
end
println("Cells of fixed values over nonuniform mesh created successfully!")
for i=1:N_mesh
    c_rand = rand(tuple(mesh_nonuniform[i].dims...))
    c_r[i]= createCellVariable(mesh_nonuniform[i], c_rand)
end
println("Cells of random values over nonuniform mesh created successfully!")
## Part III: create face variables
f_val= 0.5
# nonuniform
f_n=cell(N_mesh)
for i=1:N_mesh
    f_n[i]= createFaceVariable(mesh_nonuniform[i], f_val)
end
println("Face variable over nonuniform mesh created successfully!")
## Part IV: Test boundary conditions
BC_n=cell(N_mesh)
for i=1:N_mesh
    BC_n[i]=createBC(mesh_nonuniform[i])
    BC_n[i].left.a[:]=0.0
    BC_n[i].left.b[:]=1.0
    BC_n[i].left.c[:]=1.0
    BC_n[i].right.a[:]=0.0
    BC_n[i].right.b[:]=1.0
    BC_n[i].right.c[:]=0.0
end
println("Boundary condition over a nonuniform mesh created successfully!")
## Part V: solve a steady-state diffusion equation
c_dif=cell(N_mesh)
M_bc=cell(N_mesh)
M_dif=cell(N_mesh)
RHS_bc=cell(N_mesh)
for i=1:N_mesh
    M_bc[i], RHS_bc[i]= boundaryConditionTerm(BC_n[i])
    M_dif[i]=diffusionTerm(f_n[i])
    c_dif[i]=solveLinearPDE(mesh_nonuniform[i], M_dif[i]+M_bc[i], RHS_bc[i])
end
# visualize
# figure(1)
# for i=1:N_mesh
#     subplot(3, 3, i)
#     visualizeCells(c_dif[i])
# end
# println("Diffusion equation solved and visualized successfully")
## Part VI: solve convection diffucion equation
# nonuniform
c_conv=cell(N_mesh)
M_bc=cell(N_mesh)
M_dif=cell(N_mesh)
M_conv=cell(N_mesh)
RHS_bc=cell(N_mesh)
for i=1:N_mesh
    M_bc[i], RHS_bc[i]= boundaryConditionTerm(BC_n[i])
    M_dif[i]=diffusionTerm(f_n[i])
    M_conv[i]=convectionTerm(0.1*f_n[i])
    c_conv[i]=solveLinearPDE(mesh_nonuniform[i], M_conv[i]-M_dif[i]+M_bc[i], RHS_bc[i])
end
# visualize
# figure(2)
# for i=1:N_mesh
#     subplot(3, 3, i)
#     visualizeCells(c_conv[i])
# end
# println("Convection-Diffusion equation solved and visualized successfully")
## Part VII: test the calculus fanctions
grad_c=cell(N_mesh)
for i=1:N_mesh
    grad_c[i]=gradientTerm(c_dif[i])
end
div_c=cell(N_mesh)
for i=1:N_mesh
    div_c[i]=divergenceTerm(grad_c[i])
end
println("Gradient and divergence functions work fine!")
## Solve a dynamic equation
dt=0.1
c_old=c_n
c_trans=cell(N_mesh)
M_bc=cell(N_mesh)
M_dif=cell(N_mesh)
M_conv=cell(N_mesh)
RHS_bc=cell(N_mesh)
M_ls=cell(N_mesh)
RHS_s=cell(N_mesh)
for i=1:N_mesh
    M_bc[i], RHS_bc[i]= boundaryConditionTerm(BC_n[i])
    M_dif[i]=diffusionTerm(0.1*f_n[i])
    M_conv[i]=convectionUpwindTerm(0.01*f_n[i])
    M_ls[i]=linearSourceTerm(0.1*c_n[i])
    RHS_s[i]=constantSourceTerm(0.2*c_n[i])
end

for i=1:N_mesh
    for j=1:10
        (M_t, RHS_t)=transientTerm(c_old[i], dt)
        c_trans[i]=solveLinearPDE(mesh_nonuniform[i],
            M_t+M_ls[i]+M_conv[i]-M_dif[i]+M_bc[i], RHS_t+RHS_s[i]+RHS_bc[i])
        c_old[i]=copyCell(c_trans[i])
    end
end
# visualize
# figure(3)
#for i=1:N_mesh
#     visualizeCells(c_trans[i])
#     pause(1.5)
#end
println("Transient convection-diffucion-reaction solved successfully!")
## Part VIII: test the utilities
# only test the averaging, don"t save the result
FL=fluxLimiter("Superbee")
for i=1:N_mesh
    linearMean(c_trans[i])
    arithmeticMean(c_trans[i])
    geometricMean(D_n[i])
    harmonicMean(D_n[i])
    upwindMean(c_trans[i], f_n[i])
    tvdMean(c_trans[i], f_n[i], FL)
end
println("Averaging functions run smoothly!")
## Part IX: test the classes and operators

end # end function


text:
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{St, D},
    ::Type{S},
) where {St <: PSY.Storage, D <: AbstractStorageFormulation, S <: PM.AbstractPowerModel}
    devices = get_available_components(St, sys)

    add_variables!(container, ActivePowerInVariable, devices, D())
    add_variables!(container, ActivePowerOutVariable, devices, D())
    add_variables!(container, ReactivePowerVariable, devices, D())
    add_variables!(container, EnergyVariable, devices, D())
    if get_attribute(model, "reservation")
        add_variables!(container, ReservationVariable, devices, D())
    end

    initial_conditions!(container, devices, D())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_feedforward_arguments!(container, model, devices)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{St, D},
    ::Type{S},
) where {St <: PSY.Storage, D <: AbstractStorageFormulation, S <: PM.AbstractPowerModel}
    devices = get_available_components(St, sys)

    add_constraints!(
        container,
        OutputActivePowerVariableLimitsConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        InputActivePowerVariableLimitsConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(container, EnergyCapacityConstraint, EnergyVariable, devices, model, S)

    # Energy Balanace limits
    add_constraints!(container, EnergyBalanceConstraint, devices, model, S)

    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{St, D},
    ::Type{S},
) where {
    St <: PSY.Storage,
    D <: AbstractStorageFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(St, sys)

    add_variables!(container, ActivePowerInVariable, devices, D())
    add_variables!(container, ActivePowerOutVariable, devices, D())
    add_variables!(container, EnergyVariable, devices, D())
    if get_attribute(model, "reservation")
        add_variables!(container, ReservationVariable, devices, D())
    end

    initial_conditions!(container, devices, D())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )
    add_feedforward_arguments!(container, model, devices)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{St, D},
    ::Type{S},
) where {
    St <: PSY.Storage,
    D <: AbstractStorageFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(St, sys)

    add_constraints!(
        container,
        OutputActivePowerVariableLimitsConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        InputActivePowerVariableLimitsConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_constraints!(container, EnergyCapacityConstraint, EnergyVariable, devices, model, S)

    # Energy Balanace limits
    add_constraints!(container, EnergyBalanceConstraint, devices, model, S)

    add_feedforward_constraints!(container, model, devices)

    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{St, EnergyTarget},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractPowerModel}
    devices = get_available_components(St, sys)

    add_variables!(container, ActivePowerInVariable, devices, EnergyTarget())
    add_variables!(container, ActivePowerOutVariable, devices, EnergyTarget())
    add_variables!(container, ReactivePowerVariable, devices, EnergyTarget())
    add_variables!(container, EnergyVariable, devices, EnergyTarget())
    add_variables!(container, EnergyShortageVariable, devices, EnergyTarget())
    add_variables!(container, EnergySurplusVariable, devices, EnergyTarget())
    if get_attribute(model, "reservation")
        add_variables!(container, ReservationVariable, devices, EnergyTarget())
    end

    add_parameters!(container, EnergyTargetTimeSeriesParameter, devices, model)

    initial_conditions!(container, devices, EnergyTarget())

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_feedforward_arguments!(container, model, devices)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{St, EnergyTarget},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractPowerModel}
    devices = get_available_components(St, sys)

    add_constraints!(
        container,
        OutputActivePowerVariableLimitsConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        InputActivePowerVariableLimitsConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(container, EnergyCapacityConstraint, EnergyVariable, devices, model, S)

    # Energy Balanace limits
    add_constraints!(container, EnergyBalanceConstraint, devices, model, S)
    add_constraints!(container, EnergyTargetConstraint, devices, model, S)
    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{St, EnergyTarget},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(St, sys)

    add_variables!(container, ActivePowerInVariable, devices, EnergyTarget())
    add_variables!(container, ActivePowerOutVariable, devices, EnergyTarget())
    add_variables!(container, EnergyVariable, devices, EnergyTarget())
    add_variables!(container, EnergyShortageVariable, devices, EnergyTarget())
    add_variables!(container, EnergySurplusVariable, devices, EnergyTarget())
    if get_attribute(model, "reservation")
        add_variables!(container, ReservationVariable, devices, EnergyTarget())
    end

    add_parameters!(container, EnergyTargetTimeSeriesParameter, devices, model)

    initial_conditions!(container, devices, EnergyTarget())

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )
    add_feedforward_arguments!(container, model, devices)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{St, EnergyTarget},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(St, sys)

    add_constraints!(
        container,
        OutputActivePowerVariableLimitsConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        InputActivePowerVariableLimitsConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_constraints!(container, EnergyCapacityConstraint, EnergyVariable, devices, model, S)

    # Energy Balanace limits
    add_constraints!(container, EnergyBalanceConstraint, devices, model, S)
    add_constraints!(container, EnergyTargetConstraint, devices, model, S)
    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{St, BatteryAncillaryServices},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractPowerModel}
    devices = get_available_components(St, sys)

    add_variables!(container, ActivePowerInVariable, devices, BatteryAncillaryServices())
    add_variables!(container, ActivePowerOutVariable, devices, BatteryAncillaryServices())
    add_variables!(container, ReactivePowerVariable, devices, BatteryAncillaryServices())
    add_variables!(container, EnergyVariable, devices, BatteryAncillaryServices())
    if get_attribute(model, "reservation")
        add_variables!(container, ReservationVariable, devices, BatteryAncillaryServices())
    end

    initial_conditions!(container, devices, BatteryAncillaryServices())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    if has_service_model(model)
        add_expressions!(container, ReserveRangeExpressionLB, devices, model)
        add_expressions!(container, ReserveRangeExpressionUB, devices, model)
    end
    add_feedforward_arguments!(container, model, devices)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{St, BatteryAncillaryServices},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractPowerModel}
    devices = get_available_components(St, sys)

    add_constraints!(
        container,
        OutputActivePowerVariableLimitsConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        InputActivePowerVariableLimitsConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(container, EnergyCapacityConstraint, EnergyVariable, devices, model, S)

    # Energy Balanace limits
    add_constraints!(container, EnergyBalanceConstraint, devices, model, S)
    if has_service_model(model)
        add_constraints!(container, ReserveEnergyConstraint, devices, model, S)
        add_constraints!(container, RangeLimitConstraint, devices, model, S)
    end
    add_feedforward_constraints!(container, model, devices)

    add_constraint_dual!(container, sys, model)

    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{St, BatteryAncillaryServices},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(St, sys)

    add_variables!(container, ActivePowerInVariable, devices, BatteryAncillaryServices())
    add_variables!(container, ActivePowerOutVariable, devices, BatteryAncillaryServices())
    add_variables!(container, EnergyVariable, devices, BatteryAncillaryServices())
    if get_attribute(model, "reservation")
        add_variables!(container, ReservationVariable, devices, BatteryAncillaryServices())
    end

    initial_conditions!(container, devices, BatteryAncillaryServices())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )
    if has_service_model(model)
        add_expressions!(container, ReserveRangeExpressionLB, devices, model)
        add_expressions!(container, ReserveRangeExpressionUB, devices, model)
    end
    add_feedforward_arguments!(container, model, devices)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{St, BatteryAncillaryServices},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(St, sys)

    add_constraints!(
        container,
        OutputActivePowerVariableLimitsConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        InputActivePowerVariableLimitsConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_constraints!(container, EnergyCapacityConstraint, EnergyVariable, devices, model, S)

    # Energy Balanace limits
    add_constraints!(container, EnergyBalanceConstraint, devices, model, S)
    if has_service_model(model)
        add_constraints!(container, ReserveEnergyConstraint, devices, model, S)
        add_constraints!(container, RangeLimitConstraint, devices, model, S)
    end

    add_feedforward_constraints!(container, model, devices)

    add_constraint_dual!(container, sys, model)

    return
end


text:
using Clp: ClpSolver
using SparseArrays: spzeros

@testset "Minimum-cost flow" begin

    # bipartite oriented + source & sink
    # source 5, sink 6, v1 {1, 2} v2 {3, 4}
    g = lg.DiGraph(6)
    lg.add_edge!(g, 5, 1)
    lg.add_edge!(g, 5, 2)
    lg.add_edge!(g, 3, 6)
    lg.add_edge!(g, 4, 6)
    lg.add_edge!(g, 1, 3)
    lg.add_edge!(g, 1, 4)
    lg.add_edge!(g, 2, 3)
    lg.add_edge!(g, 2, 4)
    w = zeros(6,6)
    w[1,3] = 10.
    w[1,4] = 5.
    w[2,3] = 2.
    w[2,4] = 2.
    # v2 -> sink have demand of one
    demand = spzeros(6,6)
    demand[3,6] = 1
    demand[4,6] = 1
    capacity = ones(6,6)

    flow = mincost_flow(g, spzeros(lg.nv(g)),capacity, w, ClpSolver(), edge_demand=demand,source_nodes=[5],sink_nodes= [6])
    @test flow[5,1] == 1
    @test flow[5,2] == 1
    @test flow[3,6] == 1
    @test flow[4,6] == 1
    @test flow[1,4] == 1
    @test flow[2,4] == 0
    @test flow[2,3] == 1
    @test flow[1,3] == 0
    @test sum(diag(flow)) == 0

    # flow conservation property
    for n1 = 1:4
        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
    end

    #same test with exact flows
    flow = mincost_flow(g, spzeros(lg.nv(g)),capacity, w, ClpSolver(),edge_demand_exact=true, edge_demand=demand,source_nodes=[5],sink_nodes= [6])
    @test flow[5,1] == 1
    @test flow[5,2] == 1
    @test flow[3,6] == 1
    @test flow[4,6] == 1
    @test flow[1,4] == 1
    @test flow[2,4] == 0
    @test flow[2,3] == 1
    @test flow[1,3] == 0
    @test sum(diag(flow)) == 0

    # flow conservation property
    for n1 = 1:4
        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
    end

    #--- Same test without edge_demands & positive cost should result in null flow
    flow = mincost_flow(g, spzeros(lg.nv(g)),capacity, w, ClpSolver(),source_nodes=[5],sink_nodes= [6])
    @test flow[5,1] == 0
    @test flow[5,2] == 0
    @test flow[3,6] == 0
    @test flow[4,6] == 0
    @test flow[1,4] == 0
    @test flow[2,4] == 0
    @test flow[2,3] == 0
    @test flow[1,3] == 0
    @test sum(diag(flow)) == 0

    # flow conservation property
    for n1 = 1:4
        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
    end

    #run sam test again with nodal demands
    node_demand=spzeros(6)
    node_demand[5]=2
    node_demand[6]=-2
    #--- Same test without edge_demands
    flow = mincost_flow(g, node_demand,capacity, w, ClpSolver())
    @test flow[5,1] == 1
    @test flow[5,2] == 1
    @test flow[3,6] == 1
    @test flow[4,6] == 1
    @test flow[1,4] == 1
    @test flow[2,4] == 0
    @test flow[2,3] == 1
    @test flow[1,3] == 0
    @test sum(diag(flow)) == 0

    # flow conservation property
    for n1 = 1:4
        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
    end

    # no demand => null flow
    d2 = spzeros(6,6)
    flow = mincost_flow(g, spzeros(lg.nv(g)),capacity, w, ClpSolver(), edge_demand=d2,source_nodes=[5],sink_nodes= [6])
    for idx in 1:6
        for jdx in 1:6
            @test flow[idx,jdx] ≈ 0.0
        end
    end

    # graph without sink or source
    g = lg.DiGraph(6)
    # create circuit
    for s in 1:5
        lg.add_edge!(g, s, s+1)
    end
    lg.add_edge!(g, 6, 1)
    lg.add_edge!(g, 2, 5) # shortcut

    capacity = ones(6,6)
    demand = spzeros(6,6)
    demand[1,2] = 1
    costs = ones(6,6)
    flow = mincost_flow(g,spzeros(lg.nv(g)), capacity, costs, ClpSolver(), edge_demand=demand)
    active_flows = [(1,2), (2,5), (5,6),(6,1)]
    for s in 1:6
        for t in 1:6
            if (s,t) in active_flows
                @test flow[s,t] ≈ 1.
            else
                @test flow[s,t] ≈ 0.
            end
        end
    end
    # flow conservation property
    for n1 = 1:6
        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
    end
    # higher short-circuit cost
    costs[2,5] = 10.
    flow = mincost_flow(g,spzeros(lg.nv(g)), capacity, costs, ClpSolver(), edge_demand=demand)
    active_flows = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,1)]
    for s in 1:6
        for t in 1:6
            if (s,t) in active_flows
                @test flow[s,t] ≈ 1.
            else
                @test flow[s,t] ≈ 0.
            end
        end
    end
        # flow conservation property
    for n1 = 1:6
        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
    end
end


text:
include("rocketparts.jl")

function load_engines(path::AbstractString)
    # since we have them we shoud use both the name and company as dictionary key
    rocket_engines = Dict{Tuple{String, String}, SingleEngine}()

    # As Engheim explains, leveraging the fact that open can take a function closure,
    # and then automatically close the file, belows is best way to open file, slurp it in,
    # and then have it automatically closed
    lines = open(readlines, path)
    for line in lines[2:end]
        parts = split(line, ',')
        name, company = parts[1:2]
        # we have to rename the variables, because they are defined as functions
        # in the included file. Engheim is propbaly using a different version of SingleEngine
        basemass, basethrust, throttle, baseIsp =
                  parse.(Float64, parts[3:end])
        # the multipliers are not used or defined anywhere so give errors (in
        # function explanation Engheim leaves them out :))
        # engine = Engine(basethrust * thrust_multiplier,
        #                 baseIsp,
        #                 basemass * mass_multiplier)
        engine = Engine(basethrust,
                        baseIsp,
                        basemass)

        rocket_engines[(name,company)] = engine
    end

    rocket_engines
end

# Example
#
# julia> include("rocketsread.jl")
# load_engines (generic function with 1 method)

# julia> load_engines("data/rocket-engines.csv")
# Dict{Tuple{String, String}, SingleEngine} with 8 entries:
#   ("LV-T30", "Kerbal")            => SingleEngine(205.16, 265.0, 1.25)
#   ("BE-3PM ", "Blue Origin")      => SingleEngine(490.0, 310.0, 0.25)
#   ("LV-T45 ", "Kerbal")           => SingleEngine(167.97, 250.0, 1.5)
#   ("48-7S", "Kerbal")             => SingleEngine(16.88, 270.0, 0.1)
#   ("Merlin 1D", "SpaceX")         => SingleEngine(845.0, 282.0, 0.47)
#   ("RD-180", "NPO Energomash")    => SingleEngine(3830.0, 311.0, 5.48)
#   ("Rutherford", "Rocket Lab")    => SingleEngine(18.0, 303.0, 0.035)
#   ("RS-25", "Aerojet Rocketdyne") => SingleEngine(1860.0, 366.0, 3.527)

function load_engines_alt(path::AbstractString)
    # this is the single line read option. I kept the dict since it seems more useful
    rocket_engines = Dict{Tuple{String, String}, SingleEngine}()

    open(path) do io
        # read in the header line to skip it
        readline(io)
        # iterate with eachline over the rest of the file until hit EOF (using a map closure)
        map(eachline(io)) do line
            parts = split(line, ',')
            name, company = parts[1:2]
            basemass, basethrust, _, baseIsp =
                parse.(Float64, parts[3:end])
            rocket_engines[(name,company)] = Engine(basethrust, baseIsp, basemass)
        end
    end

    rocket_engines
end

# Examplejulia> include("rocketsread.jl")
# load_engines_alt (generic function with 1 method)

# julia> load_engines_alt("data/rocket-engines.csv")
# Dict{Tuple{String, String}, SingleEngine} with 8 entries:
#   ("LV-T30", "Kerbal")            => SingleEngine(205.16, 265.0, 1.25)
#   ("BE-3PM ", "Blue Origin")      => SingleEngine(490.0, 310.0, 0.25)
#   ("LV-T45 ", "Kerbal")           => SingleEngine(167.97, 250.0, 1.5)
#   ("48-7S", "Kerbal")             => SingleEngine(16.88, 270.0, 0.1)
#   ("Merlin 1D", "SpaceX")         => SingleEngine(845.0, 282.0, 0.47)
#   ("RD-180", "NPO Energomash")    => SingleEngine(3830.0, 311.0, 5.48)
#   ("Rutherford", "Rocket Lab")    => SingleEngine(18.0, 303.0, 0.035)
#   ("RS-25", "Aerojet Rocketdyne") => SingleEngine(1860.0, 366.0, 3.527)

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
using DynamicalSystems, Neighborhood

function correlationsum_tree(X, ε; q = 2, metric = Euclidean(), w = 0)
    tree = KDTree(X, metric)
    qf = q == 2 ? identity : x -> float(x)^(q-1)
    _correlationsum_new(tree, X, ε, w, qf)
end

function _correlationsum_tree(tree, X, ε, w, qf)
    theiler = Neighborhood.Theiler(w)
    searchtype = Neighborhood.WithinRange(maximum(ε))
    N = length(X)
    Cq = similar(ε, Float64)
    @inbounds for i in 1:N
        idxs = Neighborhood.isearch(tree, X[i], searchtype, theiler(i))
        isempty(idxs) && continue
        dist = [Neighborhood.evaluate(tree.metric, X[i], X[j]) for j in idxs]
        sort!(dist)
        for k in 1:length(ε)-1
            j = searchsortedfirst(dist, ε[k])
            Cq[k] += qf(j-1)
        end
        Cq[end] += qf(length(dist))
    end
    return Cq ./ (N*qf(N-w-1))
end

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
# ---
# title: 383. Ransom Note
# id: problem383
# author: zhwang
# date: 2022-02-15
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/ransom-note/description/>
# hidden: true
# ---
#
# Given an arbitrary ransom note string and another string containing letters
# from all the magazines, write a function that will return true if the ransom
# note can be constructed from the magazines ; otherwise, it will return false.
#
# Each letter in the magazine string can only be used once in your ransom note.
#
#
#
# **Example 1:**
#
#
#
#     Input: ransomNote = "a", magazine = "b"
#     Output: false
#
#
# **Example 2:**
#
#
#
#     Input: ransomNote = "aa", magazine = "ab"
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: ransomNote = "aa", magazine = "aab"
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * You may assume that both strings contain only lowercase letters.
#
#
## @lc code=start
using LeetCode

function  can_construct(ransom_note::String, magazine::String)::Bool
    words = Dict{Char,Int}(s => 0 for s in 'a':'z')
    for i in magazine
        words[i] += 1
    end
    for i in ransom_note
        words[i] == 0 && return false
        words[i] -= 1
    end
    true
end
## @lc code=end

## @lc test=start
@testset "383.ransom-note.jl" begin
    @test can_construct("a", "b") == false
    @test can_construct("aa", "ab") == false
    @test can_construct("aa", "aba") == true
end
## @lc test=end

## @lc info=start
# link: [solution to 383](https://leetcode-cn.com/problems/ransom-note/solution/383-shu-jin-xin-ji-shu-pan-du-by-rex-rs-o2el/)
## @lc info=end


text:
using Base: PkgId, loaded_modules, package_callbacks
using Base.Meta: isexpr
if isdefined(Base, :mapany)
  const mapany = Base.mapany
else
  mapany(f, A::AbstractVector) = map!(f, Vector{Any}(undef, length(A)), A)
end

export @require

isprecompiling() = ccall(:jl_generating_output, Cint, ()) == 1

loaded(pkg::PkgId) = haskey(Base.loaded_modules, pkg)

const notified_pkgs = [Base.PkgId(UUID(0x295af30fe4ad537b898300126c2a3abe), "Revise")]

const _callbacks = Dict{PkgId, Vector{Function}}()
callbacks(pkg::PkgId) = get!(Vector{Function}, _callbacks, pkg)

listenpkg(@nospecialize(f), pkg::PkgId) =
  loaded(pkg) ? f() : push!(callbacks(pkg), f)

function loadpkg(pkg::PkgId)
  if haskey(_callbacks, pkg)
    fs = _callbacks[pkg]
    delete!(_callbacks, pkg)
    foreach(Base.invokelatest, fs)
  end
end

function withpath(@nospecialize(f), path::String)
  tls = task_local_storage()
  hassource = haskey(tls, :SOURCE_PATH)
  hassource && (path′ = tls[:SOURCE_PATH])
  tls[:SOURCE_PATH] = path
  try
    return f()
  finally
    hassource ?
      (tls[:SOURCE_PATH] = path′) :
      delete!(tls, :SOURCE_PATH)
  end
end

function err(@nospecialize(f), listener::Module, modname::String, file::String, line)
  try
    t = @elapsed ret = f()
    @debug "Requires conditionally ran code in $t seconds: `$listener` detected `$modname`" _file = file _line = line
    ret
  catch exc
    @warn "Error requiring `$modname` from `$listener`" exception=(exc,catch_backtrace())
  end
end

function parsepkg(ex::Expr)
  isexpr(ex, :(=)) || @goto fail
  mod, id = ex.args
  (mod isa Symbol && id isa String) || @goto fail
  return id::String, String(mod::Symbol)
  @label fail
  error("Requires syntax is: `@require Pkg=\"uuid\"`")
end

function withnotifications(@nospecialize(args...))
  for id in notified_pkgs
    if loaded(id)
      mod = Base.root_module(id)
      if isdefined(mod, :add_require)
        add_require = getfield(mod, :add_require)::Function
        add_require(args...)
      end
    end
  end
  return nothing
end

function replace_include(ex::Expr, source::LineNumberNode)
  if ex.head == :call && ex.args[1] === :include && ex.args[2] isa String
    return Expr(:macrocall, :($Requires.$(Symbol("@include"))), source, ex.args[2]::String)
  end
  return Expr(ex.head, (mapany(ex.args) do arg
    isa(arg, Expr) ? replace_include(arg, source) : arg
  end)...)
end

macro require(pkg::Union{Symbol,Expr}, expr)
  pkg isa Symbol &&
    return Expr(:macrocall, Symbol("@warn"), __source__,
                "Requires now needs a UUID; please see the readme for changes in 0.7.")
  idstr, modname = parsepkg(pkg)
  pkg = :(Base.PkgId(Base.UUID($idstr), $modname))
  expr = isa(expr, Expr) ? replace_include(expr, __source__) : expr
  expr = macroexpand(__module__, expr)
  srcfile = string(__source__.file)
  srcline = __source__.line
  quote
    if !isprecompiling()
      listenpkg($pkg) do
        withpath($srcfile) do
          err($__module__, $modname, $srcfile, $srcline) do
            $(esc(:(eval($(Expr(:quote, Expr(:block,
                                            :(const $(Symbol(modname)) = Base.require($pkg)),
                                            expr)))))))
          end
        end
        $withnotifications($srcfile, $__module__, $idstr, $modname, $(esc(Expr(:quote, expr))))
      end
    end
  end
end


text:
@inline cld_fast(x, y) = cld(x, y)
@inline function cld_fast(x::I, y) where {I <: Base.BitInteger}
    d = div_fast(x, y)
    (d + (d * unsigned(y) != unsigned(x))) % I
end
@inline cld_fast(::StaticInt{N}, y) where {N} = cld_fast(N, y)
cld_fast(::StaticInt{N}, ::StaticInt{M}) where {N,M}= (StaticInt{N}() + StaticInt{M}() - One()) ÷ StaticInt{M}()
@inline function divrem_fast(x::I, y) where {I <: Integer}
    ux = unsigned(x); uy = unsigned(y)
    d = Base.udiv_int(ux, uy)
    r = ux - d * uy
    d % I, r % I
end
@inline divrem_fast(::StaticInt{x}, y::I) where {x, I <: Integer} = divrem_fast(x % I, y)
@inline div_fast(x::I, y::Integer) where {I <: Integer} = Base.udiv_int(unsigned(x), unsigned(y)) % I
@inline div_fast(::StaticInt{x}, y::I) where {x, I <: Integer} = Base.udiv_int(unsigned(x), unsigned(y)) % I
divrem_fast(::StaticInt{N}, ::StaticInt{M}) where {N,M}= divrem(StaticInt{N}(), StaticInt{M}())
div_fast(::StaticInt{N}, ::StaticInt{M}) where {N,M}= StaticInt{N}() ÷ StaticInt{M}()
@generated function div_fast(x::I, ::StaticInt{M}) where {I<:Integer,M}
    if VectorizationBase.ispow2(M)
        lm = VectorizationBase.intlog2(M)
        Expr(:block, Expr(:meta,:inline), :(x >>> $lm))
    else
        Expr(:block, Expr(:meta,:inline), :(div_fast(x, $(I(M)))))
    end
end


text:
################################################################################
# --------------- Constructor for counterfactual state:
################################################################################
struct CounterfactualState
    x::AbstractArray
    target_encoded::Union{Number, AbstractVector}
    x′::AbstractArray
    M::AbstractFittedModel
    params::Dict
    search::Union{Dict,Nothing}
end

################################################################################
# --------------- Base type for generator:
################################################################################
"""
    AbstractGenerator

An abstract type that serves as the base type for counterfactual generators.
"""
abstract type AbstractGenerator end

# Loss:
using Flux
"""
    ℓ(generator::AbstractGenerator, counterfactual_state::CounterfactualState)

The default method to apply the generator loss function to the current counterfactual state for any generator.
"""
function ℓ(generator::AbstractGenerator, counterfactual_state::CounterfactualState)

    output = :logits # currently counterfactual loss is always computed with respect to logits

    loss = getfield(Losses, generator.loss)(
        getfield(Models, output)(counterfactual_state.M, counterfactual_state.x′),
        counterfactual_state.target_encoded
    )

    return loss
end

# Complexity:
"""
    h(generator::AbstractGenerator, counterfactual_state::CounterfactualState)

The default method to apply the generator complexity penalty to the current counterfactual state for any generator.
"""
h(generator::AbstractGenerator, counterfactual_state::CounterfactualState) = generator.complexity(counterfactual_state.x-counterfactual_state.x′)

################################################################################
# Subtypes
################################################################################

include("gradient_based/functions.jl")


text:
using Reactive

# Stop the runner task

try
    println("Killing ", Reactive.runner_task) # the task switch caused here is required!
    Reactive.stop()
catch
end

include("basics.jl")
#include("gc.jl")
include("call_count.jl")
include("flatten.jl")
include("time.jl")
include("async.jl")
FactCheck.exitstatus()


text:
module DivConformingTests
using SparseMatricesCSR
import Gridap: ∇, divergence, DIV
using Gridap
using Gridap.Algebra
using Gridap.FESpaces
using GridapDistributed
using PartitionedArrays
using Test
using FillArrays

function setup_p1_model()
  ptr  = [ 1, 5, 9 ]
  data = [ 1,2,3,4, 2,5,4,6  ]
  cell_vertex_lids = Gridap.Arrays.Table(data,ptr)
  node_coordinates = Vector{Point{2,Float64}}(undef,6)

  node_coordinates[1]=Point{2,Float64}(0.0,0.0)
  node_coordinates[2]=Point{2,Float64}(0.5,0.0)
  node_coordinates[3]=Point{2,Float64}(0.0,1.0)
  node_coordinates[4]=Point{2,Float64}(0.5,1.0)
  node_coordinates[5]=Point{2,Float64}(1.0,0.0)
  node_coordinates[6]=Point{2,Float64}(1.0,1.0)

  polytope=QUAD
  scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1)
  cell_types=collect(Fill(1,length(cell_vertex_lids)))
  cell_reffes=[scalar_reffe]
  grid = Gridap.Geometry.UnstructuredGrid(node_coordinates,
                                          cell_vertex_lids,
                                          cell_reffes,
                                          cell_types,
                                          Gridap.Geometry.NonOriented())
  Gridap.Geometry.UnstructuredDiscreteModel(grid)
end

function setup_p2_model()
  ptr  = [ 1, 5, 9 ]
  data = [ 1,2,3,4, 5,1,6,3  ]
  cell_vertex_lids = Gridap.Arrays.Table(data,ptr)
  node_coordinates = Vector{Point{2,Float64}}(undef,6)

  node_coordinates[1]=Point{2,Float64}(0.5,0.0)
  node_coordinates[2]=Point{2,Float64}(1.0,0.0)
  node_coordinates[3]=Point{2,Float64}(0.5,1.0)
  node_coordinates[4]=Point{2,Float64}(1.0,1.0)
  node_coordinates[5]=Point{2,Float64}(0.0,0.0)
  node_coordinates[6]=Point{2,Float64}(0.0,1.0)

  polytope=QUAD
  scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1)
  cell_types=collect(Fill(1,length(cell_vertex_lids)))
  cell_reffes=[scalar_reffe]
  grid = Gridap.Geometry.UnstructuredGrid(node_coordinates,
                                          cell_vertex_lids,
                                          cell_reffes,
                                          cell_types,
                                          Gridap.Geometry.NonOriented())
  Gridap.Geometry.UnstructuredDiscreteModel(grid)
end

function f(model,reffe)
    V = FESpace(model,reffe,conformity=:Hdiv)
    U = TrialFESpace(V)

    das = FullyAssembledRows()
    trian = Triangulation(das,model)
    degree = 2
    dΩ = Measure(trian,degree)
    a(u,v) = ∫( u⋅v )*dΩ

    u  = get_trial_fe_basis(U)
    v  = get_fe_basis(V)
    dc = a(u,v)

    assem = SparseMatrixAssembler(U,V,das)
    data = collect_cell_matrix(U,V,a(u,v))
    A = assemble_matrix(assem,data)
    t1  = trian.trians.parts[1]
    t2  = trian.trians.parts[2]
    dc1 = dc.contribs.parts[1]
    dc2 = dc.contribs.parts[2]
    c1  = Gridap.CellData.get_contribution(dc1,t1)
    c2  = Gridap.CellData.get_contribution(dc2,t2)
    tol = 1.0e-12
    @test norm(c1[1]-c2[2]) < tol
    @test norm(c1[2]-c2[1]) < tol
end

function main(parts)
    @assert isa(parts,SequentialData)

    models=map_parts(parts) do part
      if (part==1)
        setup_p1_model()
      else
        setup_p2_model()
      end
    end

    noids,firstgid,hid_to_gid,hid_to_part=map_parts(parts) do part
      if (part==1)
        1,1,[2],Int32[2]
      else
        1,2,[1],Int32[1]
      end
    end

    gids=PRange(
      parts,
      2,
      noids,
      firstgid,
      hid_to_gid,
      hid_to_part)

    model = GridapDistributed.DistributedDiscreteModel(models,gids)

    reffe=ReferenceFE(raviart_thomas,Float64,0)
    f(model,reffe)
    reffe=ReferenceFE(QUAD, raviart_thomas, 0)
    f(model,reffe)
  end

end # module


text:
using SparseArrays

arr1 = spzeros(Int64, 5, 5)
arr2 = spzeros(Int64, 5, 1)

arr2[4] = 1
arr1[3,:] = arr2[:,1]
arr2[4] = 2

# display(arr1)
display(Matrix(arr1))


text:
module TriangulationsTests

using Gridap
using Test
import Gridap: ∇

# CartesianTriangulation

partition = (2,4)
trian = CartesianTriangulation(partition)

test_triangulation(trian)

@test string(trian) == "CartesianTriangulation object"

s = "CartesianTriangulation object:\n physdim: 2\n refdim: 2\n ncells: 8"
@test sprint(show,"text/plain",trian) == s

@test ncells(trian) == 8

coords = CellPoints(trian)

x4 = Point{2,Float64}[[0.0, -0.5], [1.0, -0.5], [0.0, 0.0], [1.0, 0.0]]

@test coords[4] ≈ x4

# CellQuadrature

quad = CellQuadrature(trian,degree=2)
q = coordinates(quad)
w = weights(quad)

# CellField

ufun(x::Point{2}) = 2*x[1]+x[2]

u = CellField(trian,ufun)
@test isa(u,CellField{2,Float64})

ufun_grad(x) = VectorValue(2.0,1.0)
∇(::typeof(ufun)) = ufun_grad

u = CellField(trian,ufun)
@test isa(u,CellField{2,Float64})

νfun(x::Point{2},u::Float64) = TensorValue(x[1], u*x[2], 0.0, u)

ν = CellField(trian,νfun,u)

@test isa(u,CellField{2,Float64})
@test isa(ν,CellField{2,TensorValue{2,Float64,4}})

#ν(u) = CellField(trian,νfun,u)
#w = ν(u) # julia nightly builds get stuck here

function μfun(x,u,v,i)
  if i == 4
    return TensorValue(x[1], u*x[2], 0.0, u)
  else
    return TensorValue(0.1, x[2], 0.0, u)
  end
end

i = zeros(Int,length(u))
i[1:2] .= 4

cf = CellField(trian,μfun,u,u,i)

_ = collect(evaluate(cf,q))

cf = CellBasis(trian,μfun,u,u,i)

b = CellBasis(trian)

cb = CellBasis(trian,μfun,b,u,i)

_ = collect(evaluate(cb,q))

end # module


text:
@symbol_func function W_M(cur_reactor::AbstractReactor)
  cur_energy = K_WM(cur_reactor)

  cur_energy *= cur_reactor.R_0 ^ 3

  cur_energy *= cur_reactor.B_0 ^ 2

  cur_energy
end


text:
module MaxHelpingHandFlagExitBlock

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/FlagExitBlock" FlagExitBlock(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight,
    tileType::String="3", flag::String="flag_exit_block", inverted::Bool=false, playSound::Bool=true, instant::Bool=false)

const placements = Ahorn.PlacementDict(
    "Flag Exit Block (max480's Helping Hand)" => Ahorn.EntityPlacement(
        FlagExitBlock,
        "rectangle",
        Dict{String, Any}(),
        Ahorn.tileEntityFinalizer
    ),
)

Ahorn.materialTileTypeKey(entity::FlagExitBlock) = "tileType"

Ahorn.editingOptions(entity::FlagExitBlock) = Dict{String, Any}(
    "tileType" => Ahorn.tiletypeEditingOptions()
)

Ahorn.minimumSize(entity::FlagExitBlock) = 8, 8
Ahorn.resizable(entity::FlagExitBlock) = true, true

Ahorn.selection(entity::FlagExitBlock) = Ahorn.getEntityRectangle(entity)

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::FlagExitBlock, room::Maple.Room) = Ahorn.drawTileEntity(ctx, room, entity, alpha=0.7)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
using HorizonSideRobots

function obstacle_cross_marking(r::Robot)
    for side in (Nord, West, Sud, Ost)
        n = direction_marking_return_steps_count(r, side) # Записывает пройденный путь, маркируя при этом клетки
        for i in 1:n # Вернуться по n шагов через препятствия
            move_through_obstacle(r, HorizonSide(mod(Int(side) + 2, 4)))
        end
    end
    putmarker!(r)
end

function direction_marking_return_steps_count(r::Robot, side::HorizonSide)
    n = 0
    while move_through_obstacle(r, side) == true
        putmarker!(r)
        n += 1
    end
    return n
end

function move_through_obstacle(r::Robot, forward_side::HorizonSide)::Bool
    side_side = HorizonSide(mod(Int(forward_side) - 1, 4)) # Направление в сторону обхода
    otherside_side = HorizonSide(mod(Int(side_side) + 2, 4)) # Обратное направление в сторону
    n = 0
    if isborder(r,forward_side)==false
        move!(r,forward_side)
        not_hit_border=true
    else
        while isborder(r,forward_side) == true
            if isborder(r, side_side) == false
                move!(r, side_side)
                n += 1
            else
                break
            end
        end
        if isborder(r,forward_side) == false
            move!(r,forward_side)
            while isborder(r, otherside_side) == true
                move!(r,forward_side)
            end
            not_hit_border = true
        else
            not_hit_border = false
        end
        while n > 0
            n -= 1
            move!(r, otherside_side)
        end
    end
    return not_hit_border
end