text:
using ImageFiltering
using ImageTransformations
using ImageCore

"""
    resizekern(source_size::Tuple{Int,Int}, dest_size::Tuple{Int,Int})

Create an image resize blur kernel, for use before reducing image size to avoid aliasing.
"""
function resizekern(source_size, dest_size)
    source_size[1:2] == dest_size[1:2] && return nothing
    # Blur before resizing to prevent aliasing (kernel size dependent on both source and target image size)
    σ = 0.1 .* (source_size[1:2] ./ dest_size[1:2])

    if first(σ) < 1
        return ImageFiltering.KernelFactors.gaussian(σ)
    else
        return ImageFiltering.KernelFactors.IIRGaussian(σ)
    end
end

"""
    sizethatfits(original_size::Tuple{Int,Int}, target_shape::Tuple{Int,Int})

Takes an original image size, and fits it to a target shape for image resizing. Maintains aspect ratio.
"""
function sizethatfits(src_size, dst_size)
    scale = minimum([dst_size[1]/src_size[1], dst_size[2]/src_size[2]])
    result = floor.(Int, src_size[[1,2]] .* scale)
    return result
end

"""
    getpadding(size_inner, size_outer)

Get padding in pixels, with size_inner placed at the center of size_outer.
size_inner must be equal or smaller than size_outer
"""
function getpadding(size_inner, size_outer)
    size_diff = size_outer[1:2] .- size_inner[1:2]
    @assert all(size_diff .>= 0) "`size_inner` should be equal or smaller in both dimensions than `size_outer`"
    padding = [floor(Int, size_diff[1]/2), floor(Int, size_diff[2]/2), ceil(Int, size_diff[1]/2), ceil(Int, size_diff[2]/2)]
end

"""
    prepareImage(img::AbstractArray{T}, model::U) where {T<:ImageCore.Colorant, U<:AbstractModel}
    prepareImage(img::AbstractArray{T}, img_resized_size::Tuple, kern) where {T<:ImageCore.Colorant}
    prepareImage!(dest_arr::AbstractArray{Float32}, img::AbstractArray{T}, kern) where {T<:Real}
    prepareImage!(dest_arr::AbstractArray{Float32}, img::AbstractArray{T}, kern) where {T<:ImageCore.Colorant}

Loads and prepares (resizes + pads) an image to fit within a given shape.
Input images should be column-major (julia default), and will be converted to row-major (darknet).
"""
function prepareImage(img::AbstractArray{T}, model::AbstractModel) where {T<:ImageCore.Colorant}
    modelInputSize = getModelInputSize(model)
    if ndims(img) == 3 && size(img)[[3,2,1]] == modelInputSize[1:3]
        return (gpu(permuteddimsview(Float32.(channelview(img)), [3,2,1])), [0,0,0,0])
    elseif size(img)[[2,1]] == modelInputSize[1:2] && modelInputSize[3] == 1
        return (gpu(Float32.(reshape(permuteddimsview(channelview(img), [2,1]), modelInputSize...))), [0,0,0,0])
    elseif size(img)[[2,1]] == modelInputSize[1:2] && modelInputSize[3] == 3
        img_chv = channelview(img)
        if ndims(img_chv) == 2
            return (gpu(Float32.(repeat(reshape(permuteddimsview(img_chv, [2,1]), size(img,2), size(img,1), 1), outer=[1,1,modelInputSize[3]]))), [0,0,0,0])
        elseif size(img_chv, 1) == 1
            return (gpu(Float32.(repeat(permuteddimsview(img_chv, [3,2,1]), outer=[1,1,modelInputSize[3]]))), [0,0,0,0])
        elseif size(img_chv, 1) == 3
            return (gpu(Float32.(permuteddimsview(img_chv, [3,2,1]))), [0,0,0,0])
        elseif size(img_chv, 1) == 4
            return (gpu(Float32.(permuteddimsview(view(img_chv,1:3,:,:), [3,2,1]))), [0,0,0,0])

        else
            error("Image element type $(eltype(img)) not supported")
        end
    else
        img_size = size(img)[[2,1]]
        img_resized_size = sizethatfits(img_size, modelInputSize)
        kern = resizekern(img_size, img_resized_size)
        return prepareImage(img, modelInputSize, kern)
    end
end
function prepareImage(img::AbstractArray{Float32}, model::AbstractModel)
    modelInputSize = getModelInputSize(model)
    if ndims(img) == 3 && size(img) == modelInputSize[1:3]
        return (gpu(img), [0,0,0,0])
    elseif ndims(img) == 3 && size(img)[[3,2,1]] == modelInputSize[1:3]
        return (gpu(permuteddimsview(img, [3,2,1])), [0,0,0,0])
    elseif size(img)[[2,1]] == modelInputSize[1:2] && modelInputSize[3] == 1
        return (gpu(Float32.(reshape(permuteddimsview(img, [2,1]), modelInputSize...))), [0,0,0,0])
    elseif ndims(img) == 2 && size(img)[[2,1]] == modelInputSize[1:2] && modelInputSize[3] == 3
        return (gpu(Float32.(repeat(reshape(permuteddimsview(img, [2,1]), size(img,2), size(img,1), 1), outer=[1,1,modelInputSize[3]]))), [0,0,0,0])
    else
        img_size = size(img)[[2,1]]
        img_resized_size = sizethatfits(img_size, modelInputSize)
        kern = resizekern(img_size, img_resized_size)
        return prepareImage(img, modelInputSize, kern)
    end
end

function prepareImage(img::AbstractArray{Float32}, modelInputSize::Tuple)
    img_size = size(img)[[2,1]]
    img_resized_size = sizethatfits(img_size, modelInputSize)
    kern = resizekern(img_size, img_resized_size)
    return prepareImage!(gpu(zeros(Float32, modelInputSize[1:3])), img, kern)
end

prepareImage(img::AbstractArray{Float32}, modelInputSize::Tuple, kern) =
    prepareImage!(gpu(zeros(Float32, modelInputSize[1:3])), img, kern)

prepareImage(img::AbstractArray{T}, modelInputSize::Tuple, kern) where {T<:ImageCore.Colorant} =
    prepareImage!(gpu(zeros(Float32, modelInputSize[1:3])), img, kern)

function prepareImage!(dest_arr::AbstractArray{Float32}, img::AbstractArray{Float32}, kern)
    #TODO: Make this multiple-dispatchy
    if ndims(img) == 3 && size(img) == size(dest_arr)
        return (gpu(img), [0,0,0,0])
    elseif ndims(img) == 3 && size(img)[[2,1,3]] == size(dest_arr)
        return (gpu(permuteddimsview(img, [2,1,3])), [0,0,0,0])
    elseif ndims(img) == 2 && size(img)[[2,1]] == size(dest_arr)[1:2]
        return (gpu(reshape(permuteddimsview(img, [2,1]), size(img,2), size(img, 1), 1)))
    elseif ndims(img) == 2
        return prepareImage!(dest_arr, colorview(Gray, img), kern)
    elseif size(img, 1) == 1
        return prepareImage!(dest_arr, colorview(Gray, img)[1,:,:], kern)
    elseif size(img, 3) == 1
        return prepareImage!(dest_arr, colorview(Gray, img)[:,:,1], kern)
    elseif size(img, 3) == 3
        return prepareImage!(dest_arr, colorview(RGB, permuteddimsview(img, [3,1,2])), kern)
    else
        error("Array needs to match dimensions exactly, or 3rd dim should be of length 1 or 3 to allow colortype transformations")
    end
