text:
include("../RLAlgorithms.jl")
using .RLAlgorithms
using Plots
using Flux
using Flux: Optimise.update!
using Statistics: mean

function act(actor, s, noisescale)
    a = actor(s)
    return clamp.(a .+ noisescale .* randn(Float32, size(a)), -1f0, 1f0)
end

function run!(actor, aopt, critic, copt, env;
        epochs=100, steps=300, maxt=100, batchsize=128, noisescale=1f-1, γ=99f-2, τ=1f-3)
    rewards = zeros(Float32, epochs)
    actortar = deepcopy(actor)
    critictar = deepcopy(critic)
    memory = ReplayMemory{Float32, Float32, Float32}(
        length(env.observationspace), length(env.actionspace), 10*steps, batchsize
    )
    actorps = params(actor)
    criticps = params(critic)

    for i in 1:epochs
        episodes!(memory, env, steps, maxt) do s
            act(actor, s, noisescale)
        end

        for _ in 1:steps
            s1, a1, r, s2, done = sample(memory)

            atar2 = actortar(s2)
            Qtars2 = vec(critictar(vcat(s2, atar2)))
            Qtar = r .+ (1 .- done) .* γ .* Qtars2
            cgs = gradient(criticps) do
                Q = vec(critic(vcat(s1, a1)))
                return Flux.mse(Q, Qtar)
            end
            update!(copt, criticps, cgs)

            ags = gradient(actorps) do
                as = actor(s1)
                return -mean(critic(vcat(s1, as)))
            end
            update!(aopt, actorps, ags)

            softupdate!(actortar, actor, τ)
            softupdate!(critictar, critic, τ)
        end
        rewards[i] = test(s -> actor(s), env)
    end
    return rewards
end

env = CartPoleContinuous{Float32}()
aDim = length(env.actionspace)
oDim = length(env.observationspace)

actor = Chain(
    Dense(oDim, 32, relu),
    Dense(32, aDim, tanh)
)
critic = Chain(
    Dense(oDim+aDim, 32, relu),
    Dense(32, 16, relu),
    Dense(16, 1)
)

r = run!(actor, ADAM(), critic, ADAM(), env)
plt = plot(r, labels="Reward")
display(plt)
#savefig(plt, "ddpg")


text:
"""

    simulatedata(; N=50, SN=100.0, Temperature=5000.0, offset=5.0, scale=11.0, seed=1)

Simulates observations from total flux, i.e. AGN + Galaxy. The observations are simulated at the LSST wavelengths. Briefly, the data are simulated as follows:

- The AGN is simulated using the `lsststd` function.

- The galaxy is simulated as follows: given temperature, we use a planck model to simulate emission at LSST wavelengths. These emissions form a vector which is then normalised and then scaled by the input argument `scale`.

The function returns two types of observations: **co-occuring** (or simultaneous) in time and **non co-occuring**.

When called, it will produce the following plots:

- One plot ...
- The other plot ...

The function is useful for verifiying the code and is not meant to necessarily present a physically realistic scenario.

# Arguments

- `N`           : integer number of observations to be simulated per filter
- `SN`          : Signal to noise
- `Temperature` : temperature of the galaxy simulated as a black body between 1000 and 12000 Kelvin
- `offset`      : positive offset from minimum observed AGN activity
- `scale`       : positive scalar mulitplied with the unit galaxy vector (see above)
- `seed`        : controls the random number generator, useful for reproducing data

# Returns

- `obstime`   : vector of observation times
- `obsflux`   : vector of non co-occuring total flux observations
- `σobs`      : observed errors of measurements
- `filteridx` : vector of integer entries that express filter membership
- `obsfluxco` : matrix of dimensions N×(number of filters) that holds co-occuring total flux observations
- `σobsco`    : matrix of dimensions N×(number of filters) that holds errors of co-occuring observations
- `lineeq`    : true total flux line
- `obstimeco` : times of co-occuring flux observations

# Example

```julia-repl
obstime, obsflux, σobs, filteridx, obsfluxco, σobsco = simulatedata(N=50, σ=0.1, Temperature=3000, offset=15.5, scale=10.0)
```
"""
function simulatedata(; N=50, SN=50.0, logs=-8.0, Temperature = 5000.0, offset=2.0, scale=20.0, seed=1, cooccur=false)

    @assert(scale > 0.0)
    @assert(Temperature > 1000.0)
    @assert(Temperature < 15_000.0)
    @assert(offset >= 0.0)
    @assert(SN > 0.0)

    rg = MersenneTwister(seed)

    numwaves = length(lsstwaves())

    # specify parameters for data generation

    a = [lsststd(zpwave = λ, z = 0.5, lum = 5e43) for λ in lsstwaves()].*5
    b = unitgalaxyvector(Temperature) * scale
    θ = [0.0; logs]

    # These are the filter memberships of the data

    filteridx = kron(collect(1:numwaves), ones(Int, N))

    # Initialise arrays to store data

    obstime   = Array{Array{Float64,1},1}(undef, numwaves)
    agnflux   = Array{Array{Float64,1},1}(undef, numwaves)
    obsflux   = Array{Array{Float64,1},1}(undef, numwaves)
    σobs      = Array{Array{Float64,1},1}(undef, numwaves)

    agnfluxco = zeros(N, numwaves)
    obsfluxco = zeros(N, numwaves)
    σobsco    = zeros(N, numwaves)

    # draw latent function from Gaussian process

    T         = 2000
    timegrid  = collect(LinRange(0.0, 1000, T))
    K         = calculatekernelmatrix(timegrid',timegrid', rbf, θ)
    z         = rand(rg, MvNormal(zeros(T), K))

    @printf("Statistics of latent signal\n")
    @show mean(z) std(z)

    sametime = sort(randperm(rg, length(timegrid))[1:N])

    for f = 1:numwaves

        # Non co-occuring observations

        timeidx    = cooccur ? sametime : sort(randperm(rg, length(timegrid))[1:N])

        obstime[f] = timegrid[timeidx] .+ offset

        agnflux[f] = a[f]*(z[timeidx] .+ offset)

        flux       = a[f]*(z[timeidx] .+ offset) .+ b[f]

        σobs[f]    = flux / SN

        obsflux[f] = flux .+ σobs[f] .* randn(rg, N)

        # Co-occuring observations

        agnfluxco[:, f] = a[f]*(z[sametime] .+ offset)

        fluxco          = a[f]*(z[sametime] .+ offset) .+ b[f]

        σobsco[:, f]    = fluxco / SN

        obsfluxco[:, f] = fluxco .+ σobsco[:, f].*randn(rg, N)

    end

    figure(1)
    cla()
    title("latent signal simulated as gp draw")
    plot(timegrid, z, "k")

    figure(2)
    cla()
    title("plotting irregular total flux time series")
    for f = 1:numwaves
        PyPlot.plot(obstime[f], obsflux[f], "o", label=@sprintf("%.2f", lsstwaves()[f]))
    end

    PyPlot.legend()

    figure(3)
    cla()
    title("plotting regular total flux time series")
    for f = 1:numwaves
        PyPlot.plot(timegrid[sametime].+offset, obsfluxco[:,f], "o", label=@sprintf("%.2f", lsstwaves()[f]))
    end

    PyPlot.legend()

    figure(4)
    cla()
    title("plotting first three filters")
    PyPlot.plot3D(obsfluxco[:, 1], obsfluxco[:, 2], obsfluxco[:, 3], "bo", alpha=0.2, label="obsmix")
    PyPlot.plot3D(agnfluxco[:, 1], agnfluxco[:, 2], agnfluxco[:, 3], "k.", alpha=1.0, label="agn")
    PyPlot.plot3D([0;b[1]], [0;b[2]], [0;b[3]], "c-", alpha=1.0, label="galaxy")

    PyPlot.legend()

    lineeq(t) = a*(t + offset) + b

    obstimeco = timegrid[sametime] .+ offset

    return reduce(vcat, obstime), reduce(vcat, obsflux), reduce(vcat, σobs), filteridx, obsfluxco, σobsco, lineeq, obstimeco

end

ext: jl


text:
using RecipesBase

struct CartPoleParams{T} <:Any where {T<:Real}
    m::T # mass of pole
    l::T # length of pole
    mc::T # mass of cart
    muc::T # some constant?
    mup::T # some constant?
    fmag::T # magnitude of force applied
	g::T  # gravity (not signed)
    CartPoleParams() = new{Float64}(0.1, 0.5, 1., 0.0005, 0.000002, 10., 9.8)
    CartPoleParams(T::Type) = new{T}(T(0.1), T(0.5), T(1.), T(0.0005), T(0.000002), T(10.), T(9.8))
	CartPoleParams(T::Type, m, l, mc, g) = new{T}(m, l, mc, T(0.0005), T(0.000002), T(10.), g)
end

mutable struct CartPole{T,TA} <: AbstractEnvironment where {T<:Real,TA<:Real}
    state::Array{T,1}
    reward::T
    done::Bool
    params::CartPoleParams{T}
    dt::T
    simSteps::Int
    t::T # time

    function CartPole(T::Type=Float64, TA::Type=Int)
        new{T,TA}(zeros(T,4), T(0.0), false,
           CartPoleParams(T), T(0.02), 1, T(0.)
       )
    end

	function CartPole(::Type{T}, ::Type{TA}, params::CartPoleParams{T}) where {T,TA}
		new{T,TA}(zeros(T,4), T(0.0), false,
		   params, T(0.02), 1, T(0.)
	   )
	end

