text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
struct Scenario <: AbstractScenario
    settings::Settings
    site::Site
    pvs::Array{PV, 1}
    wind::Wind
    storage::Storage
    electric_tariff::ElectricTariff
    electric_load::ElectricLoad
    electric_utility::ElectricUtility
    financial::Financial
    generator::Generator
    dhw_load::DomesticHotWaterLoad
    space_heating_load::SpaceHeatingLoad
    cooling_load::CoolingLoad
    existing_boiler::Union{ExistingBoiler, Nothing}
    chp::Union{CHP, Nothing}  # use nothing for more items when they are not modeled?
    flexible_hvac::Union{FlexibleHVAC, Nothing}
    existing_chiller::Union{ExistingChiller, Nothing}
    absorption_chiller::Union{AbsorptionChiller, Nothing}
end

"""
    Scenario(d::Dict; flex_hvac_from_json=false)

Constructor for Scenario struct, where `d` has upper-case keys:
- [Site](@ref) (required)
- [ElectricTariff](@ref) (required when off_grid_flag is False)
- [ElectricLoad](@ref) (required)
- [PV](@ref) (optional, can be Array)
- [Wind](@ref) (optional)
- [ElectricStorage](@ref) (optional)
- [ElectricUtility](@ref) (optional)
- [Financial](@ref) (optional)
- [Generator](@ref) (optional)
- [DomesticHotWaterLoad](@ref) (optional)
- [SpaceHeatingLoad](@ref) (optional)
- [ExistingBoiler](@ref) (optional)
- [CHP](@ref) (optional)
- FlexibleHVAC (optional)
- ExistingChiller (optional)
- AbsorptionChiller (optional)

All values of `d` are expected to be `Dicts` except for `PV`, which can be either a `Dict` or `Dict[]`.

Set `flex_hvac_from_json=true` if `FlexibleHVAC` values were loaded in from JSON (necessary to
handle conversion of Vector of Vectors from JSON to a Matrix in Julia).
"""
function Scenario(d::Dict; flex_hvac_from_json=false)
    if haskey(d, "Settings")
        settings = Settings(;dictkeys_tosymbols(d["Settings"])...)
    else
        settings = Settings()
    end

    site = Site(;dictkeys_tosymbols(d["Site"])...)

    # Check that only PV, electric storage, and generator are modeled for off-grid
    if settings.off_grid_flag
        offgrid_allowed_keys = ["PV", "ElectricStorage", "Generator", "Settings", "Site", "Financial", "ElectricLoad", "ElectricTariff", "ElectricUtility"]
        unallowed_keys = setdiff(keys(d), offgrid_allowed_keys)
        if !isempty(unallowed_keys)
            throw(@error "Currently, only PV, ElectricStorage, and Generator can be modeled when off_grid_flag is true. Cannot model $unallowed_keys.")
        end
    end

    pvs = PV[]
    if haskey(d, "PV")
        if typeof(d["PV"]) <: AbstractArray
            for (i, pv) in enumerate(d["PV"])
                check_pv_tilt!(pv, site)
                if !(haskey(pv, "name"))
                    pv["name"] = string("PV", i)
                end
                push!(pvs, PV(;dictkeys_tosymbols(pv)..., off_grid_flag = settings.off_grid_flag))
            end
        elseif typeof(d["PV"]) <: AbstractDict
            check_pv_tilt!(d["PV"], site)
            push!(pvs, PV(;dictkeys_tosymbols(d["PV"])..., off_grid_flag = settings.off_grid_flag))
        else
            error("PV input must be Dict or Dict[].")
        end
    end

    if haskey(d, "Financial")
        financial = Financial(; dictkeys_tosymbols(d["Financial"])..., off_grid_flag = settings.off_grid_flag )
    else
        financial = Financial(; off_grid_flag = settings.off_grid_flag)
    end

    if haskey(d, "ElectricUtility") && !(settings.off_grid_flag)
        electric_utility = ElectricUtility(; dictkeys_tosymbols(d["ElectricUtility"])...)
    elseif !(settings.off_grid_flag)
        electric_utility = ElectricUtility()
    elseif settings.off_grid_flag
        if haskey(d, "ElectricUtility")
            @warn "ElectricUtility inputs are not applicable when off_grid_flag is true and any ElectricUtility inputs will be ignored. For off-grid scenarios, a year-long outage will always be modeled."
        end
        electric_utility = ElectricUtility(; outage_start_time_step = 1, outage_end_time_step = settings.time_steps_per_hour * 8760)
    end

    storage_structs = Dict{String, AbstractStorage}()
    if haskey(d,  "ElectricStorage")
        storage_dict = dictkeys_tosymbols(d["ElectricStorage"])
        storage_dict[:off_grid_flag] = settings.off_grid_flag
    else
        storage_dict = Dict(:max_kw => 0.0)
    end
    storage_structs["ElectricStorage"] = ElectricStorage(storage_dict, financial)
    # TODO stop building ElectricStorage when it is not modeled by user
    #       (requires significant changes to constraints, variables)
    if haskey(d, "HotThermalStorage")
        params = HotThermalStorageDefaults(; dictkeys_tosymbols(d["HotThermalStorage"])...)
        storage_structs["HotThermalStorage"] = ThermalStorage(params, financial, settings.time_steps_per_hour)
    end
    if haskey(d, "ColdThermalStorage")
        params = ColdThermalStorageDefaults(; dictkeys_tosymbols(d["ColdThermalStorage"])...)
        storage_structs["ColdThermalStorage"] = ThermalStorage(params, financial, settings.time_steps_per_hour)
    end
    storage = Storage(storage_structs)

    electric_load = ElectricLoad(; dictkeys_tosymbols(d["ElectricLoad"])...,
                                   latitude=site.latitude, longitude=site.longitude,
                                   time_steps_per_hour=settings.time_steps_per_hour,
                                   off_grid_flag = settings.off_grid_flag
                                )

    if !(settings.off_grid_flag) # ElectricTariff only required for on-grid
        electric_tariff = ElectricTariff(; dictkeys_tosymbols(d["ElectricTariff"])...,
                                        year=electric_load.year,
                                        NEM=electric_utility.net_metering_limit_kw > 0,
                                        time_steps_per_hour=settings.time_steps_per_hour
                                        )
    else # if ElectricTariff inputs supplied for off-grid, will not be applied.
        if haskey(d, "ElectricTariff")
            @warn "ElectricTariff inputs are not applicable when off_grid_flag is true, and will be ignored."
        end
        electric_tariff = ElectricTariff(;  blended_annual_energy_rate = 0.0,
                                            blended_annual_demand_rate = 0.0,
                                            year=electric_load.year,
                                            time_steps_per_hour=settings.time_steps_per_hour
        )
    end

    if haskey(d, "Wind")
        wind = Wind(; dictkeys_tosymbols(d["Wind"])...,
                    average_elec_load=sum(electric_load.loads_kw) / length(electric_load.loads_kw))
    else
        wind = Wind(; max_kw=0)
    end

    if haskey(d, "Generator")
        generator = Generator(; dictkeys_tosymbols(d["Generator"])..., off_grid_flag=settings.off_grid_flag, analysis_years=financial.analysis_years)
    else
        generator = Generator(; max_kw=0)
    end

    max_heat_demand_kw = 0.0
    if haskey(d, "DomesticHotWaterLoad") && !haskey(d, "FlexibleHVAC")
        add_doe_reference_names_from_elec_to_thermal_loads(d["ElectricLoad"], d["DomesticHotWaterLoad"])
        dhw_load = DomesticHotWaterLoad(; dictkeys_tosymbols(d["DomesticHotWaterLoad"])...,
                                          latitude=site.latitude, longitude=site.longitude,
                                          time_steps_per_hour=settings.time_steps_per_hour
                                        )
        max_heat_demand_kw = maximum(dhw_load.loads_kw)
    else
        dhw_load = DomesticHotWaterLoad(;
            fuel_loads_mmbtu_per_hour=zeros(8760*settings.time_steps_per_hour),
            time_steps_per_hour=settings.time_steps_per_hour
        )
    end

    if haskey(d, "SpaceHeatingLoad") && !haskey(d, "FlexibleHVAC")
        add_doe_reference_names_from_elec_to_thermal_loads(d["ElectricLoad"], d["SpaceHeatingLoad"])
        space_heating_load = SpaceHeatingLoad(; dictkeys_tosymbols(d["SpaceHeatingLoad"])...,
                                                latitude=site.latitude, longitude=site.longitude,
                                                time_steps_per_hour=settings.time_steps_per_hour
                                              )

        max_heat_demand_kw = maximum(space_heating_load.loads_kw .+ max_heat_demand_kw)
    else
        space_heating_load = SpaceHeatingLoad(;
            fuel_loads_mmbtu_per_hour=zeros(8760*settings.time_steps_per_hour),
            time_steps_per_hour=settings.time_steps_per_hour
        )
    end

    flexible_hvac = nothing
    existing_boiler = nothing
    existing_chiller = nothing

    if haskey(d, "FlexibleHVAC")
        # TODO how to handle Matrix from JSON (to Dict) ?
        if flex_hvac_from_json
            flexible_hvac = FlexibleHVAC(d["FlexibleHVAC"])
        else
            flexible_hvac = FlexibleHVAC(; dictkeys_tosymbols(d["FlexibleHVAC"])...)
        end

        if sum(flexible_hvac.bau_hvac.existing_boiler_kw_thermal) ≈ 0.0 &&
            sum(flexible_hvac.bau_hvac.existing_chiller_kw_thermal) ≈ 0.0
            @warn "The FlexibleHVAC inputs indicate that no heating nor cooling is required. Not creating FlexibleHVAC model."
            flexible_hvac = nothing
        else
            # ExistingChiller and/or ExistingBoiler are added based on BAU_HVAC energy required to keep temperature within bounds
            if sum(flexible_hvac.bau_hvac.existing_boiler_kw_thermal) > 0
                boiler_inputs = Dict{Symbol, Any}()
                boiler_inputs[:max_heat_demand_kw] = maximum(flexible_hvac.bau_hvac.existing_boiler_kw_thermal)
                boiler_inputs[:time_steps_per_hour] = settings.time_steps_per_hour
                if haskey(d, "ExistingBoiler")
                    boiler_inputs = merge(boiler_inputs, dictkeys_tosymbols(d["ExistingBoiler"]))
                end
                existing_boiler = ExistingBoiler(; boiler_inputs...)
                # TODO automatically add CHP or other heating techs?
                # TODO increase max_thermal_factor_on_peak_load to allow more heating flexibility?
            end

            if sum(flexible_hvac.bau_hvac.existing_chiller_kw_thermal) > 0
                chiller_inputs = Dict{Symbol, Any}()
                chiller_inputs[:loads_kw_thermal] = flexible_hvac.bau_hvac.existing_chiller_kw_thermal
                if haskey(d, "ExistingChiller")
                    if !haskey(d["ExistingChiller"], "cop")
                        d["ExistingChiller"]["cop"] = get_existing_chiller_default_cop(; existing_chiller_max_thermal_factor_on_peak_load=1.25,
                                                                                loads_kw=nothing,
                                                                                loads_kw_thermal=chiller_inputs[:loads_kw_thermal])
                    end
                    chiller_inputs = merge(chiller_inputs, dictkeys_tosymbols(d["ExistingChiller"]))
                else
                    chiller_inputs[:cop] = get_existing_chiller_default_cop(; existing_chiller_max_thermal_factor_on_peak_load=1.25,
                                                                                loads_kw=nothing,
                                                                                loads_kw_thermal=chiller_inputs[:loads_kw_thermal])
                end
                existing_chiller = ExistingChiller(; chiller_inputs...)
            end

            if haskey(d, "SpaceHeatingLoad")
                @warn "Not using SpaceHeatingLoad because FlexibleHVAC was provided."
            end

            if haskey(d, "DomesticHotWaterLoad")
                @warn "Not using DomesticHotWaterLoad because FlexibleHVAC was provided."
            end

            if haskey(d, "CoolingLoad")
                @warn "Not using CoolingLoad because FlexibleHVAC was provided."
            end
        end
    end

    if max_heat_demand_kw > 0 && !haskey(d, "FlexibleHVAC")  # create ExistingBoler
        boiler_inputs = Dict{Symbol, Any}()
        boiler_inputs[:max_heat_demand_kw] = max_heat_demand_kw
        boiler_inputs[:time_steps_per_hour] = settings.time_steps_per_hour
        # If CHP is considered, prime_mover may inform the default boiler efficiency
        if haskey(d, "CHP")
            if haskey(d["CHP"], "prime_mover")
                boiler_inputs[:chp_prime_mover] = d["CHP"]["prime_mover"]
            end
        end
        if haskey(d, "ExistingBoiler")
            boiler_inputs = merge(boiler_inputs, dictkeys_tosymbols(d["ExistingBoiler"]))
        end
        existing_boiler = ExistingBoiler(; boiler_inputs...)
    end

    chp = nothing
    if haskey(d, "CHP")
        chp = CHP(d["CHP"])
    end

    max_cooling_demand_kw = 0
    if haskey(d, "CoolingLoad") && !haskey(d, "FlexibleHVAC")
        # Note, if thermal_loads_ton or one of the "...fraction(s)_of_electric_load" inputs is used for CoolingLoad, doe_reference_name is ignored
        add_doe_reference_names_from_elec_to_thermal_loads(d["ElectricLoad"], d["CoolingLoad"])
        d["CoolingLoad"]["site_electric_load_profile"] = electric_load.loads_kw
        # Pass ExistingChiller inputs which are used in CoolingLoad processing, if they exist
        ec_empty = ExistingChiller(; loads_kw_thermal=zeros(8760*settings.time_steps_per_hour))
        if !haskey(d, "ExistingChiller")
            d["CoolingLoad"]["existing_chiller_max_thermal_factor_on_peak_load"] = ec_empty.max_thermal_factor_on_peak_load
        else
            if haskey(d["ExistingChiller"], "cop")
                d["CoolingLoad"]["existing_chiller_cop"] = d["ExistingChiller"]["cop"]
            end
            if haskey(d["ExistingChiller"], "max_thermal_factor_on_peak_load")
                d["CoolingLoad"]["existing_chiller_max_thermal_factor_on_peak_load"] = d["ExistingChiller"]["max_thermal_factor_on_peak_load"]
            else
                d["CoolingLoad"]["existing_chiller_max_thermal_factor_on_peak_load"] = ec_empty.max_thermal_factor_on_peak_load
            end
        end
        cooling_load = CoolingLoad(; dictkeys_tosymbols(d["CoolingLoad"])...,
                                    latitude=site.latitude, longitude=site.longitude,
                                    time_steps_per_hour=settings.time_steps_per_hour
                                    )
        max_cooling_demand_kw = maximum(cooling_load.loads_kw_thermal)
    else
        cooling_load = CoolingLoad(;
            thermal_loads_ton=zeros(8760*settings.time_steps_per_hour),
            time_steps_per_hour=settings.time_steps_per_hour
        )
    end

    absorption_chiller = nothing
    if max_cooling_demand_kw > 0 && !haskey(d, "FlexibleHVAC")  # create ExistingChiller
        chiller_inputs = Dict{Symbol, Any}()
        chiller_inputs[:loads_kw_thermal] = cooling_load.loads_kw_thermal
        if haskey(d, "ExistingChiller")
            if !haskey(d["ExistingChiller"], "cop")
                d["ExistingChiller"]["cop"] = cooling_load.existing_chiller_cop
            end
            chiller_inputs = merge(chiller_inputs, dictkeys_tosymbols(d["ExistingChiller"]))
        else
            chiller_inputs[:cop] = 1.0
        end
        existing_chiller = ExistingChiller(; chiller_inputs...)

        if haskey(d, "AbsorptionChiller")
            absorption_chiller = AbsorptionChiller(; dictkeys_tosymbols(d["AbsorptionChiller"])...)
        end
    end

    return Scenario(
        settings,
        site,
        pvs,
        wind,
        storage,
        electric_tariff,
        electric_load,
        electric_utility,
        financial,
        generator,
        dhw_load,
        space_heating_load,
        cooling_load,
        existing_boiler,
        chp,
        flexible_hvac,
        existing_chiller,
        absorption_chiller
    )
