text:
mutable struct Scene <: AbstractScene
    parent
    events::Events

    px_area::Node{IRect2D}
    camera::Camera
    camera_controls::RefValue
    limits::Node{FRect3D}

    transformation::Transformation

    plots::Vector{AbstractPlot}
    theme::Attributes
    attributes::Attributes
    children::Vector{Scene}
    current_screens::Vector{AbstractScreen}
    updated::Node{Bool}
end

Base.haskey(scene::Scene, key::Symbol) = haskey(scene.attributes, key)
Base.getindex(scene::Scene, key::Symbol) = scene.attributes[key]
Base.setindex!(scene::Scene, value, key::Symbol) = (scene.attributes[key] = value)

function Scene(
        events::Events,
        px_area::Node{IRect2D},
        camera::Camera,
        camera_controls::RefValue,
        limits::Node,
        transformation::Transformation,
        plots::Vector{AbstractPlot},
        theme::Attributes,
        children::Vector{Scene},
        current_screens::Vector{AbstractScreen},
        parent = nothing,
    )
    # indicates whether we can start updating the plot
    # needs to be set in display
    updated = Node(false)

    scene = Scene(
        parent, events, px_area, camera, camera_controls, limits,
        transformation, plots, theme, Attributes(),
        children, current_screens, updated
    )
    finalizer(scene) do scene
        # save_print("Freeing scene")
        close_all_nodes(scene.events)
        close_all_nodes(scene.transformation)
        for field in (:px_area, :limits)
            close(getfield(scene, field))
        end
        disconnect!(scene.camera)
        empty!(scene.theme)
        empty!(scene.attributes)
        empty!(scene.children)
        empty!(scene.current_screens)
        return
    end
    onany(updated, px_area) do update, px_area
        if update && !(scene.camera_controls[] isa PixelCamera)
            a = scene.attributes
            if to_value(get(a, :scale_plot, false)) && !a[:raw][]
                scale_scene!(scene)
            end
            if to_value(get(a, :center, false)) && !a[:raw][]
                center!(scene)
            end
        end
        nothing
    end
    scene
end

Base.parent(scene::Scene) = scene.parent
isroot(scene::Scene) = parent(scene) === nothing
function root(scene::Scene)
    while !isroot(scene)
        scene = parent(scene)
    end
    scene
end
parent_or_self(scene::Scene) = isroot(scene) ? scene : parent(scene)

Base.size(x::Scene) = pixelarea(x) |> to_value |> widths |> Tuple
Base.size(x::Scene, i) = size(x)[i]
function Base.resize!(scene::Scene, xy::Tuple{Number, Number})
    Base.resize!(scene, IRect(0, 0, xy))
end
Base.resize!(scene::Scene, x::Number, y::Number) = resize!(scene, (x, y))
function Base.resize!(scene::Scene, rect::Rect2D)
    pixelarea(scene)[] = rect
end

"""
    `update!(p::Scene)`

Updates a `Scene` and all its children.
"""
function update!(p::Scene)
    p.updated[] = true
    for c in p.children
        update!(c)
    end
end

# Just indexing into a scene gets you plot 1, plot 2 etc
Base.iterate(scene::Scene, idx = 1) = idx <= length(scene) ? (scene[idx], idx + 1) : nothing
Base.length(scene::Scene) = length(scene.plots)
Base.lastindex(scene::Scene) = length(scene.plots)
getindex(scene::Scene, idx::Integer) = scene.plots[idx]
GeometryTypes.widths(scene::Scene) = widths(to_value(pixelarea(scene)))
struct Axis end

zero_origin(area) = IRect(0, 0, widths(area))

child(scene::Scene) = Scene(scene, lift(zero_origin, pixelarea(scene)))

"""
Creates a subscene with a pixel camera
"""
function cam2d(scene::Scene)
    sub = child(scene)
    cam2d!(sub)
    sub
end
function campixel(scene::Scene)
    sub = child(scene)
    campixel!(sub)
    sub
end

function getindex(scene::Scene, ::Type{Axis})
    for plot in scene
        isaxis(plot) && return plot
    end
    nothing
end

"""
Each argument can be named for a certain plot type `P`. Falls back to `arg1`, `arg2`, etc.
"""
function argument_names(plot::P) where P <: AbstractPlot
    argument_names(P, length(plot.converted))
end

function argument_names(::Type{<: AbstractPlot}, num_args::Integer)
    # this is called in the indexing function, so let's be a bit efficient
    ntuple(i-> Symbol("arg$i"), num_args)
end

function Base.empty!(scene::Scene)
    empty!(scene.plots)
    disconnect!(scene.camera)
    scene.limits[] = FRect3D()
    scene.camera_controls[] = EmptyCamera()
    empty!(scene.theme)
    merge!(scene.theme, _current_default_theme)
    empty!(scene.children)
    empty!(scene.current_screens)
end

limits(scene::Scene) = scene.limits
limits(scene::SceneLike) = scene.parent.limits

# Since we can use Combined like a scene in some circumstances, we define this alias
theme(x::SceneLike, args...) = theme(x.parent, args...)
theme(x::Scene) = x.theme
theme(x::Scene, key) = x.theme[key]
theme(x::AbstractPlot, key) = x.attributes[key]
theme(::Nothing, key::Symbol) = current_default_theme()[key]

Base.push!(scene::Combined, subscene) = nothing # Combined plots add themselves uppon creation
function Base.push!(scene::Scene, plot::AbstractPlot)
    push!(scene.plots, plot)
    plot.parent[] = scene
    for screen in scene.current_screens
        insert!(screen, scene, plot)
    end
end
function Base.push!(scene::Scene, plot::Combined)
    push!(scene.plots, plot)
    for screen in scene.current_screens
        insert!(screen, scene, plot)
    end
end

function Observables.connect!(scene::Scene, child::Scene)

end

function Base.push!(scene::Scene, child::Scene)
    push!(scene.children, child)
    disconnect!(child.camera)
    nodes = map([:view, :projection, :projectionview, :resolution, :eyeposition]) do field
        lift(getfield(scene.camera, field)) do val
            push!(getfield(child.camera, field), val)
        end
    end
    cameracontrols!(child, nodes)
end

events(scene::Scene) = scene.events
events(scene::SceneLike) = events(scene.parent)

camera(scene::Scene) = scene.camera
camera(scene::SceneLike) = camera(scene.parent)

cameracontrols(scene::Scene) = scene.camera_controls[]
cameracontrols(scene::SceneLike) = cameracontrols(scene.parent)

cameracontrols!(scene::Scene, cam) = (scene.camera_controls[] = cam)
cameracontrols!(scene::SceneLike, cam) = cameracontrols!(parent(scene), cam)

pixelarea(scene::Scene) = scene.px_area
pixelarea(scene::SceneLike) = pixelarea(scene.parent)

plots(scene::SceneLike) = scene.plots

const _forced_update_scheduled = Ref(false)

"""
Returns whether a scene needs to be updated
"""
function must_update()
    val = _forced_update_scheduled[]
    _forced_update_scheduled[] = false
    val
end

"""
Forces the scene to be re-rendered
"""
function force_update!()
    _forced_update_scheduled[] = true
end

const current_global_scene = Ref{Any}()

if Sys.iswindows()
    function _primary_resolution()
        # ccall((:GetSystemMetricsForDpi, :user32), Cint, (Cint, Cuint), 0, ccall((:GetDpiForSystem, :user32), Cuint, ()))
        # ccall((:GetSystemMetrics, :user32), Cint, (Cint,), 17)
        dc = ccall((:GetDC, :user32), Ptr{Cvoid}, (Ptr{Cvoid},), C_NULL)
        ntuple(2) do i
            Int(ccall((:GetDeviceCaps, :gdi32), Cint, (Ptr{Cvoid}, Cint), dc, (2 - i) + 117))
        end
    end
else
    # TODO implement osx + linux
    _primary_resolution() = (1920, 1080) # everyone should have at least a hd monitor :D
end