end

function get_gamma(env::CartPole{T})::T where {T}
	return T(1.0)
end

function randomize_env(::Type{CartPole{T,TA}}, rng::AbstractRNG) where {T, TA}
	m = rand(rng, Uniform(0.025, .25))  # mass of pole
	l = rand(rng, Uniform(0.1, 1.0))    # length of pole
	mc = rand(rng, Uniform(0.1, 5.))    # mass of cart
	g = rand(rng, Uniform(8., 10.))     # gravity
	params = CartPoleParams(T, m, l, mc, g)
	return CartPole(T, TA, params)
end

function get_state_desc(env::CartPole{T})::Array{T,2} where {T<:Real}
    ranges = zeros(T, (4,2))
    ranges[1,:] .= [-2.4, 2.4]       # x range
	ranges[2,:] .= [-10., 10.]#[-6, 6.]       # xDot range
	ranges[3,:] .= [-π / 12.0, π / 12.0] # theta range
	ranges[4,:] .= [-π, π]           # thetaDot range
	# ranges[5,:] .= [0., 20.]
    return ranges
end

function get_action_desc(env::CartPole{T,TA})::Int where{T<:Real,TA<:Integer}
    return 2
end

function get_action_desc(env::CartPole{T,TA}) where{T<:Real,TA<:AbstractFloat}
	ranges = zeros(TA, (1,2))
	ranges[1, 1] = -1.
	ranges[1, 2] = 1.
	return ranges
end

function cartpole_sim(state::Array{T,1}, constants::CartPoleParams{T}, u::T, steps::Int, dt::T) where {T<:Real}
	x, xDot, theta, thetaDot = view(state, 1:4)
	omegaDot::T = 0.
	vDot::T = 0.
	m = constants.m     # mass of pole
    l = constants.l     # length of pole
    mc = constants.mc   # mass of cart
    muc = constants.muc # some constant?
    mup = constants.mup # some constant?
	g = constants.g
	# for i in 1:steps
	omegaDot = (g * sin(theta) + cos(theta) * (muc * sign(xDot) - u - m * l * thetaDot^2 * sin(theta)) / (m + mc) - mup * thetaDot / (m * l)) / (l * (4.0 / 3.0 - m / (m + mc) * cos(theta)^2))
    vDot = (u + m * l * (thetaDot^2 * sin(theta) - omegaDot*cos(theta)) - muc*sign(xDot)) / (m + mc)
    theta += dt * thetaDot
    thetaDot += dt * omegaDot
    x += dt * xDot
    xDot += dt * vDot

    theta = mod(theta + π, 2 * π) - π
	# end
	state[1:4] .= x, xDot, theta, thetaDot
end

function step!(env::CartPole{T}, action::Int, rng::AbstractRNG)::T where {T<:Real}
    if action <= 0 || action > 2
        error("Action needs to be an integer in [1, 2]")
    end
    u::T = 0.
    if action == 1
        u = -env.params.fmag
    else
        u =  env.params.fmag
    end

	cartpole_sim(env.state, env.params, u, env.simSteps, env.dt/env.simSteps)
	env.t += env.dt
	# env.state[5] = env.t
    # env.state[1] = clamp(env.state[1], -2.4, 2.4)

    env.reward = 1.0
	env.done = terminalCheck(env)
	# if done
	# 	println(env.t, env.dt)
	# end

    return env.reward
end

function terminalCheck(env::CartPole)::Bool
	polecond = abs(env.state[3]) > (π / 15.0)
	cartcond = abs(env.state[1]) ≥ 2.4
	timecond = env.t ≥ (20. - 1e-8)
	done = polecond | cartcond | timecond
	# if done
	# 	println([polecond, cartcond, timecond])
	# end
	return done
	# return (abs(env.state[3]) > (π / 15.0)) | (abs(env.state[1]) ≥ 2.4) | (env.t ≥ (20. - 1e-8))
end

function step!(env::CartPole{T}, action::T, rng::AbstractRNG)::T where {T<:AbstractFloat}
	u = clamp(action, -1., 1.)
	u *= env.params.fmag

	cartpole_sim(env.state, env.params, u, env.simSteps, env.dt/env.simSteps)
	env.t += env.dt
    # env.state[1] = clamp(env.state[1], -2.4, 2.4)
	# env.state[5] = env.t
    env.reward = 1.0
	env.done = terminalCheck(env)

    return env.reward
end

function step!(env::CartPole{T}, action::Array{T}, rng::AbstractRNG)::T where {T<:AbstractFloat}
	u = clamp(action[1], -1., 1.)
	u *= env.params.fmag

	cartpole_sim(env.state, env.params, u, env.simSteps, env.dt/env.simSteps)
	env.t += env.dt
    # env.state[1] = clamp(env.state[1], -2.4, 2.4)
	# env.state[5] = env.t
    env.reward = 1.0
	env.done = terminalCheck(env)

    return env.reward
end

function get_state(env::CartPole{T})::Array{T,1} where {T<:Real}
    return env.state # VectorState([env.x, env.xDot, env.theta, env.thetaDot])
end

function get_reward(env::CartPole)
    return env.reward
end

function is_terminal(env::CartPole)::Bool
    return env.done
end

function env_reset!(env::CartPole, rng)
	fill!(env.state, 0.)
    env.t = 0.
	env.done = false
end

function clone(env::CartPole{T,TA}) where {T<:Real,TA<:Union{Int,<:AbstractFloat}}
	env2 = CartPole(T,TA)
	env2.state .= deepcopy(env.state)
	env2.reward = deepcopy(env.reward)
	env2.done = deepcopy(env.done)
	env2.t = deepcopy(env.t)
	env2.dt = deepcopy(env.dt)
	env2.simSteps = deepcopy(env.simSteps)
	env2.params = deepcopy(env.params)
	return env2
end

# CartPole Rendering originally from JuliaML/Reinforce.jl and has only be modified slighlty
# https://github.com/JuliaML/Reinforce.jl/blob/master/src/envs/cartpole.jl
@recipe function f(env::CartPole)
	x, xvel, θ, θvel = env.state[1:4]
	legend := false
	xlims := (-2.4, 2.4)
	l = env.params.l
	ylims := (-Inf, 2l)
	grid := false
	ticks := nothing

	# pole
	@series begin
		linecolor := :red
		linewidth := 10
		[x, x + 2l * sin(θ)], [0.0, 2l * cos(θ)]
	end

	# cart
	@series begin
		seriescolor := :black
		seriestype := :shape
		hw = 0.5
		l, r = x-hw, x+hw
		t, b = 0.0, -0.1
		[l, r, r, l], [t, t, b, b]
	end
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
using Test
using LinearAlgebra
import sbp
using sbp.Config: Tv
using sbp.Metrics: build_jacobian, build_contravariant_basis,
                   build_covariant_metric_tensor,
                   build_contravariant_metric_tensor
using sbp.Staggered: build_all_operators_2d, build_covariant_basis_mm
using sbp.Grid: grid_xp, grid_xm, grid_2d_x, grid_2d_y

sop = sbp.OP2019
nx=10
ny=20

desc = sop.load_description()
builder = n -> sop.build_operators(desc, n)
pp, mm, pm, mp = build_all_operators_2d(builder, nx, ny)

x1, h = grid_xp(nx)
y1, h = grid_xp(ny)

x = grid_2d_x(x1, ny)
y = grid_2d_y(y1, nx)

fx = x
fy = y

a = build_covariant_basis_mm(fx, fy, mm, pm)
Ga = build_covariant_metric_tensor(a)

@testset "Jacobian" begin
        Jans = ones(mm.nx * mm.ny)
        J = build_jacobian(a)
        @test isapprox(J, Jans)
end

@testset "Covariant Metric tensor" begin
         Ga = build_covariant_metric_tensor(a)
         # Check symmetry
         @test isapprox(Ga.a12, Ga.a21)

end

@testset "Contravariant basis" begin
        J = build_jacobian(a)
        b = build_contravariant_basis(J, a)

        one = ones(mm.nx * mm.ny)
        zero = zeros(mm.nx * mm.ny)

        @test isapprox(a.x_r1 .* b.r1_x + a.y_r1 .* b.r1_y, one)
        @test isapprox(a.x_r1 .* b.r2_x + a.y_r1 .* b.r2_y, zero, atol=1.0)
        @test isapprox(a.x_r2 .* b.r2_x + a.y_r2 .* b.r2_y, one)
        @test isapprox(a.x_r2 .* b.r1_x + a.y_r2 .* b.r1_y, zero, atol=1.0)

        @testset "Contravariant Metric tensor" begin
                Gb = build_contravariant_metric_tensor(b)
                # Check symmetry
                @test isapprox(Gb.b12, Gb.b21)
                # Check orthogonality
                @test isapprox(Ga.a11 .* Gb.b11 + Ga.a12 .* Gb.b12, one)
                @test isapprox(Ga.a11 .* Gb.b12 + Ga.a12 .* Gb.b22, zero,
                               atol=1.0)
                @test isapprox(Ga.a12 .* Gb.b11 + Ga.a22 .* Gb.b12, zero,
                               atol=1.0)
                @test isapprox(Ga.a21 .* Gb.b12 + Ga.a22 .* Gb.b22, one)

        end

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
### Empirical Monte Carlo tuner