end

function prepareImage!(dest_arr::AbstractArray{Float32}, img::AbstractArray{T}, kern) where {T<:ImageCore.Colorant}
    # @show typeof(img), size(img)
    imgPerm = permuteddimsview(img, [2,1])  # Convert from column-major (julia default) to row-major (darknet)

    img_resized_size = sizethatfits(size(imgPerm), size(dest_arr))
    padding = getpadding(img_resized_size, size(dest_arr))
    xidx = (padding[1]+1):(size(dest_arr, 1) - padding[3])
    yidx = (padding[2]+1):(size(dest_arr, 2) - padding[4])
    target_img_subregion = view(dest_arr, xidx, yidx, :)
    if any(size(imgPerm)[1:2] .> size(dest_arr)[1:2]) #Apply blur first if reducing size, to avoid aliasing
        imgblur = ImageFiltering.imfilter(imgPerm, kern, NA())
        imgr = ImageTransformations.imresize(imgblur, img_resized_size[1:2])
        img_chv = channelview(imgr)
        if ndims(img_chv) == 2
            img_ready = repeat(Float32.(reshape(img_chv,(size(img_chv)...,1))),
                        outer=[1,1,size(dest_arr, 3)])
        elseif (size(img_chv,1) != size(dest_arr, 3)) && (size(img_chv, 1) == 1)
            img_ready = repeat(Float32.(permuteddimsview(img_chv, [2,3,1])),
                        outer=[1,1,size(dest_arr, 3)])
        else
            img_ready = permuteddimsview(img_chv, [2,3,1])[:, :, 1:size(dest_arr,3)]
        end
    elseif any(size(imgPerm)[1:2] .< size(dest_arr)[1:2])
        imgr = ImageTransformations.imresize(imgPerm, img_resized_size[1:2])
        img_chv = channelview(imgr)
        if ndims(img_chv) == 2
            img_ready = Float32.(reshape(img_chv, (size(img_chv)...,1)))
        else
            img_ready = view(permuteddimsview(Float32.(img_chv), [2,3,1]), :, :, 1:size(dest_arr,3))
        end
    else
        img_chv = channelview(imgPerm)
        if ndims(img_chv) == 3
            img_ready = view(permuteddimsview(Float32.(img_chv), [2,3,1]), :, :, 1:size(dest_arr,3))
        else
            img_ready = Float32.(repeat(permuteddimsview(img_chv, [1,2]), outer = [1, 1, size(dest_arr,3)]))
        end
    end
    target_img_subregion .= gpu(img_ready)
    scaled_padding = padding ./ size(dest_arr)[[1,2,1,2]]
    return (dest_arr, scaled_padding)
end


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

JLLWrappers.@generate_main_file_header("rocRAND")
JLLWrappers.@generate_main_file("rocRAND", UUID("a6151927-a32b-54c0-bc8c-bbd7b3f1a996"))
end  # module rocRAND_jll

ext: jl


text:
# # Model of predator-prey dynamics

# The predator-prey model emulates the population dynamics of predator and prey animals who
# live in a common ecosystem and compete over limited resources. This model is an
# agent-based analog to the classic
# [Lotka-Volterra](https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations)
# differential equation model. This example illustrates how to develop models with
# heterogeneous agents (sometimes referred to as a *mixed agent based model*).

# The environment is a two dimensional grid containing sheep, wolves and grass. In the
# model, wolves eat sheep and sheep eat grass. Their populations will oscillate over time
# if the correct balance of resources is achieved. Without this balance however, a
# population may become extinct. For example, if wolf population becomes too large,
# they will deplete the sheep and subsequently die of starvation.

# We will begin by loading the required packages and defining three subtypes of
# `AbstractAgent`: `Sheep`, Wolf, and `Grass`. All three agent types have `id` and `pos`
# properties, which is a requirement for all subtypes of `AbstractAgent` when they exist
# upon a `GridSpace`. Sheep and wolves have identical properties, but different behaviors
# as explained below. The property `energy` represents an animals current energy level.
# If the level drops below zero, the agent will die. Sheep and wolves reproduce asexually
# in this model, with a probability given by `reproduction_prob`. The property `Δenergy`
# controls how much energy is acquired after consuming a food source.

# Grass is a replenishing resource that occupies every position in the grid space. Grass can be
# consumed only if it is `fully_grown`. Once the grass has been consumed, it replenishes
# after a delay specified by the property `regrowth_time`. The property `countdown` tracks
# the delay between being consumed and the regrowth time.

# It is also available from the `Models` module as [`Models.predator_prey`](@ref).

using Agents
using Random # hide

mutable struct Sheep <: AbstractAgent
    id::Int
    pos::Dims{2}
    energy::Float64
    reproduction_prob::Float64
    Δenergy::Float64
end

mutable struct Wolf <: AbstractAgent
    id::Int
    pos::Dims{2}
    energy::Float64
    reproduction_prob::Float64
    Δenergy::Float64
end

mutable struct Grass <: AbstractAgent
    id::Int
    pos::Dims{2}
    fully_grown::Bool
    regrowth_time::Int
    countdown::Int
end
nothing # hide

# The function `initialize_model` returns a new model containing sheep, wolves, and grass
# using a set of pre-defined values (which can be overwritten). The environment is a two
# dimensional grid space, which enables animals to walk in all
# directions. Heterogeneous agents are specified in the model as a `Union`. Agents are
# scheduled `schedule_by_type`, which randomizes the order of agents with the constraint that agents
# of a particular type are scheduled consecutively.

function initialize_model(;
    n_sheep = 100,
    n_wolves = 50,
    dims = (20, 20),
    regrowth_time = 30,
    Δenergy_sheep = 4,
    Δenergy_wolf = 20,
    sheep_reproduce = 0.04,
    wolf_reproduce = 0.05,
)
    space = GridSpace(dims, periodic = false)
    model =
        ABM(Union{Sheep,Wolf,Grass}, space, scheduler = schedule_by_type(true, true), warn = false)
    id = 0
    for _ in 1:n_sheep
        id += 1
        energy = rand(1:(Δenergy_sheep*2)) - 1
        ## Note that we must instantiate agents before adding them in a mixed-ABM
        ## to confirm their type.
        sheep = Sheep(id, (0, 0), energy, sheep_reproduce, Δenergy_sheep)
        add_agent!(sheep, model)
    end
    for _ in 1:n_wolves
        id += 1
        energy = rand(1:(Δenergy_wolf*2)) - 1
        wolf = Wolf(id, (0, 0), energy, wolf_reproduce, Δenergy_wolf)
        add_agent!(wolf, model)
    end
    for p in positions(model)
        id += 1
        fully_grown = rand(model.rng, Bool)
        countdown = fully_grown ? regrowth_time : rand(model.rng, 1:regrowth_time) - 1
        grass = Grass(id, (0, 0), fully_grown, regrowth_time, countdown)
        add_agent!(grass, p, model)
    end
    return model
end
nothing # hide

# The function `agent_step!` is dispatched on each subtype in order to produce
# type-specific behavior. The `agent_step!` is similar for sheep and wolves: both lose 1
# energy unit by moving to an adjacent position and both consume a food source if available.
# If their energy level is below zero, an agent dies. Otherwise, the agent lives and
# reproduces with some probability.