end

"""
    Scenario(fp::String)

Consruct Scenario from filepath `fp` to JSON with keys aligned with the `Scenario(d::Dict)` method.
"""
function Scenario(fp::String)
    Scenario(JSON.parsefile(fp); flex_hvac_from_json=true)
end

function check_pv_tilt!(pv::Dict, site::Site)
    if !(haskey(pv, "tilt"))
        pv["tilt"] = site.latitude
    end
end

function add_doe_reference_names_from_elec_to_thermal_loads(elec::Dict, thermal::Dict)
    string_keys = [
        "doe_reference_name",
        "blended_doe_reference_names",
        "blended_doe_reference_percents",
    ]
    for k in string_keys
        if k in keys(elec)
            if !(k in keys(thermal)) || isempty(thermal[k])
                thermal[k] = elec[k]
            end
        end
    end
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
"""
    KernelScale(k::Kernel) <: Kernel

TODO
"""
struct KernelScale{K<:Kernel} <: Kernel
end

ext: jl


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.0
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # 主成分分析をする

# - PCA の例として MNIST 画像 dataset を主成分分析を行う

# ## Load modules

using MLDatasets
using Images
using Statistics
using LinearAlgebra
using Plots

# ## Download MNIST dataset

MNIST.download(;i_accept_the_terms_of_use=true)

train_x, train_y = MNIST.traindata()
test_x, test_y = MNIST.testdata()

# ## Study data

?MNIST.convert2image

MNIST.traintensor(1)|>MNIST.convert2image

train_x |> size # W, H, Num_data

N = size(train_x)[end]

train_x_nhw = permuteddimsview(train_x, (3, 2, 1)); # transpose so that it has NxHxW layout

id = 4
colorview(Gray, train_x_nhw[id, :, :]) |> display
println("label = $(train_y[id])")

# # Apply PCA to 7

sevens = train_x[:, :, train_y .== 7];

μ_sevens = dropdims(mean(sevens, dims=3),dims=3);