immutable EmpiricalMCTuner <: MCTuner
  targetrate::Float64 # Target acceptance rate
  period::Int # Adaptation period, which determines when to attempt tuning
  adaptnsteps::Bool # Determine whether to adapt the number of steps along with stepsize
  maxnsteps::Int # Maximum number of steps allowed (for ex maximum number of leapfrog steps)
  targetlen::Float64 # Target length (this can be thought of as adaptation step size times number of steps)
  verbose::Bool # Specify whether the tuner will be in verbose or silent mode

  function EmpiricalMCTuner(targetrate::Float64, period::Int, adaptnsteps::Bool, maxnsteps::Int, targetlen::Float64,
    verbose::Bool)
    @assert 0 < targetrate < 1 "Target acceptance rate should be between 0 and 1."
    @assert period > 0 "Adaptation period should be positive."
    @assert maxnsteps > 0 "Maximum number of steps should be positive."
    @assert targetlen > 0 "Target length should be positive."
    new(targetrate, period, adaptnsteps, maxnsteps, targetlen, verbose)
  end
end

EmpiricalMCTuner(targetrate::Float64; period::Int=100, adaptnsteps::Bool=true, maxnsteps::Int=200,
  targetlen::Float64=1., verbose::Bool=false) =
  EmpiricalMCTuner(targetrate, period, adaptnsteps, maxnsteps, targetlen, verbose)

### EmpiricalMCTune holds the temporary output of a Monte Carlo sampler that uses the EmpiricalMCMCTuner

type EmpiricalMCTune <: MCTune
  step::Float64 # Stepsize of current Monte Carlo iteration (for ex leapfrog or drift stepsize)
  nsteps::Int # Number of steps of current Monte Carlo iteration (for ex number of leapfrog steps)
  accepted::Int # Number of accepted Monte Carlo samples within current tuner period
  proposed::Int # Number of proposed Monte Carlo samples within current tuner period
  rate::Float64 # Acceptance rate over current tuner period

  function EmpiricalMCTune(step::Float64, nsteps::Int, accepted::Int, proposed::Int, rate::Float64)
    @assert step > 0 "Stepsize of current MCMC iteration should be positive."
    @assert nsteps > 0 "Number of steps of current MCMC iteration should be positive."
    @assert accepted >= 0 "Number of accepted MCMC samples should be non-negative."
    @assert proposed >= 0 "Number of proposed MCMC samples should be non-negative."
    new(step, nsteps, accepted, proposed, rate)
  end
end

EmpiricalMCTune(step::Float64, nsteps::Int, accepted::Int, proposed::Int) =
  EmpiricalMCTune(step::Float64, nsteps::Int, accepted::Int, proposed::Int, NaN)
EmpiricalMCTune(step::Float64, nsteps::Int) = EmpiricalMCTune(step::Float64, nsteps::Int, 0, 0, NaN)
EmpiricalMCTune(step::Float64) = EmpiricalMCTune(step::Float64, 10, 0, 0, NaN)

reset!(tune::EmpiricalMCTune) = ((tune.accepted, tune.proposed) = (0, 0))

count!(tune::EmpiricalMCTune) = (tune.accepted += 1)

function adapt!(tune::EmpiricalMCTune, tuner::EmpiricalMCTuner)
  tune.rate = tune.accepted/tune.proposed
  tune.step *= (1/(1+exp(-11*(tune.rate-tuner.targetrate)))+0.5)
  if tuner.adaptnsteps
    tune.nsteps = min(tuner.maxnsteps, ceil(tuner.targetlen/tune.step))
  end
end


text:
# This file was generated, do not modify it. # hide
mach = machine(model, X_iris, y_iris)

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
abstract AbstractFiniteElement

abstract AbstractLagrangeFE

#type LagrangeFE{T <: AbstractReferenceElement}
#    ref_shape::T
#    basis::PolyBasis
#end

#get_coeff(LFE::LagrangeFE, i::Int) = LFE.basis[:, i]
#get_coeffs(LFE::LagrangeFE, i::Int) = LFE.basis[:, i]

typealias NodePoly MPoly{Float64}

type LagrangeFE
    shape_functions::Vector{NodePoly}
    order::Int
    points::Vector{Vector{Float64}}
end

function LagrangeFE(shape::AbstractReferenceElement, order::Int)
    points = make_lattice(shape, LagrangeFE, order)

   db = generate_dual_basis(shape, LagrangeFE, order)

    V = zeros(length(points), length(points))
    for (i, p) in enumerate(points)
        V[:, i] = evaluate_basis(db, points[i]...)
    end

    C = inv(V)

    shape_functions = NodePoly[]
    for i in 1:size(C, 1)
        p1 = copy(db)
        coeffs = C[:, i]
        i = 1
        for (m, c) in db
            p1[m] = coeffs[i]
            i += 1
        end
        push!(shape_functions, p1)
    end

    LagrangeFE(shape_functions, order, points)

end

function tabulate(ele::LagrangeFE, pts::Vector{Vector{Float64}})
    Ns = Array(Vector{Float64}, 0)

    for point in pts
        N = tabulate(ele, point)
        push!(Ns, N)
    end
    return Ns
end
function tabulate(ele::LagrangeFE, point::Vector{Float64})
    N = Float64[]
    for sf in ele.shape_functions
        push!(N, evaluate(sf, point...))
    end
    return N
end

get_points_row(::Union(ReferenceTriangle, ReferenceTetrahedron), order::Int, col::Int) = order - col
get_points_row(::ReferenceSquare, order::Int, ::Int) = order

function generate_dual_basis(shape::AbstractLagrangeFE)
    collect(0:order)
end

function generate_dual_basis(shape::ReferenceSquare, ::Type{LagrangeFE}, order::Int)
    # Lagrange elements in 2D have this type from the Pascals triangle
    #=        1
            x   y
          .  .  .  .
     x^m .  .  .  . y^m
        x^m*y    y^m*x
           .     .
             (xy)^m
    =#

    # This is an unreadable mess but it basically
    # calculates the exponents for x and y row by row
    # left to right...

    x, y = generators(MPoly{Float64}, :x, :y)

    P = zero(x)

    # Start exponent for x and y
    x_m = [0:order-1; (order-1) * ones(Int, order-1)]
    y_m = [zeros(Int, order-1); 0:order-1]

    # Length of row
    l = [1:order;order-1:-1:1]

    for i = 1:length(l), j in 0:l[i]-1
            P += x^(x_m[i] - j) * y^(y_m[i] + j)
    end

    return P

end

# These could be made much nicer and more general:

make_lattice(::ReferenceLine, ::Type{LagrangeFE}, order::Int) = linspace(-1,1, order)

function make_lattice(refele::ReferenceTriangle, ::Type{LagrangeFE}, order::Int)
    vs = vertices(refele)
    intervals = [ (vs[i] - vs[1]) / order for i in 2:length(vs)]

    edges = edges(refele)
    for edge in edges
        push!(points, vs[edge])
    end

    points = Array(Vector{Float64},0)
    for i in 0:order, j in 0:get_points_row(refele, order, i)
            push!(points, intervals[:,1] * j + i * intervals[:,2])
    end
    return points
end

function make_lattice(refele::ReferenceSquare, ::Type{LagrangeFE}, order::Int)
    order -= 1
    vs = vertices(refele)
    intervals = [(vs[2] - vs[1])/order  (vs[4] - vs[1])/order]
    x0 = [vs[1][1], vs[1][2]]

    points = Array(Vector{Float64},0)
    for i in 0:order, j in 0:order
            push!(points, x0 + intervals[:,1] * j + i * intervals[:,2])
    end
    return points
end

function make_lattice(refele::ReferenceTetrahedron, ::Type{LagrangeFE}, order::Int)
    vs = vertices(refele)
    intervals = [ (vs[i] - vs[1]) / order for i in 2:length(vs)]
    points = Array(Vector{Float64},0)
    for m in 0:order
        for i in 0:get_points_row(refele, order, m)
            for j in 0:get_points_row(refele, get_points_row(refele, order, m), i)
                push!(points, intervals[1] * j + i * intervals[2] + m * intervals[3])
            end
        end
    end

    return points
end

function make_lattice(refele::ReferenceCube, ::Type{LagrangeFE}, order::Int)
    order -= 1
    vs = vertices(refele)
    intervals = [(vs[2] - vs[1])/order  (vs[4] - vs[1])/order (vs[5] - vs[1])/order]
    x0 = [vs[1][1], vs[1][2], vs[1][3]]
    points = Array(Vector{Float64},0)
    for m in 0:order, i = 0:order, j=0:order
        push!(points, x0 + intervals[1] * j + i * intervals[2] + m * intervals[3])
    end
    return points
end

max_stars_repo_name: JuliaTagBot/JFET

max_issues_repo_name: JuliaTagBot/JFET

max_forks_repo_name: JuliaTagBot/JFET


text:
using DelimitedFiles
using StatsFuns
using TFIMSampler
using LinearAlgebra
using ProgressMeter