function agent_step!(sheep::Sheep, model)
    move!(sheep, model)
    sheep.energy -= 1
    agents = collect(agents_in_position(sheep.pos, model))
    dinner = filter!(x -> isa(x, Grass), agents)
    eat!(sheep, dinner, model)
    if sheep.energy < 0
        kill_agent!(sheep, model)
        return
    end
    if rand(model.rng) <= sheep.reproduction_prob
        reproduce!(sheep, model)
    end
end

function agent_step!(wolf::Wolf, model)
    move!(wolf, model)
    wolf.energy -= 1
    agents = collect(agents_in_position(wolf.pos, model))
    dinner = filter!(x -> isa(x, Sheep), agents)
    eat!(wolf, dinner, model)
    if wolf.energy < 0
        kill_agent!(wolf, model)
        return
    end
    if rand(model.rng) <= wolf.reproduction_prob
        reproduce!(wolf, model)
    end
end
nothing # hide

# The behavior of grass functions differently. If it is fully grown, it is consumable.
# Otherwise, it cannot be consumed until it regrows after a delay specified by
# `regrowth_time`.

function agent_step!(grass::Grass, model)
    if !grass.fully_grown
        if grass.countdown <= 0
            grass.fully_grown = true
            grass.countdown = grass.regrowth_time
        else
            grass.countdown -= 1
        end
    end
end
nothing # hide

# Sheep and wolves move to a random adjacent position with the `move!` function.
function move!(agent, model)
    neighbors = nearby_positions(agent, model)
    position = rand(model.rng, collect(neighbors))
    move_agent!(agent, position, model)
end
nothing # hide

# Sheep and wolves have separate `eat!` functions. If a sheep eats grass, it will acquire
# additional energy and the grass will not be available for consumption until regrowth time
# has elapsed. If a wolf eats a sheep, the sheep dies and the wolf acquires more energy.

function eat!(sheep::Sheep, grass_array, model)
    isempty(grass_array) && return
    grass = grass_array[1]
    if grass.fully_grown
        sheep.energy += sheep.Δenergy
        grass.fully_grown = false
    end
end

function eat!(wolf::Wolf, sheep, model)
    if !isempty(sheep)
        dinner = rand(model.rng, sheep)
        kill_agent!(dinner, model)
        wolf.energy += wolf.Δenergy
    end
end
nothing # hide

# Sheep and wolves share a common reproduction method. Reproduction has a cost of 1/2 the
# current energy level of the parent. The offspring is an exact copy of the parent, with
# exception of `id`.

function reproduce!(agent, model)
    agent.energy /= 2
    id = nextid(model)
    A = typeof(agent)
    offspring = A(id, agent.pos, agent.energy, agent.reproduction_prob, agent.Δenergy)
    add_agent_pos!(offspring, model)
    return
end
nothing # hide

# ## Running the model
# We will run the model for 500 steps and record the number of sheep, wolves and consumable
# grass patches after each step. First: initialize the model.
using InteractiveDynamics
using CairoMakie
Random.seed!(23182) # hide
n_steps = 500
model = initialize_model()

# To view our starting population, we can build an overview plot:

offset(a::Sheep) = (0.2, 0.0)
offset(a::Wolf) = (-0.2, 0.0)
offset(a::Grass) = (0.0, 0.0)
mshape(a::Sheep) = '⚫'
mshape(a::Wolf) = '▲'
mshape(a::Grass) = '■'
mcolor(a::Sheep) = RGBAf0(1.0, 1.0, 1.0, 0.8)
mcolor(a::Wolf) = RGBAf0(0.2, 0.2, 0.2, 0.8)
mcolor(a::Grass) = cgrad([:brown, :green])[a.countdown/a.regrowth_time]
figure, = abm_plot(
    model;
    resolution = (800, 600),
    offset = offset,
    am = mshape,
    as = 22,
    ac = mcolor,
    scheduler = schedule_by_type((Grass, Sheep, Wolf), false),
)
figure
# Now, lets run the simulation and collect some data.

sheep(a) = typeof(a) == Sheep
wolves(a) = typeof(a) == Wolf
grass(a) = typeof(a) == Grass && a.fully_grown
adata = [(sheep, count), (wolves, count), (grass, count)]
results, _ = run!(model, agent_step!, n_steps; adata)

# The plot shows the population dynamics over time. Initially, wolves become extinct because they
# consume the sheep too quickly. The few remaining sheep reproduce and gradually reach an
# equilibrium that can be supported by the amount of available grass.
figure = Figure(resolution = (600, 400))
ax = figure[1, 1] = Axis(figure; xlabel = "Step", ylabel = "Population")
sheepl = lines!(ax, results.step, results.count_sheep, color = :blue)
wolfl = lines!(ax, results.step, results.count_wolves, color = :orange)
grassl = lines!(ax, results.step, results.count_grass, color = :green)
figure[1, 2] = Legend(figure, [sheepl, wolfl, grassl], ["Sheep", "Wolves", "Grass"])
figure

# Altering the input conditions, we now see a landscape where all three agents find an
# equilibrium.

Random.seed!(7756) # hide
model = initialize_model(
    n_wolves = 20,
    dims = (25, 25),
    Δenergy_sheep = 5,
    sheep_reproduce = 0.2,
    wolf_reproduce = 0.08,
)
results, _ = run!(model, agent_step!, n_steps; adata)

figure = Figure(resolution = (600, 400))
ax = figure[1, 1] = Axis(figure, xlabel = "Step", ylabel = "Population")
sheepl = lines!(ax, results.step, results.count_sheep, color = :blue)
wolfl = lines!(ax, results.step, results.count_wolves, color = :orange)
grassl = lines!(ax, results.step, results.count_grass, color = :green)
figure[1, 2] = Legend(figure, [sheepl, wolfl, grassl], ["Sheep", "Wolves", "Grass"])
figure


text:
# # Modeling with non-analytic Fourier transforms

using Comrade
using Plots

# While most of the models implemented in `Comrade` have an analytic
# Fourier transform this is not required. In this notebook we will
# describe how a user can do Bayesian model fitting with a
# non-analytic model.

# The `ExtendedRing` model is an example of a non-analytic model. The
# image structure is given by
# ```math
# I(r) = \frac{\beta^\alpha}{2\pi \Gamma(\alpha)} r^{-\alpha-2}e^{-\beta/r}
# ```

# This can be created as follows

m = ExtendedRing(10.0, 8.0)
# The first argument is the radius ``r_0`` of the ring and the second is ``\alpha`` in the above
# equation. The ``beta`` in the above equation is then given by ``r_0(1+\alpha)``.

# This is an example of a ring model that has a substantially different flux profile.

plot(m, xlims=(-40.0, 40.0), ylims=(-40.0, 40.0), uvscale=identity)

# This function does not have a simple analytic Fourier transform, e.g.

Comrade.visanalytic(ExtendedRing)

# Therefore, to find the Fourier transform of the image we need to revert to numerical methods.
# For this notebook we will use the *fast Fourier transform* or FFT. Specifically we will
# use FFTW. To compute a numerical Fourier transform we first need to specify the image.

image = IntensityMap(zeros(256, 256), 50.0, 50.0)

# This will serve as our cache to store the image going forward. The next step is to create
# a model wrapper that holds the model and the image. `Comrade` provides the `modelimage`
# function to do exactly that