"""
Returns the resolution of the primary monitor.
If the primary monitor can't be accessed, returns (1920, 1080) (full hd)
"""
function primary_resolution()
    # Since this is pretty low level and os specific + we can't test on all possible
    # computers, I assume we'll have bugs here. Let's not sweat about it too much,
    # we just use primary_resolution to have a good estimate for a default window resolution
    # if this fails, only thing happening will be a too small/big window when the user doesn't give any resolution.
    try
        _primary_resolution()
    catch e
        @warn("Could not retrieve primary monitor resolution. A default resolution of (1920, 1080) is assumed!
        Error: $(sprint(io->showerror(io, e))).")
        (1920, 1080)
    end
end

"""
Returns a reasonable resolution for the main monitor.
(right now just half the resolution of the main monitor)
"""
reasonable_resolution() = primary_resolution() .÷ 2

"""
Returns the current active scene (the last scene that got created)
"""
function current_scene()
    if isassigned(current_global_scene)
        current_global_scene[]
    else
        Scene()
    end
end

# dirty hack used for Scene(@resolution)
# could possibly be removed, since the macro isn't really used anymore
Scene(::Nothing) = Scene()

function Scene(;
        kw_args...
    )
    events = Events()
    theme = current_default_theme(; kw_args...)
    resolution = theme[:resolution][]
    px_area = Observable(IRect(0, 0, resolution))
    on(events.window_area) do w_area
        if !any(x-> x ≈ 0.0, widths(w_area)) && px_area[] != w_area
            px_area[] = w_area
        end
    end
    scene = Scene(
        events,
        px_area,
        Camera(px_area),
        RefValue{Any}(EmptyCamera()),
        node(:scene_limits, FRect3D(Vec3f0(0), Vec3f0(1))),
        Transformation(),
        AbstractPlot[],
        theme,
        Scene[],
        AbstractScreen[]
    )
    current_global_scene[] = scene
    scene
end

function Scene(
        scene::Scene;
        events = scene.events,
        px_area = scene.px_area,
        cam = scene.camera,
        camera_controls = scene.camera_controls,
        boundingbox = Node(AABB(Vec3f0(0), Vec3f0(1))),
        transformation = Transformation(scene),
        theme = Theme(),
        current_screens = scene.current_screens
    )
    child = Scene(
        events,
        px_area,
        cam,
        camera_controls,
        boundingbox,
        transformation,
        AbstractPlot[],
        merge(current_default_theme(), theme),
        Scene[],
        current_screens,
        scene
    )
    push!(scene.children, child)
    child
end

function Scene(parent::Scene, area; theme...)
    events = parent.events
    px_area = lift(pixelarea(parent), to_node(area)) do p, a
        # make coordinates relative to parent
        IRect2D(minimum(p) .+ minimum(a), widths(a))
    end
    child = Scene(
        events,
        px_area,
        Camera(px_area),
        RefValue{Any}(EmptyCamera()),
        node(:scene_limits, FRect3D(Vec3f0(0), Vec3f0(1))),
        Transformation(),
        AbstractPlot[],
        current_default_theme(; theme...),
        Scene[],
        parent.current_screens,
        parent
    )
    push!(parent.children, child)
    child
end

"""
Fetches all plots sharing the same camera
"""
plots_from_camera(scene::Scene) = plots_from_camera(scene, scene.camera)
function plots_from_camera(scene::Scene, camera::Camera, list = AbstractPlot[])
    append!(list, scene.plots)
    for child in scene.children
        child.camera == camera && plots_from_camera(child, camera, list)
    end
    list
end

"""
Flattens all the combined plots and returns a Vector of Atomic plots
"""
function flatten_combined(plots::Vector, flat = AbstractPlot[])
    for elem in plots
        if (elem isa Combined)
            flatten_combined(elem.plots, flat)
        else
            push!(flat, elem)
        end
    end
    flat
end

function insertplots!(screen::AbstractDisplay, scene::Scene)
    for elem in scene.plots
        insert!(screen, scene, elem)
    end
    foreach(s-> insertplots!(screen, s), scene.children)
end
update_cam!(scene::Scene, bb::AbstractCamera, rect) = nothing

function scale_scene!(scene::Scene)
    if is2d(scene)
        area = pixelarea(scene)[]
        lims = limits(scene)[]
        # not really sure how to scale 3D scenes in a reasonable way
        mini, maxi = minimum(lims), maximum(lims)
        l = ((mini[1], maxi[1]), (mini[2], maxi[2]))
        xyzfit = fit_ratio(area, l)
        s = to_ndim(Vec3f0, xyzfit, 1f0)
        scale!(scene, s)
    end
    return scene
end

function center!(scene::Scene, padding = 0.01)
    bb = boundingbox(scene)
    bb = transformationmatrix(scene)[] * bb
    w = widths(bb)
    padd = w .* padding
    bb = FRect3D(minimum(bb) .- padd, w .+ 2padd)
    update_cam!(scene, bb)
    scene
end
parent_scene(x::Combined) = parent_scene(parent(x))
parent_scene(x::Scene) = x

Base.isopen(x::SceneLike) = events(x).window_open[]


text:
module CVRPInstanceModule

using Distributions
using Random
using LinearAlgebra

export CVRPInstance, plotInstance, randInstance, readInstance

mutable struct CVRPInstance
  # n: number of customers
  # node 1 to n are customers. node n+1 is the depot
  #nn (=n+1): number of nodes
  n::Int32
  nn::Int32
  Q::Int32
  coords::Array{Int64,2}
  demand::Array{Int32,1}
  # cost matrix
  c::Array{Float64,2}
  name::String
end
# default constructor
CVRPInstance() = CVRPInstance(1,2,100,[ 1 1 ], [1], [0 0; 0 0], "unknown" )

# n is the number of customers. Number of nodes is n+1 (called nn in CVRPInstance)
function randInstance(instance::CVRPInstance, n, Q, seed)
  instance.n = n
  nn =n+1
  instance.nn = nn
  instance.Q = Q
  instance.coords = zeros(Int64, nn,2)
  instance.demand = zeros(Int32, nn)
  Random.seed!(seed)
  # Generate random coordinates in the box (0,0) to (1000,1000)
  rand!(instance.coords, collect(0:1000))
  rand!(instance.demand, collect(1:30))
  # set demand of the depot to zero
  instance.demand[nn] = 0;

  # this is copy-paste from TSP example by ???
  instance.c = zeros(nn, nn)
  for i = 1:nn
      for j = i:nn
          d = norm(instance.coords[i,1:2] - instance.coords[j,1:2])
          instance.c[i,j] = d
          instance.c[j,i] = d
      end
  end
end

# n is the number of customers. Number of nodes is n+1 (called nn in CVRPInstance)
# Q: vehicle capacity
# maxDemand: customer demand is selected as a random integer in the interval [1,maxDemand]
# boxSize: customers are generated within a [0,boxSize] x [0,boxSize] box

# nCustClustered: How many customers are placed in clusters
# nClusters: How many clusters are there?
# For each cluster we generate a random  center (x',y') in the box. We also draw
# number "stdDev" as a rational number in the interval [minStdDevClust, maxStdDevClust]
# that indicates how "spread out the cluster is"
# coordinates of clusters are generated by drawing x from a truncated normal
# distribution with mean x' and std. deviation <stdDev>. The distribution is truncated to the
# interval [0,boxSize]. The same happens for the y coordinate.

# nCustConCircles: How many customers are drawn from concentric circles around the depot
# We generate <nConCircles> concentric circles. The circles are centered in midle of the
# coordinate box (independent of the depot location),
# the i'th (i=1..nConCircles) circle has radius (i/nConCircles)*(boxSize/2)
# implying that the larges has radius boxSize/2 and the smallest has radius
# (1/nConCircles)*(boxSize/2)

# depotMode: 1 = center of box, 2=one of the corners, 3=random position in the box

# @@@ to do: make demand generation "more interesting"
function randInstanceAdvanced(instance::CVRPInstance, n, Q, maxDemand, boxSize,
      nCustClustered, nClusters, minStdDevClust, maxStdDevClust,
      nCustConCircles, nConCircles,
      depotMode, seed)
  # some initial error checking:
  if nCustClustered+nCustConCircles > n
    throw("randInstanceAdvanced(...): nCustClustered + nCustConCircles > n")
  end
  if maxDemand > Q
    throw("randInstanceAdvanced(...): maxDemand > Q")
  end

  srand(seed)

  instance.n = n
  nn = n+1
  instance.nn = nn
  instance.Q = Q
  #instance.coords = zeros(Int64, nn,2)
  # generate demands:
  instance.demand = rand(1:maxDemand, nn)
  # set demand of the depot to zero
  instance.demand[nn] = 0;

  coordsUniform = []
  coordsClusters = []
  coordsCircles = []
  # Generate random coordinates in the box
  if nCustClustered+nCustConCircles < n
    nCustUniform = n - (nCustClustered + nCustConCircles)
    coordsUniform = rand(0:boxSize, nCustUniform, 2)
  end

  # Generate coordinates in clusters
  if nCustClustered > 0
    coordsClusters = zeros(Int64, nCustClustered,2)
    clusterCenters = rand(0:boxSize, nClusters, 2)
    clusterStdDev = rand(Uniform(minStdDevClust,maxStdDevClust), nClusters)

    distribX = []
    distribY = []
    for i=1:nClusters
      push!(distribX, TruncatedNormal(clusterCenters[i,1], clusterStdDev[i], 0, boxSize))
      push!(distribY, TruncatedNormal(clusterCenters[i,2], clusterStdDev[i], 0, boxSize))
    end
    for i=1:nCustClustered
      clusterId=rand(1:nClusters)
      x = round(rand(distribX[clusterId]))
      y = round(rand(distribY[clusterId]))
      coordsClusters[i,1] = x
      coordsClusters[i,2] = y
    end
  end

  # Generate coordinates in concentric circles
  if nCustConCircles > 0
    coordsCircles = zeros(Int64, nCustConCircles,2)
    radiuses = zeros(Float64, nConCircles)
    for i=1:nConCircles
      radiuses[i] = (i/nConCircles)*(boxSize/2)
    end
    for i=1:nCustConCircles
      # choose a random circle
      circleId = rand(1:nConCircles)
      # choose a random angle
      angle = rand(Uniform(0,2*pi))
      coordsCircles[i,1] = round(boxSize/2+radiuses[circleId]*cos(angle))
      coordsCircles[i,1] = round(boxSize/2+radiuses[circleId]*sin(angle))
    end
  end

  # depotMode: 1 = center of box, 2=one of the corners, 3=random position in the box
  depotCoord = zeros(Int64, 1,2)
  if depotMode == 1
    depotCoord = [boxSize/2 boxSize/2]
  elseif depotMode == 2
    depotCoord = [rand(0:1)*boxSize rand(0:1)*boxSize]
  else
    depotCoord = rand(0:boxSize, 1, 2)
  end
  instance.coords = vcat(coordsUniform, coordsClusters, coordsCircles, depotCoord)

  # this is copy-paste from TSP example by ???
  instance.c = zeros(nn, nn)
  for i = 1:nn
      for j = i:nn
          d = norm(instance.coords[i,1:2] - instance.coords[j,1:2])
          instance.c[i,j] = d
          instance.c[j,i] = d
      end
  end
end

function findKey(key::String, splitFile::Vector{SubString{String}})
    lcaseKey = lowercase(key)
    for i=1:length(splitFile)
      if lowercase(splitFile[i]) == lcaseKey
        return i
      end
    end
    return -1
end

function readInstance(fileName::String, roundDist::Bool)
    instance = CVRPInstance()
    open(fileName) do f
        #allFile = readstring(f)
        allFile = read(f, String)
        splitFile = split(allFile)

        instance.name = basename(fileName)
        # read DIMENSION
        i = findKey("DIMENSION", splitFile)
        if (i < 1)
            throw("readInstance(...): DIMENSION not found")
        end
        dim = parse(Int64, splitFile[i+2])
        instance.n = dim-1
        instance.nn = dim
        instance.coords = zeros(Int64, dim,2)
        instance.demand = zeros(Int64, dim)
        instance.c = zeros(Float64, dim, dim)
        # read VEHICLES
        # We ignore the number of vehicles specified in the file so this part is commented out
        #i = findKey("VEHICLES", splitFile)
        #=if (i >= 1)
            instance.nVeh = parse(Int32, splitFile[i+2])
        else
            instance.nVeh = 0
        end=#
        # read CAPACITY
        i = findKey("CAPACITY", splitFile)
        if (i < 1)
            throw("readInstance(...): CAPACITY not found")
        end
        instance.Q = parse(Int64, splitFile[i+2])
        # read NODE_COORD_SECTION
        i = findKey("NODE_COORD_SECTION", splitFile)
        if (i < 1)
            throw("readInstance(...): NODE_COORD_SECTION not found")
        end
        coords = zeros(Int64,dim,2)
        for j=1:dim
            index = 2+i+(j-1)*3
            x = parse(Int32, splitFile[index])
            y = parse(Int32, splitFile[index+1])
            coords[j,1] = x
            coords[j,2] = y
        end
        # move coordinates to the place we prefer (node 1:n are customers, node n+1 is the depot)
        instance.coords[dim,:] = coords[1,:]
        for j=1:dim-1
            instance.coords[j,:] = coords[j+1,:]
        end
        # compute distances
        for i = 1:dim
            for j = i:dim
                if roundDist
                    d = round(norm(instance.coords[i,1:2] - instance.coords[j,1:2]))
                else
                    d = norm(instance.coords[i,1:2] - instance.coords[j,1:2])
                end
                #println("Dist computed: $d")
                instance.c[i,j] = d
                instance.c[j,i] = d
            end
        end
        # read DEMAND_SECTION
        i = findKey("DEMAND_SECTION", splitFile)
        if (i < 1)
            throw("readInstance(...): DEMAND_SECTION not found")
        end
        demands = zeros(Int64,dim)
        for j=1:dim
            index = 2+i+(j-1)*2
            demands[j] = parse(Int32, splitFile[index])
        end
        # move demands to the place we prefer (node 1:n are customers, node n+1 is the depot)
        instance.demand[dim] = demands[1]
        for j=1:dim-1
            instance.demand[j] = demands[j+1]
        end
    end
    return instance
end

end

ext: jl


text:
# Data structures for model parameters
# Alex Tam, 13/10/2020

# Parameters
"Global numerical parameters"
@with_kw struct Numerical_Parameters
    nA::Int = 50 # [-] Number of actin filaments
    nM::Int = 10 # [-] Number of myosin motors
    nT::Int = 1201 # [-] Number of time steps
    dt::Float64 = 0.05 # [s] Time step size
    lxx::Float64 = 2.5 # [μm] Reference domain width (x)
    lyy::Float64 = 2.5 # [μm] Reference domain width (y)
    xTol::Float64 = 1e-8 # [-] DOF tolerance for optimisation
    fTol::Float64 = 1e-8 # [-] Objective function tolerance for optimisation
    gTol::Float64 = 1e-8 # [-] Gradient tolerance for optimisation
end

"Actin filament properties"
@with_kw struct Actin_Properties
    nSeg::Int = 5; # [-] Number of segments
    LSeg::Float64 = 0.2 # [μm] Equilibrium segment length
    k::Float64 = 1000 # [pN/μm] Spring constant
    lambda_a::Float64 = 0.05 # [pN/(μm^2)*s] Actin-background drag coefficient
    kappa::Float64 = 0.073 # [pN*μm^2] Flexural rigidity
    lambda_pf::Float64 = 30.0 # [pN/μm*s] Protein friction drag coefficient
    k_off::Float64 = 0.04 # [/filament/s] Turnover rate
    k_p::Float64 = 0.0 # [/filament/s] Polymerisation rate
    kb::Float64 = 1.380649e-5 # [μm*pN/K] Boltzmann constant
    T::Float64 = 298.15 # [K] Temperature
end

"Myosin motor properties"
@with_kw struct Myosin_Properties
    k::Float64 = 1000 # [pN/μm] Spring constant
    Fs::Float64 = 5 # [pN] Stall force
    Vm::Float64 = 0.5 # [μm/s] Load-free velocity
    k_off::Float64 = 0.35 # [/motor/s] Reference off-rate in the absence of force
    F_ref::Float64 = 12.6 # [pN] Reference force for computing force-dependent off-rate
end


text:
# ------------------------------------------------------------------------------------------
# # Introduction to DataFrames
# **[Bogumił Kamiński](http://bogumilkaminski.pl/about/), May 13, 2018**
# ------------------------------------------------------------------------------------------

using DataFrames

# ------------------------------------------------------------------------------------------
# ## Extras - selected functionalities of selected packages
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### FreqTables: creating cross tabulations
# ------------------------------------------------------------------------------------------

using FreqTables
df = DataFrame(a=rand('a':'d', 1000), b=rand(["x", "y", "z"], 1000))
ft = freqtable(df, :a, :b) # observe that dimensions are sorted if possible

ft[1,1], ft['b', "z"] # you can index the result using numbers or names

prop(ft, 1) # getting proportions - 1 means we want to calculate them in rows (first dimension)

prop(ft, 2) # and columns are normalized to 1.0 now

x = categorical(rand(1:3, 10))
levels!(x, [3, 1, 2, 4]) # reordering levels and adding an extra level
freqtable(x) # order is preserved and not-used level is shown

freqtable([1,1,2,3,missing]) # by default missings are listed

freqtable([1,1,2,3,missing], skipmissing=true) # but we can skip them

# ------------------------------------------------------------------------------------------
# ### DataFramesMeta - working on `DataFrame`
# ------------------------------------------------------------------------------------------

using DataFramesMeta
df = DataFrame(x=1:8, y='a':'h', z=repeat([true,false], outer=4))

@with(df, :x+:z) # expressions with columns of DataFrame

@with df begin # you can define code blocks
    a = :x[:z]
    b = :x[.!:z]
    :y + [a; b]
end

a # @with creates hard scope so variables do not leak out

df2 = DataFrame(a = [:a, :b, :c])
@with(df2, :a .== ^(:a)) # sometimes we want to work on raw Symbol, ^() escapes it

df2 = DataFrame(x=1:3, y=4:6, z=7:9)
@with(df2, _I_(2:3)) # _I_(expression) is translated to df2[expression]

@where(df, :x .< 4, :z .== true) # very useful macro for filtering

@select(df, :x, y = 2*:x, z=:y) # create a new DataFrame based on the old one

@transform(df, a=1, x = 2*:x, y=:x) # create a new DataFrame adding columns based on the old one

@transform(df, a=1, b=:a) # old DataFrame is used and :a is not present there

@orderby(df, :z, -:x) # sorting into a new data frame, less powerful than sort, but lightweight

@linq df |> # chaining of operations on DataFrame
    where(:x .< 5) |>
    orderby(:z) |>
    transform(x²=:x.^2) |>
    select(:z, :x, :x²)

f(df, col) = df[col] # you can define your own functions and put them in the chain
@linq df |> where(:x .<= 4) |> f(:x)

# ------------------------------------------------------------------------------------------
# ### DataFramesMeta - working on grouped `DataFrame`
# ------------------------------------------------------------------------------------------

df = DataFrame(a = 1:12, b = repeat('a':'d', outer=3))
g = groupby(df, :b)

@by(df, :b, first=first(:a), last=last(:a), mean=mean(:a)) # more convinient than by from DataFrames

@based_on(g, first=first(:a), last=last(:a), mean=mean(:a)) # the same as by but on grouped DataFrame

@where(g, mean(:a) > 6.5) # filter gropus on aggregate conditions

@orderby(g, -sum(:a)) # order groups on aggregate conditions

@transform(g, center = mean(:a), centered = :a - mean(:a)) # perform operations within a group and return ungroped DataFrame

DataFrame(g) # a nice convinience function not defined in DataFrames

@transform(g) # actually this is the same

@linq df |> groupby(:b) |> where(mean(:a) > 6.5) |> DataFrame # you can do chaining on grouped DataFrames as well

# ------------------------------------------------------------------------------------------
# ### DataFramesMeta - rowwise operations on `DataFrame`
# ------------------------------------------------------------------------------------------

df = DataFrame(a = 1:12, b = repeat(1:4, outer=3))

# such conditions are often needed but are complex to write
@transform(df, x = ifelse.((:a .> 6) .& (:b .== 4), "yes", "no"))

# one option is to use a function that works on a single observation and broadcast it
myfun(a, b) = a > 6 && b == 4 ? "yes" : "no"
@transform(df, x = myfun.(:a, :b))

# or you can use @byrow! macro that allows you to process DataFrame rowwise
@byrow! df begin
    @newcol x::Vector{String}
    :x = :a > 6 && :b == 4 ? "yes" : "no"
end

# ------------------------------------------------------------------------------------------
# ### Visualizing data with StatPlots
# ------------------------------------------------------------------------------------------

using StatPlots # you might need to setup Plots package and some plotting backend first

# we present only a minimal functionality of the package

srand(1)
df = DataFrame(x = sort(randn(1000)), y=randn(1000), z = [fill("b", 500); fill("a", 500)])

@df df plot(:x, :y, legend=:topleft, label="y(x)") # a most basic plot

@df df density(:x, label="") # density plot

@df df histogram(:y, label="y") # and a histogram

@df df boxplot(:z, :x, label="x")

@df df violin(:z, :y, label="y")

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# When we are a child, we do the following:
# 1. Set up a temp directory that we can write to.
# 2. Copy the contents of our package from `/app` (which is read-only) to our temp
#    directory (which is read-and-write).
# 3. Delete any manifest files in the temp directory,.
# 4. Instantiate the environment inside our temp directory.
original_app_directory = "/app" # this is read-only; we cannot write to this directory
new_app_directory = mktempdir(; cleanup = true) # we have write access to this directory
for element in readdir(original_app_directory)
    old_path = joinpath(original_app_directory, element)
    new_path = joinpath(new_app_directory,      element)
    cp(old_path, new_path; force = true)
end
rm.(joinpath.(Ref(new_app_directory), Base.manifest_names); force = true)
using Pkg
Pkg.activate(new_app_directory)
Pkg.instantiate()
Pkg.precompile()

# Load Sandbox, then try to launch a nested sandbox
using Sandbox, Test

rootfs_dir = Sandbox.alpine_rootfs()

config = SandboxConfig(
    # This rootfs was downloaded within the sandbox in the `Pkg.instantiate()` above
    Dict("/" => rootfs_dir),
    # Propagate our readwrite mounting into the nested sandbox
    Dict{String,String}("/tmp/readwrite" => "/tmp/readwrite"),
    persist=true,
    verbose=true,
)

# Prove that we can write into the `readwrite` location
open("/tmp/readwrite/single_nested.txt", "w") do io
    println(io, "aperture")
end

# This should always default to the unprivileged executor, since if we're nested, `FORCE_SANDBOX_MODE` should be set
with_executor() do exe
    @test success(exe, config, `/bin/sh -c "echo science > /tmp/readwrite/double_nested.txt"`)
end
@test String(read("/tmp/readwrite/double_nested.txt")) == "science\n"


text:
module SPIRV

using CEnum
using MLStyle
using UnPack
using Graphs
using Reexport
using Dictionaries
using AutoHashEquals

const Optional{T} = Union{Nothing,T}

import Base: write, show, showerror, ==

const magic_number = 0x07230203
const ID = UInt32

# generated SPIR-V wrapper
include("generated/enums.jl")
include("generated/instructions.jl")

include("utils.jl")
include("parse.jl")
include("disassemble.jl")
include("ssa.jl")
include("functions.jl")
include("spir_types.jl")
include("ir.jl")
include("assemble.jl")

# include("StructuredCFG/StructuredCFG.jl")
# @reexport using .StructuredCFG

export
        # parse
        PhysicalInstruction, PhysicalModule,
        Instruction,

        # disassemble
        disassemble,

        # assemble
        assemble,

        # SPIR-V types
        SPIRType,
        VoidType,
        ScalarType, BooleanType, IntegerType, FloatType,
        VectorType, MatrixType,
        ImageType,
        SamplerType, SampledImageType,
        ArrayType,
        OpaqueType,
        StructType,
        PointerType,

        # IR
        IR,
        ID,
        Metadata,
        SSADict,
        @inst
end


text:
using Random
using Soss
using DataFrames
using StatsPlots
using BenchmarkTools

@inline function rate_to_proportion(r::Float64,t::Float64)
    1-exp(-r*t)
end;

sir_markov = @model state,p begin
    # Unpack parameters
    β = p.β
    c = p.c
    γ = p.γ
    δt = p.δt

    # Unpack starting counts
    t0 = state.t
    S0 = state.S
    I0 = state.I
    R0 = state.R
    N = S0 + I0 + R0

    # Transitions between states
    S_I ~ Binomial(S0, rate_to_proportion(β*c*I0/N,δt))
    I_R ~ Binomial(I0, rate_to_proportion(γ,δt))

    # Updated state
    t = t0 + δt
    S = S0 - S_I
    I = I0 + S_I - I_R
    R = R0 + I_R

    next = (p=p, state=(t=t,S=S,I=I,R=R))
end;

sir_model = @model u0,p begin
    x ~ MarkovChain(p, sir_markov(state=u0,p=p))
end;

δt = 0.1
nsteps = 400
tmax = nsteps*δt;

u0 = (t=0.0, S=990, I=10, R=0); # t,S,I,R

p = (β=0.05, c=10.0, γ=0.25, δt=δt);

Random.seed!(1234);

r = rand(sir_model(u0=u0,p=p));
data = [u0]
for (n,s) in enumerate(r.x)
    n>nsteps && break
    push!(data,s)
end;

df_markov = DataFrame(data);

@df df_markov plot(:t,
    [:S :I :R],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark begin
    r = rand(sir_model(u0=u0,p=p));
    data = [u0]
    for (n,s) in enumerate(r.x)
        n>nsteps && break
        push!(data,s)
    end
end


text:
#############
# Exp & Log #
#############

struct Exp{N} <: Bijector{N} end
struct Log{N} <: Bijector{N} end

Exp() = Exp{0}()
Log() = Log{0}()

(b::Exp)(y::Real) = exp(y)
(b::Log)(x::Real) = log(x)

(b::Exp)(y) = exp.(y)
(b::Log)(x) = log.(x)

inv(b::Exp{N}) where {N} = Log{N}()
inv(b::Log{N}) where {N} = Exp{N}()

logabsdetjac(b::Exp{0}, x::Real) = x
logabsdetjac(b::Exp{0}, x::AbstractVector) = x
logabsdetjac(b::Exp{1}, x::AbstractVector) = sum(x)
logabsdetjac(b::Exp{1}, x::AbstractMatrix) = vec(sum(x; dims = 1))

logabsdetjac(b::Log{0}, x::Real) = -log(x)
logabsdetjac(b::Log{0}, x::AbstractVector) = -log.(x)
logabsdetjac(b::Log{1}, x::AbstractVector) = - sum(log.(x))
logabsdetjac(b::Log{1}, x::AbstractMatrix) = - vec(sum(log.(x); dims = 1))


text:
using BackwardAutodiff
using Test: @test, @testset

@testset "FowardAutodiff.jl tests" begin
    @testset "Basic arithmetic" begin
        x = randn()
        y = randn()
        z = randn()

        @test isapprox(D(x->2*x)(x), 2)
    end

    @testset "Transcendental Functions" begin
        x = randn()
        @test isapprox(D(exp)(x), exp(x))
        @test isapprox(D(sin)(x), cos(x))
        @test isapprox(D(cos)(x), -sin(x))
        @test isapprox(D(tan)(x), sec(x)*sec(x))
    end

    @testset "sqrt" begin
        x = randn()

        # We need to be in the complex domain explicitly, or sqrt(x < 0) will error out
        @test isapprox(D(sqrt)(x + 0im), 1/(2*sqrt(x + 0im)))
    end

    @testset "compound function bits" begin
        x = randn()
        y = randn()
        @test isapprox(D(x -> x*y)(x), y)
        @test isapprox(D(x -> x/y)(x), 1/y)
        @test isapprox(D(y -> x/y)(y), -x/(y*y))
        @test isapprox(D(x -> exp(x/y))(x), exp(x/y)/y)
        @test isapprox(D(x -> sin(x*y))(x), y*cos(x*y))
        @test isapprox(D(x -> cos(x*y))(x), -y*sin(x*y))
        @test isapprox(D(x -> x-y)(x), 1)
        @test isapprox(D(y -> x-y)(y), -1)
        @test isapprox(D(y -> x*y)(y), x)
        @test isapprox(D(x -> -x)(x), -1)
        @test isapprox(D(x -> exp(-x/pi))(x), -exp(-x/pi)/pi)
        @test isapprox(D(x -> sin(2*pi*sqrt(2)*x))(x), 2*pi*sqrt(2)*cos(2*pi*sqrt(2)*x))
        @test isapprox(D(x -> exp(-x/pi))(3.0), -exp(-3/pi)/pi)
        @test isapprox(D(x -> exp(-x)*sin(x))(x), -exp(-x)*sin(x) + exp(-x)*cos(x))
        @test isapprox(D(x -> exp(-x))(x), -exp(-x))
        @test isapprox(D(x -> x*x)(x), 2*x)
    end

    @testset "compound function" begin
        freq = sqrt(2)
        A = randn()
        B = randn()

        a = A/sqrt(A*A + B*B)
        b = B/sqrt(A*A + B*B)
        phi = atan(b, a)

        x = randn()

        function f(x)
            return exp(-x/pi)*(b*cos(2*pi*freq*x) + a*sin(2*pi*freq*x))
        end
        fprime = D(f)

        function laborious_fprime(x)
            return exp(-x/pi)*(2*pi*freq*cos(2*pi*freq*x + phi) - sin(2*pi*freq*x + phi)/pi)
        end

        @test isapprox(fprime(x), laborious_fprime(x))
    end

    @testset "partials of simple function" begin
        function linear(x, m, b)
            return m*x+b
        end
        dldx = D(1, linear)
        dldslope = D(2, linear)
        dldintercept = D(3, linear)

        x = randn()
        m = randn()
        b = randn()
        @test isapprox(dldx(x,m,b), m)
        @test isapprox(dldslope(x,m,b), x)
        @test isapprox(dldintercept(x,m,b), 1)
    end

    @testset "gradient of many-arg and vector function" begin
        function linear(x, m, b)
            m*x + b
        end

        function vlinear(xmb)
            linear(xmb...)
        end

        g = D(linear)
        vg = D(vlinear)

        x = randn()
        m = randn()
        b = randn()

        @test isapprox(g(x, m, b), [m, x, 1])
        @test isapprox(vg([x,m,b]), [m, x, 1])
    end

    @testset "multivariate polynomial" begin
        coeffs = randn(5)
        function p(x,y,z)
            coeffs[1] + x*(coeffs[2] + y*(coeffs[3] + z*(coeffs[4] + x*coeffs[5])))
        end

        function gp(x,y,z)
            [coeffs[2] + y*(coeffs[3] + z*(coeffs[4] + x*coeffs[5])) + x*y*z*coeffs[5],
             x*(coeffs[3] + z*(coeffs[4] + x*coeffs[5])),
             x*y*(coeffs[4] + x*coeffs[5])]
         end

         x,y,z = randn(3)

         @test isapprox(D(p)(x,y,z), gp(x,y,z))
    end

    @testset "log derivative" begin
        x = rand()
        @test isapprox(D(log)(x), 1/x)

        @test isapprox(D(x -> log(exp(x)))(x), 1.0)
    end
end


text:
function fullInteraction(data, numVars, bernO)
m = bernO+1;
MATRIX = data;
designMATRIX = zeros(size(data,1),1);
if numVars>1
    for j = 1:1:numVars-1
        starti = size(designMATRIX,2);
        for i =1:1:m^j
            designMATRIX = hcat(designMATRIX, MATRIX[:,i].*data[:,j*m+1:(j+1)*m]);
        end
        designMATRIX = designMATRIX[:,2:end]
        if j<numVars-1
            MATRIX=designMATRIX;
            designMATRIX= zeros(size(data,1),1);
        end
    end
else
    designMATRIX = data;
end

return designMATRIX
end


text:
using Pkg

Pkg.activate("coveragetempenv", shared=true)

Pkg.add(PackageSpec(name="CoverageTools"))

using CoverageTools

directories = get(ENV, "INPUT_DIRECTORIES", "src")
dirs = filter!(!isempty, split(directories, ","))
for dir in dirs
    isdir(dir) || error("directory \"$dir\" not found!")
end

pfs = mapreduce(process_folder, vcat, dirs)
LCOV.writefile("lcov.info", pfs)

max_stars_repo_name: garrison/julia-processcoverage

max_issues_repo_name: garrison/julia-processcoverage

max_forks_repo_name: garrison/julia-processcoverage


text:
load("utils/utils.jl")

macro evalcatch(ex)
    quote; let
        println($string(ex))
        try
            value = $ex
            println("  success; returned ", value)
        catch err
            println("  threw \"", err, '"')
        end
        println()
    end; end
end

@evalcatch @expect true
@evalcatch @expect false

@evalcatch @expect is_expr(:(x+y), :call)
@evalcatch @expect is_expr(1, :call)
@evalcatch @expect is_expr(1)

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
using MAT, LaTeXStrings, Plots, ProgressBars, Statistics, StatsBase, StatsPlots, Ranges
cd(dirname(@__FILE__))
file = matopen("Q2_data.mat")
vars = matread("Q2_data.mat")
####Δt=10ms / L=20
histo_list=[]
for i in 1:100
    for j in 1:500
        if vars["trials"][i,j]==1
            push!(histo_list, j/2000 - 0.05)
        end
    end
end

histogram(histo_list,bins=25,linecolor = :transparent,xlabel=L"Time(ms)", ylabel=L"rate(Hz)", title=L"PETH\ for\ \Delta t=10ms", label=nothing, dpi=400, color=:turquoise3)
plot!([0,0], [0,700], color=:red, label="Spike Time")
savefig("../../Fundamentals-of-Neuroscience_Fall1400/PETCH_RasterPlot/Figs/PETH_10ms.pdf")
####Δt=2.5 / L=5
histogram(histo_list, bins=100,linecolor = :transparent, xlabel=L"Time(ms)", ylabel=L"rate(Hz)", title=L"PETH\ for\ \Delta t=10ms", label=nothing, dpi=400, color=:turquoise3)
plot!([0,0], [0,400], color=:red, label="Spike Time")
savefig("../../Fundamentals-of-Neuroscience_Fall1400/PETCH_RasterPlot/Figs/PETH_5L.pdf")
####Δt=17.5 / L=35
histogram(histo_list, bins=14,linecolor = :transparent, xlabel=L"Time(ms)", ylabel=L"rate(Hz)", title=L"PETH\ for\ \Delta t=10ms", label=nothing, dpi=400, color=:turquoise3)
plot!([0,0], [0,1000], color=:red, label="Spike Time")
savefig("../../Fundamentals-of-Neuroscience_Fall1400/PETCH_RasterPlot/Figs/PETH_17.5ms.pdf")
close(file)


text:
using DataStructures

ex = "inputs/11-15/13ex.txt"
function parseinput(file)
  lines = readlines(file)
  b = findall(x->x=="", lines)[1]
  nums = map(str -> Tuple(parse(Int,str[l]) for l in findall(r"[[:digit:]]+", str)), lines[1:b-1])
  fs = []
  for (i,l) in enumerate(lines[b+1:end])
    val = parse(Int,String(split(l,'=')[end]))
    d,val = split(l,'=')[end-1][end], parse(Int,String(split(l,'=')[end]))
    push!(fs, (val, d=='y' ? :y : :x))
  end
  (Set(nums), fs)
end

dist(v, split) = 2split - v

function fold(S, loc, axis=:x)
  for (x,y) in S
    if axis == :x && x > loc
      push!(S,(dist(x,loc),y))
      delete!(S,(x,y))
    elseif axis == :y && y > loc
      push!(S,(x,dist(y,loc)))
      delete!(S,(x,y))
    end
  end
  S
end

function display(G; xline=-1, yline=-1)
  w = maximum(i[1] for i in G) + 1
  h = maximum(i[2] for i in G) + 1
  p = [[' ' for _ in 1:w] for _ in 1:h]
  for (x,y) in G
    p[y+1][x+1] = '#'
  end
  if xline != -1
    p[xline+1] = ['-' for _ in 1:w]
  end
  for l in p
    println(join(l))
  end
end

function day13(S::Set{Tuple{Int64, Int64}}, instructions)
  for (val, axis) in instructions
    S = fold(S, val, axis)
  end
  S
end

main() = display(day13(parseinput("inputs/11-15/13v2.txt")...))

using BenchmarkTools

nums, instructions = parseinput("inputs/11-15/13v2.txt")
@btime day13(nums, instructions)

main()


text:
"""
Converts a dictionary parsed by PowerModels to a System.
Currently Supports MATPOWER and PSSE data files parsed by PowerModels.
Supports kwargs to supply formatters for different device types,
such as `bus_name_formatter` or `gen_name_formatter`.

# Examples
```julia
sys = PSY.pm2ps_dict(pm_data, configpath = "ACTIVSg25k_validation.json",
                    bus_name_formatter = x->string(x["name"]*"-"*string(x["index"])),
                    load_name_formatter = x->strip(join(x["source_id"], "_")))
```
"""
function pm2ps_dict(data::Dict{String, Any}; kwargs...)
    if length(data["bus"]) < 1
        throw(DataFormatError("There are no buses in this file."))
    end

    @info "Constructing System from Power Models" data["name"] data["source_type"]

    sys = System(data["baseMVA"]; kwargs...)

    bus_number_to_bus = read_bus!(sys, data; kwargs...)
    read_loads!(sys, data, bus_number_to_bus; kwargs...)
    read_loadzones!(sys, data, bus_number_to_bus; kwargs...)
    read_gen!(sys, data, bus_number_to_bus; kwargs...)
    read_branch!(sys, data, bus_number_to_bus; kwargs...)
    read_shunt!(sys, data, bus_number_to_bus; kwargs...)
    read_dcline!(sys, data, bus_number_to_bus; kwargs...)

    check!(sys)
    return sys
end

"""
Internal component name retreval from pm2ps_dict
"""
function _get_pm_dict_name(device_dict)
    return get(device_dict, "name", string(device_dict["index"]))
end

"""
Creates a PowerSystems.Bus from a PowerSystems bus dictionary
"""
function make_bus(bus_dict::Dict{String, Any})
    bus = Bus(bus_dict["number"],
                     bus_dict["name"],
                     bus_dict["bustype"],
                     bus_dict["angle"],
                     bus_dict["voltage"],
                     bus_dict["voltagelimits"],
                     bus_dict["basevoltage"]
                     )
     return bus
 end

function make_bus(bus_name, bus_number, d, bus_types)
    bus = make_bus(Dict{String, Any}("name" => bus_name ,
                            "number" => bus_number,
                            "bustype" => bus_types[d["bus_type"]],
                            "angle" => d["va"],
                            "voltage" => d["vm"],
                            "voltagelimits" => (min=d["vmin"], max=d["vmax"]),
                            "basevoltage" => d["base_kv"]
                            ))
    return bus
end

# "From http://www.pserc.cornell.edu/matpower/MATPOWER-manual.pdf Table B-1"
@enum MatpowerBusType begin
    MATPOWER_PQ = 1
    MATPOWER_PV = 2
    MATPOWER_REF = 3
    MATPOWER_ISOLATED = 4
end

function Base.convert(::Type{BusType}, x::MatpowerBusType)
    map = Dict(MATPOWER_ISOLATED => ISOLATED,
               MATPOWER_PQ => PQ,
               MATPOWER_PV => PV,
               MATPOWER_REF => REF)
    return map[x]
end

function read_bus!(sys::System, data; kwargs...)
    @info "Reading bus data"
    bus_number_to_bus = Dict{Int, Bus}()

    bus_types = instances(MatpowerBusType)
    data = sort(collect(data["bus"]), by = x->parse(Int64, x[1]))

    if length(data) == 0
        @error "No bus data found" # TODO : need for a model without a bus
    end

    _get_name = get(kwargs, :bus_name_formatter, _get_pm_dict_name)

    for (i, (d_key, d)) in enumerate(data)
        # d id the data dict for each bus
        # d_key is bus key
        d["name"] = get(d, "name", string(d["bus_i"]))
        bus_name = _get_name(d)
        bus_number = Int(d["bus_i"])
        bus = make_bus(bus_name, bus_number, d, bus_types)
        bus_number_to_bus[bus.number] = bus
        add_component!(sys, bus; skip_validation=SKIP_PM_VALIDATION)
    end

    return bus_number_to_bus
end

function make_load(d, bus; kwargs...)
    _get_name = get(kwargs, :load_name_formatter, x->strip(join(x["source_id"])))
    return PowerLoad(;
        name= _get_name(d),
        available=true,
        model = ConstantPower::LoadModel,
        bus=bus,
        activepower=d["pd"],
        reactivepower=d["qd"],
        maxactivepower=d["pd"],
        maxreactivepower=d["qd"],
    )
end

function read_loads!(sys::System, data, bus_number_to_bus::Dict{Int, Bus}; kwargs...)
    if !haskey(data, "load")
        @error "There are no loads in this file"
        return
    end

    for d_key in keys(data["load"])
        d = data["load"][d_key]
        if d["pd"] != 0.0
            bus = bus_number_to_bus[d["load_bus"]]
            load = make_load(d, bus; kwargs...)

            add_component!(sys, load; skip_validation=SKIP_PM_VALIDATION)
        end
    end
end

function make_loadzones(d, bus_l, activepower, reactivepower; kwargs...)
    _get_name = get(kwargs, :loadzone_name_formatter, _get_pm_dict_name)
    return LoadZones(;
        number=d["index"],
        name=_get_name(d),
        buses=bus_l,
        maxactivepower=sum(activepower),
        maxreactivepower=sum(reactivepower),
    )
end

function read_loadzones!(sys::System, data, bus_number_to_bus::Dict{Int, Bus}; kwargs...)
    if !haskey(data, "areas")
        @info "There are no Load Zones data in this file"
        return
    end

    for (d_key, d) in data["areas"]
        buses = [bus_number_to_bus[b["bus_i"]] for (b_key, b) in data["bus"]
                 if b["area"] == d["index"]]
        bus_names = Set{String}()
        for bus in buses
            push!(bus_names, get_name(bus))
        end

        active_power = Vector{Float64}()
        reactive_power = Vector{Float64}()

        for (key, load) in data["load"]
            load_bus = bus_number_to_bus[load["load_bus"]]
            if get_name(load_bus) in bus_names
                push!(active_power, load["pd"])
                push!(reactive_power, load["qd"])
            end
        end

        d["name"] = get(d, "name", d_key)

        load_zones = make_loadzones(d, buses, active_power, reactive_power; kwargs...)
        add_component!(sys, load_zones; skip_validation=SKIP_PM_VALIDATION)
    end
end

function make_hydro_gen(gen_name, d, bus)
    ramp_agc = get(d, "ramp_agc", get(d, "ramp_10", get(d, "ramp_30", d["pmax"])))
    tech = TechHydro(;
        rating=calculate_rating(d["pmax"], d["qmax"]),
        primemover=convert(PrimeMovers, d["type"]),
        activepowerlimits=(min=d["pmin"], max=d["pmax"]),
        reactivepowerlimits=(min=d["qmin"], max=d["qmax"]),
        ramplimits=(up=ramp_agc / d["mbase"], down=ramp_agc / d["mbase"]),
        timelimits=nothing,
    )

    curtailcost = TwoPartCost(0.0, 0.0)

    return HydroDispatch(name = gen_name,
                         available = Bool(d["gen_status"]),
                         bus = bus,
                         activepower = d["pg"],
                         reactivepower = d["qg"],
                         tech = tech,
                         op_cost = curtailcost)
end

function make_tech_renewable(d)
    tech = TechRenewable(;
        rating=float(d["pmax"]),
        primemover=convert(PrimeMovers, d["type"]),
        reactivepowerlimits=(min=d["qmin"], max=d["qmax"]),
        powerfactor=1.0,
    )

    return tech
end

function make_renewable_dispatch(gen_name, d, bus)
    tech = make_tech_renewable(d)
    cost = TwoPartCost(0.0, 0.0)
    generator = RenewableDispatch(;
        name=gen_name,
        available=Bool(d["gen_status"]),
        bus=bus,
        activepower = d["pg"],
        reactivepower = d["qg"],
        tech=tech,
        op_cost=cost,
    )

    return generator
end

function make_renewable_fix(gen_name, d, bus)
    tech = make_tech_renewable(d)
    generator = RenewableFix(;
        name=gen_name,
        available=Bool(d["gen_status"]),
        bus=bus,
        activepower = d["pg"],
        reactivepower = d["qg"],
        tech=tech,
    )

    return generator
end

function make_generic_battery(gen_name, d, bus)

    # TODO: placeholder
    #battery=GenericBattery(;
    #    name=gen_name,
    #    available=Bool(d["gen_status"]),
    #    bus=bus,
    #    energy=,
    #    capacity=,
    #    rating=,
    #    activepower=,
    #    inputactivepowerlimits=,
    #    outputactivepowerlimits=,
    #    efficiency=,
    #    reactivepower=,
    #    reactivepowerlimits=,
    #)
    #return battery
end

"""
The polynomial term follows the convention that for an n-degree polynomial, at least n + 1 components are needed.
    c(p) = c_n*p^n+...+c_1p+c_0
    c_o is stored in the  field in of the Econ Struct
"""
function make_thermal_gen(gen_name::AbstractString, d::Dict, bus::Bus)
    if haskey(d, "model")
        model = GeneratorCostModel(d["model"])
        if model == PIECEWISE_LINEAR::GeneratorCostModel
            cost_component = d["cost"]
            power_p = [i for (ix, i) in enumerate(cost_component) if isodd(ix)]
            cost_p =  [i for (ix, i) in enumerate(cost_component) if iseven(ix)]
            cost = [(p, c) for (p, c) in zip(cost_p, power_p)]
            fixed = max(0.0, cost[1][1] - (cost[2][1] - cost[1][1])/(cost[2][2] - cost[1][2])*cost[1][2])
            cost = [(c[1] - fixed, c[2]) for c in cost]
        elseif model == POLYNOMIAL::GeneratorCostModel
            if d["ncost"] == 0
                cost = (0.0, 0.0)
                fixed = 0.0
            elseif d["ncost"] == 1
                cost = (0.0, 0.0)
                fixed = d["cost"][1]
            elseif d["ncost"] == 2
                cost = (0.0, d["cost"][1])
                fixed = d["cost"][2]
            elseif d["ncost"] == 3
                cost = (d["cost"][1], d["cost"][2])
                fixed = d["cost"][3]
            else
                throw(DataFormatError("invalid value for ncost: $(d["ncost"]). PowerSystems only supports polynomials up to second degree"))
            end
        end
        startup = d["startup"]
        shutdn = d["shutdown"]
    else
        @warn "Generator cost data not included for Generator: $gen_name"
        tmpcost = ThreePartCost(nothing)
        cost = tmpcost.variable
        fixed = tmpcost.fixed
        startup = tmpcost.startup
        shutdn = tmpcost.shutdn
    end

    # TODO GitHub #148: ramp_agc isn't always present. This value may not be correct.
    ramp_agc = get(d, "ramp_agc", get(d, "ramp_10", get(d, "ramp_30", d["pmax"])))

    tech = TechThermal(;
        rating = sqrt(d["pmax"]^2 + d["qmax"]^2),
        primemover = convert(PrimeMovers, d["type"]),
        fuel = convert(ThermalFuels, d["fuel"]),
        activepowerlimits = (min = d["pmin"], max = d["pmax"]),
        reactivepowerlimits = (min = d["qmin"], max = d["qmax"]),
        ramplimits = (up = ramp_agc / d["mbase"], down = ramp_agc / d["mbase"]),
        timelimits = nothing,
    )
    op_cost = ThreePartCost(;
        variable=cost,
        fixed=fixed,
        startup=startup,
        shutdn=shutdn,
    )

    thermal_gen = ThermalStandard(
        name=gen_name,
        available=Bool(d["gen_status"]),
        bus=bus,
        activepower = d["pg"],
        reactivepower = d["qg"],
        tech=tech,
        op_cost=op_cost,
    )

    return thermal_gen
end

"""
Transfer generators to ps_dict according to their classification
"""
function read_gen!(sys::System, data, bus_number_to_bus::Dict{Int, Bus}; kwargs...)
    @info "Reading generator data"

    if !haskey(data, "gen")
        @error "There are no Generators in this file"
        return nothing
    end

    genmap_file = get(kwargs, :genmap_file, nothing)
    genmap = get_generator_mapping(genmap_file)

    for (name, pm_gen) in data["gen"]
        if haskey(kwargs, :gen_name_formatter)
            _get_name = kwargs[:gen_name_formatter]
        elseif haskey(pm_gen, "name")
            _get_name = _get_pm_dict_name
        elseif haskey(pm_gen, "source_id")
            _get_name = d -> strip(string(d["source_id"][1]) * "-" * string(d["source_id"][2]) * "-" * string(d["index"]))
        end

        gen_name = _get_name(pm_gen)

        bus = bus_number_to_bus[pm_gen["gen_bus"]]
        pm_gen["fuel"] = get(pm_gen, "fuel", "OTHER")
        pm_gen["type"] = get(pm_gen, "type", "OT")
        @debug "Found generator" gen_name bus pm_gen["fuel"] pm_gen["type"]

        gen_type = get_generator_type(pm_gen["fuel"], pm_gen["type"], genmap)
        if gen_type == ThermalStandard
            generator = make_thermal_gen(gen_name, pm_gen, bus)
        elseif gen_type == HydroDispatch
            generator = make_hydro_gen(gen_name, pm_gen, bus)
        elseif gen_type == RenewableDispatch
            generator = make_renewable_dispatch(gen_name, pm_gen, bus)
        elseif gen_type == RenewableFix
            generator = make_renewable_fix(gen_name, pm_gen, bus)
        elseif gen_type == GenericBattery
            @warn "Skipping GenericBattery"
            continue
            # TODO
            #generator = make_generic_battery(gen_name, pm_gen, bus)
        else
            @error "Skipping unsupported generator" gen_type
            continue
        end

        add_component!(sys, generator; skip_validation=SKIP_PM_VALIDATION)
    end
end

function make_branch(name, d, bus_f, bus_t)
    primary_shunt = d["b_fr"]
    alpha = d["shift"]
    branch_type = get_branch_type(d["tap"], alpha)

    if d["transformer"]
        if branch_type == Line
            throw(DataFormatError("Data is mismatched; this cannot be a line. $d"))
        elseif branch_type == Transformer2W
            value = make_transformer_2w(name, d, bus_f, bus_t)
        elseif branch_type == TapTransformer
            value = make_tap_transformer(name, d, bus_f, bus_t)
        elseif branch_type == PhaseShiftingTransformer
            value = make_phase_shifting_transformer(name, d, bus_f, bus_t, alpha)
        else
            error("Unsupported branch type $branch_type")
        end
    else
        # The get_branch_type() logic doesn't work for this data.
        # tap can be 1.0 for this data.
        value = make_line(name, d, bus_f, bus_t)
    end

    return value
end

function make_line(name, d, bus_f, bus_t)
    pf = get(d, "pf", 0.0)
    qf = get(d, "qf", 0.0)

    return Line(;
        name=name,
        available=Bool(d["br_status"]),
        activepower_flow = pf,
        reactivepower_flow = qf,
        arc=Arc(bus_f, bus_t),
        r=d["br_r"],
        x=d["br_x"],
        b=(from=d["b_fr"], to=d["b_to"]),
        rate=d["rate_a"],
        anglelimits=(min=d["angmin"], max=d["angmax"]),
    )
end

function make_transformer_2w(name, d, bus_f, bus_t)
    pf = get(d, "pf", 0.0)
    qf = get(d, "qf", 0.0)
    return Transformer2W(;
        name=name,
        available=Bool(d["br_status"]),
        activepower_flow = pf,
        reactivepower_flow = qf,
        arc=Arc(bus_f, bus_t),
        r=d["br_r"],
        x=d["br_x"],
        primaryshunt=d["b_fr"],  # TODO: which b ??
        rate=d["rate_a"],
    )
end

function make_tap_transformer(name, d, bus_f, bus_t)
    pf = get(d, "pf", 0.0)
    qf = get(d, "qf", 0.0)
    return TapTransformer(;
        name=name,
        available=Bool(d["br_status"]),
        activepower_flow = pf,
        reactivepower_flow = qf,
        arc=Arc(bus_f, bus_t),
        r=d["br_r"],
        x=d["br_x"],
        tap=d["tap"],
        primaryshunt=d["b_fr"],  # TODO: which b ??
        rate=d["rate_a"],
    )
end

function make_phase_shifting_transformer(name, d, bus_f, bus_t, alpha)
    pf = get(d, "pf", 0.0)
    qf = get(d, "qf", 0.0)
    return PhaseShiftingTransformer(;
        name=name,
        available=Bool(d["br_status"]),
        activepower_flow = pf,
        reactivepower_flow = qf,
        arc=Arc(bus_f, bus_t),
        r=d["br_r"],
        x=d["br_x"],
        tap=d["tap"],
        primaryshunt=d["b_fr"],  # TODO: which b ??
        α=alpha,
        rate=d["rate_a"],
    )
end

function read_branch!(sys::System, data, bus_number_to_bus::Dict{Int, Bus}; kwargs...)
    @info "Reading branch data"
    if !haskey(data, "branch")
        @info "There is no Branch data in this file"
        return
    end

    _get_name = get(kwargs, :branch_name_formatter, _get_pm_dict_name)

    for (d_key, d) in data["branch"]
        d["name"] = get(d, "name", d_key)
        name = _get_name(d)
        bus_f = bus_number_to_bus[d["f_bus"]]
        bus_t = bus_number_to_bus[d["t_bus"]]
        value = make_branch(name, d, bus_f, bus_t)

        add_component!(sys, value; skip_validation=SKIP_PM_VALIDATION)
    end
end

function make_dcline(name, d, bus_f, bus_t)
    return HVDCLine(;
        name=name,
        available=Bool(d["br_status"]),
        activepower_flow = get(d, "pf", 0.0),
        arc=Arc(bus_f, bus_t),
        activepowerlimits_from=(min=d["pminf"] , max=d["pmaxf"]),
        activepowerlimits_to=(min=d["pmint"], max=d["pmaxt"]),
        reactivepowerlimits_from=(min=d["qminf"], max=d["qmaxf"]),
        reactivepowerlimits_to=(min=d["qmint"], max=d["qmaxt"]),
        loss=(l0=d["loss0"], l1 =d["loss1"]),
    )
end

function read_dcline!(sys::System, data, bus_number_to_bus::Dict{Int, Bus}; kwargs...)
    @info "Reading DC Line data"
    if !haskey(data, "dcline")
        @info "There is no DClines data in this file"
        return
    end

    _get_name = get(kwargs, :branch_name_formatter, _get_pm_dict_name)

    for (d_key, d) in data["dcline"]
        d["name"] = get(d, "name", d_key)
        name = _get_name(d)
        bus_f = bus_number_to_bus[d["f_bus"]]
        bus_t = bus_number_to_bus[d["t_bus"]]

        dcline = make_dcline(name, d, bus_f, bus_t)
        add_component!(sys, dcline, skip_validation=SKIP_PM_VALIDATION)
    end
end

function make_shunt(name, d, bus)
    return FixedAdmittance(;
        name=name,
        available=Bool(d["status"]),
        bus=bus,
        Y=(-d["gs"] + d["bs"]im),
    )
end

function read_shunt!(sys::System, data, bus_number_to_bus::Dict{Int, Bus}; kwargs...)
    @info "Reading branch data"
    if !haskey(data, "shunt")
        @info "There is no shunt data in this file"
        return
    end

    _get_name = get(kwargs, :shunt_name_formatter, _get_pm_dict_name)

    for (d_key, d) in data["shunt"]
        d["name"] = get(d, "name", d_key)
        name = _get_name(d)
        bus = bus_number_to_bus[d["shunt_bus"]]
        shunt = make_shunt(name, d, bus)

        add_component!(sys, shunt; skip_validation=SKIP_PM_VALIDATION)
    end
end


text:
using Snappy
using Random
using Test

@testset "Low-level Interfaces" begin
    SnappyOK = Snappy.SnappyOK
    SnappyInvalidInput = Snappy.SnappyInvalidInput
    SnappyBufferTooSmall = Snappy.SnappyBufferTooSmall

    function originals()
        ["", "\x00", "\x00\x00", "foo", "foobarbaz", "x"^50]
    end

    # Compress & uncompress small data using pre-allocated buffer
    buffer_size = 100
    for original in originals()
        # compress
        compressed = Array{UInt8}(undef, buffer_size)
        olen, st = Snappy.snappy_compress(Vector{UInt8}(original), compressed)
        @test st == SnappyOK
        resize!(compressed, olen)

        # uncompress
        uncompressed = Array{UInt8}(undef, buffer_size)
        olen, st = Snappy.snappy_uncompress(compressed, uncompressed)
        @test st == SnappyOK
        resize!(uncompressed, olen)

        restored = String(uncompressed)
        @test original == restored
    end

    # Prepare too small buffer and try to compress and uncompress data
    original = "orig"

    # prepare too small buffer to compress
    compressed = Array{UInt8}(undef, 1)
    _, st = Snappy.snappy_compress(Vector{UInt8}(original), compressed)
    @test st == SnappyBufferTooSmall

    # now enough buffer size
    compressed = Array{UInt8}(undef, 100)
    olen, st = Snappy.snappy_compress(Vector{UInt8}(original), compressed)
    @test st == SnappyOK
    resize!(compressed, olen)

    # again the prepared buffer is too small to uncompress
    uncompressed = Array{UInt8}(undef, 1)
    olen, st = Snappy.snappy_uncompress(compressed, uncompressed)
    @test st == SnappyBufferTooSmall

    # Break compressed data and detect it
    original = "orig"
    compressed = Array{UInt8}(undef, 100)
    olen, st = Snappy.snappy_compress(Vector{UInt8}(original), compressed)
    @test st == SnappyOK
    resize!(compressed, olen)
    @test Snappy.snappy_validate_compressed_buffer(compressed) == SnappyOK
    # perturb the compressed buffer
    compressed[1] = 0x00
    @test Snappy.snappy_validate_compressed_buffer(compressed) == SnappyInvalidInput

    # Estimate compressed size
    for original in originals()
        maxlen = Snappy.snappy_max_compressed_length(UInt(length(original)))
        compressed = Array{UInt8}(undef, 100)
        olen, st = Snappy.snappy_compress(Vector{UInt8}(original), compressed)
        @test st == SnappyOK
        #@show Int(olen), Int(maxlen)
        @test olen <= maxlen
    end

    # Estimate uncompressed size
    for original in originals()
        compressed = Array{UInt8}(undef, 100)
        olen, st = Snappy.snappy_compress(Vector{UInt8}(original), compressed)
        @test st == SnappyOK
        resize!(compressed, olen)
        olen, st = Snappy.snappy_uncompressed_length(compressed)
        @test st == SnappyOK
        @test olen == length(original)
    end
end

@testset "High-level Interfaces" begin
    # QuickCheck-like property satisfaction tests (compress ○ uncompress = id)
    Random.seed!(2014)

    # byte arrays
    randbytes(n) = rand(UInt8, n)
    for original in map(randbytes, 0:100:10000)
        @test uncompress(compress(original)) == original
    end

    # strings
    for original in map(randstring, 0:100:1000)
        @test String(uncompress(compress(Vector{UInt8}(original)))) == original
    end

    # Large data
    # 128MiB
    Random.seed!(2014)
    original = randstring(128 * 1024^2)
    original_bytes = Vector{UInt8}(original)
    @test uncompress(compress(original_bytes)) == original_bytes
end


text:
module AMLPipelineBase

include("abstracttypes.jl")
using .AbsTypes
export Machine, Computer, Workflow, Learner, Transformer
export fit!, transform!, fit_transform!

include("utils.jl")
using .Utils
using AMLPipelineBase.Utils
export holdout, kfold, score, infer_eltype, nested_dict_to_tuples,
       nested_dict_set!, nested_dict_merge, create_transformer,
       mergedict, getiris,
       skipmean,skipmedian,skipstd,
       aggregatorclskipmissing,
       getiris, getprofb

include("baselinemodels.jl")
using .BaselineModels
export Baseline, Identity

include("basefilters.jl")
using .BaseFilters
export Imputer, OneHotEncoder, Wrapper

include("featureselector.jl")
using .FeatureSelectors
export FeatureSelector, CatFeatureSelector, NumFeatureSelector, CatNumDiscriminator

include("decisiontree.jl")
using .DecisionTreeLearners
export PrunedTree, RandomForest, Adaboost

include("ensemble.jl")
using .EnsembleMethods
export VoteEnsemble, StackEnsemble, BestLearner

include("crossvalidator.jl")
using .CrossValidators
export crossvalidate

include("naremover.jl")
using .NARemovers
export NARemover

include("pipelines.jl")
using .Pipelines
export @pipeline, @pipelinex, @pipelinez
export Pipeline, ComboPipeline

end # module


text:
"""
genetrap.jl

Fit GRSM models to live cell and smFISH data
"""

"""
Gene information
"""
const genes = ["CANX";"DNAJC5";"ERRFI1";"KPNB1";"MYH9";"Rab7a";"RHOA";"RPAP3";"Sec16A";"SLC2A1"]
const genelength = Dict([("Sec16A", 42960);("SLC2A1", 33802);("ERRFI1", 14615);("RHOA", 52948);("KPNB1", 33730);("MYH9", 106741);("DNAJC5", 40930);("CANX", 32710);("Rab7a", 88663);("RPAP3", 44130)])
const MS2end = Dict([("Sec16A", 5220);("SLC2A1", 26001);("ERRFI1", 5324);("RHOA", 51109);("KPNB1", 24000);("MYH9", 71998);("DNAJC5", 14857);("CANX", 4861);("Rab7a", 83257);("RPAP3", 38610)])
const halflife = Dict([("CANX", 50.),("DNAJC5", 5.),("ERRFI1", 1.35),("KPNB1", 9.),("MYH9", 10.),("Rab7a", 50.),("RHOA", 50.),("RPAP3", 7.5),("Sec16A", 8.),("SLC2A1", 5.)])

"""
    fit_genetrap()

"""

function makeswarm_genetrap(sfile,juliafile,nchains,genes::Vector,G::Int,R::Int,nalleles::Int,root::String,folder::String,maxtime::Float64,samplesteps::Int)
write_swarmfile(sfile,nchains,juliafile,genes)
write_fitfile(juliafile,nchains,G,R,nalleles,root,folder,maxtime,samplesteps)
end

function write_fitfile(fitfile,nchains,G::Int,R::Int,nalleles::Int,root::String,folder::String,maxtime::Float64,samplesteps::Int)
        f = open(fitfile,"w")
        s =   '"'
        write(f,"@everywhere using StochasticGene\n")
        write(f,"@time fit_genetrap($nchains,ARGS[1],$G,$R,$nalleles,$s$root$s,$s$folder$s,$maxtime,$samplesteps)\n")
        close(f)
end

function fit_genetrap(nchains,gene::String,G::Int,R::Int,nalleles::Int,root::String,folder::String,maxtime::Float64,samplesteps::Int)
    println(now())
    data,model,options = genetrap(root,folder,"mean","gt",gene,G,R,nalleles,"",maxtime,samplesteps)
    print_ll(data,model)
    fit,stats,waic = run_mh(data,model,options,nchains)
    finalize(data,model,fit,stats,waic,1.,folder,root)
    println(now())
end

"""
genetrap()

Load data, model and option structures for metropolis-hastings fit
root is the folder containing data and results

FISH counts is divided by tempfish to adjust relative weights of data
set tempfish = 0 to equalize FISH and live cell counts
"""
function genetrap(root::String,infolder::String,rtype::String,label::String,gene::String,G::Int,R::Int,nalleles::Int,type::String,maxtime::Float64,samplesteps::Int,cyclesteps=0,annealsteps=0,warmupsteps=0,method=1,temp=1.,tempanneal=1.,tempfish=1.)
    # r = readrates_genetrap(root,infolder,rtype,gene,"$G$R",type)
    r = readrates_genetrap(joinpath(root,infolder),rtype,gene,label,G,R,nalleles,type)
    println(r)
    genetrap(root,r,label,gene,G,R,nalleles,type,maxtime,samplesteps,cyclesteps,annealsteps,warmupsteps,method,temp,tempanneal,tempfish)
end

function genetrap(root,r,label::String,gene::String,G::Int,R::Int,nalleles::Int,type::String,maxtime::Float64,samplesteps::Int,cyclesteps=0,annealsteps=0,warmupsteps=0,method=1,temp=1.,tempanneal=1.,tempfish=1.)
    data = iszero(R) || iszero(tempfish) ? data_genetrap_FISH(root,label,gene) : data_genetrap(root,label,gene,tempfish)
    model = model_genetrap(r,gene,G,R,nalleles,type,method)
    options = MHOptions(samplesteps,cyclesteps,warmupsteps,annealsteps,maxtime,temp,tempanneal)
    return data,model,options
end

function data_genetrap(root,label,gene,tempfish=1.)
    LC = readLCPDF_genetrap(root,gene)
    if tempfish == 0
        counts = Int(div(sum(LC[:,2]+LC[:,3]),2))
        println(counts)
        # set FISH counts to mean of live cell counts
        histFISH = readFISH_genetrap(root,gene,counts)
    else
        histFISH = readFISH_genetrap(root,gene,tempfish)
    end
    RNALiveCellData(label,gene,length(histFISH),histFISH,LC[:,1],LC[:,3],LC[:,2])
end

function data_genetrap_FISH(root,label,gene)
    histFISH = readFISH_genetrap(root,gene,1.)
    RNAData(label,gene,length(histFISH),histFISH)
end

"""
model_genetrap

load model structure
"""
function model_genetrap(r,gene::String,G::Int,R::Int,nalleles::Int,type::String,method,)
    if R == 0
        fittedparam = collect(1:2*G-1)
        decayrate = get_decay(halflife[gene])
        if r == 0
            r = setrate(G,1,decayrate,.1)[1]
        end
        d = prior_rna(r,G,1,fittedparam,decayrate,1.)
        return  GMmodel{typeof(r),typeof(d),Float64,typeof(fittedparam),Int}(G,nalleles,r,d,0.02,fittedparam,1)
    else
        propcv=[.02*ones(2*(G-1));1e-4;1e-4*ones(R);.02*ones(R)]
        fittedparam=collect(1:num_rates(G,R)-1)
        Gprior=(.01,.5)
        Sprior=(.1,.5)
        Rcv=1e-2
        rm,rcv = prior_rates_genetrap(G,R,gene,Gprior,Sprior,Rcv)
        if r == 0
            r = rm
        end
        # d = priordistributionLogNormal_genetrap(rm,rcv,G,R)
        d = priordistributionGamma_genetrap(rm,rcv,G,R)
        return GRSMmodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,R,nalleles,type,r,d,propcv,fittedparam,method)
    end
end

"""
setpriorrate(G,R,gene,Gprior::Tuple,Sprior::Tuple,Rcv::Float64)
Set prior distribution for mean and cv of rates
"""
function prior_rates_genetrap(G,R,gene,Gprior::Tuple,Sprior::Tuple,Rcv::Float64)
    n = G-1
    rm = [fill(Gprior[1],2*n);2000/(genelength[gene]-MS2end[gene]);fill(2000/MS2end[gene]*R,R);fill(Sprior[1], R);log(2.)/(60 .* halflife[gene])]
    rcv = [fill(Gprior[2],2*n);Rcv;Rcv;fill(Sprior[2], R);Rcv]
    return rm,rcv
end

"""
priordistributionLogNormal_genetrap(r,cv,G,R)

LogNormal Prior distribution
"""
function priordistributionLogNormal_genetrap(r,cv,G,R)
    sigma = sigmalognormal(cv)
    d = []
    j = 1
    #  G rates
    for i in Grange(G)
        push!(d,Distributions.LogNormal(log(r[i]),sigma[j]))
        j += 1
    end
    # initiation rate
    i = initiation(G)
    push!(d,Distributions.LogNormal(log(r[i]),sigma[j]))
    j += 1
    # Step rates are bounded by length of insertion site to end of gene, i.e sum of reciprocal rates is bounded
    t = sum(1 ./ r[Rrange(G,R)])
    push!(d,Distributions.LogNormal(-log(t),sigma[j]))
    j += 1
    # priors for splice rates
    rs = 0.
    for i in Srange(G,R)
        rs += r[i]
        push!(d,Distributions.LogNormal(log(rs),sigma[j]))
        j += 1
    end
    return d
end

"""
priordistributionGamma_genetrap(r,cv,G,R)

Gamma Prior distribution
"""
function priordistributionGamma_genetrap(rm,cv,G,R)
    n = G-1
    zet = R
    d = []
    rsig = cv .^ 2
    j = 1
    # priors for G rates
    for i in Grange(G)
        theta = rsig[j]*rm[i]
        alpha = 1/rsig[j]
        push!(d,Distributions.Gamma(alpha,theta))
        j += 1
    end
    # prior for nu1 is upper bounded by rm, i.e. distance is bounded by distance to insertion site
    theta = rsig[j]*rm[initiation(G)]
    alpha = 1/rsig[j]
    push!(d,Distributions.Gamma(alpha,theta))
    j += 1
    # priors for nu rates are bounded by length of insertion site to end of gene, i.e sum of reciprocal rates is bounded
    tm = sum(1 ./ rm[Rrange(G,R)])
    # for i in 2*n + 2 : 2*n + zet + 1
    #     tm += 1/rm[i]
    # end
    theta = rsig[j]/tm
    alpha = 1/rsig[j]
    push!(d,Distributions.Gamma(alpha,theta))
    j += 1
    # priors for ejection rates
    rs = 0
    for i in Srange(G,R)
        rs += rm[i]
        theta = rsig[j]*rs
        alpha = 1/rsig[j]
        push!(d,Distributions.Gamma(alpha,theta))
        j += 1
    end
    return d
end

"""
readLCPDF_genetrap(root,gene)

Read in dwell time PDF
"""
function readLCPDF_genetrap(root,gene)
    infile = joinpath(root,"DwellTimePDF/$(gene)_PDF.csv")
    if isfile(infile)
        LC = readdlm(infile,',')
        x = truncate_histogram(LC[:,2],.999,1000)
        LC[1:length(x),:]
    else
        println("No data for gene: ", gene)
    end
end

"""
readFISH_genetrap(root,gene,temp=1.,clone=true)

Read and assemble smFISH data
"""
function readFISH_genetrap(root::String,gene::String,temp::Float64=1.,clone=true)
    counts = countsFISH_genetrap(root,gene,clone)
    readFISH_genetrap(root,gene,Int(div(counts,temp)),clone)
end

function readFISH_genetrap(root::String,gene::String,counts::Int,clone=true)
    fishfile=joinpath(root,"Fish_4_Genes/$(gene)_steady_state_mRNA.csv")
    col = clone ? 3 : 2
    # Get smFISH RNA histograms
    if isfile(fishfile)
        x = readdlm(fishfile,',')[3:end,col]
        x = x[x .!= ""]
        x = truncate_histogram(x,.99,1000)
        x /= sum(x)
        x *= counts
        return x
    else
        println("No smFISH data for gene: ", gene)
    end
    nothing
end

"""
countsFISH_genetrap(root,gene,x,counts,clone=true)

read in smFISH counts
"""
function countsFISH_genetrap(root,gene,clone=true)
    countfile = joinpath(root,"counts/total_counts.csv")
    wttext = "WT_"
    clonetext = "_clone"
    if isfile(countfile)
        countsdata = readdlm(countfile,',')[:,:]
        if clone
            counts=(countsdata[findall(uppercase("$(gene)$clonetext") .== uppercase.(countsdata[:,1])),2])[1]
        else
            counts=(countsdata[findall("$wttext$(gene)" .== countsdata[:,1]),2])[1]
        end
    else
        counts = 1000
    end
    return counts
end

"""
readrates_genetrap(infolder::String,rtype::String,gene::String,model::String,type::String)

Read in initial rates from previous runs
"""
# function readrates_genetrap(infolder::String,rtype::String,gene::String,label,G,R,nalleles,type::String)
#     if rtype == "rml" || rtype == "rlast"
#         rskip = 1
#         rstart = 1
#         row = 1
#     elseif rtype == "rmean"
#         rskip = 2
#         rstart = 1
#     elseif rtype == "rquant"
#         rskip = 3
#         rstart = 2
#     end
#     if type == "offeject" || type == "on"
#         type = ""
#     end
#     infile = getratefile_genetrap(infolder,rtype,gene,label,G,R,nalleles,type)
#     print(gene," ","$G$R"," ",type)
#     readrates_genetrap(infile,rstart,rskip)
# end

function readrates_genetrap(infolder::String,rtype::String,gene::String,label,G,R,nalleles,type::String)
    if rtype == "ml"
        row = 1
    elseif rtype == "mean"
        row = 2
    elseif rtype == "median"
        row = 3
    elseif rtype == "last"
        row = 4
    else
        row = 2
    end
    if type == "offeject" || type == "on"
        type = ""
    end
    infile = getratefile_genetrap(infolder,rtype,gene,label,G,R,nalleles,type)
    println(gene," ","$G$R"," ",type)
    readrates_genetrap(infile,row)
end

function readrates_genetrap(infile::String,row::Int)
    if isfile(infile) && ~isempty(read(infile))
        return readrates(infile,row)
    else
        println(" no prior")
        return 0
    end
end

# function readrates_genetrap(infile::String,rstart::Int,rskip::Int)
#     if isfile(infile) && ~isempty(read(infile))
#         r = readdlm(infile)
#         r = r[:,rstart:rskip:end]
#         println(r)
#         return r
#     else
#         println(" no prior")
#         return 0
#     end
#     nothing
# end

function getratefile_genetrap(infolder::String,rtype::String,gene::String,label,G,R,nalleles,type::String)
    model = R == 0 ? "$G" : "$G$R"
    file = "rates" * "_" * label * "_" * gene * "_" * model * "_"  * "$(nalleles)" * ".txt"
    joinpath(infolder,file)
end

# function getratefolder_genetrap(root::String,infolder::String,rtype::String,gene::String,model::String,type::String)
#     results = joinpath(root,"results")
#     infolder = joinpath(infolder,gene)
#     infolder = joinpath(infolder,rtype * model * type * txtstr)
#     joinpath(results,infolder)
# end

max_stars_repo_name: nih-niddk-mbs/Transcription

max_issues_repo_name: nih-niddk-mbs/Transcription

max_forks_repo_name: nih-niddk-mbs/Transcription


text:
module MLJTuringInterface

export TuringModel

import MLJModelInterface

using MLJModelInterface:
    Probabilistic,
    Table,
    Continuous,
    metadata_model,
    metadata_pkg
using Statistics: mean
using Tables
using Turing:
    MCMCThreads,
    DynamicPPL,
    chainscat,
    group,
    predict,
    replacenames,
    sample

const MMI = MLJModelInterface
const PKG = "MLJTuringInterface"

include("interface.jl")
include("metadata.jl")

end # module


text:
@testset "Dinic" begin
    # Construct DiGraph
    flow_graph = lg.DiGraph(8)

    # Load custom dataset
    flow_edges = [
        (1, 2, 10), (1, 3, 5), (1, 4, 15), (2, 3, 4), (2, 5, 9),
        (2, 6, 15), (3, 4, 4), (3, 6, 8), (4, 7, 16), (5, 6, 15),
        (5, 8, 10), (6, 7, 15), (6, 8, 10), (7, 3, 6), (7, 8, 10)
    ]

    capacity_matrix = zeros(Int, lg.nv(flow_graph), lg.nv(flow_graph))

    for e in flow_edges
        u, v, f = e
        lg.add_edge!(flow_graph, u, v)
        capacity_matrix[u, v] = f
    end

    # Construct the residual graph
    for fg in (flow_graph, lg.DiGraph{UInt8}(flow_graph), lg.DiGraph{Int16}(flow_graph))
      residual_graph = @inferred(LightGraphsFlows.residual(fg))

      # Test with default distances
      @test @inferred(LightGraphsFlows.dinic_impl(residual_graph, 1, 8, LightGraphsFlows.DefaultCapacity(residual_graph)))[1] == 3

      # Test with capacity matrix
      @test @inferred(LightGraphsFlows.dinic_impl(residual_graph, 1, 8, capacity_matrix))[1] == 28

      # Test on disconnected graphs
      function test_blocking_flow(residual_graph, source, target, capacity_matrix, flow_matrix)
          #disconnect source
          h = copy(residual_graph)
          for dst in collect(lg.neighbors(residual_graph, source))
              lg.rem_edge!(h, source, dst)
          end
          @test @inferred(LightGraphsFlows.blocking_flow(h, source, target, capacity_matrix, flow_matrix)) == 0

          #disconnect target and add unreachable vertex
          h = copy(residual_graph)
          for src in collect(lg.inneighbors(residual_graph, target))
              lg.rem_edge!(h, src, target)
          end
          @test @inferred(LightGraphsFlows.blocking_flow(h, source, target, capacity_matrix, flow_matrix)) == 0

          # unreachable vertex (covers the case where a vertex isn't reachable from the source)
          h = copy(residual_graph)
          lg.add_vertex!(h)
          lg.add_edge!(h, lg.nv(residual_graph) + 1, target)
          capacity_matrix_ = vcat(hcat(capacity_matrix, zeros(Int, lg.nv(residual_graph))), zeros(Int, 1, lg.nv(residual_graph) + 1))
          flow_graph_  = vcat(hcat(flow_matrix, zeros(Int, lg.nv(residual_graph))), zeros(Int, 1, lg.nv(residual_graph) + 1))

          @test @inferred(LightGraphsFlows.blocking_flow(h, source, target, capacity_matrix_, flow_graph_)) > 0

          #test with connected graph
          @test @inferred(LightGraphsFlows.blocking_flow(residual_graph, source, target, capacity_matrix, flow_matrix)) > 0
      end

      flow_matrix = zeros(Int, lg.nv(residual_graph), lg.nv(residual_graph))
      test_blocking_flow(residual_graph, 1, 8, capacity_matrix, flow_matrix)
    end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to configure the MGEO.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

################################################################################
#                              Public Functions
################################################################################

"""
    conf_design_vars(T, Nv::Integer, bits::Integer, min::Number, max::Number)

This function configures `Nv` design variables in which all of them will have
the same number of bits `bits`, and maximum and minimum values `min` and `max`.

The parameter `T` selects which MGEO algorithm will be used and it can be
`MGEO_Canonical()` (**Default**) or `MGEO_Var()`.

    conf_design_vars(T::Val{:MGEO_Real}, Nv::Integer, min::Number, max::Number, σ::Number, perc_perturbation::Bool = true)

This function configures `Nv` design variables in which all of them will have
the same maximum and minimum values `min` and `max`, and the same standard
deviation `σ` used when modifying the variables. This is used for MGEO Real
algorithm.

If `perc_perturbation` is `true`, then the design variable will be perturbed
using the following equation:

    x = x + N(0, σ⋅x),

in which `x` is the design variable and `N(a,b)` is a sample from a Normal
distribution with mean `a` and standard deviation `b`.

# Returns

An array with the configured design variables.

"""
conf_design_vars(Nv::Integer, bits::Integer, min::Number, max::Number) =
    conf_design_vars(Val(:MGEO_Canonical), Nv, bits, min, max)

function conf_design_vars(T::Val{:MGEO_Canonical},
                          Nv::Integer,
                          bits::Integer,
                          min::Number,
                          max::Number)

    # Check if the number of design variables is greater than 1.
    ( Nv <= 0 ) && error("The number of design variables must be higher than 0.")

    # Check if the number of bits are greater than 0.
    ( bits <= 0 ) && error("The number of bits must be higher than 0.")

    # Check if the minimum value is smaller than maximum value.
    ( min >= max ) && error("The minimum value must be less than the maximum value for a variable.")

    # Create the array of design variables.
    design_vars = Vector{Design_Variable_MGEO_Canonical}(undef, Nv)

    # Full scale for each variable.
    fullscale = (1 << bits) - 1

    @inbounds for i=1:Nv
        design_vars[i] = Design_Variable_MGEO_Canonical(bits,
                                                        min,
                                                        max,
                                                        fullscale,
                                                        (i-1)*bits+1,
                                                        "Var. $i")
    end

    return SVector{Nv, Design_Variable_MGEO_Canonical}(design_vars)
end

function conf_design_vars(T::Val{:MGEO_Var},
                          Nv::Integer,
                          bits::Integer,
                          min::Number,
                          max::Number)

    # Check if the number of design variables is greater than 1.
    ( Nv <= 0 ) && error("The number of design variables must be higher than 0.")

    # Check if the number of bits are greater than 0.
    ( bits <= 0 ) && error("The number of bits must be higher than 0.")

    # Check if the minimum value is smaller than maximum value.
    ( min >= max ) && error("The minimum value must be less than the maximum value for a variable.")

    # Create the array of design variables.
    design_vars = Vector{Design_Variable_MGEO_Var}(undef, Nv)

    # Full scale for each variable.
    fullscale = (1 << bits) - 1

    @inbounds for i=1:Nv
        design_vars[i] = Design_Variable_MGEO_Var(bits,
                                                  min,
                                                  max,
                                                  fullscale,
                                                  (i-1)*bits+1,
                                                  "Var. $i")
    end

    return SVector{Nv, Design_Variable_MGEO_Var}(design_vars)
end

function conf_design_vars(T::Val{:MGEO_Real},
                          Nv::Integer,
                          min::Number,
                          max::Number,
                          σ::Number,
                          perc_perturbation::Bool = true)

    # Check if the number of design variables is greater than 1.
    ( Nv <= 0 ) && error("The number of design variables must be higher than 0.")

    # Check if the minimum value is smaller than maximum value.
    ( min >= max ) && error("The minimum value must be less than the maximum value for a variable.")

    # Check if standard deviation σ is higher than 0.
    ( σ <= 0 ) && error("The standard deviation must be higher than 0.")

    # Create the array of design variables.
    design_vars = Vector{Design_Variable_MGEO_Real}(undef, Nv)

    @inbounds for i=1:Nv
        design_vars[i] = Design_Variable_MGEO_Real(min, max, σ, "Var. $i",
                                                   perc_perturbation)
    end

    return SVector{Nv, Design_Variable_MGEO_Real}(design_vars)
end

"""
    conf_design_vars(bits::AbstractVector{T1}, min::AbstractVector{T2}, max::AbstractVector{T3}, var_names::AbstractVector{String}) where {T1<:Integer, T2<:Number, T3<:Number} =

This function configures the design variables specifying for each one the number
of bits `bits`, the minimum `min` and maximum `max` values, and the variable
names `var_names`.

The parameter `T` selects which MGEO algorithm will be used and it can be
`MGEO_Canonical()` (**Default**) or `MGEO_Var()`.

    conf_design_vars(T::Val{:MGEO_Real}, min::AbstractVector{T1}, max::AbstractVector{T2}, σ::AbstractVector{T3}, var_names::AbstractVector{String}, perc_perturbation::T4 = nothing) where {T1<:Number, T2<:Number, T3<:Number, T4<:Union{Nothing, Vector{Bool}}}

This function configures the design variables specifying for each one the
minimum `min` and maximum `max` values, the standard deviations `σ` used when
modifying the variable, and the variable names `var_names`.

If, for each variable, the element in the vector `perc_perturbation` is `true`,
then the design variable will be perturbed using the following equation:

    x = x + N(0, σ⋅x),

in which `x` is the design variable value and `N(a,b)` is a sample from a Normal
distribution with mean `a` and standard deviation `b`. If `perc_perturbation` is
`nothing`, then it set to `true` for all variables.

# Returns

An array with the configured design variables.

"""
conf_design_vars(bits::AbstractVector{T1},
                 min::AbstractVector{T2},
                 max::AbstractVector{T3},
                 var_names::AbstractVector{String}) where
    {T1<:Integer, T2<:Number, T3<:Number} =
    conf_design_vars(Val(:MGEO_Canonical), bits, min, max, var_names)

function conf_design_vars(T::Val{:MGEO_Canonical},
                          bits::AbstractVector{T1},
                          min::AbstractVector{T2},
                          max::AbstractVector{T3},
                          var_names::AbstractVector{String}) where
 	{T1<:Integer, T2<:Number, T3<:Number}

    # Check if the size of arrays is correct.
    Nv = length(bits)

    ( ( length(min)       != Nv ) ||
      ( length(max)       != Nv ) ||
      ( length(var_names) != Nv ) ) &&
      error("The size of the vectors does not match.")

    # Number of bits.
    num_bits = 0

    # Create the array of design variables.
    design_vars = Vector{Design_Variable_MGEO_Canonical}(undef, Nv)

    for i = 1:Nv
        # Check if the number of bits are larger than 0.
        ( bits[i] <= 0 ) && @error("The number of bits must be higher than 0.")

        # Check if the min value is smaller than max.
        ( min[i] >= max[i] ) && @error("The minimum value must be less than the maximum value for a variable.")

        design_vars[i] = Design_Variable_MGEO_Canonical(bits[i],
                                                        min[i],
                                                        max[i],
                                                        (UInt64(1) << bits[i]) - 1,
                                                        num_bits+1,
                                                        var_names[i])
        num_bits += bits[i]
    end

    return SVector{Nv, Design_Variable_MGEO_Canonical}(design_vars)
end

function conf_design_vars(T::Val{:MGEO_Var},
                          bits::AbstractVector{T1},
                          min::AbstractVector{T2},
                          max::AbstractVector{T3},
                          var_names::AbstractVector{String}) where
    {T1<:Integer, T2<:Number, T3<:Number}

    # Check if the size of arrays is correct.
    Nv = length(bits)

    ( ( length(min)       != Nv ) ||
      ( length(max)       != Nv ) ||
      ( length(var_names) != Nv ) ) &&
      error("The size of the vectors does not match.")

    # Number of bits.
    num_bits = 0

    # Create the array of design variables.
    design_vars = Vector{Design_Variable_MGEO_Var}(undef, Nv)

    for i = 1:Nv
        # Check if the number of bits are larger than 0.
        ( bits[i] <= 0 ) && @error("The number of bits must be higher than 0.")

        # Check if the min value is smaller than max.
        ( min[i] >= max[i] ) && @error("The minimum value must be less than the maximum value for a variable.")

        design_vars[i] = Design_Variable_MGEO_Var(bits[i],
                                                  min[i],
                                                  max[i],
                                                  (UInt64(1) << bits[i]) - 1,
                                                  num_bits+1,
                                                  var_names[i])
        num_bits += bits[i]
    end

    return SVector{Nv, Design_Variable_MGEO_Var}(design_vars)
end

function conf_design_vars(T::Val{:MGEO_Real},
                          min::AbstractVector{T1},
                          max::AbstractVector{T2},
                          σ::AbstractVector{T3},
                          var_names::AbstractVector{String},
                          perc_perturbation::T4 = nothing) where
    {T1<:Number, T2<:Number, T3<:Number, T4<:Union{Nothing, Vector{Bool}}}

    # Check if the size of arrays is correct.
    Nv = length(min)

    ( ( length(min)       != Nv ) ||
      ( length(max)       != Nv ) ||
      ( length(σ)         != Nv ) ||
      ( length(var_names) != Nv ) ||
      ( (perc_perturbation != nothing) && (length(perc_perturbation) != Nv) ) ) &&
      error("The size of the vectors does not match.")

    # Create the array of design variables.
    design_vars = Vector{Design_Variable_MGEO_Real}(undef, Nv)

    for i = 1:Nv
        # Check if the min value is smaller than max.
        ( min[i] >= max[i] ) && error("The minimum value must be less than the maximum value for a variable.")

        # Check if standard deviation σ is higher than 0.
        ( σ[i] <= 0 ) && error("The standard deviation must be higher than 0.")

        perc_perturbation_i =
            (perc_perturbation != nothing) ? perc_perturbation[i] : true

        design_vars[i] = Design_Variable_MGEO_Real(min[i],
                                                   max[i],
                                                   σ[i],
                                                   var_names[i],
                                                   perc_perturbation_i)
    end

    return SVector{Nv, Design_Variable_MGEO_Real}(design_vars)
end

"""
    conf_design_vars(T, bits::AbstractVector{T1}, min::AbstractVector{T2}, max::AbstractVector{T3}) where {T1<:Integer, T2<:Number, T3<:Number}

This function configures the design variables specifying for each one the number
of bits `bits`, and the minimum `min` and maximum `max` values. The variable
names will be selected as `Var. 1`, `Var. 2`, etc.

The parameter `T` selects which MGEO algorithm will be used and it can be
`MGEO_Canonical()` (**Default**) or `MGEO_Var()`.

    conf_design_vars(T::Val{:MGEO_Real}, min::AbstractVector{T1}, max::AbstractVector{T2}, σ::AbstractVector{T3}, perc_perturbation::T4 = nothing)

This function configures the design variables specifying for each one the
minimum `min` and maximum `max` values, and the standard deviations `σ` used
when modifying the variable. The variable names will be selected as `Var. 1`,
`Var. 2`, etc.

If `perc_perturbation` is `true`, then the design variable will be perturbed
using the following equation:

    x = x + N(0, σ⋅x),

in which `x` is the design variable and `N(a,b)` is a sample from a Normal
distribution with mean `a` and standard deviation `b`.

# Returns

An array with the configured design variables.

"""
conf_design_vars(bits::AbstractVector{T1},
                 min::AbstractVector{T2},
                 max::AbstractVector{T3}) where
    {T1<:Integer, T2<:Number, T3<:Number} =
    conf_design_vars(Val(:MGEO_Canonical), bits, min, max)

function conf_design_vars(T::Union{Val{:MGEO_Canonical},
                                   Val{:MGEO_Var}},
                          bits::AbstractVector{T1},
                          min::AbstractVector{T2},
                          max::AbstractVector{T3}) where
    {T1<:Integer, T2<:Number, T3<:Number}

    # Create an array with variable names.
    var_names = Array{String}(undef, length(bits))

    @inbounds for i=1:length(var_names)
        var_names[i] = "Var. $i"
    end

    return conf_design_vars(T, bits, min, max, var_names)
end

function conf_design_vars(T::Val{:MGEO_Real},
                          min::AbstractVector{T1},
                          max::AbstractVector{T2},
                          σ::AbstractVector{T3},
                          perc_perturbation::T4 = nothing) where
    {T1<:Number, T2<:Number, T3<:Number, T4<:Union{Nothing, Vector{Bool}}}

    # Create an array with variable names.
    var_names = Array{String}(undef, length(min))

    @inbounds for i=1:length(var_names)
        var_names[i] = "Var. $i"
    end

    return conf_design_vars(T, min, max, σ, var_names, perc_perturbation)
end

"""
    conf_mgeo(Nf::Int64, τ::Float64, ngen_max::Int64, run_max::Int64, design_vars::SVector{Nv, T}, mgeo_eps::Float64 = 1e-10) where {Nv,T}

Configure the MGEO.

The MGEO algorithm that will be used will be the same configured in
`design_vars`.

# Args

* `Nf`: Number of objective functions.
* `τ`: Parameter τ of MGEO.
* `ngen_max`: Maximum number of generations.
* `run_max`: Maximum number of independent runs.
* `design_vars`: Configuration of the design variables (see `conf_design_vars`).
* `mgeo_eps`: (OPTIONAL) Epsilon to compare objective functions (**Default**:
              `1e-10`)

# Returns

An instance of the structure `MGEO_Structure` with the MGEO configuration.

"""
function conf_mgeo(Nf::Integer,
                   τ::Number,
                   ngen_max::Integer,
                   run_max::Integer,
                   design_vars::SVector{Nv, T},
                   mgeo_eps::Number = 1e-10) where {Nv,T}

    # Check if the number of objective functions is higher than 1.
    ( Nf <= 0 ) && error("The number of objective functions must be higher than 0.")

    return MGEO_Structure{Nv, Nf, T}(τ, ngen_max, run_max, design_vars, mgeo_eps)
end


text:
export NeuralNetworkApproximator, ActorCritic, GaussianNetwork, DuelingNetwork

using Flux
using Random
using Distributions: Normal, logpdf
import Functors: functor
using MacroTools: @forward

"""
    NeuralNetworkApproximator(;kwargs)

Use a DNN model for value estimation.

# Keyword arguments

- `model`, a Flux based DNN model.
- `optimizer=nothing`
"""
Base.@kwdef struct NeuralNetworkApproximator{M,O} <: AbstractApproximator
    model::M
    optimizer::O = nothing
end

# some model may accept multiple inputs
(app::NeuralNetworkApproximator)(args...; kwargs...) = app.model(args...; kwargs...)

@forward NeuralNetworkApproximator.model Flux.testmode!,
Flux.trainmode!,
Flux.params,
device

functor(x::NeuralNetworkApproximator) =
    (model = x.model,), y -> NeuralNetworkApproximator(y.model, x.optimizer)

RLBase.update!(app::NeuralNetworkApproximator, gs) =
    Flux.Optimise.update!(app.optimizer, params(app), gs)

Base.copyto!(dest::NeuralNetworkApproximator, src::NeuralNetworkApproximator) =
    Flux.loadparams!(dest.model, params(src))

#####
# ActorCritic
#####

"""
    ActorCritic(;actor, critic, optimizer=ADAM())

The `actor` part must return logits (*Do not use softmax in the last layer!*), and the `critic` part must return a state value.
"""
Base.@kwdef struct ActorCritic{A,C,O} <: AbstractApproximator
    actor::A
    critic::C
    optimizer::O = ADAM()
end

functor(x::ActorCritic) =
    (actor = x.actor, critic = x.critic), y -> ActorCritic(y.actor, y.critic, x.optimizer)

RLBase.update!(app::ActorCritic, gs) = Flux.Optimise.update!(app.optimizer, params(app), gs)

function Base.copyto!(dest::ActorCritic, src::ActorCritic)
    Flux.loadparams!(dest.actor, params(src.actor))
    Flux.loadparams!(dest.critic, params(src.critic))
end

#####
# GaussianNetwork
#####

"""
    GaussianNetwork(;pre=identity, μ, logσ)

Returns `μ` and `logσ` when called.
Create a distribution to sample from
using `Normal.(μ, exp.(logσ))`.
"""
Base.@kwdef struct GaussianNetwork{P,U,S}
    pre::P = identity
    μ::U
    logσ::S
end

Flux.@functor GaussianNetwork

"""
This function is compatible with a multidimensional action space. When outputting an action, it uses `tanh` to normalize it.

- `rng::AbstractRNG=Random.GLOBAL_RNG`
- `is_sampling::Bool=false`, whether to sample from the obtained normal distribution.
- `is_return_log_prob::Bool=false`, whether to calculate the conditional probability of getting actions in the given state.
"""
function (model::GaussianNetwork)(rng::AbstractRNG, state; is_sampling::Bool=false, is_return_log_prob::Bool=false)
    x = model.pre(state)
    μ, logσ = model.μ(x), model.logσ(x)
    if is_sampling
        π_dist = Normal.(μ, exp.(logσ))
        z = rand.(rng, π_dist)
        if is_return_log_prob
            logp_π = sum(logpdf.(π_dist, z) .- (2.0f0 .* (log(2.0f0) .- z .- softplus.(-2.0f0 .* z))), dims = 1)
            return tanh.(z), logp_π
        else
            return tanh.(z)
        end
    else
        return μ, logσ
    end
end

function (model::GaussianNetwork)(state; is_sampling::Bool=false, is_return_log_prob::Bool=false)
    model(Random.GLOBAL_RNG, state; is_sampling=is_sampling, is_return_log_prob=is_return_log_prob)
end

#####
# DuelingNetwork
#####

"""
    DuelingNetwork(;base, val, adv)

Dueling network automatically produces separate estimates of the state value function network and advantage function network. The expected output size of val is 1, and adv is the size of the action space.
"""
struct DuelingNetwork{B,V,A}
    base::B
    val::V
    adv::A
end

Flux.@functor DuelingNetwork

function (m::DuelingNetwork)(state)
    x = m.base(state)
    val = m.val(x)
    return val .+ m.adv(x) .- mean(m.adv(x), dims=1)
end