text:
if !("." in LOAD_PATH)
    push!(LOAD_PATH, ".")
end
import Coverage
import Plots

# Track allocations in ClimaCore.jl plus all _direct_ dependencies
exhaustive = "exhaustive=true" in ARGS
@show exhaustive

# Packages to monitor
import ClimaCore
import Adapt
import BlockArrays
import CUDA
import CubedSphere
import DiffEqBase
import DocStringExtensions
import ForwardDiff
import GaussQuadrature
import IntervalSets
import LinearAlgebra
import RecursiveArrayTools
import Rotations
import SparseArrays
import Static
import StaticArrays
import UnPack

mod_dir(x) = dirname(dirname(pathof(x)))
deps_to_monitor = [
    mod_dir(Adapt),
    mod_dir(BlockArrays),
    mod_dir(CUDA),
    mod_dir(CubedSphere),
    mod_dir(DiffEqBase),
    mod_dir(DocStringExtensions),
    mod_dir(ForwardDiff),
    mod_dir(GaussQuadrature),
    mod_dir(IntervalSets),
    mod_dir(LinearAlgebra),
    mod_dir(RecursiveArrayTools),
    mod_dir(Rotations),
    mod_dir(SparseArrays),
    mod_dir(Static),
    mod_dir(StaticArrays),
    mod_dir(UnPack),
]
deps_to_monitor = exhaustive ? deps_to_monitor : ()
pkg_dir = mod_dir(ClimaCore)

julia_base_dir = "/central/software/julia/1.7.0/share/julia/base/"
julia_base_dir = isdir(julia_base_dir) ? (julia_base_dir,) : ()
all_dirs_to_monitor =
    [julia_base_dir..., mod_dir(ClimaCore), deps_to_monitor...]

# https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
ENV["GKSwstype"] = "nul"

# (filename, ARGs passed to script)
all_cases = [
    (joinpath(pkg_dir, "examples", "hybrid", "bubble_2d.jl"), ""),
    (joinpath(pkg_dir, "examples", "hybrid", "bubble_3d.jl"), ""),
    (
        joinpath(pkg_dir, "examples", "3dsphere", "baroclinic_wave.jl"),
        "baroclinic_wave",
    ),
    (
        joinpath(pkg_dir, "examples", "sphere", "shallow_water.jl"),
        "barotropic_instability",
    ),
    (
        joinpath(pkg_dir, "examples", "3dsphere", "solid_body_rotation_3d.jl"),
        "",
    ),
]

# only one case for exhaustive alloc analysis
all_cases = if exhaustive
    [(joinpath(pkg_dir, "examples", "hybrid", "bubble_2d.jl"), "")]
else
    all_cases
end

allocs = Dict()
for (case, args) in all_cases
    ENV["ALLOCATION_CASE_NAME"] = case
    if exhaustive
        # run(`julia --project=examples/ --track-allocation=all perf/alloc_per_case_with_init_io.jl`)
        error("Not yet supported")
    else
        run(
            `julia --project=examples/ --track-allocation=all perf/allocs_per_case.jl $args`,
        )
    end

    allocs[case] = Coverage.analyze_malloc(all_dirs_to_monitor)

    # Clean up files
    for d in all_dirs_to_monitor
        all_files = [
            joinpath(root, f) for
            (root, dirs, files) in Base.Filesystem.walkdir(d) for f in files
        ]
        all_mem_files = filter(x -> endswith(x, ".mem"), all_files)
        for f in all_mem_files
            rm(f)
        end
    end
end

@info "Post-processing allocations"

const PKG_NAME = "ClimaCore.jl"

function plot_allocs(
    folder,
    case_name,
    allocs_per_case,
    script_args,
    n_unique_bytes,
)
    p = Plots.plot()
    case_name_args =
        script_args == "" ? case_name : case_name * "_$(script_args)"
    @info "Allocations for $case_name_args"

    function filename_only(fn)
        if occursin(".jl", fn)
            fn = join(split(fn, ".jl")[1:(end - 1)], ".jl") * ".jl"
        end
        fn = replace(fn, "\\" => "/") # for windows...
        splitby = PKG_NAME
        if occursin(splitby, fn)
            fn = PKG_NAME * last(split(fn, splitby))
        end
        splitby = "climacore-ci/"
        if occursin(splitby, fn)
            fn = PKG_NAME * last(split(fn, splitby))
        end
        return fn
    end
    function compile_pkg(fn, linenumber)
        c1 = endswith(filename_only(fn), PKG_NAME)
        c2 = linenumber == 1
        return c1 && c2
    end

    filter!(x -> x.bytes ≠ 0, allocs_per_case)
    filter!(x -> !compile_pkg(x.filename, x.linenumber), allocs_per_case)

    for alloc in allocs_per_case
        println(alloc)
    end
    println("Number of allocating sites: $(length(allocs_per_case))")
    case_bytes = getproperty.(allocs_per_case, :bytes)[end:-1:1]
    case_filename = getproperty.(allocs_per_case, :filename)[end:-1:1]
    case_linenumber = getproperty.(allocs_per_case, :linenumber)[end:-1:1]
    all_bytes = Int[]
    filenames = String[]
    linenumbers = Int[]
    loc_ids = String[]
    for (bytes, filename, linenumber) in
        zip(case_bytes, case_filename, case_linenumber)
        compile_pkg(filename, linenumber) && continue # Skip loading module
        loc_id = "$(filename_only(filename))" * "$linenumber"
        if !(bytes in all_bytes) && !(loc_id in loc_ids)
            push!(all_bytes, bytes)
            push!(filenames, filename)
            push!(linenumbers, linenumber)
            push!(loc_ids, loc_id)
            if length(all_bytes) ≥ n_unique_bytes
                break
            end
        end
    end

    all_bytes = all_bytes ./ 10^3
    max_bytes = maximum(all_bytes)
    @info "$case_name_args: $all_bytes"
    xtick_name(filename, linenumber) = "$filename, line number: $linenumber"
    markershape = (:square, :hexagon, :circle, :star, :utriangle, :dtriangle)
    for (bytes, filename, linenumber) in zip(all_bytes, filenames, linenumbers)
        label = xtick_name(filename_only(filename), linenumber)
        Plots.plot!(
            [0],
            [bytes];
            seriestype = :scatter,
            label = label,
            markershape = markershape[1],
            markersize = 1 + bytes / max_bytes * 10,
        )
        markershape = (markershape[end], markershape[1:(end - 1)]...)
    end
    p1 = Plots.plot!(
        ylabel = "Allocations (KB)",
        title = case_name_args,
        legendfontsize = 6,
    )
    n_subset = min(length(allocs_per_case) - 1, 100)
    subset_allocs_per_case = allocs_per_case[end:-1:(end - n_subset)]
    p2 = Plots.plot(
        1:length(subset_allocs_per_case),
        getproperty.(subset_allocs_per_case, :bytes) ./ 1000;
        xlabel = "i-th allocating line (truncated and sorted)",
        ylabel = "Allocations (KB)",
        markershape = :circle,
    )
    Plots.plot(p1, p2, layout = Plots.grid(2, 1))
    Plots.savefig(joinpath(folder, "allocations_$case_name_args.png"))
end

folder =
    exhaustive ? "perf/allocations_output_exhaustive" :
    "perf/allocations_output"
mkpath(folder)

@info "Allocated bytes for single tendency per case:"
for (case, args) in all_cases
    plot_allocs(
        folder,
        first(split(basename(case), ".jl")),
        allocs[case],
        args,
        10,
    )
end


text:
module ICDataServer

using ICCommon
# using Plots
using JSON
import ZMQ, ODBC
using Base.Cartesian, JLD
using AxisArrays

import AxisArrays: axes
import Base.Cartesian.inlineanonymous
import Base: show, isless, getindex, push!, length, eta
import Base.Collections: PriorityQueue, enqueue!, dequeue!, peek

export newinstrument, updateinstrument, deleteinstrument, listinstruments
export newuser, updateuser, deleteuser, listusers
export newserver, updateserver, deleteserver, listservers
export newnote, listnotes
export newjob, updatejob

const ctx = ZMQ.Context()
const dbsock = ZMQ.Socket(ctx, ZMQ.REP)

include("config.jl")

const dsn = ODBC.DSN(confd["dsn"], confd["username"], confd["password"])
ZMQ.bind(dbsock, confd["dbsock"])

function serve(;debug=false)
    while true
        msg = ZMQ.recv(dbsock)
        out = convert(IOStream, msg)
        seekstart(out)
        d = deserialize(out)
        debug && println(STDOUT, d)
        handle(d)
    end
end

function handle(jr::ICCommon.NewJobRequest)
    job_id, jobsubmit = newjob(dsn; jr.args...)
    ICCommon.ssend(dbsock, (job_id, jobsubmit))
end

function handle(jr::ICCommon.UpdateJobRequest)
    response = try
        updatejob(dsn, jr.job_id; jr.args...)
        true
    catch y
        false   # catch errors so that we can issue a proper reply regardless
    end

    # reply with `true` or `false` for success
    # (a reply must be sent to avoid zmq errors with a REQ/REP pattern)
    # maybe in the future we'd actually pass the error back
    ICCommon.ssend(dbsock, response)
end

handle(r::ICCommon.ListUsersRequest) =
    ICCommon.ssend(dbsock, Array{String}(listusers(dsn)[:username]))

handle(r::ICCommon.ListInstrumentsRequest) =
    ICCommon.ssend(dbsock, listinstruments(dsn))

function listtables(dsn)
    ODBC.query(dsn,
        """
        SELECT table_name FROM information_schema.tables
            WHERE table_type = 'BASE TABLE'
                AND table_schema NOT IN ('pg_catalog', 'information_schema');
        """
    )[:table_name]
end

"""
```
the_nuclear_option(dsn)
```

Delete all tables associated with `ICDataServer`. Probably not a good idea.
"""
function the_nuclear_option(dsn)
    ODBC.execute!(dsn, """DROP TABLE IF EXISTS users, jobs, servers,
        instruments, instrumentkinds, notes CASCADE;""")
end

include("users.jl")
include("servers.jl")
include("instruments.jl")
include("jobs.jl")
include("notes.jl")
include("setup.jl")

end

ext: jl


text:
using NetCDF
using PyPlot
using Dates

path = "/data04/jmg/fsm_simulations/netcdf/forcings_st"

for variable in ["ilwr", "iswr", "pres", "rainf", "snowf", "rhum", "tair", "wind"]

    file_low = joinpath(path, "$(variable)_50km.nc")

    data_low = ncread(file_low, variable)

    file_high = joinpath(path, "$(variable)_1km.nc")

    data_high = ncread(file_high, variable)

    time = ncread(file_low, "time_str")

    time = DateTime.(time, "yyyy-mm-dd HH:MM:SS")

    figure(figsize = (12, 7))
    fill_between(time, minimum(data_high,1)[:], maximum(data_high, dims = 1)[:], facecolor = "gray", alpha = 1.0)
    fill_between(time, minimum(data_low,1)[:], maximum(data_low, dims = 1)[:], facecolor = "red", alpha = 0.3)
    title("$(variable)")

end


text:
#function used for graph coloring in SeaPearl Zoo
function labelOfVertex(g::CPLayerGraph, d::Int64)
    cpVertex = cpVertexFromIndex(g, d)
    labelOfVertex(g, cpVertex)
end

labelOfVertex(g::CPLayerGraph, d::ConstraintVertex) = string(typeof(d.constraint)), 1
labelOfVertex(g::CPLayerGraph, d::VariableVertex) = "x"*d.variable.id, 2
labelOfVertex(g::CPLayerGraph, d::ValueVertex) = string(d.value), 3


text:
using ProjAssistant
@quickactivate "Chemostat_Rath2017"

@time begin
    import DataFrames: DataFrame
    import MAT
    import CSV

    import Chemostat
    import Chemostat.MetNets
    const Ch = Chemostat
    import Chemostat_Rath2017: Human1, RathData
    const Rd = RathData
    const H1 = Human1
    const HG = H1.HumanGEM
end

## ------------------------------------------------------------------
# load model
model = HG.load_humangem_raw_model()
biomass_idx = MetNets.rxnindex(model, HG.HUMAN_BIOMASS_IDER)

## ------------------------------------------------------------------
# Biomass equation

# I will modified the biomass equation (biomass_human) of Human1 model with data
# derived from Niklas (2013): 103–114. https://doi.org/10.1016/j.ymben.2013.01.002. Table1.
# I compute de relation between the total of each group reported in Niklas 2013 with the equivalent
# group found in the model biomass, and then rescaled each group to match the reported total.
# I do not touch the energetic part of the equation, atp + h20 -> adp + h2 + pi

biomass = Dict()
for met_idx in MetNets.rxn_mets(model, biomass_idx)
    met = model.mets[met_idx]
    biomass[met] = model.S[met_idx, biomass_idx]
end

## ------------------------------------------------------------------
# Carbohydrates

ch_ids = ["m03161c"]
exp_ch_tot = 438.3 * 1e-3 # Niklas (2013): 103–114. https://doi.org/10.1016/j.ymben.2013.01.002. Table1
println("experimental total ch: ", exp_ch_tot)
model_ch_tot = abs.(sum([biomass[met] for met in ch_ids])) # The model did not include directly any carbohydrate
println("model total ch: ", model_ch_tot)
ch_factor = exp_ch_tot/model_ch_tot
println("factor exp/model: ", ch_factor)

## ------------------------------------------------------------------
# RNA

rna_ids = ["m02847c"]
model_rna_tot = abs.(sum([biomass[met] for met in rna_ids]))
println("model total rna: ", model_rna_tot)
exp_rna_tot = 176.9 * 1e-3 # Niklas (2013): 103–114. https://doi.org/10.1016/j.ymben.2013.01.002. Table1
println("experimental total rna: ", exp_rna_tot)
rna_factor = exp_rna_tot/model_rna_tot
println("factor exp/model: ", rna_factor)

## ------------------------------------------------------------------
# DNA