mimage = modelimage(m, image; alg=Comrade.FFT())

# the `alg` keyword argument then specifies that we want to use an FFT to compute the
# Fourier transform. When `modelimage` is called, the FFT is performed and then we use
# a bicubic interpolator on the resulting visibilities to construct a continuous representation
# of the Fourier transform. Once we have this everything else is the same. Namely we can
# calculatge the VLBI data products in the usual manner i.e.

u = randn(1000)/2
v = randn(1000)/2

vis = visibilities(mimage, u, v)
amp = amplitudes(mimage, u, v)

# We can also compute derivatives of these models. However, currently only `ForwardDiff`
# is able to propogate derivatives through `modelimage`. This is due to the interpolation
# step. In the future this limitation will be removed.


text:
using Test
using SafeTestsets

@testset "CircuitConstants" begin

@safetestset "SI" begin
    using CircuitConstants: T, G, M, k, m, u, n, p, f, a
    @test T == 1e12
    @test G == 1e9
    @test M == 1e6
    @test k == 1e3
    @test m == 1e-3
    @test u == 1e-6
    @test n == 1e-9
    @test p == 1e-12
    @test f == 1e-15
    @test a == 1e-18
end
@safetestset "Dups" begin
    using CircuitConstants: μ
    @test μ == 1e-6
end
@safetestset "V" begin
    using CircuitConstants: TV, GV, MV, kV, V, mV, uV, μV, nV, pV, fV, aV
    @test TV == 1e12
    @test GV == 1e9
    @test MV == 1e6
    @test kV == 1e3
    @test V == 1
    @test mV == 1e-3
    @test uV == 1e-6
    @test μV == 1e-6
    @test nV == 1e-9
    @test pV == 1e-12
    @test fV == 1e-15
    @test aV == 1e-18
end
@safetestset "A" begin
    using CircuitConstants: TA, GA, MA, kA, A, mA, uA, μA, nA, pA, fA, aA
    @test TA == 1e12
    @test GA == 1e9
    @test MA == 1e6
    @test kA == 1e3
    @test A == 1
    @test mA == 1e-3
    @test uA == 1e-6
    @test μA == 1e-6
    @test nA == 1e-9
    @test pA == 1e-12
    @test fA == 1e-15
    @test aA == 1e-18
end
@safetestset "s" begin
    using CircuitConstants: Ts, Gs, Ms, ks, s, ms, us, μs, ns, ps, fs, as
    @test Ts == 1e12
    @test Gs == 1e9
    @test Ms == 1e6
    @test ks == 1e3
    @test s == 1
    @test ms == 1e-3
    @test us == 1e-6
    @test μs == 1e-6
    @test ns == 1e-9
    @test ps == 1e-12
    @test fs == 1e-15
    @test as == 1e-18
end
@safetestset "sec" begin
    using CircuitConstants: Tsec, Gsec, Msec, ksec, sec, msec, usec, μsec, nsec, psec, fsec, asec
    @test Tsec == 1e12
    @test Gsec == 1e9
    @test Msec == 1e6
    @test ksec == 1e3
    @test sec == 1
    @test msec == 1e-3
    @test usec == 1e-6
    @test μsec == 1e-6
    @test nsec == 1e-9
    @test psec == 1e-12
    @test fsec == 1e-15
    @test asec == 1e-18
end
@safetestset "Hz" begin
    using CircuitConstants: THz, GHz, MHz, kHz, Hz, mHz, uHz, μHz, nHz, pHz, fHz, aHz
    @test THz == 1e12
    @test GHz == 1e9
    @test MHz == 1e6
    @test kHz == 1e3
    @test Hz == 1
    @test mHz == 1e-3
    @test uHz == 1e-6
    @test μHz == 1e-6
    @test nHz == 1e-9
    @test pHz == 1e-12
    @test fHz == 1e-15
    @test aHz == 1e-18
end
@safetestset "F" begin
    using CircuitConstants: TF, GF, MF, kF, F, mF, uF, μF, nF, pF, fF, aF
    @test TF == 1e12
    @test GF == 1e9
    @test MF == 1e6
    @test kF == 1e3
    @test F == 1
    @test mF == 1e-3
    @test uF == 1e-6
    @test μF == 1e-6
    @test nF == 1e-9
    @test pF == 1e-12
    @test fF == 1e-15
    @test aF == 1e-18
end
@safetestset "H" begin
    using CircuitConstants: TH, GH, MH, kH, H, mH, uH, μH, nH, pH, fH, aH
    @test TH == 1e12
    @test GH == 1e9
    @test MH == 1e6
    @test kH == 1e3
    @test H == 1
    @test mH == 1e-3
    @test uH == 1e-6
    @test μH == 1e-6
    @test nH == 1e-9
    @test pH == 1e-12
    @test fH == 1e-15
    @test aH == 1e-18
end
@safetestset "Ohms" begin
    using CircuitConstants: TOhms, GOhms, MOhms, kOhms, Ohms, mOhms, uOhms, μOhms, nOhms, pOhms, fOhms, aOhms
    @test TOhms == 1e12
    @test GOhms == 1e9
    @test MOhms == 1e6
    @test kOhms == 1e3
    @test Ohms == 1
    @test mOhms == 1e-3
    @test uOhms == 1e-6
    @test μOhms == 1e-6
    @test nOhms == 1e-9
    @test pOhms == 1e-12
    @test fOhms == 1e-15
    @test aOhms == 1e-18
end
@safetestset "Ω" begin
    using CircuitConstants: TΩ, GΩ, MΩ, kΩ, Ω, mΩ, uΩ, μΩ, nΩ, pΩ, fΩ, aΩ
    @test TΩ == 1e12
    @test GΩ == 1e9
    @test MΩ == 1e6
    @test kΩ == 1e3
    @test Ω == 1
    @test mΩ == 1e-3
    @test uΩ == 1e-6
    @test μΩ == 1e-6
    @test nΩ == 1e-9
    @test pΩ == 1e-12
    @test fΩ == 1e-15
    @test aΩ == 1e-18
end

end

max_stars_repo_name: GlenHertz/CircuitConstants

max_issues_repo_name: GlenHertz/CircuitConstants

max_forks_repo_name: GlenHertz/CircuitConstants


text:
using TimeSeriesClustering
data_path=normpath(joinpath(dirname(@__FILE__),"..","data","TS_GER_1"))
ts_input_data = load_timeseries_data(data_path; T=24, years=[2016])

attribute_weights=Dict("solar"=>1.0, "wind"=>2.0, "el_demand"=>3.0)

clust_res=run_clust(ts_input_data;n_init=10,n_clust=4,attribute_weights=attribute_weights) # default k-means


text:
# This file contains the tests for the manually compiled wrapper

using Test

include(joinpath(@__DIR__, "..", "..", "deps", "compile.jl"))

const lib = joinpath(@__DIR__, "..", "..", "deps", "PDFHighlightsWrapper")
const pdf = joinpath(@__DIR__, "..", "pdf", "TestPDF.pdf")

@testset "get_author_title" begin

    # Prepare pointers
    author_ref = Ref{Cstring}()
    title_ref = Ref{Cstring}()

    ccall(
        (:get_author_title, lib),
        Cvoid,
        (Cstring, Ref{Cstring}, Ref{Cstring}),
        pdf, author_ref, title_ref,
    )

    @test unsafe_string(author_ref[]) == "Pavel Sobolev"
    @test unsafe_string(title_ref[]) == "A dummy PDF for tests"

