text:
# geothermal energy simulation; see http://fvt.simulkade.com/posts/2017-04-07-heat-transport-in-porous-media.html
# written by Ali A. Eftekhari, while on paternity leave!
# 14 June 2017; Lyngby public library
# output:
# pressure gradient (pressure difference between the injection and the production well)
# outlet temperature: temperature at the middle of the right boundary
# vs simulation time
#

# You can install CoolProp by:
# Pkg.clone("https://github.com/vimalaad/CoolProp.jl")
# Pkg.build("CoolProp")
using CoolProp, DataFrames, PyPlot, Polynomials, JFVM
import GR
# a table of density values
depth_reservoir = 2500         # [m]
p_res = depth_reservoir/10*1e5 # [Pa]
T_res = 80 + 273.15            # [K]
T_inj = 40 + 273.15            # [K]
T     = collect(linspace(T_inj, T_res, 10))
rho_water = zeros(length(T))
mu_water  = zeros(length(T))
for i in eachindex(T)
    rho_water[i] = PropsSI("D", "T", T[i], "P", p_res, "water")           # [kg/m^3]
    mu_water[i]  = PropsSI("viscosity", "T", T[i], "P", p_res, "water")   # [Pa.s]
end

mu_fit  = polyfit(T, mu_water, 2)
rho_fit = polyfit(T, rho_water, 2)
# figure()
# subplot(2,1,1)
# scatter(T, mu_water, label = "μ data")
# plot(T, polyval(mu_fit, T))
# xlabel("T [K]")
# ylabel("Viscosity [Pa.s]")
# legend()
# axis([minimum(T), maximum(T), minimum(mu_water), maximum(mu_water)])
#
# subplot(2,1,2)
# scatter(T, rho_water, label = "ρ data")
# plot(T, polyval(rho_fit, T))
# xlabel("T [K]")
# ylabel("Density [kg/m^3]")
# legend()

flow_m3_h = 100        # [m3/h]
flow_inj  = flow_m3_h/3600 # m3/s
perm_D  = 0.02        # Darcy
perm    = perm_D*1e-12 # [m^2]
L_res    = 500            # [m] will be changed later in a loop
L_left = 500           # [m] distance between the left boundary and injection well
L_right = 500          # [m] distance between the right boundary and injection well
L = L_res + L_left + L_right # [m] domain length
W     = 500            # [m]
thick_res = 100        # [m]
poros = 0.20           # [-]
T0    = 25.0 + 273.15  # [K] reference temperature
V_dp  = 0.2            # Dykstra-Parsons coef
clx   = 0.1            # Correlation length in x direction
cly   = 0.1            # Correlation length in y direction
T_inj = 40 + 273.15    # [K]
theta_inj = T_inj - T0 # [K]
T_res = 80 + 273.15    # [K]
theta_res = T_res - T0 # [K]
λ_water  = PropsSI("conductivity", "T", T_res, "P", p_res, "water") # W/m/K
λ_rock   = 0.16        # [W/m/K] rock conductivity (for sandstone from my thesis)
λ_eff    = λ_water^poros*λ_rock^(1-poros) # effective conductivity
cp_water = PropsSI("CPMASS", "T", T_res, "P", p_res, "water") # J/kg/K
cp_rock  = 837.0       # [J/kg/K] rock heat capacity (for sandstone)
rho_rock = 2650.0      # [kg/m^3] rock density (for sandstone)

well_diameter = 7*0.0254 # [m]; converted from a well diameter of 7 inch
nx_left = 50
nx_right = 50
nx_well           = 50
ny_well           = 50
nx                = 30
ny                = 30
x_range =
  collect([linspace(0, L_left-(nx_well+1)*well_diameter, nx_left);
           linspace(L_left-nx_well*well_diameter, L_left+nx_well*well_diameter, 2*nx_well);
           linspace(L_left+(nx_well+1)*well_diameter, L_left+L_res-(nx_well+1)*well_diameter, nx);
           linspace(L_left+L_res-nx_well*well_diameter, L_left+L_res+nx_well*well_diameter, 2*nx_well);
           linspace(L_left+L_res+(nx_well+1)*well_diameter, L, nx_right)])
y_range = collect([linspace(0, 0.5*W-(ny_well+1)*well_diameter, ny);
                   linspace(0.5*W-ny_well*well_diameter, 0.5*W+ny_well*well_diameter, 2*ny_well);
                   linspace(0.5*W+(ny_well+1)*well_diameter, W, ny)])
ind_inj  = [nx_left+nx_well+1, ny+ny_well]
ind_prod = [nx_left+nx_well+nx_well+nx+nx_well, ny+ny_well]
m   = createMesh2D(x_range, y_range)  # 2D Cartesian grid
v_cell = thick_res*cellVolume(m)
Nx  = length(x_range)-1
Ny  = length(y_range)-1

perm_val   = perm # permfieldlogrnde(Nx, Ny, perm, V_dp, clx, cly)
perm_field = createCellVariable(m, perm_val)

left_range  = ny+ny_well:ny+ny_well
right_range = ny+ny_well:ny+ny_well
u_inj = flow_inj/(m.cellsize.y[ny+ny_well+1]*thick_res)
BCp = createBC(m)                 # pressure boundary
BCp.left.a[:] = 0.0
BCp.left.b[:] = 1.0
BCp.left.c[:] = p_res
BCp.right.a[:] = 0.0
BCp.right.b[:] = 1.0
BCp.right.c[:] = p_res

# BCp.left.a[left_range]   =
#     perm_field.value[2, left_range]/polyval(mu_fit, T_inj)
# BCp.left.b[left_range]   = 0.0
# BCp.left.c[left_range]   = -u_inj
# BCp.right.a[right_range] = 0.0
# BCp.right.b[right_range] = 1.0
# BCp.right.c[right_range] = p_res

BCt = createBC(m)                 # temperature boundary
# Danckwertz with assumptions (density of water is almost constant)
rho_water_inj  = polyval(rho_fit, T_inj)
# BCt.left.a[left_range]  = λ_eff
# BCt.left.b[left_range]  = -rho_water_inj*cp_water*u_inj
# BCt.left.c[left_range]  = -rho_water_inj*cp_water*u_inj*theta_inj

# figure(figsize=(8,2))
# visualizeCells(perm_field)
# title("Permeability [m^2]", fontsize = 10)
# colorbar()

theta_init = createCellVariable(m, theta_res, BCt)
theta_val  = createCellVariable(m, theta_res, BCt)
p_init     = createCellVariable(m, p_res, BCp)
p_val      = createCellVariable(m, p_res, BCp)
rho_init   = createCellVariable(m,
    polyval(rho_fit, theta_init.value+T0))
rho_val    = copyCell(rho_init)
mu_init    = createCellVariable(m,
    polyval(mu_fit, theta_init.value+T0))
mu_val     = copyCell(mu_init)

dt_init    = (L*W*thick_res*poros)/flow_inj/50 # (L/Nx)/(u_inj/poros)   # [s] time step
final_time = 100*dt_init              # [s]
t_step     = collect([linspace(0, 5*dt_init, 30);
                      linspace(5.5*dt_init, final_time, 70)])
n_steps = length(t_step) # number of time steps
# variables to be stored:
dp_res = zeros(n_steps)
T_out  = zeros(n_steps)
p_inj  = zeros(n_steps)
p_out  = zeros(n_steps)
T_out[1] = T_res

# discretization
M_BCp, RHS_BCp = boundaryConditionTerm(BCp)
M_BCt, RHS_BCt = boundaryConditionTerm(BCt)
M_conductivity = diffusionTerm(-harmonicMean(
        createCellVariable(m, λ_eff)))
q_inj = createCellVariable(m, 0)
q_inj.value[ind_inj+1 ...] = flow_inj/v_cell.value[ind_inj+1 ...]
RHSinjt = constantSourceTerm(rho_water_inj*cp_water*(T_inj-T0)*q_inj)

q_prod = createCellVariable(m, 0)
q_prod.value[ind_prod+1 ...] = -flow_inj/v_cell.value[ind_prod+1 ...]

q_cont = createCellVariable(m, 0)

  # solve pressure equation
  water_mobil = harmonicMean(perm_field./mu_val)
  rho_face = arithmeticMean(rho_val)
  M_diff_p    = diffusionTerm(-rho_face.*water_mobil)
  q_cont.value[ind_inj+1 ...] = rho_val.value[ind_inj+1 ...]*flow_inj/v_cell.value[ind_inj+1 ...]
  q_cont.value[ind_prod+1 ...] = -rho_val.value[ind_prod+1 ...]*flow_inj/v_cell.value[ind_prod+1 ...]
  RHScont = constantSourceTerm(q_cont)
  p_val       = solveLinearPDE(m, M_diff_p + M_BCp,
      RHS_BCp + RHScont)

  # velocity vector
  u = -water_mobil.*gradientTerm(p_val)
figure()
visualizeCells(p_val)
figure()
plot(m.cellcenters.x[ind_inj[1]:ind_prod[1]], p_val.value[ind_inj[1]+1:ind_prod[1]+1, ind_inj[2]+1], ".")
p_inj  = p_val.value[ind_inj+1 ...]
p_out  = p_val.value[ind_prod+1 ...]
dp_res = p_val.value[ind_inj+1 ...]-p_val.value[ind_prod+1 ...]
println("p_inj= $p_inj, p_out= $p_out, dp = $dp_res")

df1 = DataFrame(x_m = m.cellcenters.x[ind_inj[1]:ind_prod[1]], p_pa = p_val.value[ind_inj[1]+1:ind_prod[1]+1, ind_inj[2]+1])
writetable("p_profile_local_refine", df1)

max_stars_repo_name: simulkade/geothermal

max_issues_repo_name: simulkade/geothermal

max_forks_repo_name: simulkade/geothermal


text:
using YaoBase
using Test

# mocked registers
struct TestRegister{B} <: AbstractRegister{B} end

TestRegister() = TestRegister{1}()

YaoBase.nqubits(::TestRegister) = 8
YaoBase.nactive(::TestRegister) = 2

export TestInterfaceRegister
struct TestInterfaceRegister{B} <: AbstractRegister{B} end

TestInterfaceRegister() = TestInterfaceRegister{1}()

@testset "Test general interface" begin
    @test_throws MethodError nactive(TestInterfaceRegister())
    @test_throws MethodError nqubits(TestInterfaceRegister())
    @test_throws MethodError nremain(TestInterfaceRegister())
end

@testset "Test default fallbacks" begin
    @test nbatch(TestInterfaceRegister()) == 1
end

@testset "adjoint register" begin
    @test adjoint(TestRegister()) isa AdjointRegister
    @test adjoint(adjoint(TestRegister())) isa TestRegister
    @test nqubits(adjoint(TestRegister())) == 8
    @test nactive(adjoint(TestRegister())) == 2
end

@testset "Test Printing" begin
    @test repr(TestRegister()) == """
    TestRegister{1}
        active qubits: 2/8"""
end

@testset "Test adjoint printing" begin
    @test repr(adjoint(TestRegister())) == """
        adjoint(TestRegister{1})
            active qubits: 2/8"""
end

ext: jl


text:
# License for module ModiaMath.Logging: MIT
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module ModiaMath.Logger

Log model evaluations.