function clean(L; name="train", r=1, R=50)
    h = 1.0
    h = round(h, digits=2)
    dir = "/scratch/mbeach/TFIM_samples/PBC/h-1.0/L-$L/"
    dir2 = "/scratch/mbeach/new_new_cleaned_TFIM_samples/L-$L/$name/"
    Base.Filesystem.rm(dir2, force=true, recursive=true)
    mkpath(dir2)

    @info "L is $L"
    @info "h is $h"

    M = pairing(L, h)
    logZ = TFIMSampler.logZ(M)
    @info "log Z is $logZ"
    @info "Z is $(exp(logZ))"

    file = dir * "run-$r.txt"
    train = readdlm(file)
    @info "Original training set size : $(size(train))"
    @showprogress for i in r+1:R
        file = dir * "run-$i.txt"
        train = vcat(train, readdlm(file))
    end

    @info "Training set size : $(size(train))"
    un = unique(train, dims=1)
    @info "unique configurations size : $(size(un))"

    un_ps = [TFIMSampler.get_prob(un[i,:]|> BitArray, M) for i in 1:size(un,1)]
    ps = [TFIMSampler.get_prob(train[i,:]|> BitArray, M) for i in 1:size(train,1)]

    @info "logZ is approx. $logZ"
    ps = [ps[i] - logZ for i in 1:size(train,1)]
    un_ps = [un_ps[i] - logZ for i in 1:size(un,1)]
    @info "Sum of ps is $(sum(exp.(un_ps)))"

    writedlm(dir2 * "logZ_$name.txt", logZ)
    writedlm(dir2 * "configs_$name.txt", Int.(train))
    writedlm(dir2 * "unique_configs_$name.txt", Int.(un))
    writedlm(dir2 * "logps_$name.txt", ps)
    writedlm(dir2 * "unique_logps_$name.txt", un_ps)
end

clean(4, name="train", r=1, R=50)
clean(4, name="test", r=51, R=100)

# for L in [8, 16, 32, 64, 128, 256]
# for L in [8, 16, 32, 64, 128, 256, 512]
for L in [512]
    @info L
    clean(L, name="train", r=1, R=50)
    clean(L, name="test", r=51, R=100)
end


text:
"""

# ancestral_graph

$(SIGNATURES)

Ancestral graphs after marginalization and conditioning.

### Required arguments
```julia
* `amat::NamedArray{Int, 2}`           : Adjacency matrix of a DAG
```

### Optional arguments
```julia
* `m::Vector{Symbol}`                  : Nodes in DAG that are marginalized
* `c::Vector{Symbol})`                 : Nodes in DAG there are conditioned on
```

### Returns
```julia
* `ag::NamedArray`                     : Ancestral graph remaining
```

# Extended help

### Example

### Adjacency matrix used for testing in ggm

```julia
amat_data = transpose(reshape([
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,
  0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,
  0,0,0,0,1,0,1,0,1,1,0,0,0,0,0,0,
  1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,
  0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0
], (16,16)));

vars = [Symbol("n\\\$i") for i in 1:size(amat_data, 1)]
amat = NamedArray(Int.(amat_data), (vars, vars), ("Rows", "Cols"));
m = [:n3, :n5, :n6, :n15, :n16];
c = [:n4, :n7];

ag = ancestral_graph(amat; m = m, c = c)
```

### Acknowledgements

Original author:                       Kayvan Sadeghi

Translated to Julia:                   Rob J Goedman

### References

Sadeghi, K. (2011). Stable classes of graphs containing directed acyclic
graphs.

Richardson, T.S. and Spirtes, P. (2002).  Ancestral graph Markov
models {Annals of Statistics}, 30(4), 962-1030.

### Licence

The R package ggm is licensed under License: GPL-2.

The Julia translation is licenced under: MIT.

Part of the api, exported.
"""
function ancestral_graph(amat::NamedArray{Int, 2}; m=Symbol[], c=Symbol[])

  a = amat'
  vars = names(a, 1)
  s = update_s(a, c)

  ar = copy(a)
  at = 2a

  while true

    at = copy(ar)
    ar = update_21(ar, m)
    ar .+= update_22(ar, s)
    ar .+= update_23(ar, m)
    ar .+= update_24(ar, m)
    ar .+= update_25(ar, m)
    ar .+= update_26(ar, m)
    ar .+= update_27(ar, s)
    ar .+= update_28(ar, s)
    ar .+= update_29(ar, m)
    ar .+= update_30(ar, m)

    ar == at && break

  end

  ar = update_a(ar)
  ar = update_b(ar, s)
  an = update_an(ar)

  while true
    at = copy(ar)

    a27n = update_27n(ar, an)
    an .+= a27n
    ar .+= a27n
    a22n = update_22n(ar, an)
    an .+= a22n
    ar .+= a22n

    at == ar && break
  end

  #ar_old = copy(ar)
  ar = update_i(ar)
  select = setdiff(vars, vcat(c, m))
  fr = ar[select, select]

end

"""

# ancestral_graph

$(SIGNATURES)

Ancestral graphs after marginalization and conditioning.

### Required arguments
```julia
* `d::DAG`                             : DAG onject
```

### Optional arguments
```julia
* `m::Vector{Symbol}`                  : Nodes in DAG that are marginalized
* `c::Vector{Symbol})`                 : Nodes in DAG there are conditioned on
```

### Returns
```julia
* `ag::NamedArray`                     : Ancestral graph remaining
```
"""
function ancestral_graph(d::DAG; m=Symbol[], c=Symbol[])
  ancestral_graph(d.a; m=m, c=c)
end

export
  ancestral_graph


text:
"""
    module Curls
"""
module Curls

export Curl

files = (
    "1D/Curls1D.jl",
    # "2D/Curls2D.jl",
    # "3D/Curls3D.jl"
    )

using LinearAlgebra
using SparseArrays
using ..Lattices

import ..PRINTED_COLOR_DARK
import ..PRINTED_COLOR_VARIABLE
import ..PRINTED_COLOR_INSTRUCTION

struct Curl{N} curl::SparseMatrixCSC{Float64,Int} end

foreach(include,files)

# Pretty Printing
function Base.show(io::IO,cc::Curl{N}) where N
    print(io,"$(N)D ")
    printstyled(io,"Curl",color=PRINTED_COLOR_DARK)
    printstyled(io,"(call with no args or ",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ky",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,", ",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"kz",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,")",color=PRINTED_COLOR_INSTRUCTION)
end

end


text:
using PlantBiophysics
using DataFrames
using CSV
using Statistics
using Plots
using MonteCarloMeasurements
using Makie
using CairoMakie
using Dates
unsafe_comparisons(true)
constants = Constants()
CairoMakie.activate!()

# Defining data path and choosing day and tree displayed
datapath="/Users/simon/Code/PlantBiophysics-paper/data/"
date="14/11/2001"
tree=3

# Reading file
file = datapath*"/TumbarumbaGasex_Spot_Medlyn.csv"
df = read_licor6400(file)
day = df[(df.Date.==date)&(day.Tree.==tree),:]

# Fitting paramaters (photosynthesis and stomatal conductance)
day.Asim .= day.Esim .= day.Gssim .= day.Dlsim .= day.Tlsim .= 0. ± 0.
VcMaxRef,JMaxRef,RdRef,TPURef = PlantBiophysics.fit(Fvcb,day)
g0,g1 = PlantBiophysics.fit(Medlyn,day)

# Simulating
for i in 1:size(day,1)
    meteo = Atmosphere(T = day.T[i] ± 0.1, Wind = 40. ± 10., P = day.P[i] ± 0.001*day.P[i], Rh = day.Rh[i] ± 0.01 ,Cₐ= day.Cₐ[i] ± 10.)
    leaf = LeafModels(energy = Monteith(maxiter=100),
        photosynthesis = Fvcb(VcMaxRef=VcMaxRef,JMaxRef=JMaxRef,RdRef=RdRef,TPURef=TPURef,θ=0.9),
    stomatal_conductance = Medlyn(g0, g1),
    Rₛ = (day.PPFD[i] ± 0.1*day.PPFD[i])/(4.57), sky_fraction =1. ± 0.1,
    PPFD = day.PPFD[i] ± 0.1*day.PPFD[i], d = 0.01 .. 0.10)
    energy_balance!(leaf,meteo)
    day.Asim[i] = leaf[:A]
    day.Esim[i] = leaf[:λE]/(meteo.λ * constants.Mₕ₂ₒ)*1000.
    day.Gssim[i] = leaf[:Gₛ]
    day.Dlsim[i] = leaf[:Dₗ]
    day.Tlsim[i] = leaf[:Tₗ]
end

# Plotting
noto_sans = assetpath("fonts", "NotoSans-Regular.ttf")
fig = Makie.Figure(backgroundcolor = RGBf(1, 1, 1),resolution = (1000, 500),font=noto_sans,dpi=600,size=(1000,1000),xminorgridstyle=true)
#sideinfo1 = Label(fig[2,0], "Simulations", rotation = pi/2,labelsize=25)
#sideinfo2 = Label(fig[2,0], "Observations")
ax = Axis(fig[1,1],title="Leaf temperature in day "*date*" and tree "*string(tree),xlabel="Time (hh:mm:ss)",ylabel="Leaf temperature (°C)")
ax.xticks = ([1,2,3,4,5,6],string.(day.Time))
dat = Makie.scatter!(ax,day.Tₗ,color=:black,marker='∆',markersize=20)
err = errorbars!(ax,[1,2,3,4,5,6],pmean.(day.Tlsim), pstd.(day.Tlsim), pstd.(day.Tlsim),whiskerwidth = 15)
sim = lines!(ax,[1,2,3,4,5,6],pmean.(day.Tlsim),color=:black)
Makie.ylims!(ax,(14.5,27.))
axislegend(ax, [dat,sim,err], ["data","simulations","95% confidence interval"], position = :rb,
    orientation = :vertical,labelsize=13)