end

macro unsafe_wrap(array::Symbol, len::Symbol)
    return esc(:(unsafe_wrap(Array, $(array)[], $(len)[]; own = true)))
end

macro unsafe_wrap(array::Expr, len::Union{Symbol, Expr})
    return esc(:(unsafe_wrap(Array, $(array), $(len); own = true)))
end

macro unsafe_wrap_strings(array::Union{Symbol, Expr}, len::Union{Symbol, Expr})
    return esc(:(unsafe_string.(@unsafe_wrap $(array) $(len))))
end

@testset "get_lines_comments_pages" begin

    # Prepare pointers
    lines_arrays_ref = Ref{Ptr{Ptr{Cstring}}}()
    lines_x_anchors_arrays_ref = Ref{Ptr{Ptr{Cdouble}}}()
    lines_yl_anchors_arrays_ref = Ref{Ptr{Ptr{Cdouble}}}()
    lines_yu_anchors_arrays_ref = Ref{Ptr{Ptr{Cdouble}}}()
    lines_lens_ref = Ref{Ptr{Culong}}()
    comments_ref = Ref{Ptr{Cstring}}()
    pages_ref = Ref{Ptr{Cint}}()
    num_ref = Ref{Culong}(0)

    # Call
    ccall(
        # Function and library names
        (
            :get_lines_comments_pages,
            lib,
        ),
        # Type of return value
        Cvoid,
        # Types of the arguments
        (
            Cstring,
            Ref{Ptr{Ptr{Cstring}}},
            Ref{Ptr{Ptr{Cdouble}}},
            Ref{Ptr{Ptr{Cdouble}}},
            Ref{Ptr{Ptr{Cdouble}}},
            Ref{Ptr{Culong}},
            Ref{Ptr{Cstring}},
            Ref{Ptr{Cint}},
            Ref{Culong},
        ),
        # Arguments
        pdf,
        lines_arrays_ref,
        lines_x_anchors_arrays_ref,
        lines_yl_anchors_arrays_ref,
        lines_yu_anchors_arrays_ref,
        lines_lens_ref,
        comments_ref,
        pages_ref,
        num_ref,
    )

    # Steal the upper-level arrays
    lines_arrays = @unsafe_wrap lines_arrays_ref num_ref
    lines_x_anchors_arrays = @unsafe_wrap lines_x_anchors_arrays_ref num_ref
    lines_yl_anchors_arrays = @unsafe_wrap lines_yl_anchors_arrays_ref num_ref
    lines_yu_anchors_arrays = @unsafe_wrap lines_yu_anchors_arrays_ref num_ref
    lines_lens = @unsafe_wrap lines_lens_ref num_ref
    comments = @unsafe_wrap_strings comments_ref num_ref
    pages = @unsafe_wrap pages_ref num_ref

    highlights = Vector{String}(undef, num_ref[])

    highlights_x_anchors = Vector{Float64}(undef, num_ref[])
    highlights_yl_anchors = Vector{Float64}(undef, num_ref[])
    highlights_yu_anchors = Vector{Float64}(undef, num_ref[])

    for index in 1:num_ref[]

        # Steal the lower-level array
        len = lines_lens[index]
        lines = @unsafe_wrap_strings lines_arrays[index] len
        lines_x_anchors = @unsafe_wrap lines_x_anchors_arrays[index] len
        lines_yl_anchors = @unsafe_wrap lines_yl_anchors_arrays[index] len
        lines_yu_anchors = @unsafe_wrap lines_yu_anchors_arrays[index] len

        # Save the highlights anchors
        highlights_x_anchors[index] = lines_x_anchors[1]
        highlights_yl_anchors[index] = lines_yl_anchors[end]
        highlights_yu_anchors[index] = lines_yu_anchors[1]

        # Highlight string
        highlight = ""

        hyphen_chopped = false

        # Remove word transfers
        for (index, line) in enumerate(lines)
            if endswith(line, "-")
                if hyphen_chopped
                    if index == len
                        highlight = string(highlight, line)
                    else
                        highlight = string(highlight, chop(line))
                    end
                else
                    if len == 1
                        highlight = string(highlight, line)
                    elseif index == 1
                        highlight = string(highlight, chop(line))
                    elseif index == len
                        highlight = string(highlight, ' ', line)
                    else
                        highlight = string(highlight, ' ', chop(line))
                    end
                end
                hyphen_chopped = true
            else
                if hyphen_chopped
                    highlight = string(highlight, line)
                else
                    if index == 1
                        highlight = string(highlight, line)
                    else
                        highlight = string(highlight, ' ', line)
                    end
                end
                hyphen_chopped = false
            end
        end

        # Combine all lines in a highlight and save it
        highlights[index] = highlight

    end

    @test highlights == [
        "Highlight 1",
        "Highlight 2",
        "Highlight 3",
        "High-",
        "light 4",
        "Highhighlight 5",
        "6th Highhigh light-",
        "High light 7",
        "8th Highlight-",
    ]

    @test comments == [
        "Comment 1",
        ".c1 Comment 2",
        ".c2 Comment 3",
        ".c1 Comment 4",
        ".c2",
        "",
        "",
        "",
        "",
    ]

    @test pages == Int32[1, 2, 3, 4, 5, 6, 7, 8, 9]

    @test highlights_x_anchors == repeat([42.52], 9)
    @test highlights_yl_anchors == [
        repeat([24.50399999999999], 5);
        [
            53.396000000000015,
            67.84200000000001,
            38.95000000000002,
            38.95000000000002
        ]
    ]
    @test highlights_yu_anchors == repeat([39.59200000000001], 9)

end


text:
module Reflection

using FlatBuffers: get, pos, @scopedenum, offset, Table, indirect, init, Array, bytes

@scopedenum BaseType::Int8 None=0 UType=1 Bool=2 Byte=3 UByte=4 Short=5 UShort=6 Int=7 UInt=8 Long=9 ULong=10 Float=11 Double=12 String=13 Vector=14 Obj=15 Union=16 Array=17 MaxBaseType=18

@scopedenum AdvancedFeatures::UInt64 AdvancedArrayFeatures=1 AdvancedUnionFeatures=2 OptionalScalars=4 DefaultVectorsAndStrings=8

struct Type_  <: Table
	bytes::Vector{UInt8}
	pos::Base.Int
end

Base.propertynames(::Type_) = (
	:base_type,
	:element,
	:index,
	:fixed_length,
)

function Base.getproperty(x::Type_, field::Symbol)
	if field === :base_type
		o = offset(x, 4)
		o != 0 && return get(x, o + pos(x), BaseType)
		return BaseType(0)
	elseif field === :element
		o = offset(x, 6)
		o != 0 && return get(x, o + pos(x), BaseType)
		return BaseType(0)
	elseif field === :index
		o = offset(x, 8)
		o != 0 && return get(x, o + pos(x), Int32)
		return Int32(-1)
	elseif field === :fixed_length
		o = offset(x, 10)
		o != 0 && return get(x, o + pos(x), UInt16)
		return UInt16(0)
	end
	return nothing
end

struct KeyValue_ <: Table
	bytes::Vector{UInt8}
	pos::Base.Int
end

Base.propertynames(::KeyValue_) = (
	:key,
	:value,
)