# Main developer
Martin Otter, [DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
"""
module Logging

@static if VERSION >= v"0.7.0-DEV.2005" @eval using Printf end

export Logger, setLog!, setAllLogCategories!, setLogCategories!
export isLogStatistics, isLogProgress, isLogInfos, isLogWarnings, isLogEvents
export logOn!, logOff!, setLogCategories
export SimulationStatistics, reInitializeStatistics!, set_nResultsForSimulationStatistics!

"""
    mutable struct Logger - Log model evaluations
"""
mutable struct Logger
    log::Bool           # = true, if logging

    # log categories
    statistics::Bool
    progress::Bool
    infos::Bool
    warnings::Bool
    events::Bool

    Logger() = new(false, true, true, true, true, true)
end

function setLog!(logger::Logger, log::Bool)
    logger.log = log
    return nothing
end

function setAllLogCategories!(logger::Logger; default=true)
    logger.statistics = default
    logger.progress   = default
    logger.infos      = default
    logger.warnings   = default
    logger.events     = default
end

"""
    ModiaMath.setLogCategories!(obj, categories; reinit=true)

Set log categories on obj (of type ModiaMath.SimulationState, ModiaMath.AbstractSimulationModel,
or ModiaMath.Logger) as vector of symbols, e.g. setLogCategories!(simulationModel, [:LogProgess]).
Supported categories:
- `:LogStatistics`, print statistics information at end of simulation
- `:LogProgress`, print progress information during simulation
- `:LogInfos`, print information messages of the model
- `:LogWarnings`, print warning messages of the model
- `:LogEvents`, log events of the model

If option reinit=true, all previous set categories are reinitialized to be no longer present.
If reinit=false, previously set categories are not changed.
"""
function setLogCategories!(logger::Logger, categories::Vector{Symbol}; reinit=true)
    if reinit
        setAllLogCategories!(logger;default=false)
    end

    for c in categories
        if c == :LogStatistics
            logger.statistics = true
        elseif c == :LogProgress
            logger.progress = true
        elseif c == :LogInfos
            logger.infos = true
        elseif c == :LogWarnings
            logger.warnings = true
        elseif c == :LogEvents
            logger.events = true
        else
            warning("Log categorie ", c, " not known (will be ignored)")
        end
    end
    return nothing
end

"""
    ModiaMath.logOn!(obj)

Enable logging on `obj` (of type ModiaMath.SimulationState, ModiaMath.AbstractSimulationModel,
or ModiaMath.Logger)
"""
function logOn!(logger::Logger)
    logger.log = true
    return nothing
end

"""
    ModiaMath.logOff!(obj)

Disable logging on `obj` (of type ModiaMath.SimulationState, ModiaMath.AbstractSimulationModel,
or ModiaMath.Logger)
"""
function logOff!(logger::Logger)
    logger.log = false
    return nothing
end

"""
    ModiaMath.isLogStatistics(logger::ModiaMath.Logger)

Return true, if logger settings require to print **statistics** messages of the model.
"""
isLogStatistics(logger::Logger) = logger.log && logger.statistics

"""
    ModiaMath.isLogProgress(logger::ModiaMath.Logger)

Return true, if logger settings require to print **progress** messages of the model
"""
isLogProgress(logger::Logger) = logger.log && logger.progress

"""
    ModiaMath.isLogInfos(obj)

Return true, if logger settings require to print **info** messages of the model
(obj must be of type ModiaMath.SimulationState, ModiaMath.AbstractSimulationModel,
or ModiaMath.Logger).
"""
isLogInfos(logger::Logger)      = logger.log && logger.infos

"""
    ModiaMath.isLogWarnings(obj)

Return true, if logger settings require to print **warning** messages of the model
(obj must be of type ModiaMath.SimulationState, ModiaMath.AbstractSimulationModel,
or ModiaMath.Logger).
"""
isLogWarnings(logger::Logger)   = logger.log && logger.warnings

"""
    ModiaMath.isLogEvents(obj)

Return true, if logger settings require to print **event** messages of the model
(obj must be of type ModiaMath.SimulationState, ModiaMath.AbstractSimulationModel,
or ModiaMath.Logger).
"""
isLogEvents(logger::Logger)     = logger.log && logger.events

"""
    mutable struct SimulationStatistics - Collect statistics of the last simulation run.

The following data is stored in this structure:

- `cpuTimeInitialization`: CPU-time for initialization
- `cpuTimeIntegration`: CPU-time for integration
- `startTime`: start time of the integration
- `stopTime`: stop time of the integration
- `interval`: communication interval of the integration
- `tolerance`: relative tolerance used for the integration
- `nEquations`: number of equations (length of y and of yp)
- `nResults`: number of time points stored in result data structure
- `nSteps`: number of (successful) steps
- `nResidues`: number of calls to compute residues (includes residue calls for Jacobian)
- `nZeroCrossing`: number of calls to compute zero crossings
- `nJac`: number of calls to compute Jacobian
- `nTimeEvents`: number of time events
- `nRestartEvents`: number of events with integrator restart
- `nErrTestFails`: number of fails of error tests
- `h0`: stepsize used at the first step
- `hMin`: minimum integration stepsize
- `hMax`: maximum integration stepsize
- `orderMax`: maximum integration order
- `sparseSolver` = true: if sparse solver used, otherwise dense solver
- `nGroups`: if sparseSolver, number of column groups to compute Jacobian
  (e.g. if nEquations=100, nGroups=5, then 5+1=6 model evaluations are needed
  to compute the Jacobian numerically, instead of 101 model evaluations without
  taking the sparseness structure into account).
"""
mutable struct SimulationStatistics
    cpuTimeInitialization::Float64
    cpuTimeIntegration::Float64
    startTime::Float64
    stopTime::Float64
    interval::Float64
    tolerance::Float64
    nEquations::Int
    nResults::Int
    nSteps::Int
    nResidues::Int
    nZeroCrossings::Int
    nJac::Int
    nTimeEvents::Int
    nStateEvents::Int
    nRestartEvents::Int
    nErrTestFails::Int
    h0::Float64
    hMin::Float64
    hMax::Float64
    orderMax::Int
    sparseSolver::Bool
    nGroups::Int

    @static if VERSION >= v"0.7.0-DEV.2005"
        SimulationStatistics(nEquations::Int, sparseSolver::Bool, nGroups::Int) =
        new(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, nEquations, 0, 0, 0, 0, 0, 0, 0, 0, 0, floatmax(Float64),
            floatmax(Float64), 0.0, 0, sparseSolver, nGroups)
    else
        SimulationStatistics(nEquations::Int, sparseSolver::Bool, nGroups::Int) =
        new(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, nEquations, 0, 0, 0, 0, 0, 0, 0, 0, 0, realmax(Float64),
            realmax(Float64), 0.0, 0, sparseSolver, nGroups)
    end

end

function reInitializeStatistics!(stat::SimulationStatistics,
                                 startTime::Float64, stopTime::Float64, interval::Float64, tolerance::Float64)
    stat.cpuTimeInitialization = 0.0
    stat.cpuTimeIntegration    = 0.0
    stat.startTime      = startTime
    stat.stopTime       = stopTime
    stat.interval       = interval
    stat.tolerance      = tolerance
    stat.nResults       = 0
    stat.nSteps         = 0
    stat.nResidues      = 0
    stat.nZeroCrossings = 0
    stat.nJac           = 0
    stat.nTimeEvents    = 0
    stat.nStateEvents   = 0
    stat.nRestartEvents = 0
    stat.nErrTestFails  = 0

    @static if VERSION >= v"0.7.0-DEV.2005"
        stat.h0             = floatmax(Float64)
        stat.hMin           = floatmax(Float64)
    else
        stat.h0             = realmax(Float64)
        stat.hMin           = realmax(Float64)
    end
    stat.hMax           = 0.0
    stat.orderMax       = 0
end

import Base.show
Base.print(io::IO, stat::SimulationStatistics) = show(io, stat)

function Base.show(io::IO, stat::SimulationStatistics)
    @printf(io,"        cpuTime        = %.2g s (init: %.2g s, integration: %.2g s)\n",
                                           stat.cpuTimeInitialization + stat.cpuTimeIntegration,
                                           stat.cpuTimeInitialization, stat.cpuTimeIntegration)
    println(io, "        startTime      = ", stat.startTime, " s")
    println(io, "        stopTime       = ", stat.stopTime, " s")
    println(io, "        interval       = ", stat.interval, " s")
    println(io, "        tolerance      = ", stat.tolerance)
    println(io, "        nEquations     = ", stat.nEquations)
    println(io, "        nResults       = ", stat.nResults)
    println(io, "        nSteps         = ", stat.nSteps)
    println(io, "        nResidues      = ", stat.nResidues, " (includes residue calls for Jacobian)")
    println(io, "        nZeroCrossings = ", stat.nZeroCrossings)
    println(io, "        nJac           = ", stat.nJac)
    println(io, "        nTimeEvents    = ", stat.nTimeEvents)
    println(io, "        nStateEvents   = ", stat.nStateEvents)
    println(io, "        nRestartEvents = ", stat.nRestartEvents)
    println(io, "        nErrTestFails  = ", stat.nErrTestFails)
    @printf(io, "        h0             = %.2g s\n", stat.h0)
    @printf(io, "        hMin           = %.2g s\n", stat.hMin)
    @printf(io, "        hMax           = %.2g s\n", stat.hMax)
    println(io, "        orderMax       = ", stat.orderMax)
    println(io, "        sparseSolver   = ", stat.sparseSolver)

    if stat.sparseSolver
        println(io, "        nGroups        = ", stat.nGroups)
    end
end

function set_nResultsForSimulationStatistics!(stat::SimulationStatistics, nt::Int)
    stat.nResults = nt
    return nothing
end

end


text:
using Test
using CLIMAParameters
using CLIMAParameters.Atmos.EDMF

@testset "EDMF" begin

  struct EarthParameterSet <: AbstractEarthParameterSet end
  ps = EarthParameterSet()

  @test !isnan(EDMF.c_λ(ps))
  @test !isnan(EDMF.c_ε(ps))
  @test !isnan(EDMF.c_δ(ps))
  @test !isnan(EDMF.c_t(ps))
  @test !isnan(EDMF.β(ps))
  @test !isnan(EDMF.μ_0(ps))
  @test !isnan(EDMF.χ(ps))
  @test !isnan(EDMF.w_min(ps))
  @test !isnan(EDMF.lim_ϵ(ps))
  @test !isnan(EDMF.lim_amp(ps))
  @test !isnan(EDMF.a_min(ps))
  @test !isnan(EDMF.a_surf(ps))
  @test !isnan(EDMF.κ_star²(ps))
  @test !isnan(EDMF.ψϕ_stab(ps))
  @test !isnan(EDMF.α_d(ps))
  @test !isnan(EDMF.α_a(ps))
  @test !isnan(EDMF.α_b(ps))
  @test !isnan(EDMF.H_up_min(ps))
  @test !isnan(EDMF.c_d(ps))
  @test !isnan(EDMF.c_m(ps))
  @test !isnan(EDMF.c_b(ps))
  @test !isnan(EDMF.a1(ps))
  @test !isnan(EDMF.a2(ps))
  @test !isnan(EDMF.ω_pr(ps))
  @test !isnan(EDMF.Pr_n(ps))
  @test !isnan(EDMF.Ri_c(ps))
  @test !isnan(EDMF.smin_ub(ps))
  @test !isnan(EDMF.smin_rm(ps))

end


text:
"Default inference algorithm"
defalg(args...) = FailUnsat

"Default Ω to use"
# defΩ(args...) = SimpleΩ{Vector{Int}, Any}
defΩ(args...) = LinearΩ{Vector{Int}, UnitRange{Int64}, Vector{Any}}

"Default projection"
defΩProj(args...; OT = defΩ(args...)) = ΩProj{OT, idtype(OT)}

"Default callbacks"
defcb(args...) = donothing

"Single sample form from `x`"
function Base.rand(x::RandVar;
                   alg::SamplingAlgorithm = defalg(x),
                   ΩT = defΩ(alg),
                   kwargs...)
  first(rand(Random.GLOBAL_RNG, x, 1, alg; ΩT = ΩT, kwargs...))
end

function Base.rand(rng::AbstractRNG,
                   x::RandVar;
                   alg::SamplingAlgorithm = defalg(x),
                   ΩT::Type{OT} = defΩ(alg),
                   kwargs...) where {OT <: Ω}
  first(rand(rng, x, 1, alg; ΩT = ΩT, kwargs...))
end

function Base.rand(rng::AbstractRNG,
                   x::RandVar,
                   n::Integer,
                   alg::SamplingAlgorithm;
                   ΩT::Type{OT} = defΩ(alg),
                   kwargs...) where {OT <: Ω}
  logdensity = Omega.mem(logerr(indomainₛ(x)))
  ωsamples = rand(rng, ΩT, logdensity, n, alg; kwargs...)
  map(ω -> applynotrackerr(Omega.mem(x), ω), ωsamples)
end

function Base.rand(x::RandVar,
                   n::Integer;
                   alg::SamplingAlgorithm = defalg(x),
                   ΩT::Type{OT} = defΩ(alg),
                   kwargs...) where {OT <: Ω}
  rand(Random.GLOBAL_RNG, x, n, alg; ΩT = ΩT, kwargs...)
end

function Base.rand(rng::AbstractRNG,
                   x::RandVar,
                   n::Integer;
                   alg::SamplingAlgorithm = defalg(x),
                   ΩT::Type{OT} = defΩ(alg),
                   kwargs...) where {OT <: Ω}
  rand(rng, x, n, alg; ΩT = ΩT, kwargs...)
end

function Base.rand(x::RandVar,
                   n::Integer,
                   alg::SamplingAlgorithm;
                   ΩT::Type{OT} = defΩ(alg),
                   kwargs...) where {OT <: Ω}
  rand(Random.GLOBAL_RNG, x, n, alg; ΩT = ΩT, kwargs...)
end

# Condition (deprecate this?) x | y

Base.rand(x::RandVar, y::RandVar, n; kwargs...) = rand(Random.GLOBAL_RNG, cond(x, y), n; kwargs...)
Base.rand(x::RandVar, y::RandVar; kwargs...) = rand(Random.GLOBAL_RNG, cond(x, y); kwargs...)
Base.rand(rng::AbstractRNG, x::RandVar, y::RandVar, n; kwargs...) = rand(rng, cond(x, y), n; kwargs...)
Base.rand(rng::AbstractRNG, x::RandVar, y::RandVar; kwargs...) = rand(rng, cond(x, y); kwargs...)

# Return Tuple

Base.rand(x::UTuple{RandVar}, n::Integer; kwargs...) = rand(randtuple(x), n; kwargs...)
Base.rand(x::UTuple{RandVar}; kwargs...) = rand(randtuple(x); kwargs...)

Base.rand(x::UTuple{RandVar}, y::RandVar, n::Integer; kwargs...) = rand(randtuple(x), y, n; kwargs...)
Base.rand(x::UTuple{RandVar}, y::RandVar; kwargs...) = rand(randtuple(x), y; kwargs...)

# Return Ω object

"Returns Ω object, .e.g `x = normal(0, 1); rand(Ω, x > 0)`"
function ld(rng::AbstractRNG,
                   ΩT::Type{OT},
                   x::RandVar,
                   n::Integer,
                   alg::SamplingAlgorithm;
                   kwargs...) where {OT <: Ω}
  logdensity = Omega.mem(logerr(indomainₛ(x)))
  ωsamples = rand(rng, ΩT, logdensity, n, alg; kwargs...)
end

"Returns Ω object, e.g. `x = normal(0, 1); rand(Ω, cond(x, x >ₛ 2.0), 10)`"
Base.rand(rng::AbstractRNG, ::Type{Ω}, x::RandVar, n::Integer; alg::SamplingAlgorithm = defalg(x), kwargs...) =
  ld(rng, x, defΩ(alg), n, alg; kwargs...)

Base.rand(::Type{Ω}, x::RandVar, n::Integer; alg::SamplingAlgorithm = defalg(x), kwargs...) =
  ld(Random.GLOBAL_RNG, defΩ(alg), x, n, alg; kwargs...)

# Autorand

"`autorand(x::RandVar)` Single sample from `x`.  Inference `alg` and hyperparams, chosen automatically"
autorand(x::RandVar) = rand(x, 1000; alg = SSMH)[end]
Base.getindex(x::RandVar) = autorand(x)


text:
using ReachabilityAnalysis, JLD2
using ReachabilityAnalysis: add_dimension

LazySets.set_ztol(Float64, 1e-15);

examples_dir = normpath(@__DIR__, "..", "..", "..", "examples")
ISS_path = joinpath(examples_dir, "ISS", "ISS.jld2")

@load ISS_path C;
const C3 = C[3, :]; # variable y₃
const C3_ext = vcat(C3, fill(0.0, 3));

function ISSF01()
    @load ISS_path A B

    U = Hyperrectangle(low=[0.0, 0.8, 0.9], high=[0.1, 1., 1.])
    X0 = BallInf(zeros(size(A, 1)), 0.0001)
    return @ivp(x' = A*x + B*u, x(0) ∈ X0, u ∈ U, x ∈ Universe(270))
end

dirs = CustomDirections([C3, -C3]);
prob_ISSF01 = ISSF01();
sol_ISSF01 = solve(prob_ISSF01, T=20.0, alg=LGG09(δ=6e-4, template=dirs, sparse=true, cache=false));

dim(sol_ISSF01)

πsol_ISSF01 = flatten(sol_ISSF01);

dim(πsol_ISSF01)

using Plots, Plots.PlotMeasures, LaTeXStrings

fig = Plots.plot();
Plots.plot!(fig, πsol_ISSF01[1:10:end], vars=(0, 1), linecolor=:blue, color=:blue, alpha=0.8,
     xlab=L"t",
     ylab=L"y_{3}",
     xtick=[0, 5, 10, 15, 20.], ytick=[-0.00075, -0.0005, -0.00025, 0, 0.00025, 0.0005, 0.00075],
     xlims=(0., 20.), ylims=(-0.00075, 0.00075),
     bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm);
fig

function ISSC01()
    @load ISS_path A B

    A_ext = add_dimension(A, 3)
    A_ext[1:270, 271:273] = B

    U = Hyperrectangle(low=[0.0, 0.8, 0.9], high=[0.1, 1., 1.])
    X0 = BallInf(zeros(size(A, 1)), 0.0001)
    X0 = X0 * U
    return @ivp(x' = A_ext*x, x(0) ∈ X0)
end

dirs = CustomDirections([C3_ext, -C3_ext]);
prob_ISSC01 = ISSC01();
sol_ISSC01 = solve(prob_ISSC01, T=20.0, alg=LGG09(δ=0.01, template=dirs, sparse=true, cache=false));

πsol_ISSC01 = flatten(sol_ISSC01);

fig = Plots.plot();
Plots.plot!(fig, πsol_ISSC01, vars=(0, 1), linecolor=:blue, color=:blue, alpha=0.8, lw=1.0,
    xlab=L"t",
    ylab=L"y_{3}",
    xtick=[0, 5, 10, 15, 20.], ytick=[-0.0002, -0.0001, 0.0, 0.0001, 0.0002],
    xlims=(0., 20.), ylims=(-0.0002, 0.0002),
    bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm)
fig


text:
#!/usr/bin/env bash
# -*-coding:utf-8-*-
#=
exec julia --color=yes --startup-file=no -e 'include(popfirst!(ARGS))' \
    "${BASH_SOURCE[0]}" "$@"
=#

# Copyright (c) 2021, Oliver Merkel.
# Please see the AUTHORS file for details.
# All rights reserved.
#
# Use of this code is governed by a
# MIT license that can be found in the LICENSE file.

const FILENAME = "dict_en.txt"
const MAXFALSEATTEMPTS = 5

function welcome()
    println("Guess my secret word with a maximum of $MAXFALSEATTEMPTS false attempts.")
end

function get_word()
    words = readlines(FILENAME)
    word = words[rand(1:length(words))]
    return word
end

function get_test_result( word::String, guesses::String )
    result = ""
    incorrect = 0
    for c in word
        if c in guesses
            result *= c
        else
            result *= "_"
            incorrect += 1
        end
    end
    ( result, incorrect )
end

function play_again()
    yes = [ "y", "yes", "sure" ]
    no = [ "n", "no", "nope" ]
    is_yes = false
    is_no = false
    while (!is_yes && !is_no)
        println("Play again?")
        answer = lowercase(readline())
        is_yes = answer in yes
        is_no = answer in no
        if (!is_yes && !is_no)
            println("Sorry? Please answer with yes or no...")
        end
    end
    return is_yes
end

function game_loop()
    welcome()
    word = get_word()
    end_of_game = false
    false_attempts = 0
    guesses = ""
    while !end_of_game
        ( result, incorrect ) = get_test_result( word, guesses )
        println(result)
        if incorrect == 0
            end_of_game = true
            println("You won!")
        else
            println("Guess a character: ")
            guess = lowercase(strip(readline()))
            if length(guess) != 1
                println("Type a single character only! Try again...")
            elseif !all(isletter, guess)
                println("Type a alphabetic character! Try again...")
            elseif occursin(guess, guesses)
                println("You already guessed this character! Try again...")
            else
                guesses *= guess
                if !occursin(guess, word)
                    false_attempts += 1
                    println("Wrong.")
                    if false_attempts == 1
                        println("You have one false attempt.")
                    else
                        println("You have $false_attempts false attempts.")
                    end
                    if false_attempts > MAXFALSEATTEMPTS
                        end_of_game = true
                        println("You lost!")
                        println("The secret word has been ", uppercase(word))
                    end
                end
            end
        end
    end
end

function menu_loop()
    keep_playing = true
    while keep_playing
        game_loop()
        keep_playing = play_again()
    end
    println("Bye.");
end

menu_loop()

max_stars_repo_name: OMerkel/hangman

max_issues_repo_name: OMerkel/hangman

max_forks_repo_name: OMerkel/hangman


text:
using SIAN

println("Setting up the problem")

ode = @ODEmodel(
  x1'(t) = -1 * p1 * x1(t) + x2(t) + u0(t),
  x2'(t) = p3 * x1(t) - p4 * x2(t) + x3(t),
  x3'(t) = p6 * x1(t) - p7 * x3(t),
  u0'(t) = 1,
  y(t) = x1(t),
  y2(t) = u0(t)
)

res = identifiability_ode(ode, get_parameters(ode); p=0.99, p_mod=0, nthrds=1)

println(res)

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
using CovidRt
using Test

@testset "CovidRt.jl" begin
    # Write your own tests here.
end


text:
"""
    Yeoh(; λ, μ, c₂, c₃)

Hyper elastic material
#Arguments
- `λ::Float64`: Lamé parameter
- `μ::Float64`: Lamé parameter (shear modulus)
- `c₂::Float64`: Material constant
- `c₃::Float64`: Material constant
"""
struct Yeoh <: AbstractMaterial
    λ::Float64
    μ::Float64
    c₂::Float64
    c₃::Float64
end

struct YeohState <: AbstractMaterialState
end

function initial_material_state(::Yeoh)
    return YeohState()
end

function Yeoh(; λ::T, μ::T, c₂::T, c₃::T) where T <: AbstractFloat
    return Yeoh(λ, μ, c₂, c₃)
end

strainmeasure(::Yeoh) = RightCauchyGreen()

function ψ(mp::Yeoh, C::SymmetricTensor{2,3})
    J = sqrt(det(C))
    I = tr(C)
    return mp.μ/2 * (I-3) + mp.c₂*(I-3)^2 + mp.c₃*(I-3)^3 - mp.μ*log(J) + mp.λ/2 * log(J)^2
end

function material_response(mp::Yeoh, C::SymmetricTensor{2,3}, state::YeohState = YeohState(),
                           Δt=nothing; cache=nothing, options=nothing)
    invC=inv(C)
    detC=det(C);
    J=sqrt(detC);

    dlnJdc = 1/2*invC
    d2lnJdcdc = 1/2*inv(otimesu(-C,C))

    𝐈 = one(C)
    Ic = tr(C)

    S = 2*(mp.μ/2*𝐈 + 2*mp.c₂*(Ic-3)*𝐈 + 3*mp.c₃*(Ic-3)^2*𝐈 - mp.μ*dlnJdc + mp.λ* log(J) * dlnJdc)
    ∂S∂E = 4 * (2*mp.c₂*(𝐈 ⊗ 𝐈) + 6*mp.c₃*(Ic-3)*(𝐈 ⊗ 𝐈) - mp.μ*d2lnJdcdc + mp.λ*(log(J)*d2lnJdcdc + (dlnJdc ⊗ dlnJdc)))
    ∂S∂E = symmetric(∂S∂E)
    # TODO use transform tangents from traits branch

    return S, ∂S∂E, YeohState()
end


text:
using Dataframes, StatsBase
aapl = readtable("AAPL-Short,csv");

naapl = size(aapl)[1]
m1 = int(mean((aapl[:Volume]))); # => 6306547

using Match
wts = zeros(naapl);
for i in 1:naapl
    dt = aapl[:Date][i]
    wts[i] = @match dt begin
         r"^2000" => 1.0
         r"^2001" => 2.0
         r"^2002" => 4.0
         _        => 0.0
    end
end;

wv = WeightVec(wts);
m2 = int(mean(aapl[:Volume], wv))
mean(aapl[:Close])
mean(aapl[:Close], wv)
trimmean(aapl[:Close], 0.1)

std(aapl[:Close])
skewness(aapl[:Close])
kurtosis(aapl[:Close])


text:
###
### Gibbs samplers / compositional samplers.
###

const GibbsComponent = Union{Hamiltonian,MH,ESS,PG}

"""
    Gibbs(algs...)

Compositional MCMC interface. Gibbs sampling combines one or more
sampling algorithms, each of which samples from a different set of
variables in a model.

Example:
```julia
@model gibbs_example(x) = begin
    v1 ~ Normal(0,1)
    v2 ~ Categorical(5)
end
```

# Use PG for a 'v2' variable, and use HMC for the 'v1' variable.
# Note that v2 is discrete, so the PG sampler is more appropriate
# than is HMC.
alg = Gibbs(HMC(0.2, 3, :v1), PG(20, :v2))
```

Tips:
- `HMC` and `NUTS` are fast samplers, and can throw off particle-based
methods like Particle Gibbs. You can increase the effectiveness of particle sampling by including
more particles in the particle sampler.
"""
struct Gibbs{space, A<:Tuple{Vararg{GibbsComponent}}} <: InferenceAlgorithm
    algs::A   # component sampling algorithms
end

function Gibbs(algs::GibbsComponent...)
    # obtain space of sampling algorithms
    space = Tuple(union(getspace.(algs)...))

    Gibbs{space, typeof(algs)}(algs)
end

"""
    GibbsState{V<:VarInfo, S<:Tuple{Vararg{Sampler}}}

Stores a `VarInfo` for use in sampling, and a `Tuple` of `Samplers` that
the `Gibbs` sampler iterates through for each `step!`.
"""
mutable struct GibbsState{V<:VarInfo, S<:Tuple{Vararg{Sampler}}} <: AbstractSamplerState
    vi::V
    samplers::S
end

function GibbsState(model::Model, samplers::Tuple{Vararg{Sampler}})
    return GibbsState(VarInfo(model), samplers)
end

function Sampler(alg::Gibbs, model::Model, s::Selector)
    # sanity check for space
    space = getspace(alg)
    # create tuple of samplers
    i = 0
    samplers = map(alg.algs) do _alg
        i += 1
        if i == 1
            prev_alg = alg.algs[end]
        else
            prev_alg = alg.algs[i-1]
        end
        rerun = !(_alg isa MH) || prev_alg isa PG || prev_alg isa ESS
        selector = Selector(Symbol(typeof(_alg)), rerun)
        Sampler(_alg, model, selector)
    end
    # create a state variable
    state = GibbsState(model, samplers)

    # create the sampler
    info = Dict{Symbol, Any}()
    spl = Sampler(alg, info, s, state)

    # add Gibbs to gids for all variables
    vi = spl.state.vi
    for sym in keys(vi.metadata)
        vns = getfield(vi.metadata, sym).vns

        for vn in vns
            # update the gid for the Gibbs sampler
            DynamicPPL.updategid!(vi, vn, spl)

            # try to store each subsampler's gid in the VarInfo
            for local_spl in samplers
                DynamicPPL.updategid!(vi, vn, local_spl)
            end
        end
    end

    return spl
end

# Initialize the Gibbs sampler.
function sample_init!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Gibbs},
    N::Integer;
    kwargs...
)
    # Initialize each local sampler.
    for local_spl in spl.state.samplers
        sample_init!(rng, model, local_spl, N; kwargs...)
    end
end

# Finalize the Gibbs sampler.
function sample_end!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Gibbs},
    N::Integer;
    kwargs...
)
    # Finalize each local sampler.
    for local_spl in spl.state.samplers
        sample_end!(rng, model, local_spl, N; kwargs...)
    end
end

# First step.
function step!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Gibbs},
    N::Integer;
    kwargs...
)
    Turing.DEBUG && @debug "Gibbs stepping..."

    time_elapsed = 0.0

    # Iterate through each of the samplers.
    for local_spl in spl.state.samplers
        Turing.DEBUG && @debug "$(typeof(local_spl)) stepping..."

        Turing.DEBUG && @debug "recording old θ..."

        # Update the sampler's VarInfo.
        local_spl.state.vi = spl.state.vi

        # Step through the local sampler.
        time_elapsed_thin =
            @elapsed step!(rng, model, local_spl, N; kwargs...)

        # After the step, update the master varinfo.
        spl.state.vi = local_spl.state.vi

        # Uncomment when developing thinning functionality.
        # Retrieve symbol to store this subsample.
        # symbol_id = Symbol(local_spl.selector.gid)

        # # Store the subsample.
        # spl.state.subsamples[symbol_id][] = trans

        # Record elapsed time.
        time_elapsed += time_elapsed_thin
    end

    return Transition(spl)
end

# Steps 2:N
function step!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Gibbs},
    N::Integer,
    t::AbstractTransition;
    kwargs...
)
    Turing.DEBUG && @debug "Gibbs stepping..."

    time_elapsed = 0.0

    # Iterate through each of the samplers.
    for local_spl in spl.state.samplers
        Turing.DEBUG && @debug "$(typeof(local_spl)) stepping..."

        Turing.DEBUG && @debug "recording old θ..."

        # Update the sampler's VarInfo.
        local_spl.state.vi = spl.state.vi

        # Step through the local sampler.
        time_elapsed_thin =
            @elapsed trans = step!(rng, model, local_spl, N, t; kwargs...)

        # After the step, update the master varinfo.
        spl.state.vi = local_spl.state.vi

        # Uncomment when developing thinning functionality.
        # Retrieve symbol to store this subsample.
        # symbol_id = Symbol(local_spl.selector.gid)
        #
        # # Store the subsample.
        # spl.state.subsamples[symbol_id][] = trans

        # Record elapsed time.
        time_elapsed += time_elapsed_thin
    end

    return Transition(spl)
end


text:
a = [UncertainValue(Normal, 1, 0.5) for i = 1:10]
b = UncertainValue(rand(1000))
c = UncertainValue(Uniform, rand(10000))
s = rand(-20:20, 100)
e = UncertainValue(s, rand(length(s)))
f = UncertainValue(2)
g = UncertainValue(3.1)

uidxs =  UncertainIndexDataset([a; b; c; e; f; g])
uvals = UncertainValueDataset([a; b; c; e; f; g])
iv = UncertainIndexValueDataset(uidxs, uvals)

@test resample(iv, NoConstraint()) isa Tuple{Vector{Float64}, Vector{Float64}}
@test resample(iv, TruncateLowerQuantile(0.2)) isa Tuple{Vector{Float64}, Vector{Float64}}
@test resample(iv, TruncateUpperQuantile(0.2)) isa Tuple{Vector{Float64}, Vector{Float64}}
@test resample(iv, TruncateQuantiles(0.2, 0.8)) isa Tuple{Vector{Float64}, Vector{Float64}}
@test resample(iv, TruncateMaximum(10)) isa Tuple{Vector{Float64}, Vector{Float64}}
@test resample(iv, TruncateMinimum(-20)) isa Tuple{Vector{Float64}, Vector{Float64}}
@test resample(iv, TruncateRange(-20, 20)) isa Tuple{Vector{Float64}, Vector{Float64}}

@test resample(iv, NoConstraint(), 5) isa Vector{Tuple{Vector{Float64}, Vector{Float64}}}
@test resample(iv, TruncateLowerQuantile(0.2), 5) isa Vector{Tuple{Vector{Float64}, Vector{Float64}}}
@test resample(iv, TruncateUpperQuantile(0.2), 5) isa Vector{Tuple{Vector{Float64}, Vector{Float64}}}
@test resample(iv, TruncateQuantiles(0.2, 0.8), 5) isa Vector{Tuple{Vector{Float64}, Vector{Float64}}}
@test resample(iv, TruncateMaximum(10), 5) isa Vector{Tuple{Vector{Float64}, Vector{Float64}}}
@test resample(iv, TruncateMinimum(-20), 5) isa Vector{Tuple{Vector{Float64}, Vector{Float64}}}
@test resample(iv, TruncateRange(-20, 20), 5) isa Vector{Tuple{Vector{Float64}, Vector{Float64}}}

@test length(resample(iv, NoConstraint(), 5)) == 5
@test length(resample(iv, TruncateLowerQuantile(0.2), 5)) == 5
@test length(resample(iv, TruncateUpperQuantile(0.2), 5)) == 5
@test length(resample(iv, TruncateQuantiles(0.2, 0.8), 5)) == 5
@test length(resample(iv, TruncateMaximum(10), 5)) == 5
@test length(resample(iv, TruncateMinimum(-20), 5)) == 5
@test length(resample(iv, TruncateRange(-20, 20), 5)) == 5


text:
# Linear mean function

"""
    MeanLin <: Mean

Linear mean function
```math
m(x) = xᵀβ
```
with linear coefficients ``β``.
"""
mutable struct MeanLin <: Mean
    "Linear coefficients"
    β::Vector{Float64}
    "Priors for mean parameters"
    priors::Array

    """
        MeanLin(β::Vector{Float64})

    Create `MeanLin` with linear coefficients `β`.
    """
    MeanLin(β::Vector{Float64}) = new(β, [])
end

mean(mLin::MeanLin, x::AbstractVector) = dot(mLin.β, x)
mean(mLin::MeanLin, X::AbstractMatrix) = X'mLin.β

get_params(mLin::MeanLin) = mLin.β
get_param_names(::MeanLin) = [:β]
num_params(mLin::MeanLin) = length(mLin.β)

function set_params!(mLin::MeanLin, hyp::AbstractVector)
    length(hyp) == length(mLin.β) || throw(ArgumentError("Linear mean function only has $(mLin.dim) parameters"))
    copyto!(mLin.β, hyp)
end

function grad_mean(mLin::MeanLin, x::AbstractVector)
    dM_theta = x
    return dM_theta
end


text:
using ArgParse
using DrWatson
@quickactivate
using CSV
using BSON
using Random
using FileIO
using DataFrames
using PrettyTables
using PrettyTables.Crayons

using GenerativeAD.Evaluation: _prefix_symbol, PAT_METRICS, aggregate_stats_mean_max
using GenerativeAD.Evaluation: rank_table, print_rank_table

s = ArgParseSettings()
@add_arg_table! s begin
	"filename"
		arg_type = String
		default = "evaluation/images_eval.bson"
		help = "Location of cached DataFrame."
	"-c", "--criterion-metric"
		arg_type = String
		default = "val_auc"
		help = "Criterion to sort the models."
	"-r", "--rank-metric"
		arg_type = String
		default = ""
		help = "Metric to rank models."
	"-b", "--backend"
		arg_type = String
		default = "txt"
		help = "Backend for PrettyTable print. Either of [txt (default), tex, html] is allowed."
	"-o", "--output-prefix"
		arg_type = String
		default = "evaluation/images_eval"
		help = "Output prefix for storing results."
	"-v", "--verbose"
		action = :store_true
		help = "Print all results instead of storing to files."
	"-p", "--proportional"
		action = :store_true
		help = "Overloads criterion and uses pat@x% with incresing x. Prints only ranks."
	"--best-params"
		action = :store_true
		help = "Stores CSV files for each model's best parameters."
end

function main(args)
	f = datadir(args["filename"])
	df = load(f)[:df]
	@info "Loaded $(nrow(df)) rows from $f"

	if args["proportional"]
		ranks = []
		if args["rank-metric"] != ""
			for criterion in _prefix_symbol.("val", PAT_METRICS)
				df_agg = aggregate_stats_mean_max(df, criterion)
				push!(ranks, rank_table(df_agg, args["rank-metric"])[end:end, :])
			end
		else # pat/pat scenario if no rank-metric is provided
			for (criterion, metric) in zip(_prefix_symbol.("val", PAT_METRICS), _prefix_symbol.("tst", PAT_METRICS))
				df_agg = aggregate_stats_mean_max(df, criterion)
				push!(ranks, rank_table(df_agg, metric)[end:end, :])
			end
		end

		df_ranks = vcat(ranks...)
		df_ranks[:, :dataset] .= PAT_METRICS
		hl_best_rank = Highlighter(
					f = (data, i, j) -> (data[i,j] == minimum(df_ranks[i, 2:end])),
					crayon = crayon"green bold")

		rank_metric = (args["rank-metric"] != "") ? args["rank-metric"] : "tst_pat_x"

		@info "Best models chosen by validation pat@x% with incresing x"
		@info "Ranking by $(rank_metric)"

		if ~args["verbose"]
			target_filename = datadir(args["output-prefix"],"ranks_pat_$(rank_metric).txt")
			(~isdir(dirname(target_filename))) && mkdir(dirname(target_filename))
			open(target_filename, "w") do io
				pretty_table(
					io,
					df_ranks,
					formatters = ft_round(2),
					highlighters = (hl_best_rank),
				)
			end
		else
			pretty_table(
					df_ranks,
					formatters = ft_round(2),
					highlighters = (hl_best_rank),
				)
		end
	else
		df_agg = aggregate_stats_mean_max(df, Symbol(args["criterion-metric"]))
		rt = rank_table(df_agg, args["rank-metric"])

		@info "Best models chosen by $(args["criterion-metric"])"
		@info "Ranking by $(args["rank-metric"])"
		if ~args["verbose"]
			target_filename = datadir(args["output-prefix"],"$(args["criterion-metric"])_$(args["rank-metric"]).$(args["backend"])")
			(~isdir(dirname(target_filename))) && mkdir(dirname(target_filename))
			open(target_filename, "w") do io
				print_rank_table(io, rt; backend=Symbol(args["backend"]))
			end
		else
			print_rank_table(rt; backend=Symbol(args["backend"]))
		end

		if args["best-params"]
			target_dir = datadir(args["output-prefix"], "params_"*args["criterion-metric"])
			@info "Storing best parameters as separate CSVs in $(target_dir)"
			(~isdir(target_dir)) && mkdir(target_dir)
			all_models = unique(df_agg[:modelname])
			for m in all_models
				magg = filter(x -> (x.modelname == m), df_agg)
				CSV.write(joinpath(target_dir, "$(m).csv"), magg)
			end
		end
	end
end

main(parse_args(ARGS, s))


text:
include("solver_preamble.jl")
import ECOS
factory = with_optimizer(ECOS.Optimizer, verbose=false)
config = MOI.Test.TestConfig(atol=1e-5, rtol=1e-5, query=false)
@testset "Linear" begin
    Tests.linear_test(factory, config)
end
@testset "SOC" begin
    Tests.soc_test(factory, config)
end


text:
#========================================================================================#
#	Laboratory 7
#
# Encapsulation and software design.
#
# Author: Niall Palfreyman, 06/02/2022
#========================================================================================#
[
	Activity(
		"""
		In this laboratory we look at the issue of encapsulation that we saw at the end of the
		previous laboratory. There, we found that the value of the GLOBAL variable fig had been
		overwritten by the code inside the function myheatmap() ...

			THIS IS SOMETHING THAT SHOULD _NEVER_EVER_ HAPPEN!

		We are sometimes tempted to use global variables to pass data between functions, but
		this ALWAYS carries the danger that someone might change the value of these variables
		by accident - particularly since it is usually extremely difficult for users of our
		program to notice that we are using global variables to pass the data!

		A large part of the science of informatics concerns how to pass data in ways that
		protect that data from being changed by accident, and that is what we will investigate
		in this laboratory.

		In our first experiment, enter the following code, then tell me the value of paula:

		linus = [5,4,3,2,1]; paula = 5
		""",
		"",
		x -> x==5
	),
	Activity(
		"""
		Now enter the following function:

		function change_paula()
			paula = paula + 2
			paula
		end

		Then call the function change_paula() and tell me the value you get back:
		""",
		"Your result may (not) surprise you, depending on how you think about scoping rules",
		x -> x==7
	),
	Activity(
		"""
		Your result shows us two things: a) The value of the GLOBAL variable paula is available
		inside the LOCAL scope of the function change_paula(); b) We can change the value of
		paula within this LOCAL scope.

		Now tell me the current value of paula:
		""",
		"Ask Julia for the value of paula",
		x -> x==5
	),
	Activity(
		"""
		Aha! So although we can change the value of paula locally within change_paula(), this
		does not change the GLOBAL value of paula. In fact, there exist two different variables
		named paula: the GLOBAL variable containing the value 5, and a LOCAL variable containing
		the value 7. When change_paula() ends, the variables in its local scope are thrown away,
		and the LOCAL paula disappears.

		If we REALLY want to change the global value of paula, we can do so by redefining the
		function change_paula():

		function change_paula()
			global paula = paula + 2
			paula
		end

		This is what I did in the function myheatmap of laboratory 6. It is NOT a good idea! Tell
		me the value of paula now:
		""",
		"",
		x -> x==7
	),
	Activity(
		"""
		So Julia does allow us to make use of global values inside a local scope, but it forces us
		to announce this by using the keyword "global".

		There is a further issue here. Enter this code:

		function change_linus()
			linus[3] = 7
			linus
		end

		Again, the return value tells us that we are able to change the value of a local variable
		named linus. But now tell me the value of the GLOBAL variable named linus:
		""",
		"",
		x -> x==[5,4,7,2,1]
	),
	Activity(
		"""
		The problem is that linus is a Vector that refers to its contents (5,4,3,2,1). We are not
		allowed to change the value of linus, but we ARE allowed to change the contents that it
		refers to. So global variables are still unsafe! What are we to do? The solution is this:

		ALWAYS encapsulate (i.e.: wrap/hide) EVERYTHING you do inside a MODULE!

		Let's see how to do this. Enter the following code:

		module MyModule
			function change_paula1()
				global paula = 9
				paula
			end
		end

		Now call MyModule.change_paula1(), then tell me the value of paula afterwards:
		""",
		"",
		x -> x==7
	),
	Activity(
		"""
		OK, so putting change_paula1() inside the module MyModule means it cannot interfere
		with the value of our global variable paula. But wait! We know that we can load
		modules into global scope by means of the keyword "using": will that make it possible
		for users change paula's value by accident? Load the module MyModule now:

		using .MyModule

		Repeat the previous experiment - what is the value of paula afterwards?
		""",
		"",
		x -> x==7
	),
	Activity(
		"""
		Great! Now that we know how to hide variables and functions inside a module, we can do
		some real live software development! In the next laboratory, we will develop a genetic
		algorithm (GA), and GAs need to work with arrays of random values. However, generating
		random numbers is very time-expensive, so in this laboratory we develop a Casino module
		that can generate arrays of random values very quickly.

		There is just one thing we need to do first. Our modules can get quite complex, so we
		will build them up step-by-step within the Ingolstadt filesystem. Please create now a
		subfolder of Ingolstadt\\src named "Development", then create within the Development
		folder another subfolder named "Casinos". Finally, create withïn the Casinos folder an
		empty file named "Casinos.jl" that will contain our wonderful new module ...
		""",
		"Create this file structure now, before continuing",
		x -> isfile("Development/Casinos/Casinos.jl")
	),
	Activity(
		"""
		If you study the file Ingolstadt.jl, you will see that it contains my source code. It's
		all packed into a module named Ingolstadt, and this module is prefixed by a triple-quoted
		multi-line string that describes its purpose. Use my source code as a template to create
		in the file Casinos.jl a module Casinos prefixed by an explanatory help-string, and which
		contains the following test code:

		module Casinos
		test = 5
		end
		""",
		"",
		x -> true
	),
	Activity(
		"""
		Now test your new module by entering this line at the Julia prompt:

		include("src/Development/Casinos/Casinos.jl")

		This will read and parse the Casino module, and you can test it by telling me the
		answer to the following line:

		Casinos.test
		""",
		"",
		x -> x==5
	),
	Activity(
		"""
		Now we can develop our Casinos module in the file Casinos.jl. We start by setting up a
		use-case for the module - that is, we sketch out how we will want to use the module when it
		is finished. Replace your module definition in the file Casinos.jl by the following code,
		then move on to the next activity:

		module Casinos

		#-----------------------------------------------------------------------------------------
		\"""
			unittest()

		Unit-test the Casinos module.
		\"""
		function unittest()
			println("\\n============ Unit test Casinos: ===============")
			println("Casino vault of randomness 2 for matrix withdrawals up to size (2x3):")
			casino = Casino(2,3,2)
			display( casino.vault)
			println()

			println("Draw several (2x3) matrices from the casino:")
			display( draw( casino,2,3)); println()
			display( draw( casino,2,3)); println()
			display( draw( casino,2,3)); println()

			println("Finally, reshuffle the casino and redisplay its vault:")
			shuffle!(casino)
			display(casino.vault)
		end

		end # of Casinos
		""",
		"Note: you can't compile or run Casinos yet - just enter this code and move on",
		x -> true
	),
	Activity(
		"""
		You can test your new module by reincluding Casinos.jl:

		include("src/Development/Casinos/Casinos.jl")

		If you enter Casinos.unittest() at the Julia prompt, you should see that the command
		runs, but throws various errors (exceptions). We will now start to fix those errors...

		First let's get comment out the lines of unittest() that are causing problems. Insert the
		multi-line comment marker #= at the beginning of the third line of unittest() so that
		it looks like this:

		#=	casino = Casino(2,3,2)

		Next, close this multi-line comment by inserting the marker =# to the right of the final
		line of unittest():

			display(casino.vault) =#

		Reinclude Casinos.jl. Now you should be able to run Casino.unittest() without errors.
		Tell me the last six characters that you see in the output:
		""",
		"You should get two lines of output, and the final characters specify a matrix size",
		x -> x == "(2x3):"
	),
	Activity(
		"""
		To start developing our new module, we'll define the type Casino. Insert the following
		code AFTER the "module" line and BEFORE the comment box before unittest(), and check that
		everything still includes and runs properly:

		using Random

		#-----------------------------------------------------------------------------------------
		# Module types:

		\"""
			Casino

		A Casino can return arrays of random numbers in the range [0,1), up to a maximum number of rows
		(maxrows), and a maximum number of columns (maxcols). It also contains a vault of prepared
		random numbers from which it draws the arrays.
		\"""
		struct Casino
			maxrows::Int							# Maximum number of drawable rows
			maxcols::Int							# Maximum number of drawable columns
			randomness::Int							# How randomised will our withdrawals be?
			vault::Matrix							# Repository of random numbers in [0,1)

			"The one and only constructor"
			function Casino(maxrows::Int,maxcols::Int,randomness::Int=5)
				new(
					maxrows, maxcols, randomness,
					rand((maxrows+1)*randomness,(maxcols+1)*randomness)
				)
			end
		end

		#-----------------------------------------------------------------------------------------
		# Module methods:
		""",
		"Remember your code won't do anything new yet: just getting it to run is your first goal!",
		x -> true
	),
	Activity(
		"""
		Right, now let's activate the new code. Remove the open-comment marker from the third line
		of unittest() and insert it instead at the beginning of the 7th line:

		#=	println("Draw several (2x3) matrices from the casino:")

		This reveals the lines 3-6. Test that your program can now correctly create and display a
		Casino. Now tell me the size of the Casino vault, and think about why I have designed this
		size to depend on the constructor argument randomness:
		""",
		"",
		x -> x == (6,8)
	),
	Activity(
		"""
		Did you work out why the size of the vault depends on randomness? The idea is that we
		want to draw random matricess from the Casino like drawing cards from a shuffled deck of
		cards, and that is only possible if the deck contains more cards than we actually need.

		Now we'll implement the draw() method. First reveal the next four lines of unittest() and
		insert the following dummy code immediately after the "Module methods:" comment. Now test
		this dummy version of draw to make sure it is robust before moving on:

		\"""
			draw( casino, nrows, ncols)

		Draw the required number of rows and columns from the casino vault, first ensuring that the
		vault is large enough to support the withdrawal.
		\"""
		function draw( casino::Casino, nrows::Int, ncols::Int)
			if nrows > casino.maxrows || ncols > casino.maxcols
				# Repository is too small - throw exception:
				error( "Requested withdrawal is too large")
			end

			# Choose random offsets and strides for drawing a matrix of size (nrows x ncols) from
			# the vault, assuming that it is big enough to support the withdrawal:
			ones(nrows,ncols)
		end

		""",
		"",
		x -> true
	),
	Activity(
		"""
		Now comes the cool part of the code. Notice how we have made absolutely sure that by the
		time we get to the dummy line "ones(nrows,ncols)", we can rely on the values nrows and
		ncols being small enough to be able to draw our new random matrix. Now replace this
		dummy line by the following code and test it:

			vaultrows, vaultcols = size(casino.vault)
			offset_r = rand( 1 : (vaultrows-nrows))
			stride_r = (nrows <= 1) ? 1 :
							rand( 1 : (vaultrows-offset_r) ÷ (nrows-1))
			offset_c = rand( 1 : (vaultcols-ncols))
			stride_c = (ncols <= 1) ? 1 :
							rand( 1 : (vaultcols-offset_c) ÷ (ncols-1))

			# Return a randomly chosen table of slices from the vault:
			casino.vault[
				(offset_r : stride_r : (offset_r + (nrows-1)*stride_r)),
				(offset_c : stride_c : (offset_c + (ncols-1)*stride_c))
			]
		""",
		"",
		x -> true
	),
	Activity(
		"""
		Finally, we'll now give Casinos the functionality of shuffling. We want users to be
		able to shuffle the contents of the vault to give new random values. Reveal all lines
		of code in unittest(), insert the command "using Random" at the beginning of the Casinos
		module - immediately after the line "module Casinos", then insert the following code after
		draw() in Casinos. Then test the module again:

		#-----------------------------------------------------------------------------------------
		\"""
			shuffle!( casino)

		Reassign random values in the vault.
		\"""
		function shuffle!( casino::Casino)
			rand!( casino.vault)
		end

		""",
		"",
		x -> true
	),
	Activity(
		"""
		Congratulations! You have written your first Julia module! You can test its functionality
		for yourself by doing something like this at the Julia prompt:

		include("src/Development/Casinos/Casinos.jl")
		casinos = Casinos.Casino(3,3,5)
		Casinos.draw(casino,3,3)

		Make sure you also test error cases like this:

		Casinos.draw(casino,7,9)
		""",
		"",
		x -> true
	),
	Activity(
		"""
		There is just one small thing we can do to make life easier for users of the Casinos
		module: We shall expose the new functionality. After all, it's a pain to have to write
		"Casinos." in front of every command. To avoid this, insert the following lines between
		the line "module Casinos" and the line "using Random":

		# Externally callable methods of Casinos
		export Casino, draw, shuffle!

		Now reinclude Casinos.jl and repeat your tests from the previous activity, but this time
		first enter "using .Casinos" at the Julia prompt. Can you now call all the Casinos
		functionality without typing "Casinos." in front of everything?
		""",
		"",
		x -> occursin('y',lowercase(x))
	),
]


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## symbols ##

"""
    gensym([tag])

Generates a symbol which will not conflict with other variable names.
"""
gensym() = ccall(:jl_gensym, Ref{Symbol}, ())

gensym(s::String) = ccall(:jl_tagged_gensym, Ref{Symbol}, (Ptr{UInt8}, Int32), s, sizeof(s))

gensym(ss::String...) = map(gensym, ss)
gensym(s::Symbol) =
    ccall(:jl_tagged_gensym, Ref{Symbol}, (Ptr{UInt8}, Int32), s, ccall(:strlen, Csize_t, (Ptr{UInt8},), s))

"""
    @gensym

Generates a gensym symbol for a variable. For example, `@gensym x y` is transformed into
`x = gensym("x"); y = gensym("y")`.
"""
macro gensym(names...)
    blk = Expr(:block)
    for name in names
        push!(blk.args, :($(esc(name)) = gensym($(string(name)))))
    end
    push!(blk.args, :nothing)
    return blk
end

## expressions ##

function copy(e::Expr)
    n = Expr(e.head)
    n.args = copy_exprargs(e.args)
    return n
end

# copy parts of an AST that the compiler mutates
copy_exprs(@nospecialize(x)) = x
copy_exprs(x::Expr) = copy(x)
function copy_exprs(x::PhiNode)
    new_values = Vector{Any}(undef, length(x.values))
    for i = 1:length(x.values)
        isassigned(x.values, i) || continue
        new_values[i] = copy_exprs(x.values[i])
    end
    return PhiNode(copy(x.edges), new_values)
end
function copy_exprs(x::PhiCNode)
    new_values = Vector{Any}(undef, length(x.values))
    for i = 1:length(x.values)
        isassigned(x.values, i) || continue
        new_values[i] = copy_exprs(x.values[i])
    end
    return PhiCNode(new_values)
end
copy_exprargs(x::Array{Any,1}) = Any[copy_exprs(x[i]) for i in 1:length(x)]

==(x::Expr, y::Expr) = x.head === y.head && isequal(x.args, y.args)
==(x::QuoteNode, y::QuoteNode) = isequal(x.value, y.value)

"""
    macroexpand(m::Module, x; recursive=true)

Take the expression `x` and return an equivalent expression with all macros removed (expanded)
for executing in module `m`.
The `recursive` keyword controls whether deeper levels of nested macros are also expanded.
This is demonstrated in the example below:
```julia-repl
julia> module M
           macro m1()
               42
           end
           macro m2()
               :(@m1())
           end
       end
M

julia> macroexpand(M, :(@m2()), recursive=true)
42

julia> macroexpand(M, :(@m2()), recursive=false)
:(#= REPL[16]:6 =# M.@m1)
```
"""
function macroexpand(m::Module, @nospecialize(x); recursive=true)
    if recursive
        ccall(:jl_macroexpand, Any, (Any, Any), x, m)
    else
        ccall(:jl_macroexpand1, Any, (Any, Any), x, m)
    end
end

"""
    @macroexpand

Return equivalent expression with all macros removed (expanded).

There are differences between `@macroexpand` and [`macroexpand`](@ref).

* While [`macroexpand`](@ref) takes a keyword argument `recursive`, `@macroexpand`
is always recursive. For a non recursive macro version, see [`@macroexpand1`](@ref).

* While [`macroexpand`](@ref) has an explicit `module` argument, `@macroexpand` always
expands with respect to the module in which it is called.
This is best seen in the following example:
```julia-repl
julia> module M
           macro m()
               1
           end
           function f()
               (@macroexpand(@m),
                macroexpand(M, :(@m)),
                macroexpand(Main, :(@m))
               )
           end
       end
M

julia> macro m()
           2
       end
@m (macro with 1 method)

julia> M.f()
(1, 1, 2)
```
With `@macroexpand` the expression expands where `@macroexpand` appears in the code (module `M` in the example).
With `macroexpand` the expression expands in the module given as the first argument.
"""
macro macroexpand(code)
    return :(macroexpand($__module__, $(QuoteNode(code)), recursive=true))
end

"""
    @macroexpand1

Non recursive version of [`@macroexpand`](@ref).
"""
macro macroexpand1(code)
    return :(macroexpand($__module__, $(QuoteNode(code)), recursive=false))
end

## misc syntax ##

"""
    Core.eval(m::Module, expr)

Evaluate an expression in the given module and return the result.
"""
Core.eval

"""
    @inline

Give a hint to the compiler that this function is worth inlining.

Small functions typically do not need the `@inline` annotation,
as the compiler does it automatically. By using `@inline` on bigger functions,
an extra nudge can be given to the compiler to inline it.
This is shown in the following example:

```julia
@inline function bigfunction(x)
    #=
        Function Definition
    =#
end
```
"""
macro inline(ex)
    esc(isa(ex, Expr) ? pushmeta!(ex, :inline) : ex)
end

"""
    @noinline

Prevents the compiler from inlining a function.

Small functions are typically inlined automatically.
By using `@noinline` on small functions, auto-inlining can be
prevented. This is shown in the following example:

```julia
@noinline function smallfunction(x)
    #=
        Function Definition
    =#
end
```
"""
macro noinline(ex)
    esc(isa(ex, Expr) ? pushmeta!(ex, :noinline) : ex)
end

"""
    @pure ex
    @pure(ex)

`@pure` gives the compiler a hint for the definition of a pure function,
helping for type inference.

A pure function can only depend on immutable information.
This also means a `@pure` function cannot use any global mutable state, including
generic functions. Calls to generic functions depend on method tables which are
mutable global state.
Use with caution, incorrect `@pure` annotation of a function may introduce
hard to identify bugs. Double check for calls to generic functions.
"""
macro pure(ex)
    esc(isa(ex, Expr) ? pushmeta!(ex, :pure) : ex)
end

"""
    @propagate_inbounds

Tells the compiler to inline a function while retaining the caller's inbounds context.
"""
macro propagate_inbounds(ex)
    if isa(ex, Expr)
        pushmeta!(ex, :inline)
        pushmeta!(ex, :propagate_inbounds)
        esc(ex)
    else
        esc(ex)
    end
end

"""
    @polly

Tells the compiler to apply the polyhedral optimizer Polly to a function.
"""
macro polly(ex)
    esc(isa(ex, Expr) ? pushmeta!(ex, :polly) : ex)
end

## some macro utilities ##

function pushmeta!(ex::Expr, sym::Symbol, args::Any...)
    if isempty(args)
        tag = sym
    else
        tag = Expr(sym, args...)
    end

    inner = ex
    while inner.head == :macrocall
        inner = inner.args[end]::Expr
    end

    idx, exargs = findmeta(inner)
    if idx != 0
        push!(exargs[idx].args, tag)
    else
        body::Expr = inner.args[2]
        pushfirst!(body.args, Expr(:meta, tag))
    end
    ex
end

popmeta!(body, sym) = _getmeta(body, sym, true)
peekmeta(body, sym) = _getmeta(body, sym, false)

function _getmeta(body::Expr, sym::Symbol, delete::Bool)
    body.head == :block || return false, []
    _getmeta(body.args, sym, delete)
end
_getmeta(arg, sym, delete::Bool) = (false, [])
function _getmeta(body::Array{Any,1}, sym::Symbol, delete::Bool)
    idx, blockargs = findmeta_block(body, args -> findmetaarg(args,sym)!=0)
    if idx == 0
        return false, []
    end
    metaargs = blockargs[idx].args
    i = findmetaarg(blockargs[idx].args, sym)
    if i == 0
        return false, []
    end
    ret = isa(metaargs[i], Expr) ? (metaargs[i]::Expr).args : []
    if delete
        deleteat!(metaargs, i)
        isempty(metaargs) && deleteat!(blockargs, idx)
    end
    true, ret
end

# Find index of `sym` in a meta expression argument list, or 0.
function findmetaarg(metaargs, sym)
    for i = 1:length(metaargs)
        arg = metaargs[i]
        if (isa(arg, Symbol) && (arg::Symbol)    == sym) ||
           (isa(arg, Expr)   && (arg::Expr).head == sym)
            return i
        end
    end
    return 0
end

function is_short_function_def(ex)
    ex.head == :(=) || return false
    while length(ex.args) >= 1 && isa(ex.args[1], Expr)
        (ex.args[1].head == :call) && return true
        (ex.args[1].head == :where || ex.args[1].head == :(::)) || return false
        ex = ex.args[1]
    end
    return false
end

function findmeta(ex::Expr)
    if ex.head == :function || is_short_function_def(ex)
        body::Expr = ex.args[2]
        body.head == :block || error(body, " is not a block expression")
        return findmeta_block(ex.args)
    end
    error(ex, " is not a function expression")
end

findmeta(ex::Array{Any,1}) = findmeta_block(ex)

function findmeta_block(exargs, argsmatch=args->true)
    for i = 1:length(exargs)
        a = exargs[i]
        if isa(a, Expr)
            if (a::Expr).head == :meta && argsmatch((a::Expr).args)
                return i, exargs
            elseif (a::Expr).head == :block
                idx, exa = findmeta_block(a.args, argsmatch)
                if idx != 0
                    return idx, exa
                end
            end
        end
    end
    return 0, []
end

remove_linenums!(ex) = ex
function remove_linenums!(ex::Expr)
    if ex.head === :block || ex.head === :quote
        # remove line number expressions from metadata (not argument literal or inert) position
        filter!(ex.args) do x
            isa(x, Expr) && x.head === :line && return false
            isa(x, LineNumberNode) && return false
            return true
        end
    end
    for subex in ex.args
        remove_linenums!(subex)
    end
    return ex
end

macro generated()
    return Expr(:generated)
end

"""
    @generated f
    @generated(f)
`@generated` is used to annotate a function which will be generated.
In the body of the generated function, only types of arguments can be read
(not the values). The function returns a quoted expression evaluated when the
function is called. The `@generated` macro should not be used on functions mutating
the global scope or depending on mutable elements.

See [Metaprogramming](@ref) for further details.

## Example:
```julia
julia> @generated function bar(x)
           if x <: Integer
               return :(x ^ 2)
           else
               return :(x)
           end
       end
bar (generic function with 1 method)

julia> bar(4)
16

julia> bar("baz")
"baz"
```
"""
macro generated(f)
    if isa(f, Expr) && (f.head === :function || is_short_function_def(f))
        body = f.args[2]
        lno = body.args[1]
        return Expr(:escape,
                    Expr(f.head, f.args[1],
                         Expr(:block,
                              lno,
                              Expr(:if, Expr(:generated),
                                   body,
                                   Expr(:block,
                                        Expr(:meta, :generated_only),
                                        Expr(:return, nothing))))))
    else
        error("invalid syntax; @generated must be used with a function definition")
    end
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
function LinearAlgebra.svd(a::IndicesArray{T}, args...; kwargs...) where {T}
    return svd!(
        LinearAlgebra.copy_oftype(a, LinearAlgebra.eigtype(T)),
        args...;
        kwargs...
    )
end

function LinearAlgebra.svd!(a::IndicesArray, args...; kwargs...)
    inner = svd!(parent(nda), args...; kwargs...)
    return SVD(
        IndicesArray(getfield(inner, :U), axes(a)),
        getfield(inner, :S),
        IndicesArray(getfield(inner, :Vt), axes(a))
       )
end

function Base.parent(fact::SVD{T,Tr,<:IndicesArray}) where {T,Tr}
    return SVD(parent(getfield(fact, :U)), getfield(fact, :S), parent(getfield(fact, :Vt)))
end

function Base.getproperty(fact::SVD{T,Tr,<:IndicesArray}, d::Symbol) where {T,Tr}
    inner = getproperty(parent(fact), d)
    axs = axes(fact)
    if d === :U
        return IndicesArray(inner, (unsafe_axes(a, 1), nothing), AllUnique, LengthChecked)
    elseif d === :V
        return IndicesArray(inner, (nothing, unsafe_axes(a, 2)), AllUnique, LengthChecked)
    elseif d === :Vt
        return IndicesArray(inner, (unsafe_axes(a, 2), nothing), AllUnique, LengthChecked)
    else # d === :S
        return inner
    end
end


text:
module HalideCall

export HalideBuffer

const libname = "libhalidecall.so"
const libHalideCall = find_library([libname], [Pkg.dir("HalideCall", "deps")])
isempty(libHalideCall) && error("Can't find $libname")

immutable HalideBuffer
    dev::Uint64
    host::Ptr{Uint8}
    extent0::Cint
    extent1::Cint
    extent2::Cint
    extent3::Cint
    stride0::Cint
    stride1::Cint
    stride2::Cint
    stride3::Cint
    min0::Cint
    min1::Cint
    min2::Cint
    min3::Cint
    elem_size::Cint
    host_dirty::Bool
    dev_dirty::Bool
end
HalideBuffer(A::Array) = HalideBuffer(A, zeros(Int,4))
HalideBuffer(A::Array, pad) = HalideBuffer(A, pad, pad)
function HalideBuffer(A::Array, padmin, padmax)
    sz = growto4(size(A))
    st = growto4(strides(A))
    pmn = growto4(padmin)
    pmx = growto4(padmax)
    extent = sz - pmn - pmx
    HalideBuffer(0, pointer(A),
                 extent[1], extent[2], extent[3], extent[4],
                 st[1], st[2], st[3], st[4],
                 pmn[1], pmn[2], pmn[3], pmn[4],
                 sizeof(eltype(A)), false, false)
end

growto4(v::Vector) = growto4!(copy(v))
growto4(t::Tuple) = growto4!([t...])
function growto4!(v::Vector)
    length(v) <= 4 || error("Vector is too large")
    while length(v) < 4
        push!(v, 0)
    end
    v
end

end # module


text:
const RELAXED_NLSOLVE_F_TOLERANCE = :1e-4
const STRICT_NLSOLVE_F_TOLERANCE = :1e-7
const NLSOLVE_X_TOLERANCE = :1e-9
const MINIMAL_ACCEPTABLE_NLSOLVE_F_TOLERANCE = :1e-3
const MAX_INIT_RESIDUAL = 10.0
const MAX_NLSOLVE_INTERATIONS = 10
const BOUNDS_TOLERANCE = 1e-6

const ISAPPROX_ZERO_TOLERANCE = 1e-6


text:
using Base.Test
hd = homedir()
pd = Pkg.dir()

using RCall

# https://github.com/JuliaStats/RCall.jl/issues/68
@test hd == homedir()
@test pd == Pkg.dir()

tests = ["basic",
         "conversion",
         "dataframe",
         "rstr",
         "library",
         "repl",
         ]

println("Running tests:")

for t in tests
    tfile = string(t, ".jl")
    println(" * $(tfile) ...")
    include(tfile)
end


text:
@testset "cuda/conv" begin
    T = Float32
    in_channel = 3
    out_channel = 5
    N = 4
    adj = T[0 1 0 1;
           1 0 1 0;
           0 1 0 1;
           1 0 1 0]

    fg = FeaturedGraph(adj)

    @testset "GCNConv" begin
        gc = GCNConv(fg, in_channel=>out_channel) |> gpu
        @test size(gc.weight) == (out_channel, in_channel)
        @test size(gc.bias) == (out_channel,)
        @test collect(GraphSignals.adjacency_matrix(gc.fg)) == adj

        X = rand(in_channel, N) |> gpu
        Y = gc(X)
        @test size(Y) == (out_channel, N)

        g = Zygote.gradient(() -> sum(gc(X)), Flux.params(gc))
        @test length(g.grads) == 2
    end

    @testset "ChebConv" begin
        k = 6
        cc = ChebConv(fg, in_channel=>out_channel, k) |> gpu
        @test size(cc.weight) == (out_channel, in_channel, k)
        @test size(cc.bias) == (out_channel,)
        @test collect(GraphSignals.adjacency_matrix(cc.fg)) == adj
        @test cc.k == k
        @test size(cc.weight, 2) == in_channel
        @test size(cc.weight, 1) == out_channel

        X = rand(in_channel, N) |> gpu
        Y = cc(X)
        @test size(Y) == (out_channel, N)

        g = Zygote.gradient(() -> sum(cc(X)), Flux.params(cc))
        @test length(g.grads) == 2
    end

    @testset "GraphConv" begin
        gc = GraphConv(fg, in_channel=>out_channel) |> gpu
        @test size(gc.weight1) == (out_channel, in_channel)
        @test size(gc.weight2) == (out_channel, in_channel)
        @test size(gc.bias) == (out_channel,)

        X = rand(in_channel, N) |> gpu
        Y = gc(X)
        @test size(Y) == (out_channel, N)

        g = Zygote.gradient(() -> sum(gc(X)), Flux.params(gc))
        @test length(g.grads) == 3
    end

    @testset "GATConv" begin
        adj = T[1 1 0 1;
                1 1 1 0;
                0 1 1 1;
                1 0 1 1]

        fg = FeaturedGraph(adj)

        gat = GATConv(fg, in_channel=>out_channel) |> gpu
        @test size(gat.weight) == (out_channel, in_channel)
        @test size(gat.bias) == (out_channel,)

        X = rand(in_channel, N) |> gpu
        Y = gat(X)
        @test size(Y) == (out_channel, N)

        g = Zygote.gradient(() -> sum(gat(X)), Flux.params(gat))
        @test length(g.grads) == 3
    end

    @testset "GatedGraphConv" begin
        num_layers = 3
        ggc = GatedGraphConv(fg, out_channel, num_layers) |> gpu
        @test size(ggc.weight) == (out_channel, out_channel, num_layers)

        X = rand(in_channel, N) |> gpu
        Y = ggc(X)
        @test size(Y) == (out_channel, N)

        g = Zygote.gradient(() -> sum(ggc(X)), Flux.params(ggc))
        @test length(g.grads) == 6
    end

    @testset "EdgeConv" begin
        ec = EdgeConv(fg, Dense(2*in_channel, out_channel)) |> gpu
        X = rand(in_channel, N) |> gpu
        Y = ec(X)
        @test size(Y) == (out_channel, N)

        g = Zygote.gradient(() -> sum(ec(X)), Flux.params(ec))
        @test length(g.grads) == 2
    end
end


text:
using ModelConstructors, FileIO, Random, SMC, JLD2

### Construct a generic model and populate it with parameters
reg = GenericModel()
fn = dirname(@__FILE__)
reg <= Setting(:dataroot, "$(fn)/save/input_data/")
reg <= Setting(:saveroot, "$(fn)/save/")

reg <= parameter(:α1, 0., (-1e5, 1e5), (-1e5, 1e5), Untransformed(), Normal(0, 10), fixed = false)
reg <= parameter(:β1, 0., (-1e5, 1e5), (-1e5, 1e5), Untransformed(), Normal(0, 10), fixed = false)

### Estimate with SMC

## Make data
N = 100 # number observations
M = 1 # number regressors

X = rand(N)
β = 1.
α = 1.
σ2 = 1.
y = β*X .+ α #.+ randn(N) .* sqrt(σ2)

jldopen(dataroot(reg)*"reg_data.jld2", "w") do file
    file["data"] = hcat(y, X)
end

function compute_coefficients(X::Vector, Y::Vector)
    @assert length(X) == length(Y)
    X_aug = hcat(ones(length(X)), X)
    β     = inv(X_aug'* X_aug)*X_aug'*Y
    return β
end

function compute_coefficients(X::Array, Y::Vector)
    @assert size(X)[1] == length(Y)
    X_aug = hcat(ones(size(X)[1], 1), X)
    β     = inv(X_aug'* X_aug)*X_aug'*Y
    return β
end

compute_coefficients(X, y)

## Construct log-likelihood function:
function log_likelihood_fnct(p, y)
    α = p[1]
    β = p[2]
    term1 = -(N/2)*log(2*π) - (N/2)*log(σ2)
    errors = y .- α .- β .* X
    logprob = term1 - (1/(2*σ2))*dot(errors, errors)
    return logprob
end

Random.seed!(1793)
println("Starting to estimate Regression with SMC . . .")
@everywhere using SMC
smc(log_likelihood_fnct, reg.parameters, Matrix(Matrix{Float64}((y'))'), n_parts = 100, use_fixed_schedule = true) #, tempering_target = 0.99)

cloud = load("smc_cloud.jld2", "cloud")
mean(SMC.get_vals(cloud), dims = 2)


text:
#TODO a weight for y_j based on how far away from best ord it is.

function setup_sub!(subproblems,patients::IndexedTable,visits,resources::IndexedTable,timeslots,mastercalendar::Dict{Int,Date},TimeDelta,sets,months)

    for patientgroup in [i for i in 1:length(sets.Pg)]

        #TODO length of min treatment time should be calculated and used.

        setup_sub!(subproblems,patientgroup,visits,resources,timeslots,mastercalendar,TimeDelta,sets,months)
    end
end
startofyear(date::Date) = Date(Year(date))
endofyear(date::Date) = Date(Year(date),Month(12),Day(31))
endofmonth(date::Date) = date + Month(1)-Day(date+Month(1))
startofmonth(date::Date) = date-Day(date)+Day(1)

function getTimeDelta(timeDelta,visits,v1,v2)
    k = (visits[v1].req_type,visits[v2].req_type)
    if haskey(timeDelta,k)
        timeDelta[k[1],k[2]].delta
    else
        -1
    end
end

function setup_sub!(subproblems,patientgroup::Int64,visits,resources::IndexedTable,timeslots,mastercalendar::Dict{Int,Date},timeDelta,sets,months)
    M1 = 10
    M2 = 1000

    Td(v1,v2) = getTimeDelta(timeDelta,visits,v1,v2)
    p = sets.Pg[patientgroup].patients[1]
    V = sets.Vp[p].v; Dv = sets.Dv; J = sets.J ;Jd = sets.Jd; I = sets.I
    Ts(d,j,i) = Dates.value(timeslots[d,j,i].startTime)/60000000000
    Te(d,j,i) = Dates.value(timeslots[d,j,i].endTime)/60000000000

    (bestordmin,bestordmax) = reduce((min,max),visits[V],select=:bestord)
    startdate = max(minimum(mastercalendar)[2],min(startofmonth(maximum(mastercalendar)[2]-Month(months)),startofmonth(mastercalendar[bestordmin])-Month(div(months,2))))
    enddate = startdate + Month(months)-Day(1)
    subMastercalendar = MasterCalendar(mastercalendar,startdate,enddate)
    if p== 6
        println("V = $V")
        for v in V
            println(Dv[v].d)
        end
    end

        sub = Model(with_optimizer(Gurobi.Optimizer,OutputFlag=0))
        @variable(sub,xvars[v in V,d in Dv[v].d,j in Jd[d].j,i in I[d,j].i],Bin)
        @variable(sub,yvars[j = J],Bin)
        #@variable(sub,tvars[d = D,j = Jd[d].j])  #TODO Only for consultations
        @variable(sub,kvars[V,V], Bin)
        @variable(sub,gvars[V,V], Bin)
        for v1 in V, v2 in V
            if Td(v1,v2) >= 0
                fix(kvars[v1,v2],1)
            end
        end

        @objective(sub,Min, 0)

        #    @constraint(sub,lastesttime[d =D ,j = Jd[d].j,i = I[d,j].i],
        #        sum(xvars[v,d,j,i] for v in V if d in Dv[v].d && j in Jd[d].j && i in I[d,j].i) *Te(d,j,i) <= tvars[d,j])#TODO: somethings wrong here

            @constraint(sub,demand[v in V],
                sum(xvars[v,d,j,i] for d in Dv[v].d,j in Jd[d].j,i in I[d,j].i) == 1)

            @constraint(sub,openday[j in J],
                sum(xvars[v,d,j2,i] for v in V, d in Dv[v].d, j2 in Jd[d].j, i in I[d,j2].i if j2 ==j) <= M1 * yvars[j] )

            @constraint(sub,sameday[v1 in V, v2 in V,j in J; v1 != v2 ],
                    sum(xvars[v1,d,j2,i] for d in Dv[v1].d,j2 in Jd[d].j, i in I[d,j2].i if j2 == j)+sum(xvars[v2,d,j2,i] for d in Dv[v2].d,j2 in Jd[d].j, i in I[d,j2].i if j2 == j) <= gvars[v1,v2]+ 1 )

            @constraint(sub,beginsBefore[v1 in V, v2 in V,j in J; v1 != v2] ,
                -sum(Ts(d,j2,i)*xvars[v2,d,j2,i] for d in Dv[v2].d,j2 in Jd[d].j,i in I[d,j2].i if j2 == j )
                + sum(Te(d,j2,i)*xvars[v1,d,j2,i] for d in Dv[v1].d,j2 in Jd[d].j,i in I[d,j2].i if j2 == j)
                <=  (1-gvars[v1,v2])*M2+(1-kvars[v1,v2])*M2)

            @constraint(sub, numberofrelations, sum(kvars[v2,v1] for v1 in V, v2 in V if v1 != v2 ) == length(V)*(length(V)-1)/2)

            @constraint(sub,timebetweenappointments[v1 in V, v2 in V;v1 != v2 ],
                sum(j*xvars[v2,d,j,i] for d in Dv[v2].d,j in Jd[d].j,i in I[d,j].i )
                - sum(j*xvars[v1,d,j,i] for d in Dv[v1].d,j in Jd[d].j,i in I[d,j].i )
                >= Td(v1,v2)
                -(1-kvars[v1,v2])*M2)

        addPricingProblem(subproblems,patientgroup,sub,xvars,yvars,[],kvars,gvars,p,) ##TODO stop adding all that stuff

end

function solveSub!(sub,xStar)
    @constraint(sub.model,sum(sub.xvars[sol] for sol in xStar) <= length(xStar)-1)
    optimize!(sub.model)
    status = termination_status(sub.model)
    if status != MOI.TerminationStatusCode(1)
        throw("Error: Non optimal sub-problem")
    end

    return objective_value(sub.model)
end

function solveSub!(sub,sets,ϕ,θ,κ)
    p = sets.Pg[sub.intID].patients[1]

    #@objective(sub.model, Min, sum(1000*sub.yvars[j] for j in sets.J) - sum(sub.tvars[d,j]*ϕ[d,j] for d in sets.Dp[p].d,j in sets.Jd[d].j,i in sets.I[d,j].i ) - sum(sub.xvars[v,d,j,i] * θ[d,j,i] for v in sets.Vp[p].v, d in sets.Dv[v].d,j in sets.Jd[d].j,i in sets.I[d,j].i)-κ[sub.intID] )
    @objective(sub.model, Min, sum(1000*sub.yvars[j] for j in sets.J)  - sum(sub.xvars[v,d,j,i] * θ[d,j,i] for v in sets.Vp[p].v, d in sets.Dv[v].d,j in sets.Jd[d].j,i in sets.I[d,j].i)-κ[sub.intID] )

    optimize!(sub.model)
    status = termination_status(sub.model)
    if status != MOI.TerminationStatusCode(1)

        throw("Error: Non optimal sub-problem for subproblem $(sub.intID)")

        #TODO Throw warning and extend search area
    end

    return objective_value(sub.model)
end


text:
using CSV, Distributions, HypothesisTests

data1 = CSV.read("../data/machine1.csv", header=false)[:,1]
data2 = CSV.read("../data/machine2.csv", header=false)[:,1]
xBar1, xBar2 = mean(data1), mean(data2)
s1, s2 = std(data1), std(data2)
n1, n2 = length(data1), length(data2)
alpha = 0.05

v = (s1^2/n1 + s2^2/n2)^2 / ( (s1^2/n1)^2 / (n1-1) + (s2^2/n2)^2 / (n2-1) )

t = quantile(TDist(v),1-alpha/2)

println("Calculating formula: ", (xBar1 - xBar2 - t*sqrt(s1^2/n1 + s2^2/n2),
				  xBar1 - xBar2 + t*sqrt(s1^2/n1 + s2^2/n2)))
println("Using confint(): ",   confint(UnequalVarianceTTest(data1,data2),alpha))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
module Environment

using RessourceManager
using FileManager

#using Logging
#@Logging.configure(level=DEBUG)

#export LOGGERS

#export CreateLogger
#export LogTest
#export StartMainLogger

#LogLevel = Logging.DEBUG
#LOGGERS = Dict{Any,Logger}()

#function (d::Dict, s::Symbol, f::Function, args...)
#  ok = true
#	for (k,e) in d
#     ok = f(e,s,args...)
#     if !ok break end
#  end
#	ok
#end

# EVENTS = [ :OnWindowIconify, :OnWindowResize, :OnWindowClose, :OnWindowFocus,
# :OnWindowPos, :OnWindowMove, :OnWindowRefresh, :OnFramebufferResize,
# :OnCursorPos, :OnCursorMove, :OnCursorEnter, :OnDroppedFiles,
# :OnKey, :OnMouseKey, :OnScroll, :OnCharMods ]

#StartMainLogger() = Logging.configure(output=open(string(GetPath(:LOGS),"log","_",CurrentDay(),".log"),"w"))

#function ShowLog(show)
  #lvl = if log == show LogLevel else Logging.OFF end
  #Logging.configure(level=lvl)
#end

#function CreateLogger(id::Any, path::AbstractString, lvl=DEBUG)
  #(dir, name, ext) = SplitFileName(path, GetPath(:LOGS))
  #logger = Logger(name)
  #Logging.configure(logger, level=lvl, output=open(string(dir,name,"_",CurrentDay(),ext), "w"))
  #LOGGERS[id] = logger
  #logger
	#nothing
#end

#function LogTest()
	#@debug("debug message")
	#@info("info message")
	#@warn("warning message")
	#@err("error message")
	#@critical("critical message")
#end

end # Environment

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
module PopulationModule

import Random: randperm
import ..CoreModule: Options, Dataset, RecordType, string_tree
import ..EquationUtilsModule: compute_complexity
import ..LossFunctionsModule: score_func
import ..MutationFunctionsModule: gen_random_tree
import ..PopMemberModule: PopMember
# A list of members of the population, with easy constructors,
#  which allow for random generation of new populations
mutable struct Population{T<:Real}
    members::Array{PopMember{T},1}
    n::Int
end
"""
    Population(pop::Array{PopMember{T}, 1})

Create population from list of PopMembers.
"""
Population(pop::Array{PopMember{T},1}) where {T<:Real} = Population{T}(pop, size(pop, 1))
"""
    Population(dataset::Dataset{T}, baseline::T;
               npop::Int, nlength::Int=3, options::Options,
               nfeatures::Int)

Create random population and score them on the dataset.
"""
function Population(
    dataset::Dataset{T},
    baseline::T;
    npop::Int,
    nlength::Int=3,
    options::Options,
    nfeatures::Int,
) where {T<:Real}
    return Population(
        [
            PopMember(
                dataset,
                baseline,
                gen_random_tree(nlength, options, nfeatures),
                options;
                parent=-1,
                deterministic=options.deterministic,
            ) for i in 1:npop
        ],
        npop,
    )
end
"""
    Population(X::AbstractMatrix{T}, y::AbstractVector{T},
               baseline::T; npop::Int, nlength::Int=3,
               options::Options, nfeatures::Int)

Create random population and score them on the dataset.
"""
function Population(
    X::AbstractMatrix{T},
    y::AbstractVector{T},
    baseline::T;
    npop::Int,
    nlength::Int=3,
    options::Options,
    nfeatures::Int,
) where {T<:Real}
    return Population(
        Dataset(X, y), baseline; npop=npop, options=options, nfeatures=nfeatures
    )
end

# Sample 10 random members of the population, and make a new one
function sample_pop(pop::Population, options::Options)::Population
    idx = randperm(pop.n)[1:(options.ns)]
    return Population(pop.members[idx])
end

# Sample the population, and get the best member from that sample
function best_of_sample(
    pop::Population, frequencyComplexity::AbstractVector{T}, options::Options
)::PopMember where {T<:Real}
    sample = sample_pop(pop, options)

    if options.useFrequencyInTournament
        # Score based on frequency of that size occuring.
        # In the end, all sizes should be just as common in the population.
        frequency_scaling = 20
        # e.g., for 100% occupied at one size, exp(-20*1) = 2.061153622438558e-9; which seems like a good punishment for dominating the population.

        scores = []
        for member in 1:(options.ns)
            size = compute_complexity(sample.members[member].tree, options)
            frequency = (size <= options.maxsize) ? frequencyComplexity[size] : T(0)
            score = sample.members[member].score * exp(frequency_scaling * frequency)
            push!(scores, score)
        end
    else
        scores = [sample.members[member].score for member in 1:(options.ns)]
    end

    p = options.probPickFirst

    if p == 1.0
        chosen_idx = argmin(scores)
    else
        sort_idx = sortperm(scores)
        # scores[sort_idx] would put smallest first.

        k = collect(0:(options.ns - 1))
        prob_each = p * (1 - p) .^ k
        prob_each /= sum(prob_each)
        cumprob = cumsum(prob_each)
        raw_chosen_idx = findfirst(cumprob .> rand())

        # Sometimes, due to precision issues, we might have cumprob[end] < 1,
        # so we must check for nothing returned:
        if raw_chosen_idx === nothing
            chosen_idx = sort_idx[end]
        else
            chosen_idx = sort_idx[raw_chosen_idx]
        end
    end
    return sample.members[chosen_idx]
end

function finalize_scores(
    dataset::Dataset{T}, baseline::T, pop::Population, options::Options
)::Tuple{Population,Float64} where {T<:Real}
    need_recalculate = options.batching
    num_evals = 0.0
    if need_recalculate
        @inbounds @simd for member in 1:(pop.n)
            score, loss = score_func(dataset, baseline, pop.members[member].tree, options)
            pop.members[member].score = score
            pop.members[member].loss = loss
        end
        num_evals += pop.n * (options.batchSize / dataset.n)
    end
    return (pop, num_evals)
end

# Return best 10 examples
function best_sub_pop(pop::Population; topn::Int=10)::Population
    best_idx = sortperm([pop.members[member].score for member in 1:(pop.n)])
    return Population(pop.members[best_idx[1:topn]])
end

function record_population(pop::Population{T}, options::Options)::RecordType where {T<:Real}
    return RecordType(
        "population" => [
            RecordType(
                "tree" => string_tree(member.tree, options),
                "loss" => member.loss,
                "score" => member.score,
                "complexity" => compute_complexity(member.tree, options),
                "birth" => member.birth,
                "ref" => member.ref,
                "parent" => member.parent,
            ) for member in pop.members
        ],
        "time" => time(),
    )
end

end


text:
using Measures

# 1 function to be moved to Measures.jl
Base.convert(::Type{T} , x::Measure) where T<:Real  = T(x.value)

struct DiscreteSizeScale <: Gadfly.ScaleElement
    f::Function
    levels::Union{Nothing, AbstractVector}
    order::Union{Nothing, AbstractVector}
    preserve_order::Bool
end
DiscreteSizeScale(f; levels=nothing, order=nothing, preserve_order=true) =
        DiscreteSizeScale(f, levels, order, preserve_order)

element_aesthetics(scale::DiscreteSizeScale) = [:size]

default_discrete_sizes(n::Int) = range(0.45mm, 1.8mm, length=n)

"""
    size_discrete2(f=default_discrete_sizes;
                    levels=nothing, order=nothing, preserve_order=true)

A discrete size scale that maps the categorical values in the `size`
aesthetic to x-axis units or `Measure` units (from Measures.jl).  `f` is a function that produces a vector of size units.
`levels` are the categorical levels, and level order will be respected.  `order` is
a vector of integers giving a permutation of the levels default order.  If
`preserve_order` is `true` levels are ordered as they appear in the data.
"""
size_discrete2(f::Function=Gadfly.current_theme().discrete_sizemap; levels=nothing, order=nothing, preserve_order=true) =
        DiscreteSizeScale(f, levels, order, preserve_order)

function Scale.apply_scale(scale::DiscreteSizeScale, aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)

    d = []
    for (aes, data) in zip(aess, datas)
        data.size === nothing && continue
        append!(d, skipmissing(data.size))
    end
    levelset = unique(d)

    if scale.levels == nothing
        scale_levels = [levelset...]
        scale.preserve_order || sort!(scale_levels)
    else
        scale_levels = scale.levels
    end
    scale.order == nothing || permute!(scale_levels, scale.order)

    sizes = scale.f(length(scale_levels))

    size_map = Dict(s=>string(label) for (s, label) in zip(sizes, scale_levels))
    labeler(xs) = [size_map[x] for x in xs]
    key_vals = OrderedDict(s=>i for (i,s) in enumerate(sizes))

    for (aes, data) in zip(aess, datas)
        data.size === nothing && continue
        ds = discretize([d for d in skipmissing(data.size)], scale_levels)
        vals = sizes[ds.index]
        aes.size = discretize_make_ia(vals, sizes)
        aes.size_label = labeler
        aes.size_key_vals = key_vals
    end
end


text:
LarA = Lar.Arrangement

using IntervalTrees
using SparseArrays
using NearestNeighbors
using DataStructures

#-------------------------------------------------------------------------------
#	TSAS PIPELINE
#    input_collection
#    spatial_indexing
#    pairwise_decomposition
#    equivalence_congruence
#    connection_extraction
#    topological_gift_wrapping
#    boundaries_evaluation
#    holes_cotainments
#    transitive_reduction
#    cell_adjoining
#    cell_assembling
#    colleclt_output
#-------------------------------------------------------------------------------