dna_ids = ["m01721n"]
model_dna_tot = abs.(sum([biomass[met] for met in dna_ids]))
println("model total dna: ", model_dna_tot)
exp_dna_tot = 45.3 * 1e-3 # Niklas (2013): 103–114. https://doi.org/10.1016/j.ymben.2013.01.002. Table1
println("experimental total dna: ", exp_dna_tot)
dna_factor = exp_dna_tot/model_dna_tot
println("factor exp/model: ", dna_factor)

## ------------------------------------------------------------------
# Lipids

lip_ids = ["m10014c"]
model_lip_tot = abs.(sum([biomass[met] for met in lip_ids]))
println("model total lipids:", model_lip_tot)
# Niklas (2013): 103–114. https://doi.org/10.1016/j.ymben.2013.01.002. Table1
exp_lip_tot = 202.9 * 1e-3
println("experimental total lipids:", exp_lip_tot)
lip_factor = exp_lip_tot/model_lip_tot
println("factor exp/model:", lip_factor)

# Aminoacids

aa_ids = ["m10013c"]
model_prot_tot = abs.(sum([biomass[met] for met in aa_ids]))
println("model total protein: ", model_prot_tot)
exp_prot_tot = 7462.6 * 1e-3 # Niklas (2013): 103–114. https://doi.org/10.1016/j.ymben.2013.01.002. Table1
println("experimental total protein: ", exp_prot_tot)
prot_factor = exp_prot_tot/model_prot_tot
println("factor exp/model: ", prot_factor)

## ------------------------------------------------------------------
# Rescaling

# Carbohydrates
for ch in ch_ids
    biomass[ch] = biomass[ch] * ch_factor
end

# Aminoacids
for aa in aa_ids
    biomass[aa] = biomass[aa] * prot_factor
end

# lipids
for lip in lip_ids
    biomass[lip] = biomass[lip] * lip_factor
end

# DNA
for dna in dna_ids
    biomass[dna] = biomass[dna] * dna_factor
end

# RNA
for rna in rna_ids
    biomass[rna] = biomass[rna] * rna_factor
end

## ------------------------------------------------------------------
# save
sdat(HG, biomass,
    "niklas_biomass", ".jls";
    verbose = true
)


text:
function IsotropicDrag(arg0::jdouble, arg1::jdouble)
    return IsotropicDrag((jdouble, jdouble), arg0, arg1)
end