function Base.getproperty(x::KeyValue_, field::Symbol)
	if field === :key
		o = offset(x, 4)
		o != 0 && return String(x, o + pos(x))
		return string(0)
	elseif field === :value
		o = offset(x, 6)
		o != 0 && return String(x, o + pos(x))
		return string(0)
	end
	return nothing
end

struct EnumVal_ <: Table
	bytes::Vector{UInt8}
	pos::Base.Int
end

Base.propertynames(::EnumVal_) = (
	:name,
	:value,
	:union_type,
	:documentation,
)

function Base.getproperty(x::EnumVal_, field::Symbol)
	if field === :name
		o = offset(x, 4)
		o != 0 && return String(x, o + pos(x))
		return string(0)
	elseif field === :value
		o = offset(x, 6)
		o != 0 && return get(x, o + pos(x), Int64)
		return Int64(0)
	elseif field === :union_type
		o = offset(x, 10)
		if o != 0
			y = indirect(x, o + pos(x))
			return init(Type_, bytes(x), y)
		end
	elseif field === :documentation
		o = offset(x, 12)
		o != 0 && return Array{[]byte}(x, o)
	end
	return nothing
end

struct Enum_ <: Table
	bytes::Vector{UInt8}
	pos::Base.Int
end

Base.propertynames(::Enum_) = (
	:name,
	:values,
	:is_union,
	:underlying_type,
	:attributes,
	:documentation,
	:declaration_file,
)

function Base.getproperty(x::Enum_, field::Symbol)
	if field === :name
		o = offset(x, 4)
		o != 0 && return String(x, o + pos(x))
		return string(0)
	elseif field === :values
		o = offset(x, 6)
		o != 0 && return Array{EnumVal_}(x, o)
	elseif field === :is_union
		o = offset(x, 8)
		o != 0 && return get(x, o + pos(x), Bool)
		return Bool(false)
	elseif field === :underlying_type
		o = offset(x, 10)
		if o != 0
			y = indirect(x, o + pos(x))
			return init(Type_, bytes(x), y)
		end
	elseif field === :attributes
		o = offset(x, 12)
		o != 0 && return Array{KeyValue_}(x, o)
	elseif field === :documentation
		o = offset(x, 14)
		o != 0 && return Array{[]byte}(x, o)
	elseif field === :declaration_file
		o = offset(x, 16)
		o != 0 && return String(x, o + pos(x))
		return string(0)
	end
	return nothing
end

struct Field_ <: Table
	bytes::Vector{UInt8}
	pos::Base.Int
end

Base.propertynames(::Field_) = (
	:name,
	:type,
	:id,
	:offset,
	:default_integer,
	:default_real,
	:deprecated,
	:required,
	:key,
	:attributes,
	:documentation,
	:optional,
)

function Base.getproperty(x::Field_, field::Symbol)
	if field === :name
		o = offset(x, 4)
		o != 0 && return String(x, o + pos(x))
		return string(0)
	elseif field === :type
		o = offset(x, 6)
		if o != 0
			y = indirect(x, o + pos(x))
			return init(Type_, bytes(x), y)
		end
	elseif field === :id
		o = offset(x, 8)
		o != 0 && return get(x, o + pos(x), UInt16)
		return UInt16(0)
	elseif field === :offset
		o = offset(x, 10)
		o != 0 && return get(x, o + pos(x), UInt16)
		return UInt16(0)
	elseif field === :default_integer
		o = offset(x, 12)
		o != 0 && return get(x, o + pos(x), Int64)
		return Int64(0)
	elseif field === :default_real
		o = offset(x, 14)
		o != 0 && return get(x, o + pos(x), Float64)
		return Float64(0.0)
	elseif field === :deprecated
		o = offset(x, 16)
		o != 0 && return get(x, o + pos(x), Bool)
		return Bool(false)
	elseif field === :required
		o = offset(x, 18)
		o != 0 && return get(x, o + pos(x), Bool)
		return Bool(false)
	elseif field === :key
		o = offset(x, 20)
		o != 0 && return get(x, o + pos(x), Bool)
		return Bool(false)
	elseif field === :attributes
		o = offset(x, 22)
		o != 0 && return Array{KeyValue_}(x, o)
	elseif field === :documentation
		o = offset(x, 24)
		o != 0 && return Array{[]byte}(x, o)
	elseif field === :optional
		o = offset(x, 26)
		o != 0 && return get(x, o + pos(x), Bool)
		return Bool(false)
	end
	return nothing
end

struct Object_ <: Table
	bytes::Vector{UInt8}
	pos::Base.Int
end

Base.propertynames(::Object_) = (
	:name,
	:fields,
	:is_struct,
	:minalign,
	:bytesize,
	:attributes,
	:documentation,
	:declaration_file,
)

function Base.getproperty(x::Object_, field::Symbol)
	if field === :name
		o = offset(x, 4)
		o != 0 && return String(x, o + pos(x))
		return string(0)
	elseif field === :fields
		o = offset(x, 6)
		o != 0 && return Array{Field_}(x, o)
	elseif field === :is_struct
		o = offset(x, 8)
		o != 0 && return get(x, o + pos(x), Bool)
		return Bool(false)
	elseif field === :minalign
		o = offset(x, 10)
		o != 0 && return get(x, o + pos(x), Int32)
		return Int32(0)
	elseif field === :bytesize
		o = offset(x, 12)
		o != 0 && return get(x, o + pos(x), Int32)
		return Int32(0)
	elseif field === :attributes
		o = offset(x, 14)
		o != 0 && return Array{KeyValue_}(x, o)
	elseif field === :documentation
		o = offset(x, 16)
		o != 0 && return Array{[]byte}(x, o)
	elseif field === :declaration_file
		o = offset(x, 18)
		o != 0 && return String(x, o + pos(x))
		return string(0)
	end
	return nothing
end

struct RPCCall_ <: Table
	bytes::Vector{UInt8}
	pos::Base.Int
end

Base.propertynames(::RPCCall_) = (
	:name,
	:request,
	:response,
	:attributes,
	:documentation,
)

function Base.getproperty(x::RPCCall_, field::Symbol)
	if field === :name
		o = offset(x, 4)
		o != 0 && return String(x, o + pos(x))
		return string(0)
	elseif field === :request
		o = offset(x, 6)
		if o != 0
			y = indirect(x, o + pos(x))
			return init(Object_, bytes(x), y)
		end
	elseif field === :response
		o = offset(x, 8)
		if o != 0
			y = indirect(x, o + pos(x))
			return init(Object_, bytes(x), y)
		end
	elseif field === :attributes
		o = offset(x, 10)
		o != 0 && return Array{KeyValue_}(x, o)
	elseif field === :documentation
		o = offset(x, 12)
		o != 0 && return Array{[]byte}(x, o)
	end
	return nothing
end

struct Service_ <: Table
	bytes::Vector{UInt8}
	pos::Base.Int
end

Base.propertynames(::Service_) = (
	:name,
	:calls,
	:attributes,
	:documentation,
	:declaration_file,
)

function Base.getproperty(x::Service_, field::Symbol)
	if field === :name
		o = offset(x, 4)
		o != 0 && return String(x, o + pos(x))
		return string(0)
	elseif field === :calls
		o = offset(x, 6)
		o != 0 && return Array{RPCCall_}(x, o)
	elseif field === :attributes
		o = offset(x, 8)
		o != 0 && return Array{KeyValue_}(x, o)
	elseif field === :documentation
		o = offset(x, 10)
		o != 0 && return Array{[]byte}(x, o)
	elseif field === :declaration_file
		o = offset(x, 12)
		o != 0 && return String(x, o + pos(x))
		return string(0)
	end
	return nothing