fig

save(path*"../tutorials/figure_day.pdf", fig, px_per_unit = 2) # output size = 1600 x 1200 pixels

max_stars_repo_name: VEZY/PlantBiophysics-paper

max_issues_repo_name: VEZY/PlantBiophysics-paper

max_forks_repo_name: VEZY/PlantBiophysics-paper


text:
group = addgroup!(SUITE, "latency")

base_cmd = Base.julia_cmd()
if Base.JLOptions().project != C_NULL
    base_cmd = `$base_cmd --project=$(unsafe_string(Base.JLOptions().project))`
end

# make sure all artifacts are downloaded
CUDA.version()

# time to precompile the package and its dependencies
precompile_cmd =
    `$base_cmd -e "uuid = Base.UUID(\"052768ef-5323-5732-b1bb-66c8b64840ba\")
                   id = Base.PkgId(uuid, \"CUDA\")
                   Base.compilecache(id)"`
group["precompile"] = @benchmarkable run($precompile_cmd) evals=1 seconds=60

# time to actually import the package
import_cmd =
    `$base_cmd -e "using CUDA"`
group["import"] = @benchmarkable run($import_cmd) evals=1 seconds=30

# time to initialize CUDA and all other libraries
initialize_time =
    `$base_cmd -e "using CUDA
                   CUDA.version()"`
group["initialize"] = @benchmarkable run($initialize_time) evals=1 seconds=30

# time to actually compile a kernel
ttfp_cmd =
    `$base_cmd -e "using CUDA
                   kernel() = return
                   CUDA.code_ptx(devnull, kernel, Tuple{}; kernel=true)"`
group["ttfp"] = @benchmarkable run($ttfp_cmd) evals=1 seconds=60


text:
function fakebin(x)
    buffer = IOBuffer(maxsize=sizeof(x))
    for c in x write( buffer, c < '5' ? '0' : '1' ) end
    String(take!(buffer))
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Generator
        service_on::Bool
        pg::Float64
        pgd_ub::Float64
        bus::Int
        pgd_lb::Float64
        name::String
        max_run::Float64
        qg::Float64
        reserve_agc::Bool
        status::Int
        pg_0::Float64
        down_time::Float64
        cost_start_warm::Float64
        reserve_agc_pg_ub::Float64
        cost_pg_model::Int
        time_start_cold::Float64
        qg_ub::Float64
        reserve_agc_pg_lb::Float64
        cost_pg_parameters::Vector{Float64}
        fuel_type::Int
        qg_lb::Float64
        cost_start_hot::Float64
        pg_lb::Float64
        cost_start_cold::Float64
        index::Int
        must_run::Int
        pg_ub::Float64
        min_run::Float64
        time_start_hot::Float64
        up_time::Float64
        time_start_warm::Float64
        min_down::Float64
    end