function IsotropicDrag(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble)
    return IsotropicDrag((jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function drag_acceleration(obj::IsotropicDrag, arg0::AbsoluteDate, arg1::Frame, arg2::Vector3D, arg3::Rotation, arg4::jdouble, arg5::jdouble, arg6::Vector3D, arg7::Vector{jdouble})
    return jcall(obj, "dragAcceleration", Vector3D, (AbsoluteDate, Frame, Vector3D, Rotation, jdouble, jdouble, Vector3D, Vector{jdouble}), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function drag_acceleration(obj::IsotropicDrag, arg0::FieldAbsoluteDate, arg1::Frame, arg2::FieldVector3D, arg3::FieldRotation, arg4::RealFieldElement, arg5::RealFieldElement, arg6::FieldVector3D, arg7::Vector{RealFieldElement})
    return jcall(obj, "dragAcceleration", FieldVector3D, (FieldAbsoluteDate, Frame, FieldVector3D, FieldRotation, RealFieldElement, RealFieldElement, FieldVector3D, Vector{RealFieldElement}), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function get_drag_parameters_drivers(obj::IsotropicDrag)
    return jcall(obj, "getDragParametersDrivers", Vector{ParameterDriver}, ())
end


text:
export MvNormalModel, MvNormalData, MvNormalParams

"""
Model observations as drawn from a multivariate normal distribution. See
`MvNormalData` for summarizing data for estimation and likelihood calculations.
"""
struct MvNormalModel end

"""
    $SIGNATURES

Test if `matrix` is square and conformable with `vector`.
"""
function is_conformable_square(vector::AbstractVector, matrix::AbstractMatrix)
    k = length(vector)
    size(matrix) == (k, k)
end

struct MvNormalData{TW <: Real, Tm <: AbstractVector, TS <: AbstractMatrix}
    "sum of the weights (alternatively, total number of observations)"
    W::TW
    "(weighted) sample mean"
    m::Tm
    "(weighted) sample covariance matrix"
    S::TS
    @doc """
        MvNormalData(n, m, S)

    Multivariate normal model summary statistics with `n` observations, mean `m`
    and sample covariance `S`. Only saves the summary statistics.

    !!! usage

        Use `MvNormalData(X, [wv])` to construct from data.
    """ ->                      # FIXME workaround, remove -> in v0.7
    function MvNormalData(W::TW, m::Tm, S::TS) where {TW <: Real,
                                                      Tm <: AbstractVector,
                                                      TS <: AbstractMatrix}
        @argcheck is_conformable_square(m, S)
        new{TW, Tm, TS}(W, m, S)
    end
end

size(ss::MvNormalData) = ss.W, length(ss.m)

function show(io::IO, ss::MvNormalData)
    W, k = size(ss)
    println(io, "Summary statistics for multivariate normal, $(W) × $(k) samples")
    println(io, "sample mean: $(ss.m)")
    println(io, "sample covariance: $(ss.S)")
end

"""
    $SIGNATURES

Multivariate normal summary statistics from observations (each row of `X` is an
observation).
"""
function MvNormalData(X::AbstractMatrix)
    m, S = mean_and_cov(X, 1; corrected = false)
    MvNormalData(size(X, 1), vec(m), S)
end

"""
    $SIGNATURES

Multivariate normal summary statistics from observations (each row of `X` is an
observation), with weights.
"""
function MvNormalData(X::AbstractMatrix, wv::AbstractWeights)
    m, S = mean_and_cov(X, wv, 1; corrected = false)
    MvNormalData(sum(wv), vec(m), S)
end

struct MvNormalParams{Tμ, TΣ}
    "mean"
    μ::Tμ
    "variance"
    Σ::TΣ
    @doc """
        MvNormalParams(μ, Σ)

    Parameters for the multivariate normal model ``x ∼ MvNormal(μ, Σ)``.

    !!! usage

        Construct using [`MLE(::MvNormalModel, ::MvNormalData)`](@ref).
    """ ->                      # FIXME workaround, remove -> in v0.7
    function MvNormalParams(μ::Tμ, Σ::TΣ) where {Tμ, TΣ}
        @argcheck is_conformable_square(μ, Σ)
        new{Tμ, TΣ}(μ, Σ)
    end
end

MLE(::MvNormalModel, ss::MvNormalData) = MvNormalParams(ss.m, ss.S)

function loglikelihood(::MvNormalModel, ss::MvNormalData, params::MvNormalParams)
    @unpack W, m, S = ss
    @unpack μ, Σ = params
    K = length(m)
    @argcheck length(μ) == K
    C = cholfact(Σ)
    d = m-μ
    -W/2*(K*log(2*π) + logdet(C) + dot(vec(S) + vec(d * d'), vec(inv(Σ))))
end

vec_parameters(ϕ::MvNormalParams) = vec_parameters((ϕ.μ, ϕ.Σ))


text:
using Test
using HDGElasticity

l = 2.0
m = 3.0

matrix = HDGElasticity.plane_strain_voigt_hooke_matrix_2d(l,m)
testmatrix = [(l+2m)  l     0.0
                l   (l+2m) 0.0
               0.0   0.0    m]
@test all(matrix .== testmatrix)

matrix = HDGElasticity.plane_strain_voigt_hooke_matrix(l,m,2)
@test all(matrix .== testmatrix)
@test_throws ArgumentError HDGElasticity.plane_strain_voigt_hooke_matrix(l,m,3)
@test_throws ArgumentError HDGElasticity.plane_strain_voigt_hooke_matrix(l,m,1)


text:
@testset "Test MSTM v4" begin
    using STMMRunner
    using STMMRunner.MSTM.V4

    sphere1 = Sphere(r = 1.0, x = 4.2, y = 4.3, z = 4.0, m = 1.5)
    sphere2 = Sphere(r = 1.0, x = 0.0, y = 0.3, z = 1.0, m = 1.5)

    spheres = [sphere1, sphere2]
    layers = [Layer(), Layer(d = 10.0, m = 2.0), Layer(m = 1.2)]

    base = STMMConfig(
        number_processors = Sys.CPU_THREADS ÷ 2,
        scattering_map_model = 0,
        α = 5.0,
        β = 12.0,
        spheres = spheres,
        frame = IncidentFrame,
    )

    @testset "Fixed orientation" begin
        fixed_base = STMMConfig(
            base;
            orientation = FixedOrientation,
            calculate_near_field = true,
            near_field_x_min = -5.0,
            near_field_x_max = 5.0,
            near_field_y_min = -5.0,
            near_field_y_max = 5.0,
            near_field_z_min = -5.0,
            near_field_z_max = 5.0,
            near_field_step_size = 1.0
        )

        @testset "No layer boundaries, scattering_map_model = 0, frame = incident" begin
            param = STMMConfig(fixed_base)
            @test !isnothing(run_mstm(param; keep = true))
        end

        @testset "No layer boundaries, scattering_map_model = 0, frame = incident, use_fft = true" begin
            param = STMMConfig(fixed_base; use_fft_translation = true)
            @test !isnothing(run_mstm(param; keep = true))
        end

        @testset "No layer boundaries, scattering_map_model = 0, frame = target" begin
            param = STMMConfig(fixed_base; frame = TargetFrame)
            @test !isnothing(run_mstm(param; keep = true))
        end

        @testset "No layer boundaries, periodic" begin
            param = STMMConfig(fixed_base; periodic = true)
            @test !isnothing(run_mstm(param; keep = true))
        end

        @testset "No layer boundaries, scattering_map_model = 1, frame = incident" begin
            param = STMMConfig(
                fixed_base;
                scattering_map_model = 1,
                scattering_map_dimension = 20
            )
            @test !isnothing(run_mstm(param; keep = true))
        end

        @testset "With layer boundaries, scattering_map_model = 0, frame = incident" begin
            param = STMMConfig(fixed_base; layers = layers)
            @test !isnothing(run_mstm(param; keep = true))
        end

        @testset "With layer boundaries, periodic" begin
            param = STMMConfig(fixed_base; layers = layers, periodic = true)
            @test !isnothing(run_mstm(param; keep = true))
        end

        @testset "With layer boundaries, scattering_map_model = 1, frame = incident" begin
            param = STMMConfig(
                fixed_base;
                layers = layers,
                scattering_map_model = 1,
                scattering_map_dimension = 20
            )
            @test !isnothing(run_mstm(param; keep = true))
        end
    end

    @testset "Random orientation" begin
        random_base = STMMConfig(base; orientation = RandomOrientation)

        @testset "frame = incident" begin
            param = STMMConfig(random_base)
            @test !isnothing(run_mstm(param; keep = true))
        end

        @testset "frame = target" begin
            param = STMMConfig(random_base; frame = TargetFrame)
            @test !isnothing(run_mstm(param; keep = true))
        end

        @testset "use Monte Carlo" begin
            param = STMMConfig(
                random_base;
                use_monte_carlo_integration = true,
                number_incident_directions = 10
            )
            @test !isnothing(run_mstm(param; keep = true))
        end
    end

    try
        for path in readdir()
            if isdir(path) && occursin("mstm_tmp", path)
                rm(path; force = true, recursive = true)
            end
        end
    catch _ # The unlink operation might fail on Windows so we need this try-catch
    end
end


text:
module Servers

using ..IOExtras
using ..Streams
using ..Messages
using ..Parsers
using ..ConnectionPool
import ..@info, ..@warn, ..@error, ..@debug, ..@debugshow, ..DEBUG_LEVEL
using MbedTLS: SSLConfig, SSLContext, setup!, associate!, hostname!, handshake!

if !isdefined(Base, :Nothing)
    const Nothing = Void
    const Cvoid = Void
end

import ..Dates

@static if !isdefined(Base, :Distributed)
    using Distributed
end

using ..HTTP, ..Handlers

export Server, ServerOptions, serve
#TODO:
 # add in "events" handling
 # dealing w/ cookies
 # reverse proxy?
 # auto-compression
 # health report of server
 # http authentication subrequests?
 # ip access control lists?
 # JWT authentication?
 # live activity monitoring
 # live reconfigure?
 # memory/performance profiling for thousands of concurrent requests?
 # fault tolerance?
 # handle IPv6?
 # flv & mp4 streaming?
 # URL rewriting?
 # bandwidth throttling
 # IP address-based geolocation
 # user-tracking
 # WebDAV
 # FastCGI
 # default handler:
   # handle all common http requests/etc.
   # just map straight to filesystem
 # special case OPTIONS method like go?
 # buffer re-use for server/client wire-reading
 # easter egg (response 418)
mutable struct ServerOptions
    tlsconfig::HTTP.MbedTLS.SSLConfig
    readtimeout::Float64
    ratelimit::Rational{Int}
    support100continue::Bool
    chunksize::Union{Nothing, Int}
    logbody::Bool
end

abstract type Scheme end
struct http <: Scheme end
struct https <: Scheme end

ServerOptions(; tlsconfig::HTTP.MbedTLS.SSLConfig=HTTP.MbedTLS.SSLConfig(true),
                readtimeout::Float64=180.0,
                ratelimit::Rational{Int}=Int(5)//Int(1),
                support100continue::Bool=true,
                chunksize::Union{Nothing, Int}=nothing,
                logbody::Bool=true) =
    ServerOptions(tlsconfig, readtimeout, ratelimit, support100continue, chunksize, logbody)

"""
    Server(handler, logger::IO=STDOUT; kwargs...)

An http/https server. Supports listening on a `host` and `port` via the `HTTP.serve(server, host, port)` function.
`handler` is a function of the form `f(::Request, ::Response) -> HTTP.Response`, i.e. it takes both a `Request` and pre-built `Response`
objects as inputs and returns the, potentially modified, `Response`. `logger` indicates where logging output should be directed.
When `HTTP.serve` is called, it aims to "never die", catching and recovering from all internal errors. To forcefully stop, one can obviously
kill the julia process, interrupt (ctrl/cmd+c) if main task, or send the kill signal over a server in channel like:
`put!(server.in, HTTP.Servers.KILL)`.

Supported keyword arguments include:
  * `cert`: if https, the cert file to use, as passed to `HTTP.MbedTLS.SSLConfig(cert, key)`
  * `key`: if https, the key file to use, as passed to `HTTP.MbedTLS.SSLConfig(cert, key)`
  * `tlsconfig`: pass in an already-constructed `HTTP.MbedTLS.SSLConfig` instance
  * `readtimeout`: how long a client connection will be left open without receiving any bytes
  * `ratelimit`: a `Rational{Int}` of the form `5//1` indicating how many `messages//second` should be allowed per client IP address; requests exceeding the rate limit will be dropped
  * `support100continue`: a `Bool` indicating whether `Expect: 100-continue` headers should be supported for delayed request body sending; default = `true`
  * `logbody`: whether the Response body should be logged when `verbose=true` logging is enabled; default = `true`
"""
mutable struct Server{T <: Scheme, H <: HTTP.Handler}
    handler::H
    logger::IO
    in::Channel{Any}
    out::Channel{Any}
    options::ServerOptions

    Server{T, H}(handler::H, logger::IO=STDOUT, ch=Channel(1), ch2=Channel(1),
                 options=ServerOptions()) where {T, H} =
        new{T, H}(handler, logger, ch, ch2, options)
end

mutable struct RateLimit
    allowance::Float64
    lastcheck::Dates.DateTime
end

function update!(rl::RateLimit, ratelimit)
    current = Dates.now()
    timepassed = float(Dates.value(current - rl.lastcheck)) / 1000.0
    rl.lastcheck = current
    rl.allowance += timepassed * ratelimit
    return nothing
end

function check_rate_limit(tcp;
                          ratelimits=nothing,
                          ratelimit::Rational{Int}=Int(5)//Int(1), kw...)
    ip = getsockname(tcp)[1]
    rate = Float64(ratelimit.num)
    rl = get!(ratelimits, ip, RateLimit(rate, Dates.now()))
    update!(rl, ratelimit)
    if rl.allowance > rate
        @warn "throttling $ip"
        rl.allowance = rate
    end
    if rl.allowance < 1.0
        @warn "discarding connection from $ip due to rate limiting"
        return false
    else
        rl.allowance -= 1.0
    end
    return true
end

@enum Signals KILL

function serve(server::Server{T, H}, host, port, verbose) where {T, H}

    tcpserver = Ref{Base.TCPServer}()

    @async begin
        while !isassigned(tcpserver)
            sleep(1)
        end
        while true
            val = take!(server.in)
            val == KILL && close(tcpserver[])
        end
    end

    listen(host, port;
           tcpref=tcpserver,
           ssl=(T == https),
           sslconfig=server.options.tlsconfig,
           verbose=verbose,
           tcpisvalid=server.options.ratelimit > 0 ? check_rate_limit :
                                                     (tcp; kw...) -> true,
           ratelimits=Dict{IPAddr, RateLimit}(),
           ratelimit=server.options.ratelimit) do request::HTTP.Request

        handle(server.handler, request)
    end

    return
end

Server(h::Function, l::IO=STDOUT; cert::String="", key::String="", args...) = Server(HTTP.HandlerFunction(h), l; cert=cert, key=key, args...)
function Server(handler::H=HTTP.HandlerFunction((req, rep) -> HTTP.Response(200, "Hello World!")),
               logger::IO=STDOUT;
               cert::String="",
               key::String="",
               args...) where {H <: HTTP.Handler}
    if cert != "" && key != ""
        server = Server{https, H}(handler, logger, Channel(1), Channel(1), ServerOptions(; tlsconfig=HTTP.MbedTLS.SSLConfig(cert, key), args...))
    else
        server = Server{http, H}(handler, logger, Channel(1), Channel(1), ServerOptions(; args...))
    end
    return server
end

"""
    HTTP.serve([server,] host::IPAddr, port::Int; verbose::Bool=true, kwargs...)

Start a server listening on the provided `host` and `port`. `verbose` indicates whether server activity should be logged.
Optional keyword arguments allow construction of `Server` on the fly if the `server` argument isn't provided directly.
See `?HTTP.Server` for more details on server construction and supported keyword arguments.
By default, `HTTP.serve` aims to "never die", catching and recovering from all internal errors. Two methods for stopping
`HTTP.serve` include interrupting (ctrl/cmd+c) if blocking on the main task, or sending the kill signal via the server's in channel
(`put!(server.in, HTTP.Servers.KILL)`).
"""
function serve end

serve(server::Server, host=ip"127.0.0.1", port=8081; verbose::Bool=true) = serve(server, host, port, verbose)
function serve(host::IPAddr, port::Int,
                   handler=(req, rep) -> HTTP.Response(200, "Hello World!"),
                   logger::I=STDOUT;
                   cert::String="",
                   key::String="",
                   verbose::Bool=true,
                   args...) where {I}
    server = Server(handler, logger; cert=cert, key=key, args...)
    return serve(server, host, port, verbose)
end
serve(; host::IPAddr=ip"127.0.0.1",
        port::Int=8081,
        handler=(req, rep) -> HTTP.Response(200, "Hello World!"),
        logger::IO=STDOUT,
        cert::String="",
        key::String="",
        verbose::Bool=true,
        args...) =
    serve(host, port, handler, logger; cert=cert, key=key, verbose=verbose, args...)

function getsslcontext(tcp, sslconfig)
    ssl = SSLContext()
    setup!(ssl, sslconfig)
    associate!(ssl, tcp)
    handshake!(ssl)
    return ssl
end

const nosslconfig = SSLConfig()

const nolimit = typemax(Int)

"""
    HTTP.listen([host="localhost" [, port=8081]]; <keyword arguments>) do http
        ...
    end

Listen for HTTP connections and execute the `do` function for each request.

Optional keyword arguments:
 - `ssl::Bool = false`, use https.
 - `require_ssl_verification = true`, pass `MBEDTLS_SSL_VERIFY_REQUIRED` to
   the mbed TLS library.
   ["... peer must present a valid certificate, handshake is aborted if
     verification failed."](https://tls.mbed.org/api/ssl_8h.html#a5695285c9dbfefec295012b566290f37)
 - `sslconfig = SSLConfig(require_ssl_verification)`
 - `pipeline_limit = 16`, number of concurrent requests per connection.
 - `reuse_limit = nolimit`, number of times a connection is allowed to be reused
                            after the first request.
 - `tcpisvalid::Function (::TCPSocket) -> Bool`, check accepted connection before
    processing requests. e.g. to implement source IP filtering, rate-limiting,
    etc.
 - `tcpref::Ref{Base.TCPServer}`, this reference is set to the underlying
                                  `TCPServer`. e.g. to allow closing the server.

e.g.
```
    HTTP.listen() do http::HTTP.Stream
        @show http.message
        @show HTTP.header(http, "Content-Type")
        while !eof(http)
            println("body data: ", String(readavailable(http)))
        end
        setstatus(http, 404)
        setheader(http, "Foo-Header" => "bar")
        startwrite(http)
        write(http, "response body")
        write(http, "more response body")
    end

    HTTP.listen() do request::HTTP.Request
        @show HTTP.header(request, "Content-Type")
        @show HTTP.payload(request)
        return HTTP.Response(404)
    end
```
"""

listen(f, host, port; kw...) = listen(f, string(host), Int(port); kw...)

function listen(f::Function,
                host::String="127.0.0.1", port::Int=8081;
                ssl::Bool=false,
                require_ssl_verification::Bool=true,
                sslconfig::SSLConfig=nosslconfig,
                pipeline_limit::Int=ConnectionPool.default_pipeline_limit,
                tcpisvalid::Function=(tcp; kw...)->true,
                tcpref::Ref{Base.TCPServer}=Ref{Base.TCPServer}(),
                kw...)

    if sslconfig === nosslconfig
        sslconfig = SSLConfig(require_ssl_verification)
    end

    @info "Listening on: $host:$port"
    tcpserver = Base.listen(getaddrinfo(host), port)

    tcpref[] = tcpserver

    try
        while isopen(tcpserver)
            try
                io = accept(tcpserver)
            catch e
                if e isa Base.UVError
                    @warn "$e"
                    break
                else
                    rethrow(e)
                end
            end
            if !tcpisvalid(io; kw...)
                close(io)
                continue
            end
            io = ssl ? getsslcontext(io, sslconfig) : io
            let io = Connection(host, string(port), pipeline_limit, io)
                @info "Accept:  $io"
                @async try
                    handle_connection(f, io; kw...)
                catch e
                    @error "Error:   $io" e catch_stacktrace()
                finally
                    close(io)
                    @info "Closed:  $io"
                end
            end
        end
    catch e
        if typeof(e) <: InterruptException
            @warn "Interrupted: listen($host,$port)"
        else
            rethrow(e)
        end
    finally
        close(tcpserver)
    end

    return
end

"""
Start a timeout monitor task to close the `Connection` if it is inactive.
Create a `Transaction` object for each HTTP Request received.
"""

function handle_connection(f::Function, c::Connection;
                           reuse_limit::Int=nolimit,
                           readtimeout::Int=0, kw...)

    wait_for_timeout = Ref{Bool}(true)
    if readtimeout > 0
        @async while wait_for_timeout[]
            @show inactiveseconds(c)
            if inactiveseconds(c) > readtimeout
                @warn "Timeout: $c"
                writeheaders(c.io, Response(408, ["Connection" => "close"]))
                close(c)
                break
            end
            sleep(8 + rand() * 4)
        end
    end

    try
        count = 0
        while isopen(c)
            io = Transaction(c)
            handle_transaction(f, io; final_transaction=(count == reuse_limit),
                                      kw...)
            if count == reuse_limit
                close(c)
            end
            count += 1
        end
    finally
        wait_for_timeout[] = false
    end
    return
end

"""
Create a `HTTP.Stream` and parse the Request headers from a `HTTP.Transaction`.
If there is a parse error, send an error Response.
Otherwise, execute stream processing function `f`.
"""

function handle_transaction(f::Function, t::Transaction;
                            final_transaction::Bool=false,
                            verbose::Bool=false, kw...)

    request = HTTP.Request()
    http = Streams.Stream(request, ConnectionPool.getparser(t), t)

    try
        startread(http)
    catch e
        if e isa EOFError && !messagestarted(http.parser)
            return
        elseif e isa HTTP.ParsingError
            @error e
            status = e.code == :HPE_INVALID_VERSION ? 505 :
                     e.code == :HPE_INVALID_METHOD  ? 405 : 400
            write(t, Response(status, body = HTTP.Parsers.ERROR_MESSAGES[e.code]))
            close(t)
            return
        elseif e isa HeaderSizeError
            write(t, Response(413))
            close(t)
            return
        else
            rethrow(e)
        end
    end

    if verbose
        @info http.message
    end

    response = request.response
    response.status = 200
    if final_transaction || hasheader(request, "Connection", "close")
        setheader(response, "Connection" => "close")
    end

    @async try
        handle_stream(f, http)
    catch e
        if isioerror(e)
            @warn e
        else
            @error e catch_stacktrace()
        end
        close(t)
    end
    return
end

"""
Execute stream processing function `f`.
If there is an error and the stream is still open,
send a 500 response with the error message.

Close the `Stream` for read and write (in case `f` has not already done so).
"""

function handle_stream(f::Function, http::Stream)

    try
        if applicable(f, http)
            f(http)
        else
            handle_request(f, http)
        end
    catch e
        if isopen(http) && !iswritable(http)
            @error e catch_stacktrace()
            http.message.response.status = 500
            startwrite(http)
            write(http, sprint(showerror, e))
        else
            rethrow(e)
        end
    end

    closeread(http)
    closewrite(http)
    return
end

"""
Execute Request processing function `f(::HTTP.Request) -> HTTP.Response`.
"""

function handle_request(f::Function, http::Stream)
    request::HTTP.Request = http.message
    request.body = read(http)
    request.response::HTTP.Response = f(request)
    startwrite(http)
    write(http, request.response.body)
    return
end

end # module


text:
using Schistoxpkg
using Test
using Distributions
using Random

env = Environment()
humans = env.humans

N = 1000    #population size
time_step = 10.0
N_communities = 1
community_probs = [1.0]
community_contact_rate = [1.0]
density_dependent_fecundity = 0.0007
average_worm_lifespan = 5.7
max_age = 100
initial_worms = 0
initial_miracidia = 100000*N/1000
initial_miracidia_days = 3
init_env_cercariae = 100000*N/1000
worm_stages = 1
contact_rate = 0.1
max_fecundity = 0.34
female_factor = 1
male_factor = 1
age_contact_rates = [0.032,0.61, 1,0.06]
ages_for_contacts = [4,9,15,100]
contact_rate_by_age_array = fill(0.09, trunc(Int,(max_age+1)))
mda_adherence = 0.9
mda_access = 0.9
female_factor = 1
male_factor = 1
birth_rate = 28*time_step/(1000*365)
human_cercariae_prop = 1
predis_aggregation = 0.24
cercariae_survival = 1/2
miracidia_survival = 1/2
miracidia_maturity = 24
death_prob_by_age = [0.0656, 0.0093, 0.003, 0.0023, 0.0027, 0.0038, 0.0044, 0.0048, 0.0053,
                     0.0065, 0.0088, 0.0106, 0.0144, 0.021, 0.0333, 0.0529, 0.0851, 0.1366, 0.2183, 0.2998 , 0.3698, 1]

ages_for_death = [1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60,
                   65, 70, 75, 80, 85, 90, 95, 100, 110];
r = 0.03
vaccine_effectiveness  = 0.86
drug_effectiveness = 0.86
spec_ages = [7639, 7082, 6524, 5674, 4725, 4147, 3928, 3362,
              2636, 1970, 1468, 1166, 943, 718, 455, 244]
ages_per_index = 5
record_frequency = 1/24;
predis_aggregation = 0.8
scenario = "moderate adult";
use_kato_katz = 0
pars = Parameters(N, time_step, N_communities, community_probs, community_contact_rate,
        density_dependent_fecundity, average_worm_lifespan,
        max_age, initial_worms, initial_miracidia, initial_miracidia_days, init_env_cercariae,
        worm_stages, contact_rate, max_fecundity, age_contact_rates,
        ages_for_contacts, contact_rate_by_age_array, mda_adherence, mda_access,  female_factor, male_factor, miracidia_maturity,
        birth_rate, human_cercariae_prop, predis_aggregation, cercariae_survival, miracidia_survival,
        death_prob_by_age, ages_for_death, r, vaccine_effectiveness, drug_effectiveness,
        spec_ages, ages_per_index, record_frequency, use_kato_katz)
pars = make_age_contact_rate_array(pars, scenario, [],[]);

humans, miracidia, cercariae = create_population_specified_ages(pars)

@testset "create_population" begin
    @test length(humans) == N
end

humans = generate_ages_and_deaths(20000, humans, pars)
humans = update_contact_rate(humans,  pars)
@testset "update_contact_rate1" begin
    @test humans[1].age_contact_rate == pars.contact_rate_by_age_array[end]
end
humans[1].age = 0
humans = update_contact_rate(humans,  pars)

@testset "update_contact_rate2" begin
    @test humans[1].age_contact_rate == pars.contact_rate_by_age_array[1]
end

@testset "miracidia_death" begin
    @test isapprox(miracidia_death!(miracidia, pars), [pars.initial_miracidia, pars.initial_miracidia , round(pars.initial_miracidia  * pars.miracidia_survival)])
end

@testset "cercariae_death" begin
    @test cercariae_death!(cercariae, pars)== pars.init_env_cercariae  * pars.cercariae_survival

end

@testset "worm_pairs" begin
    @test isapprox(calculate_worm_pairs([[1,2],[2,2]],[[5,1],[2,0]]), [3,2])
end

humans[1].female_worms = [1,2]
humans[1].male_worms = [5,1]

female_worms = (p->p.female_worms).(humans)
male_worms = (p->p.male_worms).(humans)

@testset "worm_pairs" begin
    @test calculate_worm_pairs(female_worms, male_worms)[1] ==3
end

@testset "make_age_contact_rate_array(max_age,scenario)" begin
    @test pars.contact_rate_by_age_array[1] == 0.032
end

@testset "make_age_contact_rate_array(max_age,scenario)" begin
    @test pars.contact_rate_by_age_array[100] == 0.06
end


text:
""" Basic utilitary functions """

# All functions are public

" True if `sol1` and `sol2` are different solutions, false otherwise. "
function differentFloatSols(sol1::Sol,sol2::Sol)
    return abs(sol1.y[1]-sol2.y[1]) > 0.000001 || abs(sol1.y[2]-sol2.y[2]) > 0.000001
end

# NOTE : for domination tests, one could overwrite the operator >= of julia
# and test dominance using sol1 >= sol2

" True if the given solution `sol1` dominates the solution `sol2`. "
function dominate(sol1::Sol, sol2::Sol)
	return (sol1.y[1] >= sol2.y[1] && sol1.y[2] > sol2.y[2]) || (sol1.y[1] > sol2.y[1] && sol1.y[2] >= sol2.y[2])
end

" True if the given solution `sol` dominates the nadirPoint `nadirPoint`. "
function dominate(sol::Sol, nadirPoint::PairOfSolution)
	return sol.y[1] >= (nadirPoint.solL.y[1] + 1) && sol.y[2] >= (nadirPoint.solR.y[2]+1)
end


text:
# Example of writing synaptic kinetics

using Conductor, IfElse, OrdinaryDiffEq, Plots, Unitful, ModelingToolkit
import Unitful: mV, mS, cm, µm, pA, nA, mA, µA, ms, nS, pS
import Conductor: Na, K

Vₘ = MembranePotential()

nav_kinetics = [
    Gate(AlphaBetaRates,
         αₘ = IfElse.ifelse(Vₘ == -40.0, 1.0, (0.1*(Vₘ + 40.0))/(1.0 - exp(-(Vₘ + 40.0)/10.0))),
         βₘ = 4.0*exp(-(Vₘ + 65.0)/18.0),
         p = 3)
    Gate(AlphaBetaRates,
         αₕ = 0.07*exp(-(Vₘ+65.0)/20.0),
         βₕ = 1.0/(1.0 + exp(-(Vₘ + 35.0)/10.0)))]

kdr_kinetics = [
    Gate(AlphaBetaRates,
         αₙ = IfElse.ifelse(Vₘ == -55.0, 0.1, (0.01*(Vₘ + 55.0))/(1.0 - exp(-(Vₘ + 55.0)/10.0))),
         βₙ = 0.125 * exp(-(Vₘ + 65.0)/80.0),
         p = 4)]

@named NaV = IonChannel(Sodium, nav_kinetics, 120mS/cm^2)
@named Kdr = IonChannel(Potassium, kdr_kinetics, 36mS/cm^2)
@named leak = PassiveChannel(Leak, 0.3mS/cm^2)

# Equilibrium potentials are a implicit description of a ion concentration gradient
gradients = Equilibria([Na   =>  50.0mV,
                        K    => -77.0mV,
                        Leak => -54.4mV])

#area = 4*pi*(10µm)^2
area = 0.629e-3cm^2

@named neuron1 = Soma([NaV,Kdr,leak], gradients, holding = 5000pA, area = ustrip(Float64, cm^2, area));
@named neuron2 = Soma([NaV,Kdr,leak], gradients, area = ustrip(Float64, cm^2, area));

# Synaptic model
syn∞ = 1/(1 + exp((-35 - Vₘ)/5))
τsyn = (1 - syn∞)/(1/40)
syn_kinetics = Gate(SteadyStateTau, :s, syn∞, τsyn, 1)
EGlut = Equilibrium{Conductor.Mixed}(0mV, :Glut)

@named Glut = Conductor.SynapticChannel(Leak, [syn_kinetics], EGlut, 30nS);
topology = [neuron1 => (neuron2, Glut)];
network = Conductor.Network([neuron1, neuron2], topology)

t = 250
sim = Simulation(network, time = t*ms)
solution = solve(sim, Rosenbrock23())

# Plot at 5kHz sampling
plot(solution; plotdensity=Int(t*5))


text:
import Term: Panel, TextBox

@testset "\e[34mPanel - no content" begin
    for style in ("default", "red", "on_blue")
        testpanel(
            Panel(;fit=true, style=style), 3, 2
        )

        testpanel(
            Panel(), 88, 2
        )

        testpanel(
            Panel(; width=12, height=4, style=style), 12, 4
        )
    end

end

@testset "\e[34mPANEL - fit - measure" begin
    for style in ("default", "red", "on_blue")
        for justify in (:left, :center, :right)
            # ----------------------------- text only content ---------------------------- #
            testpanel(
                Panel("t"; fit=true, style=style), 7, 3
            )
            testpanel(
                Panel("test"; fit=true, style=style), 10, 3
            )
            testpanel(
                Panel("1234\n123456789012"; fit=true, style=style), 18, 4
            )
            testpanel(
                Panel("나랏말싸미 듕귁에 달아"; fit=true, style=style), 28, 3
            )
            testpanel(
                Panel("나랏말싸미 듕귁에 달아\n1234567890123456789012"; fit=true, style=style), 28, 4
            )
            testpanel(
                Panel("."^500; fit=true, style=style), displaysize(stdout)[2]-4, nothing
            )

            # ------------------------------- nested panels ------------------------------ #
            testpanel(
                Panel(
                    Panel("test"; fit=true, style=style);
                fit=true, style=style),
                16, 5
            )

            testpanel(
                Panel(
                    Panel(Panel("."; fit=true, style=style); fit=true, style=style);
                fit=true, style=style), 19, 7
            )

            # @test_nothrow Panel(
            #         Panel("."^250; fit=true, style=style); fit=true, style=style
            #     )
        end
    end
end

@testset "\e[34mPANEL - nofit - measure" begin
    for style in ("default", "red", "on_blue")
        for justify in (:left, :center, :right)
            # ----------------------------- text only content ---------------------------- #
            testpanel(
                Panel("t"; style=style), 88, 3
            )
            testpanel(
                Panel("test"; style=style), 88, 3
            )
            testpanel(
                Panel("1234\n123456789012"; style=style), 88, 4
            )
            testpanel(
                Panel("나랏말싸미 듕귁에 달아"; style=style), 88, 3
            )
            testpanel(
                Panel("나랏말싸미 듕귁에 달아\n1234567890123456789012"; style=style), 88, 4
            )
            testpanel(
                Panel("."^1500; style=style), 88, 21
            )

            # ------------------------------- nested panels ------------------------------ #
            testpanel(
                Panel(
                    Panel("test");
                fit=true, style=style), 94, 5
            )

            testpanel(
                Panel(
                    Panel(Panel("."); fit=true, style=style);
                fit=true, style=style), 100, 7
            )

            testpanel(
                Panel(
                    Panel("."^250); fit=true, style=style
                ), 94, 8
            )

            testpanel(
                Panel(
                    Panel("test"; style=style);
            ), 94, 5
            )

            testpanel(
                Panel(
                    Panel(Panel("."; style=style); style=style);
            ), 100, 7
            )

            testpanel(
                Panel(
                    Panel("."^250; style=style);
                ), 94, 8
            )

            testpanel(
                Panel(
                    Panel("t1"; style=style),
                    Panel("t2"; style=style),
                ), 94, 8
            )

            testpanel(
                Panel(
                    Panel("test", width=22);  width=30, height=8
                ), 30, 8
            )

            testpanel(
                Panel(
                    Panel("test", width=42);  width=30, height=8
                ), 48, 8
            )

            testpanel(
                Panel(
                    Panel("test", width=42,height=12);  width=30, height=8
                ), 48, 14
            )
        end
    end
end

@testset "\e[34mPanel + renderables" begin
    testpanel(
        Panel(
            RenderableText("x".^5)
        ), 88, 3
    )

    testpanel(
        Panel(
            RenderableText("x".^500)
        ), 88, nothing
    )

    testpanel(
        Panel(
            RenderableText("x".^5); fit=true
        ), 11, 3
    )

    testpanel(
        Panel(
            RenderableText("x".^500); fit=true
        ), displaysize(stdout)[2]-4, nothing
    )

end

@testset "\e[34mPANEL - titles" begin
    for fit in (true, false)
        for justify in (:left, :center, :right)
            for style in ("red", "bold", "default", "on_green")

                testpanel(
                    Panel("."^50, title="test",
                            title_style=style,
                            title_justify=justify,
                            subtitle="subtest",
                            subtitle_style=style,
                            subtitle_justify=justify,
                            fit=fit
                            ),
                    fit ? nothing : 88,
                    nothing
                )

                testpanel(
                    Panel(
                        Panel("."^50, title="test",
                                title_style=style,
                                title_justify=justify,
                                subtitle="subtest",
                                subtitle_style=style,
                                subtitle_justify=justify,
                                  fit=fit,
                                )
                    ),
                    fit ? nothing : 94,
                    nothing
                )
            end
        end
    end
end

@testset "\e[34mTBOX" begin
    w = displaysize(stdout)[2]

    for justify in (:left, :center, :right)

        testpanel(
            TextBox(
                "nofit"^25;
                width=1000,
                justify=justify
            ),w - 4, nothing)

        testpanel(
            TextBox(
                "truncate"^25;
                width=100,
                fit=:truncate,
                justify=justify
        ), 100, 3)

        testpanel(
            TextBox(
                "truncate"^25;
                width=100,
                justify=justify
            ), 100, 7)

        testpanel(
            TextBox(
                "truncate"^8;
                fit=:fit,
                justify=justify
            ), 68, 4)

        testpanel(
            TextBox(
                "[red]truncate[/red]"^8;
                fit=:fit,
                justify=justify
            ), 68, 4)

        testpanel(
            TextBox(
                "[red]truncate[/red]test"^8;
                fit=:fit,
                justify=justify
            ), 100, 4)

        testpanel(TextBox(
                "[red]tru\nncate[/red]test"^1;
                fit=:fit,
                justify=justify
            ), 13, 7)
    end
end

@testset "\e[34mPanel - padding" begin

    p = Panel("."^24; padding = [4, 4, 2, 2])
    testpanel(p, 88, 7)
    # @test string(p) == "\e[22m╭──────────────────────────────────────────────────────────────────────────────────────╮\e[22m\n\e[22m│\e[22m                                                                                \e[22m│\e[22m\n\e[22m│\e[22m                                                                                \e[22m│\e[22m\n\e[22m│\e[22m    ........................                                                          \e[22m│\e[22m\n\e[22m│\e[22m                                                                                \e[22m│\e[22m\n\e[22m│\e[22m                                                                                \e[22m│\e[22m\n\e[22m╰──────────────────────────────────────────────────────────────────────────────────────╯\e[22m"

    p = Panel("."^24; padding = [4, 4, 2, 2], fit=true)
    testpanel(p, 34, 7)

end


text:
#!/usr/bin/env julia

include("../src/mesh/trees.jl")

using .Trees: Tree, refine!, size, capacity, leaf_nodes, refine_box!
using TimerOutputs
using Profile

to = TimerOutput()

x_start = -16
x_end = 16
capacity_ = 1000
initial_refinement_level = 2

domain_center = (x_start + x_end) / 2
domain_length = x_end - x_start

# Create tree object
@timeit to "create tree" t = Tree(Val{1}(), capacity_, [domain_center], domain_length)
println("Initial tree:")
println(t)

# Create initial refinement
@timeit to "initial refinement" for l = 1:initial_refinement_level
  @timeit to "refine!" refine!(t)
  println("After uniform refinement to level $l:")
  println(t)
end

# Create non-uniform refinement
@timeit to "local refinement" begin
  @timeit to "refine!" refine!(t, 4)
  println("After refining node 4:")
  println(t)
  @timeit to "refine!" refine!(t, 5)
  println("After refining node 5:")
  println(t)
end

# Refine one level further
@timeit to "refine!" refine!(t)
println("After refining everything once more:")
println(t)

refine_box!(t, -8, 8)
println("After refining center nodes:")
println(t)

# Print tree information
@show leaf_nodes(t)
@show t.coordinates[:, leaf_nodes(t)]

print_timer(to)
println()


text:
"""
Abstract type for sampling with:
 * `Random.rand(BayesNet, BayesNetSampler)`
 * `Random.rand(BayesNet, BayesNetSampler, nsamples)`
 * `Random.rand!(Assignment, BayesNet, BayesNetSampler)`
"""
abstract type BayesNetSampler end

"""
Overwrites assignment with a sample from bn using the sampler
"""
Random.rand!(a::Assignment, bn::BayesNet, sampler::BayesNetSampler) =
        error("rand! not implemented for $(typeof(sampler))")

"""
Returns an assignment sampled from the bn using the provided sampler
"""
Random.rand(bn::BayesNet, sampler::BayesNetSampler) = rand!(Assignment(), bn, sampler)

"""
Generates a DataFrame containing a dataset of variable assignments.
Always return a DataFrame with `nsamples` rows.
"""
function Random.rand(bn::BayesNet, sampler::BayesNetSampler, nsamples::Integer)

    a = rand(bn, sampler)
    df = DataFrame()
    for cpd in bn.cpds
        df[!, name(cpd)] = Array{typeof(a[name(cpd)])}(undef, nsamples)
    end

    for i in 1:nsamples
        rand!(a, bn, sampler)
        for cpd in bn.cpds
            n = name(cpd)
            df[i, n] = a[n]
        end
    end

    df
end

#
# Direct Sampling
#
"""
Straightforward sampling from a BayesNet.
The default sampler.
"""
struct DirectSampler <: BayesNetSampler end

function Random.rand!(a::Assignment, bn::BayesNet, sampler::DirectSampler)
    for cpd in bn.cpds
        a[name(cpd)] = rand(cpd, a)
    end
    a
end

Random.rand!(a::Assignment, bn::BayesNet) = rand!(a, bn, DirectSampler())

Random.rand(bn::BayesNet, nsamples::Integer) = rand(bn, DirectSampler(), nsamples)

Random.rand(bn::BayesNet) = rand(bn, DirectSampler())

#
# Rejection Sampling
#
"""
Rejection Sampling in which the assignments are forced to be consistent with the provided values.
Each sampler is attempted at most `max_nsamples` times before returning an empty assignment.
"""
struct RejectionSampler <: BayesNetSampler
    evidence::Assignment
    max_nsamples::Int
end
RejectionSampler(pair::Pair{NodeName}...; max_nsamples::Integer=100) =
        RejectionSampler(Assignment(pair), max_nsamples)

function Random.rand!(a::Assignment, bn::BayesNet, sampler::RejectionSampler)
    for sample_count in 1 : sampler.max_nsamples
        rand!(a, bn)
        if consistent(a, sampler.evidence)
            return a
        end
    end
    return empty!(a)
end

Random.rand(bn::BayesNet, nsamples::Integer, evidence::Assignment) =
        rand(bn, RejectionSampler(evidence, 100), nsamples)

Random.rand(bn::BayesNet, nsamples::Integer, pair::Pair{NodeName}...) =
        rand(bn, nsamples, Assignment(pair))

#
# Likelihood Sampling
#
"""
Draw an assignment from the Bayesian network but set any variables in the evidence accordingly.
Returns the assignment and the probability weighting associated with the evidence.
"""
function get_weighted_sample!(a::Assignment, bn::BayesNet, evidence::Assignment)
    w = 1.0
    for cpd in bn.cpds
        varname = name(cpd)
        if haskey(evidence, varname)
            a[varname] = evidence[varname]
            w *= pdf(cpd, a)
        else
            a[varname] = rand(cpd, a)
        end
    end

    return (a, w)
end

get_weighted_sample!(a::Assignment, bn::BayesNet, pair::Pair{NodeName}...) =
        get_weighted_sample!(a, bn, Assignment(pair))

get_weighted_sample(bn::BayesNet, evidence::Assignment) =
        get_weighted_sample!(Assignment(), bn, evidence)

get_weighted_sample(bn::BayesNet, pair::Pair{NodeName}...) =
        get_weighted_sample(bn, Assignment(pair))

sample_weighted_dataframe(weighted_dataframe::DataFrame) =
        sample_weighted_dataframe!(Assignment(), weighted_dataframe)

"""
Likelihood Weighted Sampling
"""
struct LikelihoodWeightedSampler <: BayesNetSampler
    evidence::Assignment
end
LikelihoodWeightedSampler(pair::Pair{NodeName}...) =
        LikelihoodWeightedSampler(Assignment(pair))

function Random.rand!(a::Assignment, bn::BayesNet, sampler::LikelihoodWeightedSampler)
    get_weighted_sample!(a, sampler.weighted_dataframe)
    return a
end

Random.rand(bn::BayesNet, sampler::LikelihoodWeightedSampler, nsamples::Integer) =
        get_weighted_dataframe(bn, nsamples, sampler.evidence)

#
# only used by the src/gibbs code, not the inference code
#
"""
A dataset of variable assignments is obtained with an additional column
of weights in accordance with the likelihood of each assignment.
"""
function get_weighted_dataframe(bn::BayesNet, nsamples::Integer, evidence::Assignment)

    t = Dict{Symbol, Vector{Any}}()
    for name in names(bn)
        t[name] = Any[]
    end

    w = ones(Float64, nsamples)
    a = Assignment()

    for i in 1:nsamples
        a, weight = get_weighted_sample!(a, bn, evidence)
        for (varname, val) in a
             push!(t[varname], val)
        end
        w[i] = weight
    end
    t[:p] = w / sum(w)

    convert(DataFrame, t)
end

get_weighted_dataframe(bn::BayesNet, nsamples::Integer, pair::Pair{NodeName}...) =
        get_weighted_dataframe(bn, nsamples, Assignment(pair))

"""
Chooses a sample at random from a weighted dataframe
"""
function sample_weighted_dataframe!(a::Assignment, weighted_dataframe::DataFrame)
    p = weighted_dataframe[:, :p]
    n = length(p)
    i = 1
    c = p[1]
    u = rand()
    while c < u && i < n
        c += p[i += 1]
    end
    for varname in names(weighted_dataframe)
        if varname != :p
            a[varname] = weighted_dataframe[i, varname]
        end
    end
    return a
end


text:
#ДАНО: Робот - рядом с горизонтальной перегородкой (под ней), бесконечно продолжающейся в обе стороны, в которой имеется проход шириной в одну клетку.
#РЕЗУЛЬТАТ: Робот - в клетке под проходом

function find_passage(r::Robot)
    n=0;
    side=Ost
    while isborder(r,Nord)==true
        n+=1
        moves!(r,side,n)
        side=inverse(side)
    end
end
moves!(r,side,num_steps)=for _ in 1:num_steps move!(r,side) end
inverse(side::HorizonSide)=HorizonSide(mod(Int(side)+2,4))


text:
using Revise
using Random, SpatialIndexing

const SI = SpatialIndexing

Random.seed!(1)
seq_tree = RTree{Float64, 2}(Int, String, leaf_capacity = 20, branch_capacity = 20)
bulk_tree = RTree{Float64, 2}(Int, String, leaf_capacity = 10, branch_capacity = 10)
mbrs = Vector{SI.mbrtype(seq_tree)}()
for i in 1:10000
    t = 10 * rand()
    u, v = (0.1*t)^0.25 * randn(2)
    x = t^1.5 * sin(2*t) + u
    y = t^1.5 * cos(2*t) + v
    rmbr = SI.Rect((x, y), (x, y))
    push!(mbrs, rmbr)
end
for (i, br) in enumerate(mbrs)
    insert!(seq_tree, br, i, string(i))
end
SI.check(seq_tree)

SI.load!(bulk_tree, enumerate(mbrs), convertel -> eltype(bulk_tree)(x[2], x[1], string(x[1])))
SI.check(bulk_tree)

include(joinpath(@__DIR__, "dataframe_utils.jl"))

seq_tree_df = convert(DataFrame, seq_tree)

include(joinpath(@__DIR__, "plot_utils.jl"))

seq_tree_plot = plot(seq_tree);
open(joinpath(@__DIR__, "spiral_rtree_seq.html"), "w") do io
    PlotlyJS.savehtml(io, seq_tree_plot, :embed)
end

using ORCA
savefig(seq_tree_plot, joinpath(@__DIR__, "spiral_rtree_seq.png"), width=1000, height=800)

bulk_tree_plot = plot(bulk_tree);
open(joinpath(@__DIR__, "spiral_rtree_bulk.html"), "w") do io
    PlotlyJS.savehtml(io, bulk_tree_plot, :embed)
end


text:
@enum TokenType ILLEGAL EOF IDENT INT ASSIGN PLUS MINUS BANG ASTERISK SLASH EQ NOT_EQ LT GT COMMA SEMICOLON LPAREN RPAREN LBRACKET RBRACKET LBRACE RBRACE COLON FUNCTION MACRO LET TRUE FALSE NULL IF ELSE WHILE BREAK CONTINUE RETURN STRING

struct Token
    type::TokenType
    literal::String
end

const KEYWORDS = Dict{String,TokenType}(
    "fn" => FUNCTION,
    "let" => LET,
    "true" => TRUE,
    "false" => FALSE,
    "null" => NULL,
    "if" => IF,
    "else" => ELSE,
    "while" => WHILE,
    "break" => BREAK,
    "continue" => CONTINUE,
    "return" => RETURN,
    "macro" => MACRO,
)

lookup_indent(ident::AbstractString) =
    if ident in keys(KEYWORDS)
        KEYWORDS[ident]
    else
        IDENT
    end


text:
using Documenter, FilePathsBase, FilePathsBase.TestPaths, LinearAlgebra

makedocs(
    modules=[FilePathsBase],
    format=Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
    ),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/rofinn/FilePathsBase.jl/blob/{commit}{path}#L{line}",
    sitename="FilePathsBase.jl",
    authors="Rory Finnegan",
    checkdocs = :exports,
    # strict = true,
)

deploydocs(
    repo = "github.com/rofinn/FilePathsBase.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
##############################################
# PARTE "1": CARGA CABIDA                    #
##############################################

# cd("C:/Users/rjvia/.julia/dev/Cabida/src/")

# include("Cabida.jl")
# using .Cabida

using Cabida, NPFinancial, .poly2D, .polyShape

##############################################
# PARTE "2": GENERACIÓN DE PARÁMETROS        #
##############################################

#EJEMPLO 1 (Córdova y Figueroa con Vargas Fontecilla, Quinta Normal)
factorCorreccion = 0.8345;
x = factorCorreccion * [-7869557.5 -7869559.7 -7869535.4 -7869535.5 -7869490.4 -7869488.3]';
y = factorCorreccion * [-3952906.8 -3952965.7 -3952966.5 -3952970.9 -3952972.4 -3952909.6]';
# Calle: [1 6]
dcp = datosCabidaPredio(x, y, [1 6], [18 20], 0, 200);

"""
#EJEMPLO 2 (Catedral llegando a General Velásquez, Quinta Normal)
factorCorreccion=0.8547;
x=factorCorreccion*[-7869442.1 -7869413.9 -7869409.5 -7869438.5]';
y=factorCorreccion*[-3953939.3 -3953935.5 -3953861.6 -3953864.5]';
# Calle: [1]
dcp = datosCabidaPredio(x, y, [1], [20], 1, 200);
"""

"""
factorCorreccion = 5;
x = factorCorreccion * [0 10 15 20 25 13 12 4 9]';
y = factorCorreccion * [10 0 12 7 20 25 22 28 14]';
# Calle: [1 5]
dcp = datosCabidaPredio(x, y, [1 5], [20 20], 1, 200);
"""

dcn = datosCabidaNormativa(5, #SEPMIN (m)
                            5, #ANTEJARDIN (m)
                            2.75, #RASANTE (m de altura / m de separación)
                            24, #ALTURAMAX (m)
                            9, #MAXPISOS (unidades)
                            .7, #COEFOCUPACION (m2 / m2 de terreno)
                            1000, #SUBPREDIALMIN (m2)
                            2500, #DENSIDADMAX (Habitantes / 10000 m2 de terreno bruto)
                            4, #COEFCONSTRUCTIBILIDAD (m2 / m2 de terreno)
                            [.5, 1, 1, 1], #ESTACIONAMIENTOSPORVIV (unidades / departamento)
                            .15, #PORCADICESTACVISITAS (unidades / estacionamiento vendible)
                            34, #SUPPORESTACIONAMIENTO (m2 / Estacionamiento)
                            .25, #ESTBICICLETAPORVIV
                            3, #ESTACPORBICICLETERO
                            10, #MAXSUBTE (unidades)
                            .8, #COEFOCUPACIONEST (m2 / superficieTerreno)
                            7, #SEPESTMIN (m)
                            0 #REDUCCIONESTPORDISTMETRO
                            );

dca = datosCabidaArquitectura(2.625, #ALTURAPISO (m / piso)
                                .2, #PORCSUPCOMUN (m2 / m2 útil)
                                .05, #PORCTERRAZA (m2 / m2 útil)
                                12, #ANCHOMAX (m)
                                MATCELDASCONF_,
                                MATCONFHOR_
                                );

dcc = datosCabidaComercial([20, 30, 38, 45], #SUPDEPTOUTIL (m2)
                            [.3*1, 0, .4, 1], #MAXPORCTIPODEPTO
                            [56, 55, 54, 52], #PRECIOVENTA (UF / m2 vendible)
                            200 #PRECIOVENTAEST (UF / unidad)
                            );

dcu = datosCabidaUnit(30, #duracionProyecto
                        20, #CostoTerreno
                        2 * 1.19 / 100, #ComisionCorredor
                        0.60, #Demolicion
                        20, #LosaSNT
                        10, #LosaBNT
                        5, #EstacVisita
                        3 / 100, #ExtrasPostVenta
                        1.75 / 100, #Arquitectura
                        2.75 / 100, #GestionAdministracion
                        30, #ITO
                        30, #Contabilidad
                        0.85 / 100, #Legales
                        0.145, #Calculo
                        250, #MecanicaSuelo
                        50, #Topografia
                        .45 / 100, #ProyectoElectrico
                        0, #ProyectoSanitario
                        0, #ProyectoSeguridad
                        50, #ProyectoPavimentacion
                        25, #ProyectoBasura
                        0, #ProyectoCalefaccion
                        1.2 / 100, #Marketing
                        0.5 / 100, #PilotSalaVentas
                        0, #NotarialesLegales
                        0, #CBR
                        0, #EmpalmesAportes
                        600, #DerechosPermisosMunicipales
                        15, #InspeccionCertificacionTecnica
                        0, #Contribuciones
                        120, #CopiaPlanosArquitctura
                        60, #CopiaPlanosObra
                        300, #HallAcceso
                        400, #JardinesPaisajismo
                        0, #HabilitacionSalaUsosMultiples
                        0, #HabilitacionPiscinaEnfermeria
                        0, #CuentasServiciosDeptos
                        0, #GastosComunes
                        0, #ContribucionesDeptos
                        0, #Seguros
                        1 / 100, #ResponsabilidadCivil
                        1 / 100, #PostVentaInmobiliaria
                        0.1 / 100, #SeguroVentaEnVerde
                        1 / 100 #Imprevistos
                        );

dcf = datosCabidaFlujo([0.05, 0.05, 0.05, 0.05, 0.10, 0.70], #IngresosVentas
                        [1, 0, 0, 0, 0, 0], #CostoTerreno
                        [0.05, 0.10, 0.35, 0.35, .15, 0], #CostoConstruccion
                        [0.40, 0.15, 0.15, 0.15, .15, 0], #CostoHonorariosProyectos
                        [0, 0.10, 0.20, 0.20, 0.25, 0.25], #CostoVenta
                        [0.1, 0.40, 0.25, 0.15, 0.05, 0.05], #CostoInmobiliarioObra
                        [0, 0, 0, 0, 0.50, 0.50], #CostosHabilitacion
                        [0, 0, 0, 0, 0.25, 0.75], #CostosPuestaEnMarcha
                        [0.05, 0.05, 0.05, 0.05, 0.05, 0.75], #CostosAtencionCliente
                        [0, 0.20, 0.20, 0.20, 0.20, 0.20], #Imprevistos
                        [0.05, 0.10, 0.35, 0.35, .15, 0], #ingresoLineaCredito
                        [0, 0, 0.10, 0.35, 0.35, .2], #pagoLineaCredito
                        .27, #TasaImpuestoRenta
                        .05, #TasaInteresLineaCredito
                        0.8 * dcu.duracionProyecto, #DuracionLineaCredito
                        15 #Tasaciones
                        );

MATCELDASCONF, MATCONFHOR, MATCONFVERT = generaMatCeldasConf(dca);
NumConfig, NumCeldas = size(MATCELDASCONF);

numTiposDepto = length(dcc.SUPDEPTOUTIL);

# Posiciona el origen del predio (esq. inferior izquierda) en el punto 0,0. Obtiene estrucura con las matrices de rotación
V_predio, R = infoPredio(dcp.x, dcp.y);
numLados = 3; # length(R);

# Calcula area del predio
superficieTerreno = polyArea(V_predio);

# Calcula matriz V_areaEdif asociada a los vértices del area de edificación
V_areaEdif = copy(V_predio);
numLadosPredio = size(V_areaEdif, 1);
conjuntoLados = 1:numLadosPredio;
conjuntoLadosCalle = dcp.ladosConCalle;
conjuntoLadosVecinos = setdiff(conjuntoLados, conjuntoLadosCalle);
sepVecinos = dcn.SEPMIN;
V_areaEdif = expandPolygonSides(V_areaEdif, -sepVecinos, conjuntoLadosVecinos);
V_areaEdif = expandPolygonSides(V_areaEdif, -dcn.ANTEJARDIN, conjuntoLadosCalle);

# Descompone el area de edificación no-convexa en la unión de varias areas convexas
aux, SP = nonconv2sumofconv_v2(PolyShape([V_areaEdif],1));
numNonConv = length(SP);
s_areaEdif = Array{PolyData,1}([]);
for i = 1:numNonConv
    s_V = V_areaEdif[SP[i].points,:];
    s_A, s_b = vert2con(s_V);
    push!(s_areaEdif, PolyData(s_V, s_A, s_b, []));
end

# Genera polihedro asociado al area bruta
for i = 1:length(conjuntoLadosCalle)
    if i == 1
        V_bruto = copy(V_predio);
    end
    V_bruto = expandPolygonSide(V_bruto, min(30, dcp.ANCHOESPACIOPUBLICO[i] / 2), conjuntoLadosCalle[i]);
    global V_bruto
end

# Descompone el area bruta no-convexa en la unión de varias areas convexas
aux, SP_bruto = nonconv2sumofconv_v2(PolyShape([V_bruto],1));
desplazamientoLadoComun = dcn.ALTURAMAX / dcn.RASANTE;
s_bruto = Array{PolyData,1}([]);
for i = 1:numNonConv
    flagLadosComunes = SP_bruto[i].ladoComun .> 0;
    ladosComunes = SP_bruto[i].ladoComun[flagLadosComunes];
    numLadosComunes = sum(flagLadosComunes);
    V_bruto_i = V_bruto[SP_bruto[i].points,:];
    for j = 1:numLadosComunes
        V_bruto_i = expandPolygonSide(V_bruto_i, desplazamientoLadoComun, ladosComunes[j]);
    end
    s_A, s_b = vert2con(V_bruto_i);
    s_proyeccion = (s_A[:,1].^2 + s_A[:,2].^2).^(.5);
    push!(s_bruto, PolyData(V_bruto_i, s_A, s_b, s_proyeccion));
end
superficieTerrenoBruta = polyArea(V_bruto);

# Calcula el volumen teórico
sepAlturaMax = dcn.ALTURAMAX / dcn.RASANTE;
alturaCorteVecinos = dcn.RASANTE * sepVecinos;
alturaCorteEspPublico = dcn.ALTURAMAX;
sepEspPublico = alturaCorteEspPublico / dcn.RASANTE;
deltaSepEspPublicoVecinos = sepEspPublico - sepVecinos;
V_corteVecinos = copy(V_areaEdif);
for i = conjuntoLadosVecinos
    if i == conjuntoLadosVecinos[1]
        V_corteEspPublico = copy(V_corteVecinos);
    end
    V_corteEspPublico = expandPolygonSide(V_corteEspPublico, -deltaSepEspPublicoVecinos, conjuntoLados[i]);
    global V_corteEspPublico
end
for i = conjuntoLadosCalle
    if i == conjuntoLadosCalle[1]
        V_corteAlturaMax = copy(V_corteEspPublico);
    end
#    V_corteAlturaMax = expandPolygonSide(V_corteAlturaMax, -(sepAlturaMax - sepEspPublico), conjuntoLados[i]);
    global V_corteAlturaMax
end
V_base = copy(V_areaEdif);
numVerticesBase = size(V_base, 1);
V_volteor = [V_base zeros(numVerticesBase, 1);
             V_corteVecinos alturaCorteVecinos * ones(numVerticesBase, 1);
             V_corteAlturaMax dcn.ALTURAMAX * ones(numVerticesBase, 1)
             ];

maxDeptos = dcn.DENSIDADMAX / 4 * superficieTerrenoBruta / 10000;

##############################################
# PARTE "3": DEFINICIÓN SOLVER               #
##############################################

using JuMP, Cbc, GAMS
#m = Model(Cbc.Optimizer)
#set_optimizer_attribute(m, "ratioGap", 0.001)
#set_optimizer_attribute(m, "threads", 3)

m = Model(GAMS.Optimizer)
set_optimizer_attribute(m, "OptCR", 0.001)
set_optimizer_attribute(m, "Threads", 3)

##############################################
# PARTE "4": VARIABLES DE DECISION           #
##############################################

@variables(m, begin
    0 <= pos_x[i = 1:NumCeldas]
    0 <= pos_y[i = 1:NumCeldas]
    0 <= largo_x[i = 1:NumCeldas, r = 1:numLados]
    0 <= largo_y[i = 1:NumCeldas, r = 1:numLados]
    0 <= largo[i = 1:NumCeldas]
    0 <= areaCelda[z = 1:dcn.MAXPISOS, i = 1:NumCeldas]
    0 <= areaBasal[z = 1:dcn.MAXPISOS]
    0 <= altura[i = 1:NumCeldas] <= dcn.ALTURAMAX
    pisoActivoCeldas[i = 1:NumCeldas, j = 1:dcn.MAXPISOS], Bin
    pisoActivo[z = 1:dcn.MAXPISOS], Bin
    0 <= numDeptosTipo[u = 1:numTiposDepto]
    configSel[i = 1:NumConfig], Bin
    d[i = 1:numLados], Bin
    q1[i = 1:NumCeldas, j = 1:numNonConv], Bin
    q2[i = 1:NumCeldas, j = 1:numNonConv], Bin
    q3[i = 1:NumCeldas, j = 1:numNonConv], Bin
    q4[i = 1:NumCeldas, j = 1:numNonConv], Bin
    w1[i = 1:NumCeldas, j = 1:numNonConv], Bin
    w2[i = 1:NumCeldas, j = 1:numNonConv], Bin
    w3[i = 1:NumCeldas, j = 1:numNonConv], Bin
    w4[i = 1:NumCeldas, j = 1:numNonConv], Bin
    0 <= CostoUnitTerreno
end)

##############################################
# PARTE "5": EXPRESIONES AUXILIARES          #
##############################################

# Cálculo Número de Estacionamientos
estacionamientosViviendas = sum(dcn.ESTACIONAMIENTOSPORVIV .* numDeptosTipo)
estacionamientosVisitas = estacionamientosViviendas * dcn.PORCADICESTACVISITAS;
estacionamientosNormales = estacionamientosViviendas + estacionamientosVisitas
estacionamientosDiscapacitados = (maxDeptos <= 20) ? 1 :
                                 ((maxDeptos <= 50) ? 2 :
                                 ((maxDeptos <= 200) ? 3 :
                                 ((maxDeptos <= 400) ? 4 :
                                 ((maxDeptos <= 500) ? 5 : (0.01*maxDeptos)))))
estacionamientosBicicletas = estacionamientosNormales * dcn.ESTBICICLETAPOREST;
descuentoEstCercaniaMetro = estacionamientosNormales * 0.5*dcn.REDUCCIONESTPORDISTMETRO
descuentoEstBicicletas = estacionamientosBicicletas/dcn.BICICLETASPOREST
estacionamientosNormales = estacionamientosNormales - descuentoEstCercaniaMetro - descuentoEstBicicletas
cambioEstBicicletas = estacionamientosNormales / 3
estacionamientosVendibles = estacionamientosViviendas - descuentoEstCercaniaMetro - descuentoEstBicicletas - cambioEstBicicletas
estacionamientosBicicletas = estacionamientosBicicletas + cambioEstBicicletas * dcn.BICICLETASPOREST
estacionamientos = estacionamientosVendibles + estacionamientosVisitas + estacionamientosDiscapacitados + estacionamientosBicicletas / dcn.BICICLETASPOREST;

# Cálculo de Superficies
superficieUtilDepto = numDeptosTipo .* dcc.SUPDEPTOUTIL;
superficieUtil = sum(superficieUtilDepto); # Superficie Util incluye terrazas.
superficieTerrazaDepto = 2 .* superficieUtilDepto .* dca.PORCTERRAZA;
superficieVendibleDepto = superficieUtilDepto .* (1-dca.PORCTERRAZA) + superficieTerrazaDepto * .5; # = superficieUtilDepto
superficieVendible = sum(superficieVendibleDepto);
superficieComun = dca.PORCSUPCOMUN * superficieUtil;
superficieLosaSNT = superficieVendible + superficieComun;
superficieLosaBNT = dcn.SUPPORESTACIONAMIENTO * estacionamientosVendibles;

# Cálculo de Ingresos por Ventas
ingresosVentaDeptos = sum(superficieVendibleDepto .* dcc.PRECIOVENTA);
ingresosVentaEstacionamientos = estacionamientosVendibles * dcc.PRECIOVENTAEST;
IngresosVentas = ingresosVentaDeptos + ingresosVentaEstacionamientos;

# Cálculo de Costos de Terreno y Construcción
CostoTerreno = CostoUnitTerreno * superficieTerreno * (1 + dcu.ComisionCorredor);
CostoConstruccionSNT = dcu.LosaSNT * superficieLosaSNT;
CostoConstruccionBNT = dcu.LosaBNT * superficieLosaBNT;
CostoEstacVisita = dcu.EstacVisita * dcn.SUPPORESTACIONAMIENTO * estacionamientosVisitas;
devolucionIVA = (CostoConstruccionSNT + CostoConstruccionBNT + CostoEstacVisita) / 1.19 * 0.19 * 0.65;
CostoConstruccion = dcu.Demolicion * superficieTerreno + CostoConstruccionSNT + CostoConstruccionBNT + CostoEstacVisita +
                    dcu.ExtrasPostVenta * (CostoConstruccionSNT + CostoConstruccionBNT) - devolucionIVA;

# Cálculo de Costos de Permisos, Habilitación, Puesta en Marcha y Atención de Clientes
CostoInmobiliarioObra = dcu.EmpalmesAportes * sum(numDeptosTipo) + dcu.DerechosPermisosMunicipales +
                    dcu.duracionProyecto * dcu.InspeccionCertificacionTecnica + dcu.Contribuciones * (dcu.duracionProyecto / 12) * 4 +
                    dcu.CopiaPlanosArquitctura + dcu.CopiaPlanosObra;
CostosHabilitacion = dcu.HallAcceso + dcu.JardinesPaisajismo + dcu.HabilitacionSalaUsosMultiples + dcu.HabilitacionPiscinaEnfermeria;
CostosPuestaEnMarcha = dcu.duracionProyecto * dcu.CuentasServiciosDeptos + dcu.duracionProyecto * dcu.GastosComunes +
                    dcu.ContribucionesDeptos + dcu.Seguros;
CostosAtencionCliente = dcu.ResponsabilidadCivil * CostoConstruccion + dcu.PostVentaInmobiliaria * CostoConstruccion +
                    dcu.SeguroVentaEnVerde * CostoConstruccion;

# Cálculo de Honorarios de Proyecto, Administración y Gestión de Venta
CostoHonorariosProyectos = dcu.Arquitectura * ingresosVentaDeptos + dcu.GestionAdministracion * ingresosVentaDeptos +
                    dcu.duracionProyecto * (dcu.ITO + dcu.Contabilidad) + dcu.Legales * CostoConstruccion +
                    dcu.Calculo * (superficieLosaSNT + superficieLosaBNT) + dcu.MecanicaSuelo + dcu.Topografia +
                    dcu.ProyectoElectrico * CostoConstruccion + dcu.ProyectoSanitario * CostoConstruccion +
                    dcu.ProyectoSeguridad * CostoConstruccion + dcu.ProyectoPavimentacion + dcu.ProyectoBasura + dcu.ProyectoCalefaccion;
CostoPromocionVenta = dcu.Marketing * IngresosVentas + dcu.PilotSalaVentas * IngresosVentas;
CostoOtrosVenta = dcu.NotarialesLegales * sum(numDeptosTipo) + dcu.CBR * sum(numDeptosTipo);
CostoVenta = CostoPromocionVenta + CostoOtrosVenta;

# Cálculo de Otros Costos
Imprevistos = dcu.Imprevistos * IngresosVentas;
DebitoIVA = (IngresosVentas - CostoTerreno) * 0.19 / 1.19;
CreditoIVA = (CostoConstruccion + CostoVenta + CostosHabilitacion) * 0.19 / 1.19;
CostoPagoIVA = DebitoIVA - CreditoIVA;

# Cálculo Costo Total
CostoTotal = CostoTerreno + CostoConstruccion + CostoHonorariosProyectos + CostoVenta + CostoInmobiliarioObra + CostosHabilitacion +
            CostosPuestaEnMarcha + CostosAtencionCliente + Imprevistos + CostoPagoIVA;

# Variables tipo flag asociadas a las celdas activas según configuración seleccionada
flagConfHor = configSel' * MATCONFHOR;
flagConfVert = configSel' * MATCONFVERT;
celdasConfig = configSel' * MATCELDASCONF;

# Variables auxiliares asociadas a la posición de los vertices de las celdas.
p1 = [(R[r].mat * ([pos_x'; pos_y'] - repeat(R[r].cr, 1, NumCeldas)) + repeat(R[r].cr, 1, NumCeldas))' for r=1:numLados];
p2 = [(R[r].mat * ([pos_x' + largo_x[:,r]';pos_y'] - repeat(R[r].cr, 1, NumCeldas)) + repeat(R[r].cr, 1, NumCeldas))' for r=1:numLados];
p3 = [(R[r].mat * ([pos_x' + largo_x[:,r]';pos_y' + largo_y[:,r]'] - repeat(R[r].cr, 1, NumCeldas)) + repeat(R[r].cr, 1, NumCeldas))' for r=1:numLados];
p4 = [(R[r].mat * ([pos_x'; pos_y' + largo_y[:,r]'] - repeat(R[r].cr, 1, NumCeldas)) + repeat(R[r].cr, 1, NumCeldas))' for r=1:numLados];

##############################################
# PARTE "6": RESTRICCIONES DEL MIP           #
##############################################

@constraints(m, begin
# Restricción de Altura Máxima y Área Basal Máxima (Coeficiente de Ocupación)
    [i = 1:NumCeldas], sum(pisoActivoCeldas[i,:]) * dca.ALTURAPISO <= altura[i]
    [i = 1:NumCeldas, z = 1:dcn.MAXPISOS - 1], pisoActivoCeldas[i, z + 1] <= pisoActivoCeldas[i, z]
    [z = 1:dcn.MAXPISOS], sum(pisoActivoCeldas[:, z]) <= pisoActivo[z] * 500
    areaPorPiso[z = 1:dcn.MAXPISOS], areaBasal[z] <= pisoActivo[z] * dcn.COEFOCUPACION * superficieTerreno
    maxPisosCeldas[i = 1:NumCeldas], sum(pisoActivoCeldas[i,:]) <= dcn.MAXPISOS
    alturaMaxCeldas[i = 1:NumCeldas], altura[i] <= celdasConfig[i] * min(dcn.ALTURAMAX, dcn.MAXPISOS*3.5)

# Restricciones que establecen relaciones entre areaBasal, numDeptos, largo y ancho
    [i = 1:NumCeldas, z = 1:dcn.MAXPISOS], areaCelda[z,i] <= 5000 * pisoActivoCeldas[i,z]
    [z = 1:dcn.MAXPISOS], areaCelda[z,:]  .<= largo .* dca.ANCHOMAX
    areaBasal .<= sum(areaCelda, dims = 2)
    superficieLosaSNT <= sum(areaBasal)
    maxCoefConstructibilidad, sum(areaBasal) <= superficieTerreno * dcn.COEFCONSTRUCTIBILIDAD * (1 + 0.3 * dcp.FUSIONTERRENOS)

# Restricciones de separación predial mínima
    restPredialPunto1[j = 1:NumCeldas, i = 1:numNonConv, r = 1:numLados], s_areaEdif[i].A * p1[r][j,:] .<= s_areaEdif[i].b .+ 500 * (1 - q1[j,i]) .+ 500 * (1 - d[r])
    restPredialPunto2[j = 1:NumCeldas, i = 1:numNonConv, r = 1:numLados], s_areaEdif[i].A * p2[r][j,:] .<= s_areaEdif[i].b .+ 500 * (1 - q2[j,i]) .+ 500 * (1 - d[r])
    restPredialPunto3[j = 1:NumCeldas, i = 1:numNonConv, r = 1:numLados], s_areaEdif[i].A * p3[r][j,:] .<= s_areaEdif[i].b .+ 500 * (1 - q3[j,i]) .+ 500 * (1 - d[r])
    restPredialPunto4[j = 1:NumCeldas, i = 1:numNonConv, r = 1:numLados], s_areaEdif[i].A * p4[r][j,:] .<= s_areaEdif[i].b .+ 500 * (1 - q4[j,i]) .+ 500 * (1 - d[r])
    [j = 1:NumCeldas], sum(q1[j,:]) == 1
    [j = 1:NumCeldas], sum(q2[j,:]) == 1
    [j = 1:NumCeldas], sum(q3[j,:]) == 1
    [j = 1:NumCeldas], sum(q4[j,:]) == 1

# Restricción de Rasante
    restRasantePunto1[j = 1:NumCeldas, i = 1:numNonConv, r = 1:numLados], s_bruto[i].A * dcn.RASANTE ./ s_bruto[i].proyeccion * p1[r][j,:] .+ altura[j] .<= s_bruto[i].b * dcn.RASANTE ./ s_bruto[i].proyeccion .+ 500 * (1 - w1[j,i]) .+ 500 * (1 - d[r])
    restRasantePunto2[j = 1:NumCeldas, i = 1:numNonConv, r = 1:numLados], s_bruto[i].A * dcn.RASANTE ./ s_bruto[i].proyeccion * p2[r][j,:] .+ altura[j] .<= s_bruto[i].b * dcn.RASANTE ./ s_bruto[i].proyeccion .+ 500 * (1 - w2[j,i]) .+ 500 * (1 - d[r])
    restRasantePunto3[j = 1:NumCeldas, i = 1:numNonConv, r = 1:numLados], s_bruto[i].A * dcn.RASANTE ./ s_bruto[i].proyeccion * p3[r][j,:] .+ altura[j] .<= s_bruto[i].b * dcn.RASANTE ./ s_bruto[i].proyeccion .+ 500 * (1 - w3[j,i]) .+ 500 * (1 - d[r])
    restRasantePunto4[j = 1:NumCeldas, i = 1:numNonConv, r = 1:numLados], s_bruto[i].A * dcn.RASANTE ./ s_bruto[i].proyeccion * p4[r][j,:] .+ altura[j] .<= s_bruto[i].b * dcn.RASANTE ./ s_bruto[i].proyeccion .+ 500 * (1 - w4[j,i]) .+ 500 * (1 - d[r])
    [j = 1:NumCeldas], sum(w1[j,:]) == 1
    [j = 1:NumCeldas], sum(w2[j,:]) == 1
    [j = 1:NumCeldas], sum(w3[j,:]) == 1
    [j = 1:NumCeldas], sum(w4[j,:]) == 1

# Restricción de Ancho máximo
    [r = 1:numLados], largo_x[:,r] .<= dca.ANCHOMAX .+ flagConfHor' * 500
    [r = 1:numLados], largo_y[:,r] .<= dca.ANCHOMAX .+ flagConfVert' * 500
    [r = 1:numLados], largo_x[:,r] .>= dca.ANCHOMAX * celdasConfig' .- 500 * (1 - d[r])
    [r = 1:numLados], largo_y[:,r] .>= dca.ANCHOMAX * celdasConfig' .- 500 * (1 - d[r])
    [r = 1:numLados], largo_x[:,r] .<= celdasConfig' * 500
    [r = 1:numLados], largo_y[:,r] .<= celdasConfig' * 500
    [r = 1:numLados], largo_x[:,r] .<= d[r] * 500
    [r = 1:numLados], largo_y[:,r] .<= d[r] * 500

# Restricciones de No Interesección de Celdas
    [r = 1:numLados], pos_x[2] + largo_x[2,r] <= pos_x[3]
    [r = 1:numLados], pos_x[3] + largo_x[3,r] <= pos_x[4]
    [r = 1:numLados], pos_y[2] + largo_y[2,r] <= pos_y[1]
    [r = 1:numLados], pos_y[3] + largo_y[3,r] <= pos_y[1]
    [r = 1:numLados], pos_y[4] + largo_y[4,r] <= pos_y[1]
    [r = 1:numLados], pos_y[5] + largo_y[5,r] <= pos_y[2]
    [r = 1:numLados], pos_y[5] + largo_y[5,r] <= pos_y[3]
    [r = 1:numLados], pos_y[5] + largo_y[5,r] <= pos_y[4]

# Restricciones de Rentabilidad Mínima
    restRentabilidadMin, IngresosVentas >= 1.20 * CostoTotal

# Restricciones que obligan a elegir una configuración y una rotación
    sum(configSel) == 1
    sum(d) == 1

# Restricción Suma de Departamentos menor a Máximo Número de Departamentos
    RestDensidadMax, sum(numDeptosTipo) <= maxDeptos

# Restricciones de Demanda
    RestDemanda, numDeptosTipo .<= maxDeptos .* dcc.MAXPORCTIPODEPTO

# Restricciones Auxiliares
    largo .<= celdasConfig' * 500
    largo .<= sum(largo_x, dims = 2) .+ flagConfVert' * 500
    largo .<= sum(largo_y, dims = 2) .+ flagConfHor' * 500

# Restricciones de Testeo
    #numDeptosTipo[1:end-1] .== 0; #Sólo deptos tipo 4
    #sum(configSel[1:2]) == 0  #Elimina configuraciones con 3 naves

end)

##############################################
# PARTE "7": FUNCIÓN OBJETIVO Y EJECUCIÓN    #
##############################################

@objective(m, Max, CostoUnitTerreno)
#@objective(m, Max, sum(numDeptosTipo))
#@objective(m, Max, IngresosVentas-CostoTotal)

# Resuelve el problema de optimización
JuMP.optimize!(m)

##############################################
# PARTE "8": PRESENTACION DE RESULTADOS      #
##############################################

# Cálculo del Flujo de Caja y la Tir
fl_IngresosVentas = dcf.IngresosVentas .* JuMP.value(IngresosVentas);
fl_CostoTerreno = dcf.CostoTerreno .* JuMP.value(CostoTerreno);
fl_CostoConstruccion = dcf.CostoConstruccion .* JuMP.value(CostoConstruccion);
fl_CostoHonorariosProyectos = dcf.CostoHonorariosProyectos .* JuMP.value(CostoHonorariosProyectos);
fl_CostoVenta = dcf.CostoVenta .* JuMP.value(CostoVenta);
fl_CostoInmobiliarioObra = dcf.CostoInmobiliarioObra .* JuMP.value(CostoInmobiliarioObra);
fl_CostosHabilitacion = dcf.CostosHabilitacion .* CostosHabilitacion;
fl_CostosPuestaEnMarcha = dcf.CostosPuestaEnMarcha .* CostosPuestaEnMarcha;
fl_CostosAtencionCliente = dcf.CostosAtencionCliente .* JuMP.value(CostosAtencionCliente);
fl_Imprevistos = dcf.Imprevistos .* JuMP.value(Imprevistos);
fl_CostoPagoIVA = fl_IngresosVentas .* (1 - JuMP.value(CostoTerreno) / sum(fl_IngresosVentas)) .* (0.19 / 1.19) -
                (fl_CostoConstruccion + fl_CostoVenta + fl_CostosHabilitacion) .* (0.19 / 1.19);
fl_FlujoCajaNetoAntesImpuesto = fl_IngresosVentas - fl_CostoTerreno - fl_CostoConstruccion - fl_CostoHonorariosProyectos - fl_CostoVenta -
                            fl_CostoInmobiliarioObra - fl_CostosHabilitacion - fl_CostosPuestaEnMarcha - fl_CostosAtencionCliente -
                            fl_Imprevistos - fl_CostoPagoIVA;
tirAntesImpuesto = (irr(fl_FlujoCajaNetoAntesImpuesto) + 1)^2 - 1;

numFlujos = length(fl_FlujoCajaNetoAntesImpuesto);
fl_FlujoCajaNetoAntesImpuestos_acum = cumsum(fl_FlujoCajaNetoAntesImpuesto);
vecCorrel = 1:numFlujos;
flag = cumsum(fl_FlujoCajaNetoAntesImpuestos_acum .> 0) .== 1;
flagCorrel = vecCorrel .* flag;
posNoNeg = findall(x->x > 0, flagCorrel)
fl_ImpuestoRenta = zeros(Float64, (numFlujos,));
fl_ImpuestoRenta[posNoNeg] = fl_FlujoCajaNetoAntesImpuestos_acum[posNoNeg] * dcf.TasaImpuestoRenta;
fl_FlujoCajaNetoDespuesImpuesto = fl_FlujoCajaNetoAntesImpuesto - fl_ImpuestoRenta;
tirDespuesImpuesto = (irr(fl_FlujoCajaNetoDespuesImpuesto) + 1)^2 - 1;

fl_LineaCredito = (dcf.ingresoLineaCredito - dcf.pagoLineaCredito) * JuMP.value(CostoConstruccion);
fl_MontoUtilizadoLineaCredito = cumsum(fl_LineaCredito);
fl_InteresLineaCredito = fl_MontoUtilizadoLineaCredito * dcf.TasaInteresLineaCredito / 2;
fl_CostoOperacionalCredito = fl_CostoConstruccion ./ sum(fl_CostoConstruccion) * dcf.DuracionLineaCredito * dcf.Tasaciones;
fl_FlujoCajaApalancadoNetoAntesImpuesto = fl_FlujoCajaNetoAntesImpuesto + fl_LineaCredito - fl_InteresLineaCredito - fl_CostoOperacionalCredito;
tirApalancadoAntesImpuestos = (irr(fl_FlujoCajaApalancadoNetoAntesImpuesto) + 1)^2 - 1;

# Construye estructura con los resultados de la optimización
sn = salidaNormativa(
    maxDeptos, # maxNumDeptos
    dcn.COEFOCUPACION * superficieTerreno, # maxOcupacion
    superficieTerreno * dcn.COEFCONSTRUCTIBILIDAD * (1 + 0.3 * dcp.FUSIONTERRENOS), # maxConstructibilidad
    dcn.MAXPISOS, # maxPisos
    dcn.ALTURAMAX, # maxAltura
    JuMP.value(estacionamientosVendibles), # minEstacionamientosVendible
    JuMP.value(estacionamientosVisitas), # minEstacionamientosVisita
    estacionamientosDiscapacitados # minEstacionamientosDiscapacitados
)

sa = salidaArquitectonica(
    JuMP.value.(numDeptosTipo), #numDeptosTipo
    sum(JuMP.value.(numDeptosTipo)), #numDeptos
    maximum(JuMP.value.(areaBasal)), #ocupacion
    sum(JuMP.value.(areaBasal)), #constructibilidad
    maximum(cumsum(JuMP.value.(pisoActivoCeldas),dims=2)), #numPisos
    maximum(JuMP.value.(altura)), #altura
    JuMP.value(superficieVendible), #superficieUtilSNT
    JuMP.value(superficieComun), #superficieComunSNT
    JuMP.value(superficieVendible)+JuMP.value(superficieComun), #superficieEdificadaSNT
    sort(JuMP.value.(areaBasal), rev=true), #superficiePorPiso
    JuMP.value(estacionamientosVendibles), #estacionamientosVendibles
    JuMP.value(estacionamientosVisitas), #estacionamientosVisita
    JuMP.value(estacionamientosVendibles)+JuMP.value(estacionamientosVisitas), #numEstacionamientos
    JuMP.value(estacionamientosBicicletas) #numBicicleteros
    )

si = salidaIndicadores(
    JuMP.value(IngresosVentas), #IngresosVentas
    JuMP.value(CostoTotal), #CostoTotal
    JuMP.value(IngresosVentas)-JuMP.value(CostoTotal), #MargenAntesImpuesto
    tirAntesImpuesto, #TirAntesImpuestos
    dcf.TasaImpuestoRenta*(JuMP.value(IngresosVentas)-JuMP.value(CostoTotal)), #ImpuestoRenta
    JuMP.value(IngresosVentas)-JuMP.value(CostoTotal) - dcf.TasaImpuestoRenta*(JuMP.value(IngresosVentas)-JuMP.value(CostoTotal)), #UtilidadDespuesImpuesto
    (JuMP.value(IngresosVentas)-JuMP.value(CostoTotal))/JuMP.value(CostoTotal), #RentabilidadTotalBruta
    (JuMP.value(IngresosVentas)-JuMP.value(CostoTotal) - dcf.TasaImpuestoRenta*(JuMP.value(IngresosVentas)-JuMP.value(CostoTotal)))/JuMP.value(CostoTotal), #RentabilidadTotalNeta
    JuMP.value(CostoTerreno)/JuMP.value(IngresosVentas), #IncidenciaTerreno
)

st = salidaTerreno(
    superficieTerreno, #superficieTerreno
    polyArea(V_bruto), #superficieBruta
    polyArea(V_areaEdif), #superficieEdificacion
    JuMP.value(CostoUnitTerreno * superficieTerreno), #costoTerrenoAntesCorr
    JuMP.value(CostoUnitTerreno), #costoUnitTerrenoAntesCorr
    JuMP.value(CostoUnitTerreno * superficieTerreno * dcu.ComisionCorredor), #CostoCorredor
    JuMP.value(CostoTerreno), #costoTerreno
    JuMP.value(CostoTerreno)/superficieTerreno #costoUnitTerreno
)

sm = salidaMonetaria(
    JuMP.value(IngresosVentas), #IngresosVentas
    JuMP.value(CostoTotal), #CostoTotal
    JuMP.value(CostoTerreno), #CostoTerreno
    JuMP.value(CostoTerreno)/superficieTerreno, #CostoUnitarioTerreno
    JuMP.value(CostoConstruccion), #CostoConstruccion
    JuMP.value(CostoHonorariosProyectos), #CostoHonorariosProyectos
    JuMP.value(CostoVenta), #CostoVenta
    JuMP.value(CostoInmobiliarioObra), #CostoInmobiliarioObra
    (CostosHabilitacion), #CostosHabilitacion
    (CostosPuestaEnMarcha), #CostosPuestaEnMarcha
    JuMP.value(CostosAtencionCliente), #CostosAtencionCliente
    JuMP.value(Imprevistos), #Imprevistos
    JuMP.value(CostoPagoIVA) #CostoPagoIVA
)

sf = salidaFlujoCaja(
    fl_IngresosVentas, #IngresosVentas
    fl_CostoTerreno, #CostoTerreno
    fl_CostoConstruccion, #CostoConstruccion
    fl_CostoHonorariosProyectos, #CostoHonorariosProyectos
    fl_CostoVenta, #CostoVenta
    fl_CostoInmobiliarioObra, #CostoInmobiliarioObra
    fl_CostosHabilitacion, #CostosHabilitacion
    fl_CostosPuestaEnMarcha, #CostosPuestaEnMarcha
    fl_CostosAtencionCliente, #CostosAtencionCliente
    fl_Imprevistos, #Imprevistos
    fl_CostoPagoIVA, #CostoPagoIVA
    fl_FlujoCajaNetoAntesImpuesto, #flujoCajaNetoAntesImpuestos
    tirAntesImpuesto, #TirAntesImpuestos
    fl_ImpuestoRenta, #ImpuestoRenta
    fl_FlujoCajaNetoDespuesImpuesto, #NetoDespuesImpuestos
    tirDespuesImpuesto, #TirDespuesImpuestos
    fl_LineaCredito, #LineaCredito
    fl_MontoUtilizadoLineaCredito, #MontoUtilizadoLineaCredito
    fl_InteresLineaCredito, #InteresLineaCredito
    fl_CostoOperacionalCredito, #CostoOperacionalCredito
    fl_FlujoCajaApalancadoNetoAntesImpuesto, #ApalancadoNetoAntesImpuesto
    tirApalancadoAntesImpuestos #TirApalancadoAntesImpuestos
)

vecPosNumLado = collect(1:numLados);
r_opt = Int(round(vecPosNumLado' * JuMP.value.(d)));
p1_vec = (R[r_opt].mat * ([pos_x'; pos_y'] - repeat(R[r_opt].cr, 1, NumCeldas)) + repeat(R[r_opt].cr, 1, NumCeldas))';
p2_vec = (R[r_opt].mat * ([pos_x' + largo_x[:,r_opt]';pos_y'] - repeat(R[r_opt].cr, 1, NumCeldas)) + repeat(R[r_opt].cr, 1, NumCeldas))';
p3_vec = (R[r_opt].mat * ([pos_x' + largo_x[:,r_opt]';pos_y' + largo_y[:,r_opt]'] - repeat(R[r_opt].cr, 1, NumCeldas)) + repeat(R[r_opt].cr, 1, NumCeldas))';
p4_vec = (R[r_opt].mat * ([pos_x'; pos_y' + largo_y[:,r_opt]'] - repeat(R[r_opt].cr, 1, NumCeldas)) + repeat(R[r_opt].cr, 1, NumCeldas))';
sg = salidaGeometria(
    JuMP.value.(configSel), #config
    vec(celdasConfig), #celdasConfig
    p1_vec, #p1
    p2_vec, #p2
    p3_vec, #p3
    p4_vec, #p4
    pos_x, #pos_x
    pos_y, #pos_y
    JuMP.value.(largo) #largo
);

println("Características Generales de la Cabida Óptima:")
println("----------------------------------------------")
println("N° Deptos: ", sum(round.(sa.numDeptosTipo)' .- round.([.6 .3 .1 0]*(sum(round.(sa.numDeptosTipo))-sum(sa.numDeptosTipo)))), " (Máx. Densidad = ", round(Int, sn.maxNumDeptos), ")")
println("N° Deptos por Tipo: ", round.(sa.numDeptosTipo)' .- round.([.6 .3 .1 0]*(sum(round.(sa.numDeptosTipo))-sum(sa.numDeptosTipo))))
println("Ocupacion: ", round(sa.ocupacion, digits=2), " m2 (Máx. Ocupación = ", round(sn.maxOcupacion, digits=2), " m2)")
println("Constructibilidad: ", round(sa.constructibilidad, digits=2), " m2 (Máx. Constructibilidad = ", round(sn.maxConstructibilidad, digits=2), " m2)")
println("N° Pisos: ", Int(round(sa.numPisos)), " (Máx. Pisos = ", Int(sn.maxPisos), ")")
println("Altura: ", round(sa.altura, digits=2), " m2 (Máx. Altura = ", round(sn.maxAltura, digits=2), " m2)")
println("Superficie Útil SNT: ", round(sa.superficieUtilSNT, digits=2), " m2")
println("Superficie Común SNT: ", round(sa.superficieComunSNT, digits=2), " m2")
println("Superficie Edificada SNT: ", round(sa.superficieEdificadaSNT, digits=2), " m2")
println("Superficie por Piso: ", round.(sa.superficiePorPiso, digits=2))
println("N° Estac. Vendibles: ", round(Int, sa.estacionamientosVendibles))
println("N° Estac. Visita: ", round(Int, sa.estacionamientosVisita))
println("N° Estac. Discapacitados: ", round(Int, sn.minEstacionamientosDiscapacitados))
println("N° Estac. Totales: ", round(Int, sa.estacionamientosVendibles) + round(Int, sa.estacionamientosVisita) + round(Int, sn.minEstacionamientosDiscapacitados))
println("N° Estac. Bicicletas: ", round(Int, sa.numBicicleteros))

println("")
println("Características del Terreno:")
println("----------------------------")
println("Superficie Terreno: ", round(st.superficieTerreno, digits=2), " m2")
println("Superficie Bruta: ", round(st.superficieBruta, digits=2), " m2")
println("Superficie Edificacion: ", round(st.superficieEdificacion, digits=2), " m2")
println("Costo Terreno antes Corredor: UF ", round(st.costoTerrenoAntesCorr, digits=2), " (", round(st.costoUnitTerrenoAntesCorr, digits=2), " UF/m2)")
println("Costo Corredor: UF ", round(st.CostoCorredor, digits=2), (" (IVA Inc.)"))
println("Costo Total Terreno: UF ", round(st.costoTerreno, digits=2), " (", round(st.costoUnitTerreno, digits=2), " UF/m2)")

println("")
println("Indicadores de Desempeño:")
println("-------------------------")
println("Ingresos por Ventas: UF ", round(si.IngresosVentas, digits=2))
println("Costo Total: UF ", round(si.CostoTotal, digits=2))
println("Margen antes de Impuesto: UF ", round(si.MargenAntesImpuesto, digits=2))
println("Rentabilidad Total Bruta: UF+ ", round(si.RentabilidadTotalBruta, digits=2)*100," %")
println("TIR antes Impuestos: UF+ ", round(si.TirAntesImpuestos, digits=2)*100," %")
println("Impuesto Renta: UF ", round(si.ImpuestoRenta, digits=2))
println("Utilidad después Impuesto: UF ", round(si.UtilidadDespuesImpuesto, digits=2))
println("Rentabilidad Total Neta: UF+ ", round(si.RentabilidadTotalNeta, digits=2)*100," %")
println("Incidencia Terreno: ", round(si.IncidenciaTerreno, digits=2))

# Grafica Predio
ps_predio = PolyShape([V_predio],1)
fig, ax, ax_mat = plotPolyshape3d_v2(ps_predio, 0, nothing, nothing, nothing, "green", 0.25)

# Grafica cabida óptima
vecPosNumCeldas = collect(1:NumCeldas);
celdasConfig = JuMP.value.(configSel)' * MATCELDASCONF;
setCeldasActivas = [i for i in Int.(round.(vecPosNumCeldas' .* celdasConfig)) if i >= 1]
fig, ax, ax_mat = plotCabidaOptima(setCeldasActivas, JuMP.value.(pisoActivoCeldas), JuMP.value.(pos_x), JuMP.value.(pos_y), JuMP.value.(largo_x), JuMP.value.(largo_y), R, r_opt, dca, fig, ax, ax_mat)

# Grafica Volumen Teórico
ps_volteor = PolyShape([V_volteor],1)
fig, ax, ax_mat = plotPolyshape3d_v2(ps_volteor, nothing, fig, ax, ax_mat, "grey", 0.25)