end

struct SchemaFile_ <: Table
	bytes::Vector{UInt8}
	pos::Base.Int
end

Base.propertynames(::SchemaFile_) = (
	:filename,
	:included_filenames,
)

function Base.getproperty(x::SchemaFile_, field::Symbol)
	if field === :filename
		o = offset(x, 4)
		o != 0 && return String(x, o + pos(x))
		return string(0)
	elseif field === :included_filenames
		o = offset(x, 6)
		o != 0 && return Array{[]byte}(x, o)
	end
	return nothing
end

struct Schema_ <: Table
	bytes::Vector{UInt8}
	pos::Base.Int
end

Base.propertynames(::Schema_) = (
	:objects,
	:enums,
	:file_ident,
	:file_ext,
	:root_table,
	:services,
	:advanced_features,
	:fbs_files,
)

function Base.getproperty(x::Schema_, field::Symbol)
	if field === :objects
		o = offset(x, 4)
		o != 0 && return Array{Object_}(x, o)
	elseif field === :enums
		o = offset(x, 6)
		o != 0 && return Array{Enum_}(x, o)
	elseif field === :file_ident
		o = offset(x, 8)
		o != 0 && return String(x, o + pos(x))
		return string(0)
	elseif field === :file_ext
		o = offset(x, 10)
		o != 0 && return String(x, o + pos(x))
		return string(0)
	elseif field === :root_table
		o = offset(x, 12)
		if o != 0
			y = indirect(x, o + pos(x))
			return init(Object_, bytes(x), y)
		end
	elseif field === :services
		o = offset(x, 14)
		o != 0 && return Array{Service_}(x, o)
	elseif field === :advanced_features
		o = offset(x, 16)
		o != 0 && return get(x, o + pos(x), AdvancedFeatures)
		return AdvancedFeatures(1)
	elseif field === :fbs_files
		o = offset(x, 18)
		o != 0 && return Array{SchemaFile_}(x, o)
	end
	return nothing
end

end #module


text:
function _convert(ctx::ConverterContext, decl::CLVarDecl)
	kind = isconst(type(decl)) ? "CBinding.Cglobalconst" : "CBinding.Cglobal"
	typ  = _convertType(ctx, type(decl))
	name = _convertName(ctx, decl)
	csym = repr(Symbol(spelling(decl)))

	if !(name in ctx.oneofs)
		push!(ctx.oneofs, name)

		_export(ctx, name)
		def = "const $(name) = $(kind){$(typ)}()"
		push!(ctx.converted, JuliaizedC(
			decl,
			def,
			:atcompile_bindings,
		))

		load = "CBinding.bind($(name), $(csym), $(_gensym(ctx, "libraries"))...)"
		push!(ctx.converted, JuliaizedC(
			decl,
			load,
			:atload,
		))
	end
end


text:
using JuliaBlackBoard
## Authoring a test using pools means creating several differnt .txt files
## This shows how they can all be generated from one script
## The basic idea is to use `POOL` defined below to open a new pool
## file to write to for each desired pool (in the directory of the script)
##
## One way to use
## Author a test without randomization. Once happy,
## edit each file to be randomized to add the `POOL(io) do iop` call
## and randomize the question within the `POOL()...end` block
##

OPEN = JuliaBlackBoard.OPEN(@__FILE__)
POOL = JuliaBlackBoard.POOL(@__FILE__)

OPEN() do io

    q = mt"""
    # Question 1

How goes?
"""
    question(io, ESS, q)

    ## Pool question 1
    ## written to script_directory/script_name-pool-i.txt where i is 1 below
    POOL(io) do iop
        q = mt"""
# Addition

${{:n}} + {{:m}} = $ ?
"""

        for n in 2:5, m in 2:5
            question(iop, NUM, q(n=n,m=m), n+m)
        end
    end

    ## Pool qustion 2
    POOL(io) do iop
        q = mt"""
# Times tables:

${{:n}} * {{:m}} =$ ?
"""
        for n in 2:7, m in 2:7
            question(iop, NUM, q(n=n,m=m), n*m)
        end
    end

end


text:
using Documenter, WingBiology

makedocs(sitename = "WingBiology")

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
module Minvar
##################Dependencies####################
using Revise

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

if "$(pwd())\\data" ∉ LOAD_PATH
  push!(LOAD_PATH,"$(pwd())\\data")
end

using DataFrames, Distributions, GZip, ForwardDiff, Dates, LinearAlgebra, Statistics,
   Measures, Formatting, StatsBase, Distributions, DataStructures, Finometrics, Distributed,
  Serialization,  CSV, Random, Distributed#, Gadfly, Fontconfig, Cairo NOTE: uncomment for graphs

#########################Exported functions##########
export testportfolio,
  testcov,
  testtrisectedportfolio,
  matlab2julia,
  mathematica2julia,
  testω,
  testposteriors,
  estimateminvariance,
  initializepar, #need this to set up the parallelization
  cleanup

################Priors and initial vlaues############
#Priors
const D_θG = 0.16^2/20 #current vix converted to daily
const D_δ²G = .76/20*2 #reasonably diffuse prior from volatiltiy of the vix multiplied by 2
const D_αG = 1.01 #Pick a diffuse prior, but want the expectation to exist
const D_βG = 7.0e-5 / 2.0 * (D_αG - 1.0) #derive from vix of vix (See documentation)
const D_αP = D_αG
const D_βP = D_βG * 2.0 * (D_αP - 1.0) #expect this is more volatile.
const D_θSG = D_θG
const D_δ²SG = D_δ²G*2.0 #more uncertain
const D_θSGP = D_θG
const D_δ²SGP = D_δ²G*4.0 #most uncertain

#Initial Parameter Values
const D_σ²G=.001
const D_ζ²G= .004
const D_ζ²P = .009
const D_SG = .002
const D_SGP = .003

##############Other constants and methedological parameters############
#const TEX_PATH = "C:\\Users\\Clinton\\Dropbox\\Apps\\Overleaf\\Endowment Project"
#const GRAPH_PATH = TEX_PATH * "\\sub\\fig"
#const WORKING_PATH = pwd() * "\\working"

const OUT_NAME = "test"
const OUT_PATH = "output"

const N_BURN = 100
const N_SAMPLES = 100

const DATA_NAME = "sampledata"
const DATA_PATH = pwd() * "\\data"

const DATE_BEGIN = Date(2015,1,1)
const DATE_END =  Date(2016,12,31)
const DATE_FORMAT = "yyyymmdd"

const MAX_WORKERS = Base.Sys.CPU_THREADS ÷ 2 # gets num of physical cores

const NUM_TEST_ASSETS = 100 #this shouldn't make a big difference in the results
const MIN_ADJUSTMENT = 0.75 #smaller values->more extreme weights allowed for test portfolio

#without this throws an error on each run if already defined
if !(@isdefined RANDOM_WEIGHT_DIST)
  const RANDOM_WEIGHT_DIST = TriangularDist(-1.0,1.5,1/NUM_TEST_ASSETS)
end