colorview(Gray, μ_sevens')

colorview(Gray,(sevens[:,:,1] .- μ_sevens)')

X_ = sum(sevens .- μ_sevens, dims=3)/√N
X = reshape(X_, 784)
X |> size

A = X * X'; # variance-covariance matrix w.r.t sample

F = eigen(A);
λ = F.values[end]
v = F.vectors[:, end];

# visualize the first principal vector
normalize_img(v) = (v .- minimum(v)) ./ (maximum(v)-minimum(v))
colorview(Gray, reshape(normalize_img(v),28,28)') |> display
v

# +
λs = F.values[end:-1:end-5]
vs = F.vectors[:, end:-1:end-5]
normalize(v) = (v .- minimum(v)) ./ (maximum(v)-minimum(v))
gs = []
for i in 1:size(vs)[2]
    v = reshape(vs[:, i] |> normalize,28,28)'
    push!(gs, colorview(Gray, v))
end

colorview.(Gray,gs)
# -

# # Apply PCA to whole dataset

train_x, train_y = MNIST.traindata()
μ_train_x = mean(train_x,dims=3)
X_ = sum(train_x .- μ_train_x, dims=3)/√N
X = reshape(X_, 784)
A = X * X';

F=eigen(A)
λ1 = F.values[end]
λ2 = F.values[end-1]
λ3 = F.values[end-2]
w1 = F.vectors[:, end]
w2 = F.vectors[:, end-1]
w3 = F.vectors[:, end-2];

# ## Project 2D space

p = plot()
data = reshape(train_x,(28^2, N))
for i in 1:10
    d = data[:,train_y .== (i-1)][:, 1:100]
    x1 = w1' * d
    x2 = w2' * d
    scatter!(p,x1,x2, label=false, color=palette(:tab10)[i])
end
p

# ## Project 3D space

# +
#=
If we load `ORCA` package or use plotting operation based on PlotlyJS backend over the sysimage generated by PackageCompiler using my method, it will cause the following error.
"""
InitError: IOError: could not spawn `'/usr/bin/xvfb-run /root/.julia/conda/3/envs/_ORCA_jl_/bin/orca' server --port=7982 --graph-only`: no such file or directory (ENOENT)
during initialization of module ORCA
"""
We will ignore the error at the moment, since there is no harm with respect to our purpose.
Here is my workaround "Just load `ORCA` and ignore the error above or plot a small dummy data in advance"
=#

using ORCA

# +
plotlyjs()

p = plot(xlabel="x", ylabel="y" , zlabel="z")
data = reshape(train_x,(28^2, N))
target = [0,1,2,3,4,5,6,7,8,9]
for i in target
    d = data[:,train_y .== i][:, 1:100]
    x1 = w1' * d
    x2 = w2' * d
    x3 = w3' * d
    color = palette(:tab10)[i+1]
    scatter!(p,
        x1, x2, x3,
        label=false, color=color,
        markersize = 2,
    )
end

display(palette(:tab10)) # show color palette
p


text:
# Try to use Wang-Landau to estimate "density of states" of "energy"=path length

using BilliardModels

function find_free_path(billiard_table, x, v)
    l = Vector2D(0, 0)
    p = ParticleOnLattice(x, v, l)
    num_collisions = 1
    xs, ls, free_paths = billiard_dynamics_on_lattice(p, billiard_table, num_collisions)

    free_paths[1]
end

unif_rand(delta) = delta*(rand()-0.5)

function propose(x, y, theta, radius, delta)
    x_new = x + unif_rand(delta)
    y_new = y + unif_rand(delta)

     # check outside disc:
    if x*x + y*y < radius^2
        x_new, y_new = x, y   # return to where it was
    else
        # periodize:
        (x_new > 0.5) && (x_new -= 1)
        (x_new < -0.5) && (x_new += 1)
        (y > 0.5) && (y -= 1)
        (y < -0.5) && (y += 1)
    end

    theta += unif_rand(delta)  # in any case change theta

    x, y, theta
end

acceptance_prob(x, y, theta) = 1.0  # standard random walk

function explore_phase_space(radius=0.1, delta=0.1, num_steps=1000)

    billiard_table = Sinai_billiard(radius, true, true)  # periodic in x and y

    xx, vv = initial_condition(billiard_table, -.5, .5, -.5, .5)

    theta = atan2(vv[2], vv[1])  # this is how I generated the velocity to start with (from an angle) -- a bit ridiculous

    xs = [xx[1]]
    ys = [xx[2]]
    thetas = [theta]
    free_paths = [find_free_path(billiard_table, xx, vv)]

    num_free_path_bins = 1000
    free_path_distribution = zeros(num_free_path_bins)

    for i in 1:num_steps
        x, y = xx[1], xx[2]

        x_new, y_new, theta_new = propose(x, y, theta, radius, delta)

        if rand() < acceptance_prob(x, y, theta)
            x, y, theta = x_new, y_new, theta_new
        end

        xx_new = Vector2D(x, y)
        vv = Vector2D(cos(theta), sin(theta))

        free_path = find_free_path(billiard_table, xx_new, vv)

        push!(xs, x)
        push!(ys, y)
        push!(thetas, theta)
        push!(free_paths, free_path)

        free_path_bin = ceil(free_path)
        (free_path_bin > num_free_path_bins) && (free_path_bin = num_free_path_bins)

        free_path_distribution[free_path_bin] += 1

        xx = xx_new

    end

    xs, ys, thetas, free_paths, free_path_distribution
end

max_stars_repo_name: dpsanders/BilliardModel

max_issues_repo_name: dpsanders/BilliardModel

max_forks_repo_name: dpsanders/BilliardModel


text:
using Pkg; Pkg.develop(url="https://github.com/timholy/QuadDIRECT.jl.git");
using GalacticQuadDIRECT, GalacticOptim
using Test

@testset "GalacticQuadDIRECT.jl" begin
    rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
    x0 = zeros(2)
    _p = [1.0, 100.0]
    l1 = rosenbrock(x0, _p)
    optprob = OptimizationFunction(rosenbrock)
    prob = OptimizationProblem(optprob, x0, _p, lb=[-1.0, -1.0], ub=[0.8, 0.8])

    sol = solve(prob, QuadDirect(); splits=([-0.5, 0.0, 0.5], [-0.5, 0.0, 0.5]))
    @test 10 * sol.minimum < l1
end


text:
using Test
using DataStructure101
const DS = DataStructure101

@testset "Test Stack" begin

    # Initial a Stack with data '0'
    T_e = Int64
    T = DS.Stack{T_e}
    stake = T(0)

    # Type test
    @test stake isa T

    # Push test: Push data into stake 20 times by 1:20
    data= T_e[0]
    passed_check_data = []
    for i = 1:20
        push!(stake, i)
        pushfirst!(data, i)
        push!(passed_check_data, stake.data == data)
    end
    @test all(passed_check_data)

    # Pop test: Pop data from stack 5 times
    passed_check_data = []
    for i = 1:5
        pop!(stake)
        deleteat!(data, 1)
        push!(passed_check_data, stake.data == data)
    end
    @test all(passed_check_data)

    # Show Test
    @test repr(stake) == "Stack([15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0])"

end


text:
using Pkg
Pkg.activate("./")
ENV["JULIA_MPI_BINARY"]="system"; using Pkg; Pkg.build("MPI"; verbose=true)

max_stars_repo_name: mcreel/ILSNM

max_issues_repo_name: mcreel/JuliaMPIMonteCarlo

max_forks_repo_name: mcreel/JuliaMPIMonteCarlo


text:
"""
    @record command_buffer [create_info] commands

Convenience macro for recording a sequence of API commands into a command buffer `cbuff`.
All calls to API commands have `cbuff` inserted as their first argument, and are wrapped inside
`begin_command_buffer(cbuff, info)` and `end_command_buffer(cbuff)`.

!!! warning
    An expression is assumed to be an API command if it begins with `cmd_`.
    Make sure that all functions that you call satisfy this assumption.

The two-argument version of this macro simply passes in a default `CommandBufferBeginInfo()`.
"""
macro record(cbuff, info, cmds)
    api_calls = postwalk(cmds) do ex
        @match ex begin
            :($f($(args...))) => startswith(string(f), "cmd_") ? :($f($cbuff, $(args...))) : ex
            :($f($(args...); $(kwargs...))) => startswith(string(f), "cmd_") ? :($f($cbuff, $(args...); $(kwargs...))) : ex
            _ => ex
        end
    end
    quote
        $(esc(:(begin_command_buffer($(esc(cbuff)), $(esc(info))))))
        $(esc(api_calls))
        $(esc(:(end_command_buffer($(esc(cbuff))))))
    end
end

macro record(cbuff, cmds)
    :(@record $cbuff CommandBufferBeginInfo() $(esc(cmds)))
end

"""
Manages one command pool per thread.
"""
struct ThreadedCommandPool
    pools::Dictionary{Int,CommandPool}
end

function ThreadedCommandPool(device, disp::QueueDispatch, usage = dictionary(1:Base.Threads.nthreads() .=> QUEUE_COMPUTE_BIT | QUEUE_GRAPHICS_BIT | QUEUE_TRANSFER_BIT))
    pools = dictionary(thread => CommandPool(device, info(queue(disp, usage[thread])).queue_family_index) for thread in keys(usage))
    ThreadedCommandPool(pools)
end

function Vulkan.CommandPool(pool::ThreadedCommandPool)
    pool.pools[Base.Threads.threadid()]
end

Base.convert(T::Type{CommandPool}, pool::ThreadedCommandPool) = T(pool)


text:
using DashBootstrapComponents, DashHtmlComponents

radioitems = dbc_formgroup([
    dbc_label("Choose one"),
    dbc_radioitems(
        options=[
            Dict("label" => "Option 1", "value" => 1),
            Dict("label" => "Option 2", "value" => 2),
            Dict("label" => "Disabled Option", "value" => 3, "disabled" => true),
        ],
        value=1,
        id="radioitems-input",
    ),
]);

checklist = dbc_formgroup([
    dbc_label("Choose a bunch"),
    dbc_checklist(
        options=[
            Dict("label" => "Option 1", "value" => 1),
            Dict("label" => "Option 2", "value" => 2),
            Dict("label" => "Disabled Option", "value" => 3, "disabled" => true),
        ],
        value=[1],
        id="checklist-input",
    ),
]);

switches = dbc_formgroup([
    dbc_label("Toggle a bunch"),
    dbc_checklist(
        options=[
            Dict("label" => "Option 1", "value" => 1),
            Dict("label" => "Option 2", "value" => 2),
            Dict("label" => "Disabled Option", "value" => 3, "disabled" => true),
        ],
        value=[1],
        id="switches-input",
        switch=true,
    ),
]);

inputs = html_div([
    dbc_form([radioitems, checklist, switches]),
    html_p(id="radioitems-checklist-output"),
]);

callback!(
    app,
    Output("radioitems-checklist-output", "children"),
    Input("radioitems-input", "value"),
    Input("checklist-input", "value"),
    Input("switches-input", "value"),
) do radio_items_value, checklist_value, switches_value

    n_checkboxes = length(checklist_value)
    n_switches = length(switches_value)

    output_string =
        """Radio button $radio_items_value, $n_checkboxes """ *
        """checklist item$(n_checkboxes != 1 ? "s" : "") and $n_switches """ *
        """switch$(n_switches != 1 ? "es" : "") selected."""

    return output_string
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
function twitter(id)
    json = fetch_twitter(id)
    return json[:html]
end

@memoize function fetch_twitter(id)
    url = "https://publish.twitter.com/oembed?url=https://twitter.com/andypiper/status/$id"
    response = HTTP.get(url)
    json = JSON3.read(String(response.body))
end

struct Twitter <: ShortCode
    id::Int64
end

function Base.show(io::IO, ::MIME"text/plain", tweet::Twitter)
    json = fetch_twitter(tweet.id)
    tweet_as_text = "$(json[:author_name]): $(json[:url])"
    print(io, tweet_as_text)
end

function Base.show(io::IO, ::MIME"text/html", tweet::Twitter)
    print(io, twitter(tweet.id))
end


text:
## extend logical operators so that we can have a notation like

#f(x,y) < 0 or f(x,y) == 0

#import Base: <, <=, ==, !=, !==, >=, >
import Base: &, |, !

abstract type Predicate end

"""

A predicate is defined in terms of a function of two variables, an
inquality, and either another function or a real number.  They are
conveniently created by the functions `Lt`, `Le`, `Eq`, `Neq`, `Ge`,
and `Gt`. The equivalent unicode operators:

* `≪` (`\\ll[tab]`),
* `≦` (`\\leqq[tab]`),
* `⩵` (`\\Equal[tab]`),
* `≶` (`\\lessgtr[tab]`)  or `≷` (`\\gtrless[tab]`),
* `≧` (`\\geqq[tab]`),
* `≫` (`\\leqq[tab]`) may also be used.

The use of Julia's usual comparison operators is no longer supported.

To combine predicates, `&` and `|` can be used.

To negate a predicate, `!` is used.

"""
mutable struct Pred <: Predicate
    f::Function
    op
    val
end

## meta these
preds = [(:Lt, :≪, :<), # \ll
         (:Le, :≦, :<=), # \leqq
         (:Eq, :⩵, :(==)), # \Equal
         (:Ge, :≧, :>=), # \gegg
         (:Gt, :≫, :>) # \gg
         ]

for (fn, uop, op) in preds
    fnname =  string(fn)
    @eval begin
        ($fn)(f::Function, x::Real) = Pred(f, $op, x)
        ($uop)(f::Function, x::Real) = ($fn)(f, x)
        ($fn)(f::Function, g::Function) = $(fn)((x,y) -> f(x,y) - g(x,y), 0)
        ($uop)(f::Function, g::Function) = ($fn)(f, g)
    end
    eval(Expr(:export, fn))
    eval(Expr(:export, uop))
end

Neq(f::Function, x::Real) = Pred(f, !== , x)
Neq(f::Function, g::Function) = Neq((x,y) -> f(x,y) - g(x,y), 0)

≶(x::Real, y::Real) = (x != y)
≶(f::Function, x::Real) = Neq(f, x)
≶(f::Function, g::Function) = Neq(f, g)

≷(x::Real, y::Real) = (x != y)
≷(f::Function, x::Real) = Neq(f, x)
≷(f::Function, g::Function) = Neq(f, g)

"""

Predicates can be joined together with either `&` or `|`. Individual
predicates can be negated with `!`. The parsing rules require the
individual predicates to be enclosed with parentheses, as in `(f==0) | (g==0)`.

"""
mutable struct Preds <: Predicate
    ps
    ops
end

## Some algebra for Pred and Preds
(&)(r1::Pred, r2::Pred) = Preds([r1,r2], Any[&])
(|)(r1::Pred, r2::Pred) = Preds([r1,r2], Any[|])

(&)(ps::Preds, r1::Pred) = Preds([ps.ps, r1], [ps.ops, &])
(&)(r1::Pred, ps::Preds) = ps & r1
(|)(ps::Preds, r1::Pred) = Preds(vcat(ps.ps, r1), vcat(ps.ops, |))
(|)(r1::Pred, ps::Preds) = ps | r1


text:
module QCompEmulator

export my_f

end


text:
# Prelim 1 CHEME 5440/7770 Cornell University
# Jack Crowley, April 2021

# load the include file -
include("Include.jl")

# main method -
function main(parameters_file_path::String, data_file_path::String)

    # initialize
    Dataset = readdlm(data_file_path, ',', Float64)
    Concentration_35AMP_data = Dataset[:,1]             # units: mM
    V_data = Dataset[:,2]                               # units: uM/hr
    Confidence_Int_data = Dataset[:,3]                  # units: uM/hr

    N = length(Concentration_35AMP_data)
    Sim = zeros(N)
    Sim2 = zeros(N)
    Data = zeros(N)

    # Part B:  Estimate fitting constants from experimental dataset
    # Choice of state weights, effector binding constant, and order parameter
    Weights = [1.0, 0.05, 5000.0]
    K_35AMP = 190.0
    n_35AMP = 10.0
    Weights2 = [4600.0, 228.0, 1.5*16240.0]
    K_35AMP2 = 160.0
    n_35AMP2 = 4.2

    try

        # load the parameter file, and create the problem dictionary -
        problem_dictionary = generate_problem_dictionary(parameters_file_path)
        print("F6P Concentration = $(problem_dictionary["Fructose_6_Phosphate_concentration"]) uM\n")

        # TODO:  Solve things
        for (index, x) ∈ enumerate(Concentration_35AMP_data)
            Returned_rate = calculate_reaction_rate(problem_dictionary, Weights, K_35AMP, n_35AMP, x*1000)
            Sim[index] = Returned_rate
            print("\n\nReturned_rate = $(Returned_rate)")
            Returned_rate2 = calculate_reaction_rate(problem_dictionary, Weights2, K_35AMP2, n_35AMP2, x*1000)
            Sim2[index] = Returned_rate2

            Data[index] = V_data[index]
            print("\nExperimental rate = $(V_data[index])")
        end

        # Part C:  Plot data and mathematical model
        plot(Concentration_35AMP_data, Data, st = :line, color = "black", title="PFK-Catalyzed Reaction Rate", label = "Experimental")
        plot!(Concentration_35AMP_data, Sim, st = :line, color = "red", legend=:right, label = "Computed")
        plot!(Concentration_35AMP_data, Sim2, st = :line, color = "blue", label = "Params Estimated")
        #title("PFK-Catalyzed Reaction Rate")
        xlabel!("[3'-5' AMP] (mM)")
        ylabel!("r̂ (uM/hr)")
        savefig(path_to_plot)

    catch error

        # grab the error message, and post -
        @error "Oooops! Hmmmmm. Something happend" exception=(error, catch_backtrace())
    end

    return nothing
end

# setup the path to the parameters file -
path_to_parameters_file = joinpath(_PATH_TO_CONFIG, "Parameters.toml")

# setup the path to the experimental data file
path_to_data_file = joinpath(_PATH_TO_CONFIG, "Data-3-5-AMP.csv")

# execute
results_array = main(path_to_parameters_file, path_to_data_file)


text:
using Test
using FvCFD

include("shockTubeFunctions.jl")
include("test.jl")
include("testMeshes.jl")

@testset "Vector Utilities" begin
    v1 = [1, 2, 3]
    @test FvCFD.mag(v1) == 3.7416573867739413

    v2 = [2, 3, 4]
    @test FvCFD.dot(v1, v2) == 20.0

    @test almostEqual(FvCFD.normalize(v1), [0.267261242, 0.534522484, 0.801783726], 1e-9)

    v1 = [ 1.0, 2.0, 3.0 ]
    v2 = [ 2.0, 2.0, 2.0 ]
    crossProduct = [ -2.0, 4.0, -2.0 ]
    crossProduct2 = FvCFD.cross(v1, v2)
    @test almostEqual(crossProduct, crossProduct2)
end;

@testset "Shock Tube Initialization" begin
    nCells = 4

    # FDM
    dx = [ 0.25, 0.25, 0.25, 0.25 ]
    P = [ 1, 1, 0.1, 0.1 ]
    T = [ 0.00348432, 0.00348432, 0.00278746, 0.00278746 ]
    U = [ 0, 0, 0, 0 ]
    dx2, P2, T2, U2 = initializeShockTubeFDM(4)
    @test almostEqual(dx2, dx)
    @test almostEqual(P2, P)
    @test almostEqual(T2, T)
    @test almostEqual(U2, U)

    # FVM
    U = [ [0.0,0.0,0.0],  [0.0,0.0,0.0],  [0.0,0.0,0.0],  [0.0,0.0,0.0] ]
    cells = [ [4,1], [1,2], [2,3], [3,5] ]
    cVols = [ 0.0025, 0.0025, 0.0025, 0.0025 ]
    cCenters = [ [0.125,0.0,0.0], [0.375,0.0,0.0], [0.625,0.0,0.0], [0.875,0.0,0.0] ]
    faces = [ [1,2], [2,3], [3,4], [-1,1], [4,-1] ]
    fAVecs = [ [0.01,0,0], [0.01,0,0], [0.01,0,0], [0.01,0,0], [0.01,0,0] ]
    fCenters = [ [0.25,0.0,0.0], [0.5,0.0,0.0], [0.75,0.0,0.0], [0.0,0.0,0.0], [1.0, 0.0, 0.0] ]
    boundaryFaces = [ [4,], [5,] ]
    mesh1 = [ cells, cVols, cCenters, faces, fAVecs, fCenters, boundaryFaces ]

    mesh2, P2, T2, U2 = initializeShockTubeFVM(4)
    @test almostEqual(P2, P)
    @test almostEqual(T2, T)
    for i in [ 1,3,4,5,6,7 ]
        @test size(mesh1[i], 1) == size(mesh2[i], 1)
        for a in 1:size(mesh1[i], 1)
            @test almostEqual(mesh1[i][a], mesh2[i][a])
        end
    end
    @test almostEqual(mesh1[2], mesh2[2])
end;

@testset "Decode Primitives" begin
    rho = 2
    xMom = 1
    eV2 = 3

    fluid = FvCFD.Fluid(1005, 287.05, 1.4)

    U = 0.5
    e = 1.5 - 0.125
    T = e/(fluid.Cp - fluid.R)
    P = rho*fluid.R*T
    expectedResult = [ P, T, U, 0, 0 ]

    primitives = zeros(5)
    state = [ rho, xMom, 0, 0, eV2] # y-momentum and z-momentum = 0

    FvCFD.decodePrimitives3D!(primitives, state, fluid)

    for i in 1:5
        @test almostEqual(expectedResult[i], primitives[i])
    end
end;

@testset "Encode Primitives" begin
    rho = 2
    xMom = 1
    eV2 = 3

    fluid = FvCFD.Fluid(1005, 287.05, 1.4)

    U = 0.5
    e = 1.5 - 0.125
    T = e/(fluid.Cp - fluid.R)
    P = rho*fluid.R*T
    expectedResult = [ rho, xMom, 0, 0, eV2 ]

    primitives = Matrix{Float64}(undef, 1, 5)
    primitives[1,:] = [ P, T, U, 0, 0 ]

    state = FvCFD.encodePrimitives3D(primitives, fluid)
    for i in 1:5
        @test almostEqual(expectedResult[i], state[1,i])
    end
end;

@testset "Green-Gauss Gradient" begin
    #### Simple uniform gradient test case ####
    mesh, cellPrimitives = initializeShockTube3DFVM(4)
    faceVals = zeros(21,1) # 4 internal faces, 2 on the ends, and 16 around the outside
    faceVals[:,1] = [ 2, 3, 4, 1, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ] # Faces are numbered such that boundary faces come last, spatially values increase from one to five from left to right
    expectedGradient = [ [4.0,0.0,0.0], [4.0,0.0,0.0], [4.0,0.0,0.0], [4.0,0.0,0.0] ] # Expect uniform gradient
    gradient = FvCFD.greenGaussGrad(mesh, faceVals, true)

    for i in 1:4
        @test almostEqual(expectedGradient[i], gradient[i,1,:])
    end

    ########## Case from Moukalled pg. 249 ##########
    cells, cVols, cCenters, faces, fAVecs, fCenters, boundaryFaces = meshMoukalled813()
    # Adjust properties to match fudged values in Moukalled
    cCenters = [ [1.75, 2.0, 0.5], [4.25, 3.5, 0.5] ]
    cells = [ [1, 2, 3, 4, 5 ], [1], [2], [3], [4], [5] ]
    cVols = [ 8.625, 1, 1, 1, 1, 1 ]
    faces = [ [1,2], [1,3], [1,4], [1,5], [1,6] ]
    fCenters = fCenters[1:5]
    cSizes = zeros(6,3)
    mesh = Mesh( cells, cVols, cCenters, cSizes, faces, fAVecs, fCenters, boundaryFaces )

    function samplePhiDistribution(x, y, z)
        return x^2 + y^2 + x^2*y^2
    end

    phiFaceVals = zeros(5,1)
    for i in 1:5
        phiFaceVals[i,1] = samplePhiDistribution(fCenters[i]...)
    end

    grad = FvCFD.greenGaussGrad(mesh, phiFaceVals, true)
    @test almostEqual(grad[1,1,:], [ 20.63111, 17.31454, 0.0 ], 1e-5)
end;

### Laplacian calculation not required for Euler calculations ####
# Tests need to be fixed up if moving towards Navier-Stokes
# @testset "Laplacian" begin
    # cells, cVols, cCenters, faces, fAVecs, fCenters, boundaryFaces = meshMoukalled813()
    # # Adjust properties to match fudged values in Moukalled
    # cCenters = [ [1.75, 2.0, 0.5], [4.25, 3.5, 0.5], [0,0,0], [0,0,0], [0,0,0], [0,0,0] ]
    # cells = [ [1, 2, 3, 4, 5 ], [1], [2], [3], [4], [5] ]
    # cVols = [ 8.625, 1, 1, 1, 1, 1 ]
    # faces = [ [1,2] ]
    # fCenters = fCenters[1:5]
    # cSizes = zeros(5,3)
    # mesh = Mesh( cells, cVols, cCenters, cSizes, faces, fAVecs, fCenters, boundaryFaces )

    # oldGradient = [ [20.63111, 17.31454, 0.0], [112.625, 133.4375, 0.0], [0,0,0], [0,0,0], [0,0,0] ]
    # gradVals = Array{Array{Float64,1},2}(undef, 5, 1)
    # for i in 1:5
    #     gradVals[i,1] = oldGradient[i]
    # end
    # faceGrad = linInterp_3D(mesh, gradVals)

    # Sf = [ 2.5, 1.0, 0.0 ]
    # n = normalize(Sf)
    # CF = (cCenters[2] - cCenters[1])
    # e = normalize(CF)

    # #### Test Calculation of Ef ####
    # Ef2 = laplacian_Ef("None", Sf, n, CF, e)
    # Ef = Sf
    # @test Ef == Ef2
    # Ef2 = laplacian_Ef("MinCorr", Sf, n, CF, e)
    # Ef = [ 2.279, 1.368, 0.0 ]
    # @test almostEqual(Ef2, Ef, 3)
    # Ef2 = laplacian_Ef("OrthogCorr", Sf, n, CF, e)
    # Ef = [ 2.309, 1.385, 0.0 ]
    # @test almostEqual(Ef2, Ef, 3)
    # Ef2 = laplacian_Ef("OverRelax", Sf, n, CF, e)
    # Ef = [ 2.339, 1.403, 0.0 ]
    # @test almostEqual(Ef2, Ef, 3)

    # #### Test Calculation of Laplacian ####
    # # Note that values here are the negative of those presented in Moukalled, since we are interested in only the Laplacian,
    # # not the diffusive flux, which is related to the negative Laplacian
    # dPhi = 251.578125 - 19.3125
    # orthoGrad = dPhi / mag(CF)

    # faceIntegral = 0.9122*dPhi - 13.014
    # faceIntegral2 = laplacian_FaceFlux(mesh, "MinCorr", 1, faceGrad[1], orthoGrad)
    # @test almostEqual(faceIntegral, faceIntegral2, 1)
    # faceIntegral = 0.924*dPhi - 16.294
    # faceIntegral2 = laplacian_FaceFlux(mesh, "OrthogCorr", 1, faceGrad[1], orthoGrad)
    # @test almostEqual(faceIntegral, faceIntegral2, 1)
    # faceIntegral = 0.936*dPhi - 19.649
    # faceIntegral2 = laplacian_FaceFlux(mesh, "OverRelax", 1, faceGrad[1], orthoGrad)
    # @test almostEqual(faceIntegral, faceIntegral2, 1)

    # mesh, P, T, U = initializeShockTubeFVM(4)
    # Vals = [ 1, 2, 4, 8 ]
    # faceGrads = [ 4, 8, 16, 0, 0 ]
    # lapl = [ 0.0, 16.0, 32.0, 0.0 ]
    # lapl2 = laplacian(mesh, "None", Vals)[1]
    # for i in 1:4
    #     @test almostEqual(lapl[i], lapl2[i])
    # end
# end;

@testset "Linear Interpolation" begin
    mesh, cellPrimitives = initializeShockTube3DFVM(4)
    cellVals = zeros(4,1)
    cellVals[:,1] = [ 1, 2, 3, 4 ]
    faceVals = [ 1.5, 2.5, 3.5, 0, 0 ]
    faceVals2 = FvCFD.linInterp_3D(mesh, cellVals)
    @test almostEqual(faceVals, faceVals2[1:5,1])

    cells, cVols, cCenters, faces, fAVecs, fCenters, boundaryFaces = meshMoukalled813()
    # Adjust properties to match fudged values in Moukalled
    cCenters = [ [1.75, 2.0, 0.5], [4.25, 3.5, 0.5], [0,0,0], [0,0,0], [0,0,0], [0,0,0] ]
    cells = [ [1, 2, 3, 4, 5 ], [1], [2], [3], [4], [5] ]
    cVols = [ 8.625, 1, 1, 1, 1, 1 ]
    faces = [ [1,2] ]
    fCenters = fCenters[1:5]
    cSizes = zeros(5,3)
    mesh = Mesh( cells, cVols, cCenters, cSizes, faces, fAVecs, fCenters, boundaryFaces )

    oldGradient = [ [20.63111, 17.31454, 0.0], [112.625, 133.4375, 0.0], [0,0,0], [0,0,0], [0,0,0] ]
    gradients = zeros(5,3)
    for i in 1:5
        gradients[i,:] = oldGradient[i]
    end
    faceGrad = FvCFD.linInterp_3D(mesh, gradients)
    @test almostEqual(faceGrad[1,:], [ 66.628055, 75.37602, 0.0 ])
end;

@testset "JST Solver" begin
    nCells = 4
    boundaryConditions = [ zeroGradientBoundary, [], zeroGradientBoundary, [], emptyBoundary, [] ]
    mesh, cellPrimitives = initializeShockTube3DFVM(nCells, Pratio=0.1)
    mesh = Mesh(mesh.cells, mesh.cVols, mesh.cCenters, mesh.cellSizes, mesh.faces, mesh.fAVecs, mesh.fCenters, mesh.boundaryFaces)
    cellPrimitives2 = solve(mesh, "fakePath", cellPrimitives, boundaryConditions, initDt=0.051, endTime=0.05, createRestartFile=false, createVTKOutput=false)

    P = [ 1, 0.99855553, 0.087046989, 0.1 ]
    U = [ 0, 0.090015665, 0.720126353, 0 ]
    T = [ 0.00348432, 0.00347929, 0.0024263, 0.00278746 ]

    @test almostEqual(P, cellPrimitives2[:,1], 6)
    @test almostEqual(U, cellPrimitives2[:,3], 6)
    @test almostEqual(T, cellPrimitives2[:,2], 6)
end;

max_stars_repo_name: henrystoldt/fvCFD


text:
__precompile__()
module Drivers

using ForSyDeIO.Models
using ForSyDeIO.Models.Traits
using EzXML

export load_model, write_model

load_model_forxml(name::String) = load_model_forxml(open(name, "r"))

load_model_forxml(instream::IOStream) = load_model_forxml(readxml(instream))

function load_model_forxml(doc::EzXML.Document)::Models.ForSyDeModel
    new_model = Models.ForSyDeModel()
    root_node = root(doc)
    model_ns = [Pair("forsyde", namespace(root_node))]
    for vnode in findall("/forsyde:graphml/forsyde:graph//forsyde:node", root_node, model_ns)
        identifier = vnode["id"]
        traits = Set{Models.VertexTrait}()
        for tname in split(vnode["traits"], ";")
            if Traits.is_trait(tname)
                push!(traits, Traits.make_trait_vertex(tname))
            end
        end
        ports = Set{String}()
        for pnode in findall("forsyde:port", vnode, model_ns)
            push!(ports, pnode["name"])
        end
        # ports = Set(pnode["name"] for pnode in findall("forsyde:port", vnode, model_ns))
		# traits = filter(!isnothing, Set(Traits.make_trait(tname) for tname in split(vnode["traits"], ";")))
        props = Dict{String,Models.PropertyElement}()
        for propnode in findall("forsyde:data", vnode, model_ns)
            props[propnode["attr.name"]] = load_model_property_forxml(propnode)
        end
        vertex = Models.Vertex(
            identifier,
            ports,
		 	props,
		 	traits
        )
        push!(new_model, vertex)
    end
    for enode in findall("/forsyde:graphml/forsyde:graph//forsyde:edge", root_node, model_ns)
        traits = Set{Models.EdgeTrait}()
        source = Ref(getindex(new_model, enode["source"]))
        target = Ref(getindex(new_model, enode["target"]))
        sourceport = haskey(enode, "sourceport") ? enode["sourceport"] : nothing
        targetport = haskey(enode, "targetport") ? enode["targetport"] : nothing
        for tname in split(enode["traits"], ";")
            if Traits.is_trait(tname)
                push!(traits, Traits.make_trait_edge(tname))
            end
        end
        push!(new_model, Models.Edge(source, target, sourceport, targetport, traits))
    end
    return new_model
end # function

function load_model_property_forxml(node)::Models.PropertyElement
        if node["attr.type"] == "object"
        new_dict = Dict{String,Models.PropertyElement}()
        for e in eachelement(node)
            new_dict[e["attr.name"]] = load_model_property_forxml(e)
        end
    return Models.PropertyDict(new_dict)
    elseif node["attr.type"] == "array"
        new_array = Vector{Models.PropertyElement}()
        for e in eachelement(node)
            push!(new_array, load_model_property_forxml(e))
        end
        return Models.PropertyArray(new_array)
    elseif node["attr.type"] == "int" || node["attr.type"] == "integer"
        return parse(Int, node.content)
    elseif node["attr.type"] == "bool" || node["attr.type"] == "boolean"
        return parse(Bool, node.content)
    elseif node["attr.type"] == "float"
        return parse(Float32, node.content)
    elseif node["attr.type"] == "double"
        return parse(Float64, node.content)
    elseif node["attr.type"] == "str" || node["attr.type"] == "string"
        return node.content
    else
    end
end

function load_model(name::AbstractString)::Models.ForSyDeModel
    if endswith(name, ".forxml")
        return load_model_forxml(name)
    else
        error("File '$name' is not in a supported readable format.")
    end
end

function write_model_property_forxml(nodeprop::Models.PropertyArray)
    newnode = ElementNode("data")
    link!(newnode, AttributeNode("attr.type", "array"))
    for (i, v) in enumerate(nodeprop)
        childnode = write_model_property_forxml(v)
        link!(childnode, AttributeNode("attr.name", string(i)))
        link!(newnode, childnode)
    end
    return newnode
end

function write_model_property_forxml(nodeprop::Models.PropertyDict)
    newnode = ElementNode("data")
    link!(newnode, AttributeNode("attr.type", "object"))
    for (k, v) in nodeprop
        childnode = write_model_property_forxml(v)
        link!(childnode, AttributeNode("attr.name", k))
        link!(newnode, childnode)
    end
    return newnode
end

function write_model_property_forxml(nodeprop::Int)
    newnode = ElementNode("data")
    link!(newnode, AttributeNode("attr.type", "integer"))
    link!(newnode, TextNode(string(nodeprop)))
    return newnode
end

function write_model_property_forxml(nodeprop::Float32)
    newnode = ElementNode("data")
    link!(newnode, AttributeNode("attr.type", "float"))
    link!(newnode, TextNode(string(nodeprop)))
    return newnode
end

function write_model_property_forxml(nodeprop::Float64)
    newnode = ElementNode("data")
    link!(newnode, AttributeNode("attr.type", "double"))
    link!(newnode, TextNode(string(nodeprop)))
    return newnode
end

function write_model_property_forxml(nodeprop::Bool)
    newnode = ElementNode("data")
    link!(newnode, AttributeNode("attr.type", "bool"))
    link!(newnode, TextNode(string(nodeprop)))
    return newnode
end

function write_model_property_forxml(nodeprop::String)
    newnode = ElementNode("data")
    link!(newnode, AttributeNode("attr.type", "string"))
    link!(newnode, TextNode(string(nodeprop)))
    return newnode
end

function build_model_forxml(model::Models.ForSyDeModel)::EzXML.Node
    graphmlnode = ElementNode("graphml")
    link!(graphmlnode, AttributeNode("xmlns", "http://graphml.graphdrawing.org/xmlns"))
    graphnode = ElementNode("graph")
    link!(graphnode, AttributeNode("id", "model"))
    link!(graphnode, AttributeNode("edgedefault", "directed"))
    link!(graphmlnode, graphnode)
    for v in model
        vnode = ElementNode("node")
        link!(vnode, AttributeNode("id", v.id))
        for port in v.ports
            portnode = ElementNode("port")
            link!(portnode, AttributeNode("name", port))
            link!(vnode, portnode)
        end
        for (propname, propelem) in v.properties
            datanode = write_model_property_forxml(propelem)
            link!(datanode, AttributeNode("attr.name", propname))
            link!(vnode, datanode)
        end
        link!(vnode, AttributeNode("traits", join([string(t) for t in v.vertex_traits], ";")))
        link!(graphnode, vnode)
    end
    for e in model.edges
        enode = ElementNode("edge")
        link!(enode, AttributeNode("source", Models.src(e).id))
        link!(enode, AttributeNode("target", Models.dst(e).id))
            if !isnothing(e.source_port)
            link!(enode, AttributeNode("sourceport", e.source_port))
        end
            if !isnothing(e.target_port)
            link!(enode, AttributeNode("targetport", e.target_port))
        end
        link!(enode, AttributeNode("traits", join([string(t) for t in e.edge_traits], ";")))
        link!(graphnode, enode)
    end
    return graphmlnode
end

function write_model_forxml(model::Models.ForSyDeModel, outstream::IO)
    doc = XMLDocument()
    setroot!(doc, build_model_forxml(model))
    print(outstream, doc)
end

function write_model_forxml(model, name::AbstractString)
    open(name, "w") do f
        write_model_forxml(model, f)
    end
end

function write_model(model::Models.ForSyDeModel, name::AbstractString)
    if endswith(name, ".forxml")
        write_model_forxml(model, name)
    else
        error("File '$name' is not in a supported writable format.")
    end
end

end # module

max_stars_repo_name: forsyde/forsyde-io

max_issues_repo_name: forsyde/forsyde-io

max_forks_repo_name: forsyde/forsyde-io


text:
module GrigoriadisKhachiyan

using ..RandUtils

# Code from Dilys Thomas <dilys@cs.stanford.edu>

function myunifskew(n)
    A = zeros(n, n)
    for i=1:n
        for j=1:i-1
            temp = samerand()
            if (temp < 0.5)
                temp = samerand()
                A[i,j]= temp
                A[j,i]= -A[i,j]
            else
                temp = samerand()
                A[j,i]= temp
                A[i,j]= -A[j,i]
            end
        end
    end
    return A
end

function perf_gk(n, myeps)
    A = myunifskew(n)
    g = length(myeps)
    iteration = zeros(g)
    times = zeros(g)

    for KK = 1:g
        eps = myeps[KK]
        e, X, U = ones(n), zeros(n), zeros(n)
        t, stop, iter = 0, 0, 0
        p = e./n
        epse = eps .* e
        csum = zeros(n)

        while stop != 1
            t = t+1
            iter = t

            for i=1:n
                csum[i] = sum(p[1:i])
            end

            marker = samerand()

            k = 1
            for i = 2:n
                if csum[i-1] <= marker && marker <= csum[i]
                    k=i
                    break
                end
            end

            X[k] += 1

            for i=1:n
                U[i] += A[i,k]
            end

            s = sum(p[1:n] .* exp.((eps/2) .* A[1:n,k]))

            for i=1:n
                p[i] = (p[i]*exp((eps/2)*A[i,k])) / s
            end

            u = U ./ t

            true_count = 0
            for i=1:n
                if u[i] <= epse[i]
                    true_count += 1
                end
            end

            if true_count == n
                stop = 1
            end
        end

        times[KK] = 0
        iteration[KK] = iter
        x = X/t
        etx=sum(x)
        AX = A*X
        Ax = A*x
        error = abs.(AX) .- abs.(U)
        Axepse = 0

        for i=1:n
            if Ax[i]<=epse[i]
                Axepse = Axepse+1
            end
        end

        errorlmt = 0

        for i=1:n
            if error[i]<1e-8
                errorlmt = errorlmt+1
            end
        end
    end
end

end # module


text:
# Autogenerated wrapper script for LHAPDF_jll for powerpc64le-linux-gnu-cxx03
export libLHAPDF

JLLWrappers.@generate_wrapper_header("LHAPDF")
JLLWrappers.@declare_library_product(libLHAPDF, "libLHAPDF.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libLHAPDF,
        "lib/libLHAPDF.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# From algorithm 3.1

function bracket_minimum(f, x=0, s=1e-2, k=2.0)

    a, ya = x, f(x)
    b, yb = a + s, f(a + s)
    if yb > ya
        a, b = b, a
        ya, yb = yb, ya
        s = -s
    end
    while true
        c, yc = b + s, f(b + s)
        if yc > yb
            return (min(a,c), max(a,c))
        end
        a, ya, b, yb = b, yb, c, yc
        s *= k
    end
end

function test_f(x)
     return x^3 + 2 * x^2 + x + 3
end

@time a, b = bracket_minimum(test_f)


text:
include("gausslegendre.jl")
include("getdigit_n.jl")
include("myprintf.jl")
using Printf
using .GaussLegendre
using .GetDigit_n
using .MyPrintf

const DIGIT = 17
const MAXN = 1000

function atan_gl(n)
    return BigFloat(4) * GaussLegendre.gauleg(x::BigFloat -> BigFloat(1) / (BigFloat(1) + x * x), BigFloat(0), BigFloat(1), n)
end

function main()
    n = GetDigit_n.getn_linear(atan_gl, MAXN, DIGIT, true)
    precision = GetDigit_n.getdigit(atan_gl, n, DIGIT, true)
    precision = GetDigit_n.getdigit_2(atan_gl, n, DIGIT, precision)

    setprecision(precision)
    MyPrintf.myprintf("π（計算値） = ", DIGIT, atan_gl(n))
    MyPrintf.myprintf("π（厳密値） = ", DIGIT, BigFloat(π))
    @printf("n = %d, 計算%d桁で、%d桁求まりました\n", n, precision, DIGIT)
end

main()


text:
vis = Visualizer()

@test vis.path == MeshCat.Path(["meshcat"])
@test vis[:foo].path == MeshCat.Path(["meshcat", "foo"])
@test vis[:foo, :bar].path == MeshCat.Path(["meshcat", "foo", "bar"])
@test vis[:foo, "bar"].path == MeshCat.Path(["meshcat", "foo", "bar"])
@test vis[:foo, "/bar"].path == MeshCat.Path(["bar"])
@test vis[:foo, "/bar", :baz].path == MeshCat.Path(["bar", "baz"])
@test vis["foo/bar/baz"].path == MeshCat.Path(["meshcat", "foo", "bar", "baz"])
@test vis["/foo/bar/baz"].path == MeshCat.Path(["foo", "bar", "baz"])


text:
function contract!(tiJaB_d_temp3, tiJaB_i, Wmnij)
    rvir = axes(tiJaB_d_temp3, 4)
    nvir = last(rvir)
    rocc = axes(tiJaB_d_temp3, 1)
    @inbounds @fastmath for b in rvir, a in b:nvir, j in rocc, i in j:last(rocc)
        temp = zero(eltype(tiJaB_i))
        for n in rocc, m in rocc
            temp += tiJaB_i[m,n,a,b]*Wmnij[m,n,i,j]
        end
        tiJaB_d_temp3[i,j,a,b] = temp
        tiJaB_d_temp3[j,i,b,a] = temp
    end
end

function contracttest1!(tiJaB_d_temp3, tiJaB_i, Wmnij)
    rvir = axes(tiJaB_d_temp3, 4)
    nvir = last(rvir)
    rocc = axes(tiJaB_d_temp3, 1)
    for b in rvir, j in rocc
        @turbo for a in b:nvir, i in j:last(rocc)
            temp = zero(eltype(tiJaB_i))
            for n in rocc, m in rocc
                temp += tiJaB_i[m,n,a,b]*Wmnij[m,n,i,j]
            end
            tiJaB_d_temp3[i,j,a,b] = temp
            tiJaB_d_temp3[j,i,b,a] = temp
        end
    end
    tiJaB_d_temp3
end
function contracttest2!(tiJaB_d_temp3, tiJaB_i, Wmnij)
    rvir = axes(tiJaB_d_temp3, 4)
    nvir = last(rvir)
    rocc = axes(tiJaB_d_temp3, 1)
    for b in rvir, a in b:nvir, j in rocc
        @turbo for i in j:last(rocc)
            temp = zero(eltype(tiJaB_i))
            for n in rocc, m in rocc
                temp += tiJaB_i[m,n,a,b]*Wmnij[m,n,i,j]
            end
            tiJaB_d_temp3[i,j,a,b] = temp
            tiJaB_d_temp3[j,i,b,a] = temp
        end
    end
    tiJaB_d_temp3
end

@testset "Tensors" begin
    @show Float64, @__LINE__
    LA, LIM = 31, 23;
    A = rand(LIM, LIM, LA, LA);
    B = rand(LIM, LIM, LIM, LIM);

    C1 = fill(-999.999, LIM, LIM, LA, LA);
    C2 = fill(-999.999, LIM, LIM, LA, LA);
    C3 = fill(-999.999, LIM, LIM, LA, LA);
    # C1 = Array{Float64}(undef, LIM, LIM, LA, LA);
    # C2 = similar(C1); C3 = similar(C1);

    contract!(C1, A, B)
    @test C1 ≈ contracttest1!(C2, A, B)
    @test C1 ≈ contracttest2!(C3, A, B)

    Apermute = PermutedDimsArray(permutedims(A, (2,4,1,3)), (3,1,4,2));
    Bpermute = PermutedDimsArray(permutedims(B, (2,1,4,3)), (2,1,4,3));
    @test C1 ≈ contracttest1!(C2, Apermute, Bpermute)
    @test C1 ≈ contracttest2!(C3, Apermute, Bpermute)
end


text:
#!/usr/bin/env julia

using Glob

# include("common.jl")

# run a benchmark once, returning the measurement data
function run_benchmark(dir)
    output_file = "nvprof.csv"
    output_pattern = "nvprof.csv.*"

    # delete old profile output
    rm.(glob(output_pattern, dir))

    # run and measure
    out = Pipe()
    cmd = ```
        nvprof
        --profile-from-start off
        --concurrent-kernels off
        --profile-child-processes
        --unified-memory-profiling off
        --print-api-trace
        --print-gpu-trace
        --normalized-time-unit us
        --csv
        --log-file $output_file.%p
        ./profile --depwarn=no --check-bounds=no --math-mode=fast
    ```
    cmd_success = cd(dir) do
        success(pipeline(ignorestatus(cmd), stdout=out, stderr=out))
    end
    close(out.in)
    output_files = glob(output_pattern, dir)

    # read all data
    if !cmd_success
        println(readstring(out))
        error("benchmark did not succeed")
    else
        # when precompiling, an additional process will be spawned,
        # but the resulting CSV should not contain any data.
        output_data = []
        for output_file in output_files
            data = read_data(output_file)
            if data != nothing
                push!(output_data, data)
            end
        end

        if length(output_data) == 0
            error("no output files")
        elseif length(output_data) > 1
            error("too many output files")
        else
            return output_data[1]
        end
    end
end

function read_data(output_path)
    output = readlines(output_path; chomp=false)
    any(line->contains(line, "No kernels were profiled."), output) && return nothing

    # skip nvprof comments
    comments = findlast(line->occursin(r"^==\d+==", line), output)

    # NOTE: the first row often contains missing values (because it is the NVTX range),
    #       which confuses CSV.jl (https://github.com/JuliaData/CSV.jl/issues/156)
    CSV.read(output_path; header=comments+1, datarow=comments+3,
             types=Dict("Duration"              => Union{Missing,Float64},
                        "Grid X"                => Union{Missing,Int},
                        "Grid Y"                => Union{Missing,Int},
                        "Grid Z"                => Union{Missing,Int},
                        "Block X"               => Union{Missing,Int},
                        "Block Y"               => Union{Missing,Int},
                        "Block Z"               => Union{Missing,Int},
                        "Registers Per Thread"  => Union{Missing,Int},
                        "Static SMem"           => Union{Missing,Float64},
                        "Dynamic SMem"          => Union{Missing,Float64},
                        "Size"                  => Union{Missing,Float64},
                        "Throughput"            => Union{Missing,Float64},
                        "SrcMemType"            => Union{Missing,String},
                        "DstMemType"            => Union{Missing,String},
                        "Device"                => Union{Missing,String},
                        "Context"               => Union{Missing,String},
                        "Stream"                => Union{Missing,String},
                        "Correlation_ID"        => Union{Missing,Int}))
end

function process_data(raw_data, suite, benchmark)
    # extract kernel timings
    kernel_data = filter(entry -> begin
            !(startswith(entry[:Name], "cu") && ismissing(entry[:Device])) &&
            !occursin(r"^\[CUDA .+\]$", entry[:Name]) &&
            !occursin(r"^\[Range .+\] .+", entry[:Name])
        end, raw_data)

    # demangle kernel names
    for kernel in eachrow(kernel_data)
        jl_match = match(r"ptxcall_(.*)_[0-9]+", kernel[:Name])
        if jl_match != nothing
            kernel[:Name] = jl_match.captures[1]
            continue
        end

        cu_match = match(r"(.*)\(.*", kernel[:Name])
        if cu_match != nothing
            kernel[:Name] = cu_match.captures[1]
            continue
        end

        error("could not match kernel name $(kernel[:Name])")
    end

    # extract NVTX range timings
    if any(name->startswith(name, "[Range"), raw_data[:Name])
        range_start = filter(row -> startswith(row[:Name], "[Range start]"), raw_data)
        range_end   = filter(row -> startswith(row[:Name], "[Range end]"),   raw_data)
        @assert size(range_start) == size(range_end)

        ranges = range_start
        ranges[:Name] = map(name->'#' * match(r"\[Range start\] (.+) \(Domain: .+\)",
                                              name).captures[1],
                            ranges[:Name])
        ranges[:Duration] = range_end[:Start] - range_start[:Start]

        kernel_data = vcat(kernel_data, ranges)
    end

    # generate a nicer table
    rows = size(kernel_data, 1)
    data = DataFrame(suite = suite,
                     benchmark = benchmark,
                     target = kernel_data[:Name],
                     time = kernel_data[:Duration])

    # extract additional host-side compilation time measurements from compiles.dat:
    # -    main(ARGS)
    # +    open("compiles.dat", "w") do io
    # +        ccall(:jl_dump_compiles, Nothing, (Ptr{Cvoid},), io.handle)
    # +        main(ARGS)
    # +        ccall(:jl_dump_compiles, Nothing, (Ptr{Cvoid},), C_NULL)
    # +    end
    if suite == "julia_cuda"
        compiles = joinpath(root, suite, benchmark, "compiles.dat")
        # measure everything that's defined in Main (modules can be precompiled)
        # except for kernels (already contained in the CUDAnative timings)
        kernels = filter(target->!startswith(target, '#'), unique(data[:target]))
        measurements = 0
        open(compiles) do io
            for line in eachline(io)
                if occursin(r"\bMain\.", line) &&
                   !any(occursin(Regex("\\b$kernel\\b"), line) for kernel in kernels)
                    duration = parse(Int, first(split(line)))
                    push!(data, [suite, benchmark, "#julia", duration/1000])
                    measurements += 1
                end
            end
        end

        # if we don't have any measurement (eg. all functions interpreted),
        # add a dummy timing
        if measurements == 0
            push!(data, [suite, benchmark, "#julia", 0])
        end
    end

    # pull apart iterations of irregular kernels
    if haskey(irregular_kernels, benchmark)
        counters = Dict{String,Int64}()
        bad_kernels = irregular_kernels[benchmark]
        for i in 1:size(data, 1)
            target = data[:target][i]
            if target in bad_kernels
                j = get(counters, target, 0) + 1
                data[:target][i] *= "#$j"
                counters[target] = j
            end
        end
    end

    return data
end

# check if measurements are accurate enough
function is_accurate(data)
    grouped = summarize(data)

    # calculate relative error
    grouped[:ε] = map(t->t.err / abs(t.val), grouped[:time])

    return all(i->i>=MIN_KERNEL_ITERATIONS, grouped[:target_iterations]) &&
           all(val->val<MAX_KERNEL_ERROR, grouped[:ε])
end

function measure(host=gethostname())
    # find benchmarks common to all suites
    benchmarks = Dict()
    for suite in suites
        entries = readdir(joinpath(root, suite))
        benchmarks[suite] = filter(entry->(isdir(joinpath(root,suite,entry)) &&
                                           isfile(joinpath(root,suite,entry,"profile"))), entries)
    end
    common_benchmarks = intersect(values(benchmarks)...)

    # collect measurements
    measurements = DataFrame(suite=String[], benchmark=String[], target=String[],
                             time=Float64[], execution=Int64[])
    for suite in suites, benchmark in common_benchmarks
        @info "Processing $suite/$benchmark"
        dir = joinpath(root, suite, benchmark)
        cache_path = joinpath(dir, "profile_$host.csv")

        data = if isfile(cache_path)
            try
                CSV.read(cache_path)
            catch
                @warn "Corrupt cache file at $cache_path"
                nothing
            end
        else
            nothing
        end

        if data == nothing
            iter = 1
            t0 = time()
            data = nothing
            while true
                new_data = process_data(run_benchmark(dir), suite, benchmark)
                new_data[:execution] = iter
                iter += 1

                if data == nothing
                    data = new_data
                else
                    data = vcat(data, new_data)
                end

                is_accurate(data)                    && break
                iter >= MAX_BENCHMARK_RUNS           && break
                (time()-t0) >= MAX_BENCHMARK_SECONDS && break
            end

            CSV.write(cache_path, data)
        end

        measurements = vcat(measurements, data)
    end

    CSV.write("measurements_$host.dat", measurements)

    return
end


text:
function test_method(func)
    frame = zeros(101, 101)
    frame[77, 77] = 1.0
    # central fwhm is all NaN
    @test func(frame, (51, 51), 3) === NaN
    # centered on point is infinite, because bkg stddev is 0
    @test func(frame, (77, 77), 3) ≈ Inf
end

@testset "snr" begin
    test_method(snr)
end

@testset "significance" begin
    test_method(significance)

    X = randn(rng, 101, 101) .+ 10
    # expect small sample to matter ~3fwhm out
    @test abs(significance(X, (58, 51), 3)) < abs(snr(X, (58, 51), 3))
end

@testset "noise" begin
    data = ones(101, 101)

    @test noise(data, (51, 51), 10) === NaN
end

@testset "detectionmap - $func" for func in (snr, significance, noise)
    data = ones(101, 101)
    _map = detectionmap(func, data, 10)
    # center should be fill value
    @test _map[51, 51] ≈ 0
    # there should be some values for the rest
    @test maximum(_map) > 0
    @test abs(minimum(_map)) ≥ 0
end

@testset "stim map" begin
    S = randn(rng, 100, 512, 512)
    angles = sort!(90rand(rng, 100))

    sm = stimmap(S, angles)
    @test size(sm) == (size(S, 2), size(S, 3))
    @test all(isfinite, sm)
    @test eltype(sm) == eltype(S)

    smt = stim_threshold(sm, S, angles)
    smt2 = stim_threshold(S, angles)
    @test smt ≈ smt2
    @test smt ≈ 0.5 rtol=1e-2

end

@testset "slim map" begin
    cube, angles = BetaPictoris[:cube, :pa]

    resid_cubes = subtract.(PCA.(5:10), Ref(cube))
    stim_maps = stimmap.(resid_cubes, Ref(angles))
    # 1
    stim_av, mask = slimmap(stim_maps; N=20)
    slim1 = stim_av .* mask
    #2
    stim_av2, mask2 = slimmap(resid_cubes, angles; N=20)
    slim2 = stim_av2 .* mask2

    @test slim1 ≈ slim2

    # make sure beta pic b was found
    @test slim1[62, 62] > 0.5
end


text:
# limits number of bytes written to an io stream (originally from https://github.com/JuliaDebug/Debugger.jl/blob/master/src/limitio.jl)
# useful to detect messages that would go over limit when converted to JSON.
mutable struct LimitIO{IO_t <: IO} <: IO
    io::IO_t
    maxbytes::Int
    n::Int # max bytes to write
end
LimitIO(io::IO, maxbytes) = LimitIO(io, maxbytes, 0)

struct LimitIOException <: Exception end

function Base.write(io::LimitIO, v::UInt8)
    io.n > io.maxbytes && throw(LimitIOException())
    io.n += write(io.io, v)
end


text:
################################################################################
### Lattice Iterator
################################################################################

struct SuperfluidStiffness{IT, N, T} <: AbstractLatticeIteratorTemplate
    template::IT
    dq::Vector{Float64}
end
function SuperfluidStiffness(template::DeferredLatticeIteratorTemplate, dq)
    SuperfluidStiffness(template, dq)
end

struct _SuperfluidStiffness{LI <: DeferredLatticeIterator} <: LatticeIterationWrapper{LI}
    iter::LI
    weights::Matrix{ComplexF64}
end
function (x::SuperfluidStiffness)(cache::LatticeIteratorCache, lattice::AbstractLattice)
    # We need to calculate 0.25 (Λₓₓ(dq) - Λₓₓ(0)) where
    # Λₓₓ(q) = CCS[i, j, k] * dot(eₓ, dir[j]) * dot(eₓ, dir[k]) *
    #           cis(- dot(dq, dir[i], + 0.5 (dir[j] - dir[k])))
    # we collect the 0.25, the cis from q=dq and the cis from q=0 (1) into weights
    # Idk why the dot(eₓ, dir[j/k]) is there but https://arxiv.org/pdf/1912.08848.pdf
    # has them so I'll have them too (in J_x as prefactor 1 or 2)
    # The 0.5 (dir[j] - dir[k]) comes from the same paper, eq 17-19
    all_dirs = directions(lattice)
    hop_idxs = x.template.directions
    dr = normalize(x.dq)
    fs = map(idx -> dot(dr, all_dirs[idx]), hop_idxs)
    weights = Matrix{ComplexF64}(undef, length(all_dirs), length(hop_idxs), length(hop_idxs))
    for (i, dr12) in enumerate(all_dirs)
        for (j, dj) in enumerate(hop_idxs), (k, dk) in enumerate(hop_idxs)
            weights[i,j,k] = 0.25 * fs[j] * fs[k] *
                (cis(-dot(x.dq, dr12 + 0.5 * (dirs[dj] - dirs[dk]))) - 1)
        end
    end

    _SuperfluidStiffness(x.template(cache, lattice), weights)
end

Base.iterate(s::_SuperfluidStiffness) = iterate(s.iter)
Base.iterate(s::_SuperfluidStiffness, state) = iterate(s.iter, state)
Base.length(s::_SuperfluidStiffness) = length(s.iter)
Base.eltype(s::_SuperfluidStiffness) = eltype(s.iter)

################################################################################
### Adjustments to measurement call stack
################################################################################

# binner contains Complex numbers
_binner_zero_element(mc, model, ::SuperfluidStiffness, eltype) = zero(ComplexF64)

# compute 0.25 (Λₓₓ(dq) - Λₓₓ(0)) = sum_{i,j,k} weight[i,j,k] * CCS[i,j,k]
function commit!(s::_SuperfluidStiffness, m)
    @assert size(s.weight) == size(m.temp)
    final = zero(ComplexF64)
    for i in eachindex(m.temp)
        final += s.weight[i] * m.temp[i]
    end
    push!(m.observable, final)
end

################################################################################
### Measurement constructor
################################################################################

"""
    superfluid_stiffness(dqmc, model[; Ls, kwargs...])

Returns a measurement that calculates the superfluid stiffness as
`0.25 (Λₓₓ((dq, 0)) - Λₓₓ(0))`.

Note that this might be slightly different from `0.25 (-K_x - Λₓₓ(0))` and that
some sources use `0.25 (Λₓₓ((dq, 0)) - Λₓₓ((0, dq)))` instead. (I.e. a
longitudinal Λₓₓ minus a transversal Λₓₓ.)
"""
function superfluid_stiffness(
        dqmc::DQMC, model::Model; Ls = size(lattice(model)), kwargs...
    )
    # TODO
    # I'm still not sure how exactly shift_dir/long should be picked.

    # 0.25 [Λₓₓ(q = (Δq, 0) - Λₓₓ(0))]
    shift_dir = Float64[1, 0]

    # find all hopping directions (skipping on-site)
    dir2srctrg = mc[Dir2SrcTrg()]
    T = hopping_matrix(mc, model)
    valid_directions = Int64[]

    for i in 2:length(dir2srctrg)
        for (src, trg) in dir2srctrg[i]
            if T[trg, src] != 0
                push!(valid_directions, i)
                break
            end
        end
    end

    # reduce to directions that have a positive component in shift_dir
    dirs = directions(lattice(dqmc))
    filter!(idx -> dot(dirs[idx], shift_dir) > 0, valid_directions)

    # create lattice iterator
    long = 2pi / Ls[1] * shift_dir
    li = SuperfluidStiffness(EachLocalQuadByDistance(valid_directions), long)

    current_current_susceptibility(dqmc, model, lattice_iterator = li)
end


text:
import DataFrames: DataFrameRow
import DataArrays: NA
import Compat: @compat, parse

function test_dbi()
    @test Postgres <: DBI.DatabaseSystem

    conn = connect(Postgres, "localhost", "postgres", "", "julia_test")

    @test isa(conn, DBI.DatabaseHandle)
    @test isdefined(conn, :status)

    stmt = prepare(conn, "SELECT 1::bigint, 2.0::double precision, 'foo'::character varying, " *
                         "'foo'::character(10), NULL;")
    result = execute(stmt)
    testdberror(result, PostgreSQL.PGRES_TUPLES_OK)

    iterresults = Vector{Any}[]
    for row in result
        @test row[1] === @compat Int64(1)
        @test_approx_eq row[2] 2.0
        @test typeof(row[2]) == Float64
        @test row[3] == "foo"
        @test typeof(row[3]) <: AbstractString
        @test row[4] == "foo       "
        @test typeof(row[4]) <: AbstractString
        @test row[5] === Union{}
        push!(iterresults, row)
    end

    allresults = fetchall(result)

    @test iterresults == allresults

    dfresults = fetchdf(result)

    dfrow = Any[x[2] for x in DataFrameRow(dfresults, 1)]
    dfrow[5] = Union{}

    @test dfrow == allresults[1]

    finish(stmt)

    create_str = """CREATE TEMPORARY TABLE testdbi (
            id serial PRIMARY KEY,
            combo double precision,
            quant double precision,
            name varchar,
            color text,
            bin bytea,
            is_planet bool,
            num_int numeric(80,0),
            num_float numeric(80,10)
        );"""

    run(conn, create_str)

    data = Vector[
        Any[1, 4, "Spam spam eggs and spam", "red", (UInt8)[0x01, 0x02, 0x03, 0x04], Union{}, BigInt(123), parse(BigFloat, "123.4567")],
        Any[5, 8, "Michael Spam Palin", "blue", (UInt8)[], true, -3, parse(BigFloat, "-3.141592653")],
        Any[3, 16, Union{}, Union{}, Union{}, false, Union{}, Union{}],
        Any[NA, 32, "Foo", "green", (UInt8)[0xfe, 0xdc, 0xba, 0x98, 0x76], true, 9876, parse(BigFloat, "9876.54321")]
    ]

    insert_str = "INSERT INTO testdbi (combo, quant, name, color, bin, is_planet, num_int, num_float) " *
                 "VALUES(\$1, \$2, \$3, \$4, \$5, \$6, \$7, \$8);"

    stmt = prepare(conn, insert_str)
    for row in data
        execute(stmt, row)
        testdberror(stmt, PostgreSQL.PGRES_COMMAND_OK)
    end
    finish(stmt)

    stmt = prepare(conn, "SELECT combo, quant, name, color, bin, is_planet, num_int, num_float FROM testdbi ORDER BY id;")
    result = execute(stmt)
    testdberror(stmt, PostgreSQL.PGRES_TUPLES_OK)
    rows = fetchall(result)
    @test rows[1] == data[1]
    @test rows[2] == data[2]
    @test rows[3] == data[3]
    @test rows[4][1] == Union{}
    @test rows[4][2] == data[4][2]
    @test rows[4][3] == data[4][3]
    @test rows[4][4] == data[4][4]
    @test rows[4][5] == data[4][5]
    @test rows[4][6] == data[4][6]
    @test rows[4][7] == data[4][7]
    @test rows[4][8] == data[4][8]

    finish(stmt)

    disconnect(conn)

    conn = connect(Postgres, "localhost", "postgres", "", "julia_test")
    run(conn, create_str)
    stmt = prepare(conn, insert_str)
    executemany(stmt, data)
    testdberror(stmt, PostgreSQL.PGRES_COMMAND_OK)
    finish(stmt)

    stmt = prepare(conn, "SELECT combo, quant, name, color, bin, is_planet, num_int, num_float FROM testdbi ORDER BY id;")
    result = execute(stmt)
    testdberror(stmt, PostgreSQL.PGRES_TUPLES_OK)
    rows = fetchall(result)
    @test rows[1] == data[1]
    @test rows[2] == data[2]
    @test rows[3] == data[3]
    @test rows[4][1] == Union{}
    @test rows[4][2] == data[4][2]
    @test rows[4][3] == data[4][3]
    @test rows[4][4] == data[4][4]
    @test rows[4][5] == data[4][5]
    @test rows[4][6] == data[4][6]
    @test rows[4][7] == data[4][7]
    @test rows[4][8] == data[4][8]

    finish(stmt)

    @test escapeliteral(conn, 3) == 3
    @test escapeliteral(conn, 3.3) == 3.3
    @test escapeliteral(conn, "foo") == "'foo'"
    @test escapeliteral(conn, "fo\u2202") == "'fo\u2202'"
    @test escapeliteral(conn, SubString("myfood", 3, 5)) == "'foo'"

    disconnect(conn)

    # Test copy
    conn = connect(Postgres, "localhost", "postgres", "", "julia_test")
    create_str = """CREATE TEMPORARY TABLE testcopycsv(
      intdata INTEGER,
      otherint INTEGER,
      textdata TEXT,
      color VARCHAR,
      floatdata DOUBLE PRECISION
    );"""

    run(conn, create_str)
    PostgreSQL.copy_from(conn, "testcopycsv", "test_data.csv", "csv")

    stmt = prepare(conn, "SELECT * FROM testcopycsv WHERE intdata > 1;")
    rs = execute(stmt)
    rows = fetchall(rs)
    @test rows[1][1] == 2
    @test rows[1][4] == "green"
    @test rows[2][3] == "i'm sick of eggs now"
    testdberror(stmt, PostgreSQL.PGRES_TUPLES_OK)
    finish(stmt)

    @test_throws ErrorException PostgreSQL.copy_from(conn, "testcopycsv", "test_malformed_data.csv", "binary")

    disconnect(conn)

    # Test arrays
    conn = connect(Postgres, "localhost", "postgres", "", "julia_test")
    create_str = """CREATE TEMPORARY TABLE testarrays(
      intdata INTEGER,
      floatdata DOUBLE PRECISION,
      textdata TEXT,
      varchardata VARCHAR
    );"""

    run(conn, create_str)

    data = Vector[
        Any[1, 1.0, "1", "one"],
        Any[2, 2.0, "2", "two"],
        Any[3, 3.0, "3", "three"],
    ]

    insert_str = "INSERT INTO testarrays(intdata, floatdata, textdata, varchardata) " *
                 "VALUES(\$1, \$2, \$3, \$4);"

    stmt = prepare(conn, insert_str)
    for row in data
        execute(stmt, row)
        testdberror(stmt, PostgreSQL.PGRES_COMMAND_OK)
    end
    finish(stmt)

    stmt = prepare(conn, "SELECT * FROM testarrays WHERE textdata = ANY(\$1);")
    rs = execute(stmt, Any[["1", "2"]])
    for (i, row) in enumerate(rs)
      @test row[1] == data[i][1]
      @test row[2] == data[i][2]
      @test row[3] == data[i][3]
      @test row[4] == data[i][4]
    end
    testdberror(stmt, PostgreSQL.PGRES_TUPLES_OK)
    finish(stmt)

    stmt = prepare(conn, "SELECT * FROM testarrays WHERE textdata = ANY(\$1::varchar[]);")
    rs = execute(stmt, Any[["one", "two"]])
    for (i, row) in enumerate(rs)
      @test row[1] == data[i][1]
      @test row[2] == data[i][2]
      @test row[3] == data[i][3]
      @test row[4] == data[i][4]
    end
    testdberror(stmt, PostgreSQL.PGRES_TUPLES_OK)
    finish(stmt)

    stmt = prepare(conn, "SELECT * FROM testarrays WHERE intdata = ANY(\$1);")
    rs = execute(stmt, Any[[1, 2]])
    for (i, row) in enumerate(rs)
      @test row[1] == data[i][1]
      @test row[2] == data[i][2]
      @test row[3] == data[i][3]
      @test row[4] == data[i][4]
    end
    testdberror(stmt, PostgreSQL.PGRES_TUPLES_OK)
    finish(stmt)

    stmt = prepare(conn, "SELECT * FROM testarrays WHERE floatdata = ANY(\$1);")
    rs = execute(stmt, Any[[1.0, 2.0]])
    for (i, row) in enumerate(rs)
      @test row[1] == data[i][1]
      @test row[2] == data[i][2]
      @test row[3] == data[i][3]
      @test row[4] == data[i][4]
    end
    testdberror(stmt, PostgreSQL.PGRES_TUPLES_OK)
    finish(stmt)

    stmt = prepare(conn, """
    SELECT
    '{1,2}'::integer[],
    '{1.0,2.0}'::double precision[],
    '{"aaa", "bbb"}'::text[],
    '{"aaa", "bbb"}'::varchar[];
    """)
    rs = execute(stmt)
    for row in rs
      @test row[1] == [1,2]
      @test row[2] == [1.0,2.0]
      @test row[3] == ["aaa","bbb"]
      @test row[4] == ["aaa","bbb"]
    end
    testdberror(stmt, PostgreSQL.PGRES_TUPLES_OK)
    finish(stmt)

    disconnect(conn)
end

test_dbi()

max_stars_repo_name: beffaCo/updatePostgresql

max_forks_repo_name: beffaCo/updatePostgresql


text:
@testset "Diagonal" begin
    let rng = MersenneTwister(123456), N = 10
        λ_2 = x->diagm(2 => x)
        λ_m3 = x->diagm(-3 => x)
        λ_0 = x->diagm(0 => x)
        λ_false = x->diagm(false => x)
        λ_true = x->diagm(true => x)
        for _ in 1:10
            x, vx = randn.(Ref(rng), [N, N])
            @test check_errs(x->diagm(0 => x), diagm(0 => randn(rng, N)), x, vx)

            @test check_errs(λ_2, λ_2(randn(rng, N)), x, vx)
            @test check_errs(λ_m3, λ_m3(randn(rng, N)), x, vx)
            @test check_errs(λ_0, λ_0(randn(rng, N)), x, vx)
            @test check_errs(λ_false, λ_false(randn(rng, N)), x, vx)
            @test check_errs(λ_true, λ_true(randn(rng, N)), x, vx)
        end
    end
    let rng = MersenneTwister(123456), N = 10
        for _ in 1:10
            A = randn(rng, N, N)
            VA = randn(rng, N, N)
            @test check_errs(x -> diag(x), randn(rng, N), A, VA)

            # Check various diagonals.
            for k = -3:3
                @test check_errs(x -> diag(x, k), randn(rng, N - abs(k)), A, VA)
            end
        end
    end
    let rng = MersenneTwister(123456), N = 10
        for _ in 1:10
            A = randn(rng, N)
            VA = randn(rng, N)
            @test check_errs(Diagonal, Diagonal(randn(rng, N)), A, VA)
        end
    end
    let rng = MersenneTwister(123456), N = 10
        for _ in 1:10
            A = randn(rng, N, N)
            VA = randn(rng, N, N)
            @test check_errs(Diagonal, Diagonal(randn(rng, N)), A, VA)
        end
    end
    let rng = MersenneTwister(123456), N = 10
        for _ in 1:10
            A = Diagonal(randn(rng, N))
            VA = Diagonal(randn(rng, N))
            @test check_errs(det, 10.0, A, VA)
        end
    end
    let rng = MersenneTwister(123456), N = 10
        for _ in 1:10
            A = Diagonal(exp.(randn(rng, N)))
            VA = Diagonal(randn(rng, N))
            @test check_errs(logdet, 10.0, A, VA)
        end
    end

    # test that both diagonal implementations are correct and the same
    let
        f1(x) = sum(Diagonal(x))
        f2(x) = sum(diagm(0 => x))

        @test ∇(f1)(ones(4))[1] == ones(4)
        @test ∇(f2)(ones(4))[1] == ones(4)
    end
end


text:
# Julia wrapper for header: include/blockmat.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# TODO detect size and use bitstype because if the DLL changes and gets
#      compiled with NOSHORTS we are screwed with the following code...
@static if Sys.iswindows()
    const csdpshort = Cushort
  else
    const csdpshort = CSDP_INT
end

# begin enum blockcat
const blockcat = Cuint
const DIAG = (blockcat)(0)
const MATRIX = (blockcat)(1)
const PACKEDMATRIX = (blockcat)(2)
# end enum blockcat

struct blockdatarec
    _blockdatarec::Ptr{Cdouble}
end

struct blockrec <: AbstractMatrix{Cdouble}
    data::blockdatarec
    blockcategory::blockcat
    blocksize::csdpshort
end

mutable struct blockmatrix <: AbstractBlockMatrix{Cdouble}
    nblocks::CSDP_INT
    blocks::Ptr{blockrec}
end
blockmatrix() = blockmatrix(CSDP_INT(0), C_NULL)

mutable struct sparseblock
    next::Ptr{sparseblock}
    nextbyblock::Ptr{sparseblock}
    entries::Ptr{Cdouble}
    iindices::Ptr{csdpshort}
    jindices::Ptr{csdpshort}
    numentries::CSDP_INT
    blocknum::csdpshort
    blocksize::csdpshort
    constraintnum::csdpshort
    issparse::csdpshort
end
function sparseblock(next, jblock, blocknum, constr)
    # See easysdp.c
    blocksize = jblock.n
    numentries = length(jblock.i)
    issparse = numentries <= blocksize / 4 || numentries <= 15 # FIXME also if category (which is in C...) is DIAG
end

# If I add mutable here I get : ReadOnlyMemoryError() in initsoln (I know it is counter-intuitive)
struct constraintmatrix
    blocks::Ptr{sparseblock}
end

# Skipping MacroDefinition: ijtok ( iiii , jjjj , lda ) ( ( jjjj - 1 ) * lda + iiii - 1 )
# Skipping MacroDefinition: ijtokp ( iii , jjj , lda ) ( ( iii + jjj * ( jjj - 1 ) / 2 ) - 1 )
# Skipping MacroDefinition: ktoi ( k , lda ) ( ( k % lda ) + 1 )
# Skipping MacroDefinition: ktoj ( k , lda ) ( ( k / lda ) + 1 )

mutable struct paramstruc
    axtol::Cdouble
    atytol::Cdouble
    objtol::Cdouble
    pinftol::Cdouble
    dinftol::Cdouble
    maxiter::CSDP_INT
    minstepfrac::Cdouble
    maxstepfrac::Cdouble
    minstepp::Cdouble
    minstepd::Cdouble
    usexzgap::CSDP_INT
    tweakgap::CSDP_INT
    affine::CSDP_INT
    perturbobj::Cdouble
    fastmode::CSDP_INT
end

function paramstruc(options::Dict)
    paramstruc(get(options, :axtol      , 1.0e-8),
               get(options, :atytol     , 1.0e-8),
               get(options, :objtol     , 1.0e-8),
               get(options, :pinftol    , 1.0e8),
               get(options, :dinftol    , 1.0e8),
               get(options, :maxiter    , 100),
               get(options, :minstepfrac, 0.90),
               get(options, :maxstepfrac, 0.97),
               get(options, :minstepp   , 1.0e-8),
               get(options, :minstepd   , 1.0e-8),
               get(options, :usexzgap   , 1),
               get(options, :tweakgap   , 0),
               get(options, :affine     , 0),
               get(options, :perturbobj , 1),
               get(options, :fastmode   , 0))
end


text:
using Test

include("example.jl")
# add here the file with your unit tests

include("beesalgorithm.jl")

max_stars_repo_name: kirstvh/STMOZOO

max_issues_repo_name: kirstvh/STMOZOO

max_forks_repo_name: kirstvh/STMOZOO


text:
using DataFrames
using Query
using PyPlot
using LaTeXStrings

include("schema.jl")

# number of metainference runs to a solid or dashed linestyle
const linestyles = Dict([1 => "-", 100 => "--"])

# proposal name to color
const colors = Dict([PRIOR_PROPOSAL_NAME => "darkcyan", OPTIMAL_PROPOSAL_NAME => "red"])

# short names for the legend for the different types of proposals
proposal_short_names = Dict([PRIOR_PROPOSAL_NAME => "prior proposal", OPTIMAL_PROPOSAL_NAME => "optimal proposal"])

# where to read data from
const DATA_DIR = "data"

# where to write plots to
const PLOT_DIR = "plots"
if !Base.Filesystem.isdir(PLOT_DIR)
    Base.Filesystem.mkdir(PLOT_DIR)
end

const PLOT_FIGURE_SIZE = (3, 2.5)
const LEGEND_FIGURE_SIZE = (6.8, 1.4)

function plot_aide_estimates(data::DataFrame, gold_standard_name::String, proposal_name::String,
                             color::String)

    # for each meta-inference value, plot kl estimate vs number of particles
    num_metainference_list = sort(unique(data[COL_NUM_METAINFERENCE]))
    for num_metainference in num_metainference_list
        result = @from i in data begin
            @where ((i.gold_standard_name == gold_standard_name) &&
                    (i.proposal_name == proposal_name) &&
                    (i.num_metainference == num_metainference))
            @select {i.num_particles, i.aide_estimate}
            @collect DataFrame
        end
        plt[:plot](result[COL_NUM_PARTICLES], result[COL_AIDE_ESTIMATE],
                   linestyle=linestyles[num_metainference], color=colors[proposal_name])
    end
end

function plot_markers(data::DataFrame, gold_standard_name::String)
    for (num_particles, proposal_name) in [
        (1, PRIOR_PROPOSAL_NAME),
        (10, PRIOR_PROPOSAL_NAME),
        (100, OPTIMAL_PROPOSAL_NAME)]
        num_metainference = 1
        result = @from i in data begin
            @where ((i.gold_standard_name == gold_standard_name) &&
                    (i.proposal_name == proposal_name) &&
                    (i.num_metainference == num_metainference) &&
                    (i.num_particles == num_particles))
            @select {i.aide_estimate}
            @collect DataFrame
        end
        aide_estimate = result[COL_AIDE_ESTIMATE]
        @assert length(aide_estimate) == 1
        plt[:scatter]([num_particles], aide_estimate, color="white", s=75,
                      edgecolor="black", zorder=100)
    end

    # NOTE: labels are placed at hardcoded locations
    const fontsize = 26
    ax = plt[:gca]()
    plt[:text](0.05, 0.85, "A", fontsize=fontsize, transform=ax[:transAxes])
    plt[:text](0.5, 0.30, "B", fontsize=fontsize, transform=ax[:transAxes])
    plt[:text](0.87, 0.10, "C", fontsize=fontsize, transform=ax[:transAxes])
end

function draw_zero_line(ax, num_particles_list::Vector{Int})
    xlim = [minimum(num_particles_list), maximum(num_particles_list)]
    plt[:plot](xlim, [0, 0], color="black")
    ax[:set_xlim](xlim)
end

function format_axes(num_particles_list::Vector{Int})
    ax = plt[:gca]()
    draw_zero_line(ax, num_particles_list)
    ax[:set_xscale]("log")
    ax[:set_xlabel]("Number of particles")
    ax[:set_ylabel]("AIDE estimate (nats)")
    plt[:tight_layout](pad=0)
end

function plot_aide_estimate_exact_gold_standard(data::DataFrame, fname::String)
    num_particles_list = sort(unique(convert(Vector{Int}, data[COL_NUM_PARTICLES])))
    plt[:figure](figsize=PLOT_FIGURE_SIZE)
    plot_aide_estimates(data, EXACT_GOLD_STANDARD, PRIOR_PROPOSAL_NAME, colors[PRIOR_PROPOSAL_NAME])
    plot_aide_estimates(data, EXACT_GOLD_STANDARD, OPTIMAL_PROPOSAL_NAME, colors[OPTIMAL_PROPOSAL_NAME])
    plot_markers(data, EXACT_GOLD_STANDARD)
    format_axes(num_particles_list)
    plt[:savefig](fname)
end

function plot_aide_estimate_approximate_gold_standard(data::DataFrame, fname::String)
    num_particles_list = sort(unique(convert(Vector{Int}, data[COL_NUM_PARTICLES])))
    plt[:figure](figsize=PLOT_FIGURE_SIZE)
    plot_aide_estimates(data, APPROXIMATE_GOLD_STANDARD, PRIOR_PROPOSAL_NAME, colors[PRIOR_PROPOSAL_NAME])
    plot_aide_estimates(data, APPROXIMATE_GOLD_STANDARD, OPTIMAL_PROPOSAL_NAME, colors[OPTIMAL_PROPOSAL_NAME])
    plot_markers(data, APPROXIMATE_GOLD_STANDARD)
    format_axes(num_particles_list)
    plt[:savefig](fname)
end

function plot_legend(data::DataFrame, fname::String)
    plt[:figure](figsize=LEGEND_FIGURE_SIZE)
    for proposal_name in [PRIOR_PROPOSAL_NAME, OPTIMAL_PROPOSAL_NAME]
        for num_metainference in sort(unique(data[COL_NUM_METAINFERENCE]))
            if num_metainference == 1
                metainference_label = "$num_metainference meta-inference run (\$M_t = $num_metainference\$)"
            else
                metainference_label = "$num_metainference meta-inference runs (\$M_t = $num_metainference\$)"
            end
            label=latexstring("SMC, $(proposal_short_names[proposal_name]), $metainference_label")
            plt[:plot]([], [], linestyle=linestyles[num_metainference], color=colors[proposal_name], linewidth=2, label=label)
        end
    end
    plt[:axis]("off")
    plt[:legend](fontsize=14, ncol=1, frameon=false, borderpad=0, borderaxespad=0, loc="center")
    plt[:tight_layout](pad=0)
    plt[:savefig](fname)
end

# do plotting
data = readtable("$DATA_DIR/aide_estimates.csv")
plot_aide_estimate_exact_gold_standard(data, "$PLOT_DIR/divergences_exact_gold_standard.pdf")
plot_aide_estimate_approximate_gold_standard(data, "$PLOT_DIR/divergences_approximate_gold_standard.pdf")
plot_legend(data, "$PLOT_DIR/divergences_legend.pdf")
println("done!")


text:
# Fonction qu'on veut tester
# include("../NLPAtXmod.jl")
# importall NLPAtXmod

# On vérifie que le constructeur pour problème sans contrainte fonctionne
uncons_nlp_at_x = NLPAtX(zeros(10))

@test (false in (uncons_nlp_at_x.x .== 0.0))  == false #assez bizarre comme test...
@test (false in (isnan.(uncons_nlp_at_x.fx))) == false
@test (false in (isnan.(uncons_nlp_at_x.gx))) == false
@test (false in (isnan.(uncons_nlp_at_x.Hx))) == false
@test (false in (isnan.(uncons_nlp_at_x.mu))) == false
@test (false in (isnan.(uncons_nlp_at_x.cx))) == false
@test (false in (isnan.(uncons_nlp_at_x.Jx))) == false

@test (false in (isnan.(uncons_nlp_at_x.lambda))) == false
@test (false in (isnan.(uncons_nlp_at_x.start_time))) == false
@test (!(uncons_nlp_at_x.evals == nothing))

# On vérifie que le constucteur pour problème avec contrainte fonctionne
cons_nlp_at_x = NLPAtX(zeros(10), zeros(10))

@test (false in (uncons_nlp_at_x.x .== 0.0)) == false #assez bizarre comme test...
@test (false in (isnan.(uncons_nlp_at_x.fx))) == false
@test (false in (isnan.(uncons_nlp_at_x.gx))) == false
@test (false in (isnan.(uncons_nlp_at_x.Hx))) == false
@test (false in (isnan.(uncons_nlp_at_x.mu))) == false
@test (false in (isnan.(uncons_nlp_at_x.cx))) == false
@test (false in (isnan.(uncons_nlp_at_x.Jx))) == false
@test (false in (uncons_nlp_at_x.lambda .== 0.0)) == false
@test (false in (isnan.(uncons_nlp_at_x.start_time))) == false

# On vérifie que la fonction update! fonctionne
update!(uncons_nlp_at_x, x = ones(10), fx = 1.0, gx = ones(10))
update!(uncons_nlp_at_x, Hx = ones(10,10), mu = 1.0, cx = 1.0, Jx = ones(10,10))
update!(uncons_nlp_at_x, lambda = ones(10), tmps = 1.0)

@test (false in (uncons_nlp_at_x.x .== 1.0)) == false #assez bizarre comme test...
@test uncons_nlp_at_x.fx == 1.0
@test (false in (uncons_nlp_at_x.gx .== 1.0)) == false
@test (false in (uncons_nlp_at_x.Hx .== 1.0)) == false
@test uncons_nlp_at_x.mu == 1.0
@test uncons_nlp_at_x.cx == 1.0
@test (false in (uncons_nlp_at_x.Jx .== 1.0)) == false
@test (false in (uncons_nlp_at_x.lambda .== 1.0)) == false
@test uncons_nlp_at_x.start_time == 1.0

nlp_64 = NLPAtX(ones(10))
nlp_64.x = ones(10)
nlp_64.fx = 1.0
nlp_64.gx = ones(10)

# nlp_32 = convert_nlp(Float32, nlp_64)
# @test typeof(nlp_32.x[1]) == Float32
# @test typeof(nlp_32.fx[1]) == Float32
# @test typeof(nlp_32.gx[1]) == Float32
# @test isnan(nlp_32.mu[1])
# @test isnan(nlp_32.start_time)
#
# @test typeof(nlp_64.x[1]) == Float64
# @test typeof(nlp_64.fx[1]) == Float64
# @test typeof(nlp_64.gx[1]) == Float64
# @test isnan(nlp_64.mu[1])
# @test isnan(nlp_64.start_time)


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#
"""
Conic models and utilities.
"""
module Models

using DocStringExtensions
import Hypatia.Cones

"""
$(TYPEDEF)

Conic model corresponding to Hypatia's primal general conic form.
"""
mutable struct Model{T <: Real}
    n::Int
    p::Int
    q::Int
    obj_offset::T
    c::Vector{T}
    A
    b::Vector{T}
    G
    h::Vector{T}
    cones::Vector{Cones.Cone{T}}
    cone_idxs::Vector{UnitRange{Int}}
    nu::T

    function Model{T}(
        c::Vector{T},
        A,
        b::Vector{T},
        G,
        h::Vector{T},
        cones::Vector{<:Cones.Cone{T}};
        obj_offset::T = zero(T),
        ) where {T <: Real}
        model = new{T}()

        model.n = length(c)
        model.p = length(b)
        model.q = length(h)
        model.obj_offset = obj_offset
        model.c = c
        model.A = A
        model.b = b
        model.G = G
        model.h = h
        model.cones = cones
        model.cone_idxs = build_cone_idxs(model.q, model.cones)
        model.nu = isempty(cones) ? zero(T) : sum(Cones.get_nu, cones)

        return model
    end
end

function build_cone_idxs(q::Int, cones::Vector{Cones.Cone{T}}) where {T <: Real}
    cone_idxs = Vector{UnitRange{Int}}(undef, length(cones))
    prev_idx = 0
    for (k, cone) in enumerate(cones)
        dim = Cones.dimension(cone)
        cone_idxs[k] = (prev_idx + 1):(prev_idx + dim)
        prev_idx += dim
    end
    @assert q == prev_idx
    return cone_idxs
end

get_cone_idxs(model::Model) = model.cone_idxs

# make the model's A and G matrices dense
function densify!(model::Model{T}) where {T <: Real}
    model.A = convert(Matrix{T}, model.A)
    model.G = convert(Matrix{T}, model.G)
    return model
end

end


text:
module StructArrays

using Base: tuple_type_cons, tuple_type_head, tuple_type_tail, tail

export StructArray, StructVector, LazyRow, LazyRows
export collect_structarray
export replace_storage

include("interface.jl")
include("structarray.jl")
include("utils.jl")
include("collect.jl")
include("sort.jl")
include("lazy.jl")
include("tables.jl")
include("staticarrays_support.jl")

# Implement refarray and refvalue to deal with pooled arrays and weakrefstrings effectively
import DataAPI: refarray, refvalue
using DataAPI: defaultarray

refarray(s::StructArray) = StructArray(map(refarray, components(s)))

function refvalue(s::StructArray{T}, v::Tup) where {T}
    createinstance(T, map(refvalue, components(s), v)...)
end

# Use Adapt allows for automatic conversion of CPU to GPU StructArrays
import Adapt
Adapt.adapt_structure(to, s::StructArray) = replace_storage(x->Adapt.adapt(to, x), s)

end # module


text:
@testset "LPHRep" begin
    hs = HalfSpace([1, 0], 1)
    hp = HyperPlane([0, 1], 1)
    h = hs ∩ hp
    model = Model()
    @variable(model, [1:2] in h)
    lph = hrep(model)
    @test nhalfspaces(lph) == 1
    @test first(halfspaces(lph)) == hs
    @test nhyperplanes(lph) == 1
    @test first(hyperplanes(lph)) == hp
end


text:
@connector function HeatPort(; name)
    @variables T(t), Q_flow(t) # Temperature and Heat-flow-rate
    ODESystem(Equation[], t, [T, Q_flow], [], name=name)
end

function ModelingToolkit.connect(::Type{<:HeatPort}, ps...)
    eqs = [
           0 ~ sum(p->p.Q_flow, ps)
          ]
    for i in 1:length(ps)-1
        push!(eqs, ps[i].T ~ ps[i+1].T)
    end

    return eqs
end