# Arguments
- `service_on::Bool`: in service
- `pg::Float64`: active/real power generation output
- `pgd_ub::Float64`: maximum active power increase per hour
- `bus::Int`: connecting bus id
- `pgd_lb::Float64`: maximum active power increase per hour
- `name::String`: a flexible name for components, non required to be unique
- `max_run::Float64`: unit of time (hours)
- `qg::Float64`: reactive/imaginary power generation output
- `reserve_agc::Bool`: 1/0
- `status::Int`: a 0/1 value indicating if the component should be omitted or not (0 => omitted)
- `pg_0::Float64`: previous active power output
- `down_time::Float64`: time down previously
- `cost_start_warm::Float64`: cost in dollars of starting a unit
- `reserve_agc_pg_ub::Float64`: 1/0
- `cost_pg_model::Int`: cost model type (1=pwl, 2=polynomial)
- `time_start_cold::Float64`: threshold of time (est. 16h+)
- `qg_ub::Float64`: maximum reactive power generation
- `reserve_agc_pg_lb::Float64`: 1/0
- `cost_pg_parameters::Vector{Float64}`: either points for a pwl cost or polynomial coeffents (from lowest degree first)
- `fuel_type::Int`: fuel type
- `qg_lb::Float64`: minimum active power generation
- `cost_start_hot::Float64`: cost in dollars of starting a unit
- `pg_lb::Float64`: minimum active power generation
- `cost_start_cold::Float64`: cost in dollars of starting a unit
- `index::Int`: a unique identifier for components
- `must_run::Int`: enum
- `pg_ub::Float64`: maximum active power generation
- `min_run::Float64`: unit of time (hours)
- `time_start_hot::Float64`: threshold of time (est. 9h)
- `up_time::Float64`: time up previously
- `time_start_warm::Float64`: threshold of time (est. 9-16h)
- `min_down::Float64`: unit of time (hours)
"""
mutable struct Generator
    "in service"
    service_on::Bool
    "active/real power generation output"
    pg::Float64
    "maximum active power increase per hour"
    pgd_ub::Float64
    "connecting bus id"
    bus::Int
    "maximum active power increase per hour"
    pgd_lb::Float64
    "a flexible name for components, non required to be unique"
    name::String
    "unit of time (hours)"
    max_run::Float64
    "reactive/imaginary power generation output"
    qg::Float64
    "1/0"
    reserve_agc::Bool
    "a 0/1 value indicating if the component should be omitted or not (0 => omitted)"
    status::Int
    "previous active power output"
    pg_0::Float64
    "time down previously"
    down_time::Float64
    "cost in dollars of starting a unit"
    cost_start_warm::Float64
    "1/0"
    reserve_agc_pg_ub::Float64
    "cost model type (1=pwl, 2=polynomial)"
    cost_pg_model::Int
    "threshold of time (est. 16h+)"
    time_start_cold::Float64
    "maximum reactive power generation"
    qg_ub::Float64
    "1/0"
    reserve_agc_pg_lb::Float64
    "either points for a pwl cost or polynomial coeffents (from lowest degree first)"
    cost_pg_parameters::Vector{Float64}
    "fuel type"
    fuel_type::Int
    "minimum active power generation"
    qg_lb::Float64
    "cost in dollars of starting a unit"
    cost_start_hot::Float64
    "minimum active power generation"
    pg_lb::Float64
    "cost in dollars of starting a unit"
    cost_start_cold::Float64
    "a unique identifier for components"
    index::Int
    "enum"
    must_run::Int
    "maximum active power generation"
    pg_ub::Float64
    "unit of time (hours)"
    min_run::Float64
    "threshold of time (est. 9h)"
    time_start_hot::Float64
    "time up previously"
    up_time::Float64
    "threshold of time (est. 9-16h)"
    time_start_warm::Float64
    "unit of time (hours)"
    min_down::Float64
end

function Generator(; service_on, pg, pgd_ub, bus, pgd_lb, name, max_run, qg, reserve_agc, status, pg_0, down_time, cost_start_warm, reserve_agc_pg_ub, cost_pg_model, time_start_cold, qg_ub, reserve_agc_pg_lb, cost_pg_parameters, fuel_type, qg_lb, cost_start_hot, pg_lb, cost_start_cold, index, must_run, pg_ub, min_run, time_start_hot, up_time, time_start_warm, min_down, )
    Generator(service_on, pg, pgd_ub, bus, pgd_lb, name, max_run, qg, reserve_agc, status, pg_0, down_time, cost_start_warm, reserve_agc_pg_ub, cost_pg_model, time_start_cold, qg_ub, reserve_agc_pg_lb, cost_pg_parameters, fuel_type, qg_lb, cost_start_hot, pg_lb, cost_start_cold, index, must_run, pg_ub, min_run, time_start_hot, up_time, time_start_warm, min_down, )
end

# Constructor for demo purposes; non-functional.
function Generator(::Nothing)
    Generator(;
        service_on=1,
        pg=1.2,
        pgd_ub=1.1,
        bus=1,
        pgd_lb=1.1,
        name="1",
        max_run=0,
        qg=1.1,
        reserve_agc=1,
        status=1,
        pg_0=1,
        down_time=1.1,
        cost_start_warm=1,
        reserve_agc_pg_ub=0,
        cost_pg_model=1,
        time_start_cold=1,
        qg_ub=0,
        reserve_agc_pg_lb=0,
        cost_pg_parameters=Any[1.1],
        fuel_type=1,
        qg_lb=0,
        cost_start_hot=1,
        pg_lb=1,
        cost_start_cold=1,
        index=1,
        must_run=1,
        pg_ub=2,
        min_run=0,
        time_start_hot=1,
        up_time=1.1,
        time_start_warm=1,
        min_down=0,
    )
end

"""Get Generator service_on."""
get_service_on(value::Generator) = value.service_on
"""Get Generator pg."""
get_pg(value::Generator) = value.pg
"""Get Generator pgd_ub."""
get_pgd_ub(value::Generator) = value.pgd_ub
"""Get Generator bus."""
get_bus(value::Generator) = value.bus
"""Get Generator pgd_lb."""
get_pgd_lb(value::Generator) = value.pgd_lb
"""Get Generator name."""
get_name(value::Generator) = value.name
"""Get Generator max_run."""
get_max_run(value::Generator) = value.max_run
"""Get Generator qg."""
get_qg(value::Generator) = value.qg
"""Get Generator reserve_agc."""
get_reserve_agc(value::Generator) = value.reserve_agc
"""Get Generator status."""
get_status(value::Generator) = value.status
"""Get Generator pg_0."""
get_pg_0(value::Generator) = value.pg_0
"""Get Generator down_time."""
get_down_time(value::Generator) = value.down_time
"""Get Generator cost_start_warm."""
get_cost_start_warm(value::Generator) = value.cost_start_warm
"""Get Generator reserve_agc_pg_ub."""
get_reserve_agc_pg_ub(value::Generator) = value.reserve_agc_pg_ub
"""Get Generator cost_pg_model."""
get_cost_pg_model(value::Generator) = value.cost_pg_model
"""Get Generator time_start_cold."""
get_time_start_cold(value::Generator) = value.time_start_cold
"""Get Generator qg_ub."""
get_qg_ub(value::Generator) = value.qg_ub
"""Get Generator reserve_agc_pg_lb."""
get_reserve_agc_pg_lb(value::Generator) = value.reserve_agc_pg_lb
"""Get Generator cost_pg_parameters."""
get_cost_pg_parameters(value::Generator) = value.cost_pg_parameters
"""Get Generator fuel_type."""
get_fuel_type(value::Generator) = value.fuel_type
"""Get Generator qg_lb."""
get_qg_lb(value::Generator) = value.qg_lb
"""Get Generator cost_start_hot."""
get_cost_start_hot(value::Generator) = value.cost_start_hot
"""Get Generator pg_lb."""
get_pg_lb(value::Generator) = value.pg_lb
"""Get Generator cost_start_cold."""
get_cost_start_cold(value::Generator) = value.cost_start_cold
"""Get Generator index."""
get_index(value::Generator) = value.index
"""Get Generator must_run."""
get_must_run(value::Generator) = value.must_run
"""Get Generator pg_ub."""
get_pg_ub(value::Generator) = value.pg_ub
"""Get Generator min_run."""
get_min_run(value::Generator) = value.min_run
"""Get Generator time_start_hot."""
get_time_start_hot(value::Generator) = value.time_start_hot
"""Get Generator up_time."""
get_up_time(value::Generator) = value.up_time
"""Get Generator time_start_warm."""
get_time_start_warm(value::Generator) = value.time_start_warm
"""Get Generator min_down."""
get_min_down(value::Generator) = value.min_down

"""Set Generator service_on."""
set_service_on!(value::Generator, val::Bool) = value.service_on = val
"""Set Generator pg."""
set_pg!(value::Generator, val::Float64) = value.pg = val
"""Set Generator pgd_ub."""
set_pgd_ub!(value::Generator, val::Float64) = value.pgd_ub = val
"""Set Generator bus."""
set_bus!(value::Generator, val::Int) = value.bus = val
"""Set Generator pgd_lb."""
set_pgd_lb!(value::Generator, val::Float64) = value.pgd_lb = val
"""Set Generator name."""
set_name!(value::Generator, val::String) = value.name = val
"""Set Generator max_run."""
set_max_run!(value::Generator, val::Float64) = value.max_run = val
"""Set Generator qg."""
set_qg!(value::Generator, val::Float64) = value.qg = val
"""Set Generator reserve_agc."""
set_reserve_agc!(value::Generator, val::Bool) = value.reserve_agc = val
"""Set Generator status."""
set_status!(value::Generator, val::Int) = value.status = val
"""Set Generator pg_0."""
set_pg_0!(value::Generator, val::Float64) = value.pg_0 = val
"""Set Generator down_time."""
set_down_time!(value::Generator, val::Float64) = value.down_time = val
"""Set Generator cost_start_warm."""
set_cost_start_warm!(value::Generator, val::Float64) = value.cost_start_warm = val
"""Set Generator reserve_agc_pg_ub."""
set_reserve_agc_pg_ub!(value::Generator, val::Float64) = value.reserve_agc_pg_ub = val
"""Set Generator cost_pg_model."""
set_cost_pg_model!(value::Generator, val::Int) = value.cost_pg_model = val
"""Set Generator time_start_cold."""
set_time_start_cold!(value::Generator, val::Float64) = value.time_start_cold = val
"""Set Generator qg_ub."""
set_qg_ub!(value::Generator, val::Float64) = value.qg_ub = val
"""Set Generator reserve_agc_pg_lb."""
set_reserve_agc_pg_lb!(value::Generator, val::Float64) = value.reserve_agc_pg_lb = val
"""Set Generator cost_pg_parameters."""
set_cost_pg_parameters!(value::Generator, val::Vector{Float64}) = value.cost_pg_parameters = val
"""Set Generator fuel_type."""
set_fuel_type!(value::Generator, val::Int) = value.fuel_type = val
"""Set Generator qg_lb."""
set_qg_lb!(value::Generator, val::Float64) = value.qg_lb = val
"""Set Generator cost_start_hot."""
set_cost_start_hot!(value::Generator, val::Float64) = value.cost_start_hot = val
"""Set Generator pg_lb."""
set_pg_lb!(value::Generator, val::Float64) = value.pg_lb = val
"""Set Generator cost_start_cold."""
set_cost_start_cold!(value::Generator, val::Float64) = value.cost_start_cold = val
"""Set Generator index."""
set_index!(value::Generator, val::Int) = value.index = val
"""Set Generator must_run."""
set_must_run!(value::Generator, val::Int) = value.must_run = val
"""Set Generator pg_ub."""
set_pg_ub!(value::Generator, val::Float64) = value.pg_ub = val
"""Set Generator min_run."""
set_min_run!(value::Generator, val::Float64) = value.min_run = val
"""Set Generator time_start_hot."""
set_time_start_hot!(value::Generator, val::Float64) = value.time_start_hot = val
"""Set Generator up_time."""
set_up_time!(value::Generator, val::Float64) = value.up_time = val
"""Set Generator time_start_warm."""
set_time_start_warm!(value::Generator, val::Float64) = value.time_start_warm = val
"""Set Generator min_down."""
set_min_down!(value::Generator, val::Float64) = value.min_down = val

max_stars_repo_name: NREL-SIIP/CommonTransmissionModel

max_issues_repo_name: NREL-SIIP/CommonTransmissionModel

max_forks_repo_name: NREL-SIIP/CommonTransmissionModel


text:
abstract type AbstractKernelEmbedding end

measure(ke::AbstractKernelEmbedding) = ke.measure

struct KernelEmbedding{Tk<:Kernel,Tm,Tσ<:Real,Tl} <: AbstractKernelEmbedding
    kernel::Tk # Kernel function
    σ::Tσ # Kernel variance
    l::Tl # Kernel lengthscale
    measure::Tm # Measure
end

function KernelEmbedding(bquad::AbstractBQ, prior)
    return KernelEmbedding(bquad.kernel, bquad.σ, bquad.l, prior)
end

scale(ke::KernelEmbedding) = ke.σ

@doc raw"""
    kernel_mean(ke::KernelEmbedding{Kernel,Measure}, samples::AbstractVector)

Compute the kernel mean of the kernel embedding `ke` for each one of
the `samples` $$x_i$$:
```math
    z_i = \int k(x, x_i)d\mu(x)
```
"""
kernel_mean

function kernel_mean(ke::KernelEmbedding{<:SqExponentialKernel,<:AbstractMvNormal}, samples::AbstractVector)
    z = samples .- Ref(mean(measure(ke)))
    B = Λ(ke.l)
    return scale(ke) / sqrt(det(inv(B) * cov(measure(ke)) + I)) *
           exp.(- PDMats.invquad.(Ref(PDMat(B + cov(measure(ke)))), z) / 2)
end

@doc raw"""
    kernel_variance(ke::KernelEmbedding{Kernel,Measure})