#############File linkage###############################
include("code\\MinvarPar.jl")
include("code\\Universe.jl")
include("code\\Portfolio.jl")
include("code\\Parameters.jl")
include("code\\Posteriors.jl")
include("code\\TrisectedPortfolio.jl")
include("code\\X2Julia.jl")
include("code\\Estimation.jl")
include("code\\Testsandoneoffs.jl")
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module Gen
    include("address.jl")
    include("assignment.jl")
    include("homogenous_trie.jl")
    include("generator.jl")
    include("distribution.jl")
    include("dsl_common.jl")
    include("generators/generators.jl")
    include("injective.jl")
    include("selection.jl")
    include("inference.jl")
end # module

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
## Uniform ##

@scalar_rule(
    uniformlogpdf(a, b, x),
    @setup(
        insupport = a <= x <= b,
        diff = b - a,
        c = insupport ? inv(diff) : inv(one(diff)),
        z = insupport ? zero(x) : oftype(x, NaN),
    ),
    (c, -c, z),
)

# StatsFuns: https://github.com/JuliaStats/StatsFuns.jl/pull/106

## Beta ##

@scalar_rule(
    betalogpdf(α::Real, β::Real, x::Number),
    @setup(z = digamma(α + β)),
    (
        log(x) + z - digamma(α),
        log1p(-x) + z - digamma(β),
        (α - 1) / x + (1 - β) / (1 - x),
    ),
)

## Gamma ##

@scalar_rule(
    gammalogpdf(k::Real, θ::Real, x::Number),
    @setup(
        invθ = inv(θ),
        xoθ = invθ * x,
        z = xoθ - k,
    ),
    (
        log(xoθ) - digamma(k),
        invθ * z,
        - (1 + z) / x,
    ),
)

## Chisq ##

@scalar_rule(
    chisqlogpdf(k::Real, x::Number),
    @setup(hk = k / 2),
    (
        (log(x) - logtwo - digamma(hk)) / 2,
        (hk - 1) / x - one(hk) / 2,
    ),
)

## FDist ##

@scalar_rule(
    fdistlogpdf(ν1::Real, ν2::Real, x::Number),
    @setup(
        xν1 = x * ν1,
        temp1 = xν1 + ν2,
        a = (x - 1) / temp1,
        νsum = ν1 + ν2,
        di = digamma(νsum / 2),
    ),
    (
        (-log1p(ν2 / xν1) - ν2 * a + di - digamma(ν1 / 2)) / 2,
        (-log1p(xν1 / ν2) + ν1 * a + di - digamma(ν2 / 2)) / 2,
        ((ν1 - 2) / x - ν1 * νsum / temp1) / 2,
    ),
)

## TDist ##

@scalar_rule(
    tdistlogpdf(ν::Real, x::Number),
    @setup(
        νp1 = ν + 1,
        xsq = x^2,
        invν = inv(ν),
        a = xsq * invν,
        b = νp1 / (ν + xsq),
    ),
    (
        (digamma(νp1 / 2) - digamma(ν / 2) + a * b - log1p(a) - invν) / 2,
        - x * b,
    ),
)

## Binomial ##

@scalar_rule(
    binomlogpdf(n::Real, p::Real, k::Real),
    @setup(z = digamma(n - k + 1)),
    (
        digamma(n + 2) - z + log1p(-p) - 1 / (1 + n),
        (k / p - n) / (1 - p),
        z - digamma(k + 1) + logit(p),
    ),
)

## Poisson ##

@scalar_rule(
    poislogpdf(λ::Number, x::Number),
    ((iszero(x) && iszero(λ) ? zero(x / λ) : x / λ) - 1, log(λ) - digamma(x + 1)),
)

## PoissonBinomial

function ChainRulesCore.rrule(
    ::typeof(Distributions.poissonbinomial_pdf_fft), p::AbstractVector{<:Real}
)
    y = Distributions.poissonbinomial_pdf_fft(p)
    A = poissonbinomial_partialderivatives(p)
    function poissonbinomial_pdf_fft_pullback(Δy)
        p̄ = InplaceableThunk(
            @thunk(A * Δy),
            Δ -> LinearAlgebra.mul!(Δ, A, Δy, true, true),
        )
        return (NO_FIELDS, p̄)
    end
    return y, poissonbinomial_pdf_fft_pullback
end

if isdefined(Distributions, :poissonbinomial_pdf)
    function ChainRulesCore.rrule(
        ::typeof(Distributions.poissonbinomial_pdf), p::AbstractVector{<:Real}
    )
        y = Distributions.poissonbinomial_pdf(p)
        A = poissonbinomial_partialderivatives(p)
        function poissonbinomial_pdf_pullback(Δy)
            p̄ = InplaceableThunk(
                @thunk(A * Δy),
                Δ -> LinearAlgebra.mul!(Δ, A, Δy, true, true),
            )
            return (NO_FIELDS, p̄)
        end
        return y, poissonbinomial_pdf_pullback
    end
end


text:
struct ElectricPotential{T, N, S} <: AbstractArray{T, N}
    data::Array{T, N}
    grid::Grid{T, N, S}
end

@inline size(ep::ElectricPotential{T, N, S}) where {T, N, S} = size(ep.data)
@inline length(ep::ElectricPotential{T, N, S}) where {T, N, S} = length(ep.data)
@inline getindex(ep::ElectricPotential{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = getindex(ep.data, I...)
@inline getindex(ep::ElectricPotential{T, N, S}, i::Int) where {T, N, S} = getindex(ep.data, i)
@inline getindex(ep::ElectricPotential{T, N, S}, s::Symbol) where {T, N, S} = getindex(ep.grid, s)

"""
    ElectricPotential(setup::PotentialSimulationSetup{T, 3, :cylindrical} ; kwargs...)::ElectricPotential{T, 3, :cylindrical}

Extracts the electric potential from `setup` and extrapolate it to an 2π grid.

For 2D grids (r and z) the user has to set the keyword `n_points_in_φ::Int`, e.g.: `n_points_in_φ = 36`.
"""
function ElectricPotential(setup::PotentialSimulationSetup{T, 3, :cylindrical} ; kwargs...)::ElectricPotential{T, 3, :cylindrical} where {T}
    return get_2π_potential(ElectricPotential{T, 3, :cylindrical}(setup.potential, setup.grid); kwargs...)
end

"""
    ElectricPotential(setup::PotentialSimulationSetup{T, 3, :cartesian} ; kwargs...)::ElectricPotential{T, 3, :cartesian}

Extracts the electric potential from `setup`.
"""
function ElectricPotential(setup::PotentialSimulationSetup{T, 3, :cartesian} )::ElectricPotential{T, 3, :cartesian} where {T}
    return ElectricPotential{T, 3, :cartesian}(setup.potential, setup.grid)
end

function NamedTuple(ep::ElectricPotential{T, 3}) where {T}
    return (
        grid = NamedTuple(ep.grid),
        values = ep.data * u"V",
    )
end
Base.convert(T::Type{NamedTuple}, x::ElectricPotential) = T(x)

function ElectricPotential(nt::NamedTuple)
    grid = Grid(nt.grid)
    T = typeof(ustrip(nt.values[1]))
    S = get_coordinate_system(grid)
    N = get_number_of_dimensions(grid)
    ElectricPotential{T, N, S}( ustrip.(uconvert.(u"V", nt.values)), grid)
end
Base.convert(T::Type{ElectricPotential}, x::NamedTuple) = T(x)