Compute the kernel variance of the given kernel embedding:
```math
    C = \int\int k(x,x')d\mu(x)d\mu(x')
```
"""
kernel_variance

function kernel_variance(ke::KernelEmbedding{<:SqExponentialKernel,<:AbstractMvNormal})
    B = Λ(ke.l)
    return scale(ke) / sqrt(det(2 * inv(B) * cov(measure(ke)) + I))
end

# Turn the lengthscale into a Diagonal matrix of noise
Λ(l::Real) = abs2(l) * I
Λ(l::AbstractVector) = Diagonal(abs2.(l))
Λ(l::LowerTriangular) = Cholesky(l, :L, 1)
Λ(l::AbstractMatrix) = l * l'

# Turns a transform into a lengthscale
transform_param(t::ScaleTransform) = inv(first(t.s))
transform_param(t::ARDTransform) = inv.(t.v)
transform_param(t::LinearTransform) = inv(t.A)


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct EXPIC1 <: AVR
        Tr::Float64
        Ka::Float64
        Ta::Float64
        Va_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Ta_2::Float64
        Ta_3::Float64
        Ta_4::Float64
        Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Kf::Float64
        Tf_1::Float64
        Tf_2::Float64
        Efd_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Ke::Float64
        Te::Float64
        E_sat::Tuple{Float64, Float64}
        Se::Tuple{Float64, Float64}
        Kp::Float64
        Ki::Float64
        Kc::Float64
        V_ref::Float64
        saturation_coeffs::Tuple{Float64, Float64}
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        states_types::Vector{StateTypes.StateType}
        internal::InfrastructureSystemsInternal
    end

Generic Proportional/Integral Excitation System

# Arguments
- `Tr::Float64`: Regulator input filter time constant in s, validation range: `(0, 0.5)`, action if invalid: `warn`
- `Ka::Float64`: Voltage regulator gain, validation range: `(1, 500)`, action if invalid: `warn`
- `Ta::Float64`: Voltage regulator time constant in s, validation range: `(0, 10)`, action if invalid: `warn`
- `Va_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Limits for pi controler `(Vr_min, Vr_max)`
- `Ta_2::Float64`: Voltage regulator time constant in s, validation range: `(0, nothing)`
- `Ta_3::Float64`: Voltage regulator time constant in s, validation range: `(0, nothing)`, action if invalid: `warn`
- `Ta_4::Float64`: Voltage regulator time constant in s, validation range: `(0, nothing)`, action if invalid: `warn`
- `Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Voltage regulator limits (regulator output) (Vi_min, Vi_max)
- `Kf::Float64`: Rate feedback gain, validation range: `(0, 0.3)`, action if invalid: `warn`
- `Tf_1::Float64`: Rate Feedback time constant in s, validation range: `(eps(), 15)`
- `Tf_2::Float64`: Rate Feedback time constant in s, validation range: `(0, 5)`, action if invalid: `warn`
- `Efd_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Field Voltage regulator limits (regulator output) (Efd_min, Efd_max)
- `Ke::Float64`: Exciter constant, validation range: `(0, 1)`, action if invalid: `warn`
- `Te::Float64`: Exciter time constant, validation range: `(0, 2)`, action if invalid: `warn`
- `E_sat::Tuple{Float64, Float64}`: Exciter output voltage for saturation factor: (E1, E2)
- `Se::Tuple{Float64, Float64}`: Exciter saturation factor at exciter output voltage: (Se(E1), Se(E2))
- `Kp::Float64`: Potential source gain, validation range: `(0, 5)`, action if invalid: `warn`
- `Ki::Float64`: current source gain, validation range: `(0, 1.1)`
- `Kc::Float64`: Exciter regulation factor, validation range: `(0, 2)`, action if invalid: `warn`
- `V_ref::Float64`: Reference Voltage Set-point, validation range: `(0, nothing)`
- `saturation_coeffs::Tuple{Float64, Float64}`: Coefficients (A,B) of the function: Se(V) = B(V - A)^2/V
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states are:
	Vm: Sensed terminal voltage,
	Vr1: First Lead-lag state,
	Vr2: Second regulator lead-lag state,
	Vr2: Third regulator lead-lag state
	Vf: Exciter output
	Vr3: First feedback integrator,
	Vr4: second feedback integrator
- `n_states::Int`: EXPIC1 has 6 states
- `states_types::Vector{StateTypes.StateType}`: EXPIC has 6 states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct EXPIC1 <: AVR
    "Regulator input filter time constant in s"
    Tr::Float64
    "Voltage regulator gain"
    Ka::Float64
    "Voltage regulator time constant in s"
    Ta::Float64
    "Limits for pi controler `(Vr_min, Vr_max)`"
    Va_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Voltage regulator time constant in s"
    Ta_2::Float64
    "Voltage regulator time constant in s"
    Ta_3::Float64
    "Voltage regulator time constant in s"
    Ta_4::Float64
    "Voltage regulator limits (regulator output) (Vi_min, Vi_max)"
    Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Rate feedback gain"
    Kf::Float64
    "Rate Feedback time constant in s"
    Tf_1::Float64
    "Rate Feedback time constant in s"
    Tf_2::Float64
    "Field Voltage regulator limits (regulator output) (Efd_min, Efd_max)"
    Efd_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Exciter constant"
    Ke::Float64
    "Exciter time constant"
    Te::Float64
    "Exciter output voltage for saturation factor: (E1, E2)"
    E_sat::Tuple{Float64, Float64}
    "Exciter saturation factor at exciter output voltage: (Se(E1), Se(E2))"
    Se::Tuple{Float64, Float64}
    "Potential source gain"
    Kp::Float64
    "current source gain"
    Ki::Float64
    "Exciter regulation factor"
    Kc::Float64
    "Reference Voltage Set-point"
    V_ref::Float64
    "Coefficients (A,B) of the function: Se(V) = B(V - A)^2/V"
    saturation_coeffs::Tuple{Float64, Float64}
    ext::Dict{String, Any}
    "The states are:
	Vm: Sensed terminal voltage,
	Vr1: First Lead-lag state,
	Vr2: Second regulator lead-lag state,
	Vr2: Third regulator lead-lag state
	Vf: Exciter output
	Vr3: First feedback integrator,
	Vr4: second feedback integrator"
    states::Vector{Symbol}
    "EXPIC1 has 6 states"
    n_states::Int
    "EXPIC has 6 states"
    states_types::Vector{StateTypes.StateType}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function EXPIC1(Tr, Ka, Ta, Va_lim, Ta_2, Ta_3, Ta_4, Vr_lim, Kf, Tf_1, Tf_2, Efd_lim, Ke, Te, E_sat, Se, Kp, Ki, Kc, V_ref=1.0, saturation_coeffs=PowerSystems.get_avr_saturation(E_sat, Se), ext=Dict{String, Any}(), )
    EXPIC1(Tr, Ka, Ta, Va_lim, Ta_2, Ta_3, Ta_4, Vr_lim, Kf, Tf_1, Tf_2, Efd_lim, Ke, Te, E_sat, Se, Kp, Ki, Kc, V_ref, saturation_coeffs, ext, [:Vm, :Vr1, :Vr2, :Vf, :Vr3, :Vr4], 6, [StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Differential, StateTypes.Hybrid], InfrastructureSystemsInternal(), )
end

function EXPIC1(; Tr, Ka, Ta, Va_lim, Ta_2, Ta_3, Ta_4, Vr_lim, Kf, Tf_1, Tf_2, Efd_lim, Ke, Te, E_sat, Se, Kp, Ki, Kc, V_ref=1.0, saturation_coeffs=PowerSystems.get_avr_saturation(E_sat, Se), ext=Dict{String, Any}(), states=[:Vm, :Vr1, :Vr2, :Vf, :Vr3, :Vr4], n_states=6, states_types=[StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Differential, StateTypes.Hybrid], internal=InfrastructureSystemsInternal(), )
    EXPIC1(Tr, Ka, Ta, Va_lim, Ta_2, Ta_3, Ta_4, Vr_lim, Kf, Tf_1, Tf_2, Efd_lim, Ke, Te, E_sat, Se, Kp, Ki, Kc, V_ref, saturation_coeffs, ext, states, n_states, states_types, internal, )
end

# Constructor for demo purposes; non-functional.
function EXPIC1(::Nothing)
    EXPIC1(;
        Tr=0,
        Ka=0,
        Ta=0,
        Va_lim=(min=0.0, max=0.0),
        Ta_2=0,
        Ta_3=0,
        Ta_4=0,
        Vr_lim=(min=0.0, max=0.0),
        Kf=0,
        Tf_1=0,
        Tf_2=0,
        Efd_lim=(min=0.0, max=0.0),
        Ke=0,
        Te=0,
        E_sat=(0.0, 0.0),
        Se=(0.0, 0.0),
        Kp=0,
        Ki=0,
        Kc=0,
        V_ref=0,
        saturation_coeffs=(0.0, 0.0),
        ext=Dict{String, Any}(),
    )
end

"""Get [`EXPIC1`](@ref) `Tr`."""
get_Tr(value::EXPIC1) = value.Tr
"""Get [`EXPIC1`](@ref) `Ka`."""
get_Ka(value::EXPIC1) = value.Ka
"""Get [`EXPIC1`](@ref) `Ta`."""
get_Ta(value::EXPIC1) = value.Ta
"""Get [`EXPIC1`](@ref) `Va_lim`."""
get_Va_lim(value::EXPIC1) = value.Va_lim
"""Get [`EXPIC1`](@ref) `Ta_2`."""
get_Ta_2(value::EXPIC1) = value.Ta_2
"""Get [`EXPIC1`](@ref) `Ta_3`."""
get_Ta_3(value::EXPIC1) = value.Ta_3
"""Get [`EXPIC1`](@ref) `Ta_4`."""
get_Ta_4(value::EXPIC1) = value.Ta_4
"""Get [`EXPIC1`](@ref) `Vr_lim`."""
get_Vr_lim(value::EXPIC1) = value.Vr_lim
"""Get [`EXPIC1`](@ref) `Kf`."""
get_Kf(value::EXPIC1) = value.Kf
"""Get [`EXPIC1`](@ref) `Tf_1`."""
get_Tf_1(value::EXPIC1) = value.Tf_1
"""Get [`EXPIC1`](@ref) `Tf_2`."""
get_Tf_2(value::EXPIC1) = value.Tf_2
"""Get [`EXPIC1`](@ref) `Efd_lim`."""
get_Efd_lim(value::EXPIC1) = value.Efd_lim
"""Get [`EXPIC1`](@ref) `Ke`."""
get_Ke(value::EXPIC1) = value.Ke
"""Get [`EXPIC1`](@ref) `Te`."""
get_Te(value::EXPIC1) = value.Te
"""Get [`EXPIC1`](@ref) `E_sat`."""
get_E_sat(value::EXPIC1) = value.E_sat
"""Get [`EXPIC1`](@ref) `Se`."""
get_Se(value::EXPIC1) = value.Se
"""Get [`EXPIC1`](@ref) `Kp`."""
get_Kp(value::EXPIC1) = value.Kp
"""Get [`EXPIC1`](@ref) `Ki`."""
get_Ki(value::EXPIC1) = value.Ki
"""Get [`EXPIC1`](@ref) `Kc`."""
get_Kc(value::EXPIC1) = value.Kc
"""Get [`EXPIC1`](@ref) `V_ref`."""
get_V_ref(value::EXPIC1) = value.V_ref
"""Get [`EXPIC1`](@ref) `saturation_coeffs`."""
get_saturation_coeffs(value::EXPIC1) = value.saturation_coeffs
"""Get [`EXPIC1`](@ref) `ext`."""
get_ext(value::EXPIC1) = value.ext
"""Get [`EXPIC1`](@ref) `states`."""
get_states(value::EXPIC1) = value.states
"""Get [`EXPIC1`](@ref) `n_states`."""
get_n_states(value::EXPIC1) = value.n_states
"""Get [`EXPIC1`](@ref) `states_types`."""
get_states_types(value::EXPIC1) = value.states_types
"""Get [`EXPIC1`](@ref) `internal`."""
get_internal(value::EXPIC1) = value.internal

"""Set [`EXPIC1`](@ref) `Tr`."""
set_Tr!(value::EXPIC1, val) = value.Tr = val
"""Set [`EXPIC1`](@ref) `Ka`."""
set_Ka!(value::EXPIC1, val) = value.Ka = val
"""Set [`EXPIC1`](@ref) `Ta`."""
set_Ta!(value::EXPIC1, val) = value.Ta = val
"""Set [`EXPIC1`](@ref) `Va_lim`."""
set_Va_lim!(value::EXPIC1, val) = value.Va_lim = val
"""Set [`EXPIC1`](@ref) `Ta_2`."""
set_Ta_2!(value::EXPIC1, val) = value.Ta_2 = val
"""Set [`EXPIC1`](@ref) `Ta_3`."""
set_Ta_3!(value::EXPIC1, val) = value.Ta_3 = val
"""Set [`EXPIC1`](@ref) `Ta_4`."""
set_Ta_4!(value::EXPIC1, val) = value.Ta_4 = val
"""Set [`EXPIC1`](@ref) `Vr_lim`."""
set_Vr_lim!(value::EXPIC1, val) = value.Vr_lim = val
"""Set [`EXPIC1`](@ref) `Kf`."""
set_Kf!(value::EXPIC1, val) = value.Kf = val
"""Set [`EXPIC1`](@ref) `Tf_1`."""
set_Tf_1!(value::EXPIC1, val) = value.Tf_1 = val
"""Set [`EXPIC1`](@ref) `Tf_2`."""
set_Tf_2!(value::EXPIC1, val) = value.Tf_2 = val
"""Set [`EXPIC1`](@ref) `Efd_lim`."""
set_Efd_lim!(value::EXPIC1, val) = value.Efd_lim = val
"""Set [`EXPIC1`](@ref) `Ke`."""
set_Ke!(value::EXPIC1, val) = value.Ke = val
"""Set [`EXPIC1`](@ref) `Te`."""
set_Te!(value::EXPIC1, val) = value.Te = val
"""Set [`EXPIC1`](@ref) `E_sat`."""
set_E_sat!(value::EXPIC1, val) = value.E_sat = val
"""Set [`EXPIC1`](@ref) `Se`."""
set_Se!(value::EXPIC1, val) = value.Se = val
"""Set [`EXPIC1`](@ref) `Kp`."""
set_Kp!(value::EXPIC1, val) = value.Kp = val
"""Set [`EXPIC1`](@ref) `Ki`."""
set_Ki!(value::EXPIC1, val) = value.Ki = val
"""Set [`EXPIC1`](@ref) `Kc`."""
set_Kc!(value::EXPIC1, val) = value.Kc = val
"""Set [`EXPIC1`](@ref) `V_ref`."""
set_V_ref!(value::EXPIC1, val) = value.V_ref = val
"""Set [`EXPIC1`](@ref) `saturation_coeffs`."""
set_saturation_coeffs!(value::EXPIC1, val) = value.saturation_coeffs = val
"""Set [`EXPIC1`](@ref) `ext`."""
set_ext!(value::EXPIC1, val) = value.ext = val
"""Set [`EXPIC1`](@ref) `states_types`."""
set_states_types!(value::EXPIC1, val) = value.states_types = val


text:
using Distributed
addprocs(16)
@everywhere include("../src/get_stable_sens.jl")
s = LinRange(0.,1.,25)
get_stable_sens(s)

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
using ReinforcementLearningCore

(app::NeuralNetworkApproximator)(args...; kwargs...) = app.model(args...; kwargs...)

using AbstractTrees
using TensorBoardLogger: TBLogger

AbstractTrees.children(t::StructTree{T}) where {T<:Union{TBLogger}} = ()


text:
# Geometry which displays points at given (x, y) positions.

struct PointGeometry <: Gadfly.GeometryElement
    tag::Symbol
end
PointGeometry(; tag=empty_tag) = PointGeometry(tag)

const point = PointGeometry

element_aesthetics(::PointGeometry) = [:x, :y, :size, :color, :shape]

# Generate a form for a point geometry.
#
# Args:
#   geom: point geometry.
#   theme: the plot's theme.
#   aes: aesthetics.
#
# Returns:
#   A compose Form.
#
function render(geom::PointGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.point", aes, :x, :y)
    Gadfly.assert_aesthetics_equal_length("Geom.point", aes, :x, :y)

    default_aes = Gadfly.Aesthetics()
    default_aes.shape = Function[Shape.circle]
    default_aes.color = discretize_make_ia(RGBA{Float32}[theme.default_color])
    default_aes.size = Measure[theme.point_size]
    aes = inherit(aes, default_aes)

    if eltype(aes.size) <: Int
      size_min, size_max = extrema(aes.size)
      size_range = size_max - size_min
      point_size_range = theme.point_size_max - theme.point_size_min
      interpolate_size(x) = theme.point_size_min + (x-size_min) / size_range * point_size_range
    end

    ctx = context()

    for (x, y, color, size, shape) in Compose.cyclezip(aes.x, aes.y, aes.color, aes.size, aes.shape)
        shapefun = typeof(shape) <: Function ? shape : theme.point_shapes[shape]
        sizeval = typeof(size) <: Int ? interpolate_size(size) : size
        strokecolor = aes.color_key_continuous != nothing && aes.color_key_continuous ?
                    theme.continuous_highlight_color(color) :
                    theme.discrete_highlight_color(color)
        class = svg_color_class_from_label(aes.color_label([color])[1])
        compose!(ctx, (context(),
              (context(), shapefun([x], [y], [sizeval]), svgclass("marker")),
              fill(color), stroke(strokecolor), svgclass(class)))
    end

    compose!(ctx, linewidth(theme.highlight_width))

    return compose!(context(order=4), svgclass("geometry"), ctx)
end


text:
println("PGM.jl Unit Tests")
println("==============================")

include("utils.jl")
include("lists.jl")
include("variables.jl")
include("graphs.jl")
# include("jointpmf.jl")
# include("dmrf.jl")

println()


text:
# Autogenerated wrapper script for Perl_jll for x86_64-apple-darwin
export libperl, perl

using Readline_jll
JLLWrappers.@generate_wrapper_header("Perl")
JLLWrappers.@declare_library_product(libperl, "@rpath/libperl.dylib")
JLLWrappers.@declare_executable_product(perl)
function __init__()
    JLLWrappers.@generate_init_header(Readline_jll)
    JLLWrappers.@init_library_product(
        libperl,
        "lib/libperl.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        perl,
        "bin/perl",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
abstract type AbstractOperationsModel end

mutable struct PowerOperationModel{M <: AbstractOperationsModel,
                                   T <: PM.AbstractPowerFormulation}
    op_model::Type{M}
    transmission::Type{T}
    devices::Dict{Symbol, DeviceModel}
    branches::Dict{Symbol, DeviceModel}
    services::Dict{Symbol, ServiceModel}
    system::PSY.System
    canonical_model::PSI.CanonicalModel

    function PowerOperationModel(op_model::Type{M},
                                transmission::Type{T},
                                devices::Dict{Symbol, DeviceModel},
                                branches::Dict{Symbol, DeviceModel},
                                services::Dict{Symbol, ServiceModel},
                                system::PSY.System;
                                optimizer::Union{Nothing,JuMP.OptimizerFactory}=nothing,
                                kwargs...) where {M <: AbstractOperationsModel,
                                                  T <: PM.AbstractPowerFormulation}

        ps_model = build_canonical_model(transmission,
                                        devices,
                                        branches,
                                        services,
                                        system,
                                        optimizer;
                                        kwargs...)

        new{M, T}(op_model,
                  transmission,
                  devices,
                  branches,
                  services,
                  system,
                  ps_model)

    end

end