text:
using .NetworkResilience

n = 100
ns = 50
nd = 50
q = 0.0
κ = 5.0
I = 1.0
D = 1.0
tol = 5e-4
k = 4
ensemblesize = 1000
alpha = 2.0
acs = []

println("Running cascacde simulations for q=0")
for i=1:ensemblesize
	println(i)
	ac, fmax, iter = swingcascadebisection(1, tol, n, ns, nd, q, k, κ, I, D)
	push!(acs, ac)
end

fname = "data/hist_q0_50_50_0.txt"
open(fname, "w") do f
	for i=1:length(acs)
		x = acs[i]
		write(f, "$x \n")
	end
end

q = 0.1
acs = []

println("Running cascacde simulations for q=0.1")
for i=1:ensemblesize
	println(i)
	ac, fmax, iter = swingcascadebisection(1, tol, n, ns, nd, q, k, κ, I, D)
	push!(acs, ac)
end

fname = "data/hist_q1_50_50_0.txt"
open(fname, "w") do f
	for i=1:length(acs)
		x = acs[i]
		write(f, "$x \n")
	end
end

max_stars_repo_name: o-smith/GridResilience

max_issues_repo_name: o-smith/GridResilience

max_forks_repo_name: o-smith/GridResilience


text:
@taylorize function vanderpol!(dx, x, params, t)
    local μ = 1.0
    dx[1] = x[2]
    dx[2] = (μ * x[2]) * (1 - x[1]^2) - x[1]
    return dx
end

function vanderpol()
   X0 = Hyperrectangle(low=[1.25, 2.35], high=[1.55, 2.45])
   ivp = @ivp(x' = vanderpol!(x), dim: 2, x(0) ∈ X0)
   tspan = (0.0, 5.0)
   return ivp, tspan
end

ext: jl


text:
using Test
using WordTokenizers

@testset "Words" begin
    @test split("A good day", Words) == ["A", "good", "day"]
    @test split("it cannot cannot be today", penn_tokenize) == ["it", "can", "not", "cannot", "be", "today"]
end

@testset "Sentences" begin
    @test split("Never going to give you up. Never going to let you down.", Sentences) ==
        ["Never going to give you up.", "Never going to let you down."]

    @test split.(split("Never going to give you up. Never going to let you down.", Sentences), Words) ==
        [["Never", "going", "to", "give", "you", "up", "."],
        ["Never", "going", "to", "let", "you", "down", "."]]

end


text:
using DashBootstrapComponents, DashHtmlComponents

backdrop_selector = html_div(
    [
        dbc_label("Backdrop:"),
        dbc_radioitems(
            id = "offcanvas-backdrop-selector",
            options = [
                Dict("label" => "True (default)", "value" => true),
                Dict("label" => "False", "value" => false),
                Dict("label" => "Static (no dismiss)", "value" => "static"),
            ],
            inline = true,
            value = true,
        ),
    ],
    className = "mb-2",
)

offcanvas = html_div([
    backdrop_selector,
    dbc_button("Open backdrop offcanvas", id = "open-offcanvas-backdrop", n_clicks = 0),
    dbc_offcanvas(
        html_p("Here is some content..."),
        id = "offcanvas-backdrop",
        title = "Offcanvas with/without backdrop",
        is_open = false,
    ),
])

callback!(
    app,
    Output("offcanvas-backdrop", "backdrop"),
    Input("offcanvas-backdrop-selector", "value"),
) do backdrop
    return backdrop
end;

callback!(
    app,
    Output("offcanvas-backdrop", "is_open"),
    Input("open-offcanvas-backdrop", "n_clicks"),
    State("offcanvas-backdrop", "is_open"),
) do n1, is_open
    return n1 > 0 ? is_open == 0 : is_open
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
include("main.jl");
include("SHEMS_optimizer_peer.jl");

function roll_SHEMS(market_flag, n_peers, n_market, h_start, h_end, h_predict, h_control, case)
    # Initialize technical setup__________________________________________________
    sh, pv, hp, fh, hw, b, m, p = set_SHEMS_parameters(market_flag, n_peers, n_market,
                                                h_start, h_end, h_predict, h_control, case);
    # intitial run_______________________________________________________________
    sh.soc_b, sh.soc_fh, sh.soc_hw, results_new  = SHEMS_optimizer_peer(sh, pv, hp, fh, hw, b, m, p);
    results = results_new;
    sh.h_start = m.h_start + m.h_control;

    # loop runs for rest of horizon______________________________________________
    for h = sh.h_start:m.h_control:(m.h_end-m.h_predict)
        sh.h_start = h;
        sh.soc_b, sh.soc_fh, sh.soc_hw, results_new  = SHEMS_optimizer_peer(sh, pv, hp, fh, hw, b, m, p);
        results = vcat(results, results_new);
    end

    # write to results folder____________________________________________________
    write_to_results_file(results, m, market_flag, n_peers, n_market, case)
    return nothing
end

function set_SHEMS_parameters(market_flag, n_peers, n_market, h_start, h_end, h_predict, h_control, case)
    # Initialize technical setup________________________________________________
    # Model_SHEMS(h_start, h_end, h_predict, h_control, big, rolling_flag, solver, mip_gap, output_flag, presolve_flag)
    m = Model_SHEMS(h_start, h_end,  h_predict, h_control, 60, "Gurobi", 0.005f0, 0, -1);
    # PV(eta_i)
    pv = PV(0.95f0);
    # HeatPump(rate_max)
    hp = HeatPump(3.0f0);
    # ThermalStorage(volume, loss, t_supply, soc_min, soc_max)
    fh = ThermalStorage(10.0f0, 0.045f0, 30.0f0, 20.0f0, 22.0f0);
    hw = ThermalStorage(180.0f0, 0.035f0, 45.0f0, 20.0f0, 180.0f0);
    # Battery(eta_i, eta_b, soc_min, soc_max, rate_max, loss)
    b = Battery(0.98f0, 0.0f0, 10.0f0, 4.6f0, 0.00003f0);
    # SHEMS(market_flag, n_peers, n_market, soc_b, soc_fh, soc_hw, h_start)
    sh = SHEMS(market_flag, n_peers, n_market, ones(n_peers*n_market).*0.0f0,
                ones(n_peers*n_market).*22.0f0, ones(n_peers*n_market).*180.0f0, h_start);

    if case==1 # feed-in tariff, full charge
        # Pricing(feedin, grid, peer, network, tax, eeg, others)
        p = Pricing(Matrix(CSV.read("data/tariff_matrix_case_1.csv")), 9.87, 12.19, 24.39, 7.71, 2.05, 6.756, 2.667)
    elseif case==2 # no feed-in tariff, full charge
        p = Pricing(Matrix(CSV.read("data/tariff_matrix_case_2.csv")), 9.87, 12.19, 24.39, 7.71, 2.05, 6.756, 2.667)
    elseif case==3 # feed-in tariff, support charge
        p = Pricing(Matrix(CSV.read("data/tariff_matrix_case_3.csv")), 9.87, 12.19, 24.39, 7.71, 2.05, 6.756, 2.667)
    elseif case==4 # no feed-in tariff, support charge
        p = Pricing(Matrix(CSV.read("data/tariff_matrix_case_4.csv")), 9.87, 12.19, 24.39, 7.71, 2.05, 6.756, 2.667)
    elseif case==5 # no feed-in tariff, area network
        p = Pricing(Matrix(CSV.read("data/tariff_matrix_case_5.csv")), 9.87, 12.19, 24.39, 7.71, 2.05, 6.756, 2.667)
    elseif case==6 # no feed-in tariff, area network
        p = Pricing(Matrix(CSV.read("data/tariff_matrix_case_6.csv")), 9.87, 12.19, 24.39, 7.71, 2.05, 6.756, 2.667)
    end

    return sh, pv, hp, fh, hw, b, m, p
end

function write_to_results_file(results, m, market_flag=1, n_peers=3, n_market=1, case=1)
    date=200916;
    CSV.write("results/$(date)_results_$(m.h_predict)_$(m.h_control)_$(m.h_start)-$(m.h_end)"*
                "_$(market_flag)_$(n_peers)_$(n_market)_$(case).csv",
            DataFrame(results), header=["Temp_FH", "Vol_HW",
            "SOC_B", "T_FH_plus", "T_FH_minus", "V_HW_plus", "V_HW_minus",  "profits", "sell",
            "COP_FH", "COP_HW", "Peer",
            "PV_DE", "B_DE", "GR_DE", "PV_B", "PV_GR", "PV_HP","GR_HP", "B_HP", "HP_FH", "HP_HW",
            "PV_PM", "B_PM", "PM_DE", "PM_HP",
            "month", "day", "hour"]);
    return nothing
end

function COPcalc(ts, t_outside)
    # Calculate coefficients of performance for every time period (1:h_predict)
    return cop = max.((5.8*ones(size(t_outside,1))) -(1. /14) * abs.((ts.t_supply*ones(size(t_outside,1)))-
                        t_outside), 0);
end

function run_cases(start_n, stop_n, start_c, stop_c)
    for n in start_n:stop_n
        print("$n  ");
        for c in start_c:stop_c
            print("$c  ")
            println(@elapsed roll_SHEMS(1, n, 1, 1, 8760, 36, 12, c))
        end
    end
    return nothing
end

max_stars_repo_name: lilanger/PEERS

max_issues_repo_name: lilanger/PEERS

max_forks_repo_name: lilanger/PEERS


text:
#============================================
    knn.jl

    Runs k-nearest-neighbors on review data
    to predict the average reviews of a business
    from it's location alone

    Haversine distance is used as the distance
    metric.
============================================#
using DataFrames
using ProgressMeter
include("confusion.jl")

#============================================
    Define the actual model and distance
    metrics
============================================#

"Haversin function"
@inline hsin(θ::Number) = sin(θ/2)^2

"""distance
computes the Haversin distance, or
the great circle approximation of distance on
earth
"""
@inline function distance(p1, p2)
    # convert to radians
    la1 = p1[1] * π/180
    lo1 = p1[2] * π/180
    la2 = p2[1] * π/180
    lo2 = p2[2] * π/180

    # Earth radius in METERS
    r = 6378100

    h = hsin(la2-la1) + cos(la1)*cos(la2)*hsin(lo2-lo1)

    return 2*r*asin(sqrt(h))
end

"""knn
runs knn on the dataset X with
labels y⃗, returning the prediction
for point x. The rows of X are datapoints
corresponding to the labels in the
rows of y⃗. x is a column vector.

This performs _classification_, so it
takes the _mode_ of the nearest k
points as the prediction, not the mean.
"""
function knn(x::Array, X::Array, y⃗::Array, k::Integer)
    @assert k < size(X)[1]
    @assert size(X)[1] == size(y⃗)[1]
    @assert size(x)[2] == size(X)[2]
    # create data matrix to sort
    A = [X y⃗]

    # create comparator to sort by the distance
    # to the point x
    lt = (x′,x′′) -> distance(x,x′[1:end-1]) < distance(x,x′′[1:end-1])

    # get initial neighbors by pulling from the top
    neighbors = sortrows(A[1:k,:], lt=lt)

    for i=(k+1):size(X)[1]
        for j=k:-1:1
            if lt(A[i,:],neighbors[j,:])
                # insert the new x into the array
                neighbors[end,:] = A[i,:]
                neighbors = sortrows(neighbors, lt=lt)
            else
                break
            end
        end
    end

    # predict the (rounded) mean of the nearest
    # neighbors' labels.
    return round(Integer,mean(neighbors[:,end]))
end

#============================================
    Run knn at different values of K on
    the actual data
============================================#

println("==> Loading business review data")
df = readtable("../stars.csv", separator='|')
X = Array(df[:,[:latitude,:longitude]])
y⃗ = round(Int,Array(df[:,:stars]*2)) #scale to range from 1-10

ks = [3]

println("==> Running knn for k=$(ks)")
acc = Accuracy(10)

log = open("knn.log","a")
write(log,"""
-----------------------------------
--        $(now())
-----------------------------------
""")

for k in ks
    log = open("knn.log", "a")
    start = now()
    @showprogress 1 "==> k=$(k)" for i in 1:size(X)[1]
        # remove the datapoint we're predicting on
        pred = knn(X[i,:], X[1:end .!=i,:], y⃗[1:end .!=i], k)
        predict!(acc, pred, y⃗[i])
    end
    println("==> Done with k=$(k)!\n")

    # write results to the file
    write(log, """
-----------------------------------
k=$(k)
elapsed:   $(now()-start)
mean iter: $(convert(Int,now()-start)/size(X)[1])
accuracy:  $(accuracy(acc))
precision: $(precision(acc))
recall:    $(recall(acc))\n
confusion:
$(acc.prediction)\n\n""")
    close(log)
end

max_stars_repo_name: cdipaolo/inflection

max_issues_repo_name: cdipaolo/inflection

max_forks_repo_name: cdipaolo/inflection


text:
using MIDI

cd(@__DIR__)

@testset "MIDI IO" begin
    midi = load("doxy.mid")
    @test midi.tpq == 960
    @test length(midi.tracks) == 4

    notes = getnotes(midi, 4)

    @test length(notes) == 533
    @test notes.tpq == 960

    @test notes[1].pitch == 65
    @test notes[1].velocity == 69
    @test notes[1].position == 7427
    @test notes[1].duration == 181
    @test notes[1].channel == 0

    @test notes[2].pitch == 70
    @test notes[2].velocity == 85
    @test notes[2].position == 7760
    @test notes[2].duration == 450
    @test notes[2].channel == 0

end

@testset "Deprecated MIDI IO" begin
    midi = load("doxy.mid")
    @test_deprecated readMIDIFile("doxy.mid")
    @test_deprecated writeMIDIFile("test.mid", midi)
    @test_deprecated writeMIDIFile("test.midi", getnotes(midi))
    rm("test.mid")
    rm("test.midi")
end


text:
include(joinpath("..","src","Ogs6InputFileHandler.jl"))
include(joinpath("..","src","VTUFileHandler.jl"))

file = joinpath(".","Point_injection","square_1e2_lin.prj")
ogs6md = read(Ogs6ModelDef, file)
write(ogs6md)

file1 = joinpath(".","Point_injection","expected_square_1e0_ts_10_t_1000.000000.vtu")
file2 = joinpath(".","Point_injection","expected_square_1e0_lin_ts_10_t_50000.000000.vtu")
file3 = joinpath(".","Point_injection","expected_square_1e0_ts_10_t_50000.000000.vtu")

vtu1 = VTUFile(file1);
vtu2 = VTUFile(file2);
vtu3 = VTUFile(file3);

write(vtu1)

tmp = similar(vtu1.data)

max_stars_repo_name: baxmittens/XMLParser

max_issues_repo_name: baxmittens/XMLParser

max_forks_repo_name: baxmittens/XMLParser


text:
@with_kw mutable struct Agent
    get_initial_entity::Function # Returns an entity IC
    model::DriverModel # The driver model associated with this agent
    entity_dim::Int # The dimension of the entity
    disturbance_dim::Int # The disturbance dimension of this agent
    entity_to_vec::Function # A Function that converts the agent to a vector of length o_dim
    disturbance_to_vec::Union{Function, Nothing} = nothing # A Function that converts an agent action to a vector of length a_dim
    vec_to_entity::Union{Function, Nothing} = nothing # A Function that converts a vector of length o_dim to an entity
    vec_to_disturbance::Union{Function, Nothing} = nothing # A Function that converts a vector of length a_dim to an action
    disturbance_model = nothing  # Model of the disturbances (supports logpdf, rand and actions)
end

id(a::Agent) = a.get_initial_entity().id

# Construct a regular Blinker vehicle agent
function BlinkerVehicleAgent(get_veh::Function, model::TIDM;
    entity_dim = BLINKERVEHICLE_ENTITY_DIM,
    disturbance_dim=BLINKERVEHICLE_DISTURBANCE_DIM,
    entity_to_vec = BlinkerVehicle_to_vec,
    disturbance_to_vec = BlinkerVehicleControl_to_vec,
    vec_to_entity = vec_to_BlinkerVehicle,
    vec_to_disturbance = vec_to_BlinkerVehicleControl,
    disturbance_model = get_bv_actions())
    Agent(get_veh, model, entity_dim, disturbance_dim, entity_to_vec,
          disturbance_to_vec,  vec_to_entity, vec_to_disturbance, disturbance_model)
end

# Construct a regular adversarial pedestrian agent
function NoisyPedestrianAgent(get_ped::Function, model::AdversarialPedestrian;
    entity_dim = PEDESTRIAN_ENTITY_DIM,
    disturbance_dim = PEDESTRIAN_DISTURBANCE_DIM,
    entity_to_vec = NoisyPedestrian_to_vec,
    disturbance_to_vec = PedestrianControl_to_vec,
    vec_to_entity = vec_to_NoisyPedestrian_fn(DEFAULT_CROSSWALK_LANE),
    vec_to_disturbance = vec_to_PedestrianControl,
    disturbance_model = get_ped_actions())
    Agent(get_ped, model, entity_dim, disturbance_dim, entity_to_vec,
          disturbance_to_vec,  vec_to_entity, vec_to_disturbance, disturbance_model)
end

# Definition of the adversarial driving mdp
mutable struct AdversarialDrivingMDP <: MDP{Scene, Vector{Disturbance}}
    agents::Vector{Agent} # All the agents ordered by (adversaries..., sut, others...)
    vehid2ind::Dict{Int64, Int64} # Dictionary that maps vehid to index in agent list
    num_adversaries::Int64 # The number of adversaries
    roadway::Roadway # The roadway for the simulation
    dt::Float64 # Simulation timestep
    last_observation::Array{Float64} # Last observation of the vehicle state
    disturbance_model # Model used for disturbances. supports `logpdf` and `rand` and `actions`
    γ::Float64 # discount
    ast_reward::Bool # A function that gives action log prob.
    no_collision_penalty::Float64 # penalty for not getting a collision (for ast reward)
    scale_reward::Bool #whether or not to scale the AST reward
    end_of_road::Float64 # specify an early end of the road
end

# Constructor
function AdversarialDrivingMDP(sut::Agent, adversaries::Vector{Agent}, road::Roadway, dt::Float64;
                               other_agents::Vector{Agent} = Agent[],
                               γ = 1,
                               ast_reward = false,
                               no_collision_penalty = 1e3,
                               scale_reward = true,
                               end_of_road = Inf,)
    agents = [adversaries..., sut, other_agents...]
    d = Dict(id(agents[i]) => i for i=1:length(agents))
    Na = length(adversaries)
    o = Float64[] # Last observation

    m = combine_disturbance_models(adversaries)
    AdversarialDrivingMDP(agents, d, Na, road, dt, o, m, γ,
                         ast_reward, no_collision_penalty, scale_reward, end_of_road)
end

# Returns the intial state of the mdp simulator (as a Deteriministic distribution)
function POMDPs.initialstate(mdp::MDP{Scene, A}, rng::AbstractRNG = Random.GLOBAL_RNG) where A
     POMDPModelTools.Deterministic(Scene([a.get_initial_entity(rng) for a in agents(mdp)]))
 end

# The generative interface to the POMDP
function POMDPs.gen(mdp::MDP{Scene, A}, s::Scene, a::A, rng::Random.AbstractRNG = Random.GLOBAL_RNG) where A
    mdp.last_observation = convert_s(AbstractArray, s, mdp)
    sp = step_scene(mdp, s, a, rng)
    r = reward(mdp, s, a, sp)
    (sp=sp, r=r)
end

# Get the reward from the actions taken and the next state
function POMDPs.reward(mdp::AdversarialDrivingMDP, s::Scene, a::Vector{Disturbance}, sp::Scene)
    iscollision = length(sp) > 0 && ego_collides(sutid(mdp), sp)
    if mdp.ast_reward
        isterm = isterminal(mdp, sp)
        r = logpdf(mdp, s, a)
        r += iscollision * abs(mdp.no_collision_penalty)
        # r = (isterm && !iscollision)*(-abs(mdp.no_collision_penalty)) + logpdf(mdp, s, a)
        mdp.scale_reward && (r = r / abs(mdp.no_collision_penalty))
        return Float32(r)
    else
        return Float32(iscollision)
    end
end

# Discount factor for the POMDP (Set to 1 because of the finite horizon)
POMDPs.discount(mdp::MDP{Scene, A}) where A = mdp.γ

# The simulation is terminal if there is collision with the ego vehicle or if the maximum simulation time has been reached
POMDPs.isterminal(mdp::MDP{Scene, A}, s::Scene) where A = !(sutid(mdp) in s)|| any_collides(s)

# Define the set of actions, action index and probability
POMDPs.actions(mdp::AdversarialDrivingMDP) = get_actions(mdp.disturbance_model)
POMDPs.actionindex(mdp::AdversarialDrivingMDP, a::Vector{Disturbance}) = Int32(get_actionindex(mdp.disturbance_model, a))

# The default disturbance policy according to the disturbance models
#TODO: Switch the "Function policy" to something that can get logpdf?
default_policy(mdp::AdversarialDrivingMDP, rng::AbstractRNG = Random.GLOBAL_RNG) = FunctionPolicy((s) -> rand(rng, mdp.disturbance_model))

# POMDPs.actionindex(mdp::AdversarialDrivingMDP, a::Vector{Disturbance}) = findfirst(actions(mdp) .== a)
Distributions.logpdf(mdp::AdversarialDrivingMDP, a::Vector{Disturbance}) = logpdf(mdp.disturbance_model, a, mdp)
Distributions.logpdf(mdp::AdversarialDrivingMDP, s::Scene, a::Vector{Disturbance}) = logpdf(mdp, a)
Distributions.logpdf(mdp::AdversarialDrivingMDP, h::SimHistory) = sum([logpdf(mdp, s, a) for (s,a) in eachstep(h, (:s, :a))])
Distributions.logpdf(mdp::AdversarialDrivingMDP, as::Vector) = sum([logpdf(mdp, a) for a in as])

## Helper functions

# Step the scene forward by one timestep and return the next state
function step_scene(mdp::AdversarialDrivingMDP, s::Scene, actions::Vector{Disturbance}, rng::AbstractRNG = Random.GLOBAL_RNG)
    entities = []

    # Loop through the adversaries and apply the instantaneous aspects of their disturbance
    for (adversary, action) in zip(adversaries(mdp), actions)
        update_adversary!(adversary, action, s)
    end

    # Loop through the vehicles in the scene, apply action and add to next scene
    for (i, veh) in enumerate(s)
        m = model(mdp, veh.id)
        observe!(m, s, mdp.roadway, veh.id)
        a = rand(rng, m)
        bv = Entity(propagate(veh, a, mdp.roadway, mdp.dt), veh.def, veh.id)
        !end_of_road(bv, mdp.roadway, mdp.end_of_road) && push!(entities, bv)
    end
    isempty(entities) ? Scene(typeof(sut(mdp).get_initial_entity())) : Scene([entities...])
end

# Returns the list of agents in the mdp
agents(mdp::MDP{Scene,A}) where A = mdp.agents

# Returns the list of adversaries in the mdp
adversaries(mdp::MDP{Scene,A}) where A = view(mdp.agents, 1:mdp.num_adversaries)

# Returns the model associated with the vehid
model(mdp::MDP{Scene,A}, vehid::Int) where A = mdp.agents[mdp.vehid2ind[vehid]].model

# Returns the system under test
sut(mdp::MDP{Scene,A}) where A  = mdp.agents[mdp.num_adversaries + 1]

# Returns the sut id
sutid(mdp::MDP{Scene,A}) where A = id(sut(mdp))

function update_adversary!(adversary::Agent, action::Disturbance, s::Scene)
    index = findfirst(id(adversary), s)
    isnothing(index) && return nothing # If the adversary is not in the scene then don't update
    adversary.model.next_action = action # Set the adversaries next action
    veh = s[index] # Get the actual entity
    state_type = typeof(veh.state) # Find out the type of its state
    s[index] =  Entity(state_type(veh.state, noise = action.noise), veh.def, veh.id) # replace the entity in the scene
end

## SUT driving MDP
mutable struct DrivingMDP <: MDP{Scene, BlinkerVehicleControl}
    agents::Vector{Agent} # All the agents ordered by (adversaries..., sut, others...)
    vehid2ind::Dict{Int64, Int64} # Dictionary that maps vehid to index in agent list
    num_adversaries::Int64 # The number of adversaries
    roadway::Roadway # The roadway for the simulation
    dt::Float64 # Simulation timestep
    last_observation::Array{Float64} # Last observation of the vehicle state
    γ::Float64 # discount
    end_of_road::Float64 # Early stopping of road
    per_timestep_penalty::Float64
    v_des::Float64
end

# Constructor
function DrivingMDP(sut::Agent, adversaries::Vector{Agent}, road::Roadway, dt::Float64; γ = 1, end_of_road = Inf, per_timestep_penalty = 0, v_des = 25)
    agents = [adversaries..., sut]
    d = Dict(id(agents[i]) => i for i=1:length(agents))
    DrivingMDP(agents, d, length(adversaries), road, dt, Float64[], γ, end_of_road, per_timestep_penalty, v_des)
end

# Get the reward from the actions taken and the next state
function POMDPs.reward(mdp::DrivingMDP, s::Scene, a::BlinkerVehicleControl, sp::Scene)
    id = sutid(mdp)
    v = vel(get_by_id(s, id))

    r = -abs(mdp.per_timestep_penalty)
    # If the simulation ends but the SUT is not at the end of the road, big penalty
    if ego_collides(id, sp)
        r += -1
    elseif isterminal(mdp, sp)
        r += 1
    end
    r += -.1 * (v < 0)
    r += -0.001 * abs(v  - mdp.v_des)
    r
end

function step_scene(mdp::DrivingMDP, s::Scene, action::BlinkerVehicleControl, rng::AbstractRNG = Random.GLOBAL_RNG)
    entities = []
    sid = sutid(mdp)

    # Choose random actions for the adversaries
    for adversary in adversaries(mdp)
        adv_action = rand(rng, adversary.disturbance_model, adversary.vec_to_disturbance)
        update_adversary!(adversary, adv_action, s)
    end

    # Loop through the vehicles in the scene, apply action and add to next scene
    for (i, veh) in enumerate(s)
        if veh.id == sid # for the sut, use the prescribed action
            a = action
        else # For the other vehicles use their  model
            m = model(mdp, veh.id)
            observe!(m, s, mdp.roadway, veh.id)
            a = rand(rng, m)
        end
        bv = Entity(propagate(veh, a, mdp.roadway, mdp.dt), veh.def, veh.id)
        !end_of_road(bv, mdp.roadway, mdp.end_of_road) && push!(entities, bv)
    end
    isempty(entities) ? Scene(typeof(sut(mdp).get_initial_entity())) : Scene([entities...])
end

# Define the set of actions, action index and probability
POMDPs.actions(mdp::DrivingMDP) = [BlinkerVehicleControl(a = -4.), BlinkerVehicleControl(a= -2.), BlinkerVehicleControl(a = 0.), BlinkerVehicleControl(a = 1.5), BlinkerVehicleControl(a = 3.)]
POMDPs.actionindex(mdp::DrivingMDP, a::BlinkerVehicleControl) = findfirst([a] .== actions(mdp))

max_stars_repo_name: ancorso/MAST

max_issues_repo_name: ancorso/MAST


text:
"""
lsnormal(A,b)

Solve a linear least squares problem by the normal equations.
Returns the minimizer of ||b-Ax||.
"""
function lsnormal(A,b)

N = A'*A;  z = A'*b;
R = cholesky(N).U
w = forwardsub(R',z)                   # solve R'z=c
x = backsub(R,w)                       # solve Rx=z

return x
end

"""
lsqrfact(A,b)

Solve a linear least squares problem by QR factorization. Returns
the minimizer of ||b-Ax||.
"""
function lsqrfact(A,b)

Q,R = qr(A)
c = Q'*b
x = backsub(R,c)

return x
end

"""
qrfact(A)

QR factorization by Householder reflections. Returns Q and R.
"""
function qrfact(A)

    m,n = size(A)
    Qt = Matrix(Diagonal(ones(m)))
    R = float(copy(A))
    for k in 1:n
      z = R[k:m,k]
      v = [ -sign(z[1])*norm(z) - z[1]; -z[2:end] ]
      nrmv = norm(v)
      if nrmv < eps() continue; end  # skip this iteration
      v = v / nrmv;                  # simplifies other formulas
      # Apply the reflection to each relevant column of A and Q
      for j in 1:n
        R[k:m,j] -= v*( 2*(v'*R[k:m,j]) )
      end
      for j in 1:m
        Qt[k:m,j] -= v*( 2*(v'*Qt[k:m,j]) )
      end
    end

    return Qt',triu(R)
end


text:
# !! auto-generated code, do not edit !!
module std_msgs
using RobotOSData.Messages
struct MultiArrayDimension <: Readable
    label::String
    size::UInt32
    stride::UInt32
end
struct MultiArrayLayout <: Readable
    dim::Vector{MultiArrayDimension}
    data_offset::UInt32
end
struct UInt32MultiArray <: Readable
    layout::MultiArrayLayout
    data::Vector{UInt32}
end
struct UInt16MultiArray <: Readable
    layout::MultiArrayLayout
    data::Vector{UInt16}
end
struct ByteMultiArray <: Readable
    layout::MultiArrayLayout
    data::Vector{UInt8}
end
struct UInt64MultiArray <: Readable
    layout::MultiArrayLayout
    data::Vector{UInt64}
end
struct Char <: Readable
    data::UInt8
end
struct UInt8MultiArray <: Readable
    layout::MultiArrayLayout
    data::Vector{UInt8}
end
struct Int8MultiArray <: Readable
    layout::MultiArrayLayout
    data::Vector{Int8}
end
struct ColorRGBA <: Readable
    r::Float32
    g::Float32
    b::Float32
    a::Float32
end
struct Duration <: Readable
    data::ROSDuration
end
struct Int64MultiArray <: Readable
    layout::MultiArrayLayout
    data::Vector{Int64}
end
struct Time <: Readable
    data::ROSTime
end
struct Int32MultiArray <: Readable
    layout::MultiArrayLayout
    data::Vector{Int32}
end
struct Header <: Readable
    seq::UInt32
    stamp::ROSTime
    frame_id::String
end
struct Float32MultiArray <: Readable
    layout::MultiArrayLayout
    data::Vector{Float32}
end
struct Float64MultiArray <: Readable
    layout::MultiArrayLayout
    data::Vector{Float64}
end
struct Empty <: Readable
end
struct Int16MultiArray <: Readable
    layout::MultiArrayLayout
    data::Vector{Int16}
end
struct Byte <: Readable
    data::UInt8
end
end


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

# %%
# Ensure your working directory is: https://github.com/henriquebecker91/phd/tree/master/latex/revised_PPG2KP
include("notebook_setup.jl")

# %%
# Read the data, show nothing.
raw_fr_csv_path = "./data/faithful_reimplementation.csv"
raw_fr_df = DataFrame(CSV.File(raw_fr_csv_path))
nothing
#showtable(fr_df)

# %%
# Clean the data a little.
fr_df = let fr_df = deepcopy(raw_fr_df)
    # Keep only the instance name (not the path).
    @with(fr_df, :instance_name .= basename.(:instance_name))
    @with(fr_df, :datafile .= basename.(:datafile))
    fr_df[!, :finished] = .!fr_df[!, :run_ended_by_exception]
    colnames = names(fr_df)
    for name in colnames
        column = fr_df[!, name]
        if eltype(column) <: AbstractFloat
            fr_df[!, name] = ifelse.(isnan.(column), missing, column)
        elseif eltype(column) <: Integer && occursin(r"^qt_", string(name))
            fr_df[!, name] = ifelse.(column .< 0, missing, column)
        end
    end
    if any(name -> name ∈ colnames, ["pricing_method", "disabled_redundant_cut", "disabled_cut_position"])
        args2id = Dict{NTuple{3, Bool}, String}(
            (false, true, true) => "complete",
            (false, false, true) => "only_redundant_cut",
            (false, true, false) => "only_cut_position",
            (false, false, false) => "both_reductions",
            (true, false, false) => "priced"
        )
        fr_df.model_variant = getindex.((args2id,), tuple.(
            fr_df.pricing_method .== "furini", fr_df.disabled_redundant_cut, fr_df.disabled_cut_position
        ))
        select!(fr_df, Not([:pricing_method, :disabled_redundant_cut, :disabled_cut_position]))
    end

    fr_df
end
nothing

# %%
# Shows the cleaned data.
println(names(fr_df))
showtable(fr_df)

# %%
# Shows info about the unfinished runs. Theoretically those three conditions should always
# appear together in a row, uf some row do not have all of them, then an error happened that
# was not a timeout.
let
    unfinished = @linq fr_df |>
        where(.!:finished .| :had_timeout .| (:build_stop_reason .== "NOT_REACHED")) |>
        select(:instance_name, :model_variant, :datafile, :finished, :had_timeout, :build_stop_reason)
    showtable(unfinished)
end

# %%
# Shows info about OPTIMAL_FOUND runs
@linq fr_df |>
    where(:build_stop_reason .== "FOUND_OPTIMUM") |>
    select(:instance_name, :model_variant, :datafile, :finished, :build_stop_reason) |>
    (tmp_df -> for s in sort(tmp_df[!, :instance_name]); println(s); end)()

# %%
# Check if the number of lines per model variant is correct.
@linq fr_df |>
    groupby(:model_variant) |>
    based_on(; qt = length(:model_variant))

# %%
# The data extracted from DOI 10.6092/unibo/amsdottorato/7399 by hand. The amount of variables in many
# model variants for the 59 instances used by Thomopulos in their 2016 thesis.
# The columns 'vars' and 'plates' have absolute values for rows in which 'model_variant' column is
# 'complete', and percentages for the rest, as in the thesis tables.
raw_th59_csv_path = "./data/thomopulos_2016.csv"
raw_th59 = DataFrame(CSV.File(raw_th59_csv_path))
showtable(raw_th59)
#nothing

# %%
const THOMOPULOS_THESIS_INSTANCES = vcat(
  # unweighted
  "gcut" .* string.(1:12)
  , String.(split("wang20 2s 3s A1s A2s STS2s STS4s"))
  , ["OF1", "OF2", "W", "CHL1s", "CHL2s"]
  , "A" .* string.(3:5)
  , "CHL" .* string.(5:7)
  , ["CU1", "CU2"]
  , "Hchl" .* split("3s 4s 6s 7s 8s")
  # weighted
  , "cgcut" .* string.(1:3)
  , "okp" .* string.(1:5)
  , String.(split("HH 2 3 A1 A2 STS2 STS4 CHL1 CHL2"))
  , "CW" .* string.(1:3)
  , "Hchl" .* ["2", "9"]
) :: Vector{String}
nothing

# %%
# Absolutize the vars and plates columns, i.e., now thay all refer to absolute
# quantities of variables and plates, instead of having all rows with
# 'model_variant != complete' with percentages (as it is in the CSV file).
th59 = let th59 = deepcopy(raw_th59)
    # Defines the 'origin' column.
    th59[!, :origin] .= "thomopulosthesis"
    # Fix the missing plates. The number of plates is not described in the thesis but
    # a safe interpretation is that the number stayed the same from the previous step.
    if any(ismissing, th59[!, :plates])
        reduced_plates = Dict{String, Float64}()
        restricted_plates = Dict{String, Float64}()

        @byrow! th59 begin
            if :model_variant == "both_reductions"
                reduced_plates[:instance_name] = :plates
            elseif :model_variant == "restricted"
                restricted_plates[:instance_name] = :plates
            end
        end
        th59 = @byrow! th59 begin
            if :model_variant == "priced"
                :plates = reduced_plates[:instance_name]
            elseif :model_variant == "after_restricted_pricing"
                :plates = restricted_plates[:instance_name]
            end
        end
    end
    #show(@where(th59, ismissing.(:plates)))
    th59[!, :plates] = Float64.(th59[!, :plates])
    if eltype(th59[!, :vars]) != Int
        th59_complete = @where(th59, :model_variant .== "complete")
        stats_type = NamedTuple{(:vars, :plates),Tuple{Int64,Int64}}
        th59_complete_dict = Dict{String, stats_type}()
        @byrow! th59_complete begin
            th59_complete_dict[:instance_name] = (vars = :vars, plates = :plates)
        end
        for column in [:vars, :plates]
            complete_stats = getindex.(getindex.((th59_complete_dict,), th59.instance_name), column)
            th59[!, column] = ifelse.(
                th59.model_variant .== "complete",
                complete_stats,
                round.(Int64, complete_stats .* th59[!, column] ./ 100.0)
            )
        end
    end
    th59
end
showtable(th59)
nothing

# %%
@show unique(th59[!, :model_variant])
@show names(fr_df)[.!isnothing.(match.(r"qt_cmvars.*", names(fr_df)))]
nothing

# %%
# Put the faithful reimplementation table in the same format as the data extracted from
# Thomopulos 2016 (initially, only include the rows with model_variant in "complete", "only_cut_position",
# "only_redundant_cut" and "both_reductions").
short_fr_df = let short_fr_df = deepcopy(fr_df)
    sel_cols = [
        :instance_name, :model_variant, :qt_pevars_after_preprocess, :qt_cmvars_after_preprocess,
        :qt_plates_after_preprocess
    ]
    short_fr_df = select(short_fr_df, sel_cols)
    short_fr_df[!, :origin] .= "reimplementation"
    variants = ("complete", "only_cut_position", "only_redundant_cut", "both_reductions")
    short_fr_df = @where(short_fr_df, in.(:model_variant, (variants,)))
    # basically the 'x' and 'y' of the model
    var_categories = [:qt_pevars_after_preprocess, :qt_cmvars_after_preprocess]
    select!(short_fr_df, var_categories => ((a,b) -> a .+ b) => :vars, Not(var_categories))
    rename!(short_fr_df, :qt_plates_after_preprocess => :plates)

    short_fr_df
end

# %%
# Add thee rows with model_variant equal to 'priced', 'restricted', 'after_restricted_pricing',
# and purged to 'short_fr_df', those come from the 'priced' runs (i.e., the 'priced' runs execute
# the restricted pricing as a phase of their process, so we do not need runs just to
# get such values, they are just in another column of 'priced' rows).
short_fr_df = if "restricted" in short_fr_df[!, :model_variant]
    short_fr_df
else let short_fr_df = deepcopy(short_fr_df)
    needed_cols = [
        :instance_name, :model_variant,
        # restricted
        :qt_pevars_restricted, :qt_cmvars_restricted, :qt_plates_restricted,
        # priced restricted
        :qt_pevars_priced_restricted, :qt_cmvars_priced_restricted, :qt_plates_priced_restricted,
        # after purge
        :qt_pevars_after_purge, :qt_cmvars_after_purge, :qt_plates_after_purge,
        # after pricing
        :qt_pevars_after_pricing, :qt_cmvars_after_pricing, :qt_plates_after_pricing
    ]
    # Get just the needed columns and rows.
    cut = @linq fr_df |> select(needed_cols) |> where(:model_variant .== "priced")
    # Create the origin column.
    cut[!, :origin] .= "reimplementation"

    function to_th59_format!(df, col_suffix, final_cols, variant_name)
        plate_colname = Symbol("qt_plates_" * string(col_suffix))
        var_colnames = Symbol[
            Symbol("qt_pevars_" * string(col_suffix)),
            Symbol("qt_cmvars_" * string(col_suffix))
        ]
        rename!(df, plate_colname => :plates)
        transform!(df, var_colnames => ((a,b) -> a .+ b) => :vars)
        df[!, :model_variant] .= variant_name
        select!(df, final_cols)
        df
    end

    final_columns = [:instance_name, :model_variant, :vars, :plates, :origin]
    restricted = to_th59_format!(deepcopy(cut), "restricted", final_columns, "restricted")
    priced = to_th59_format!(deepcopy(cut), "after_pricing", final_columns, "priced")
    priced_r = to_th59_format!(deepcopy(cut), "priced_restricted", final_columns, "after_restricted_pricing")
    purged = to_th59_format!(deepcopy(cut), "after_purge", final_columns, "purged")

    vcat(short_fr_df, restricted, priced_r, priced, purged; cols = :setequal)
end #= let =# end #= end =#

# %%
# Create a mock "purged" rowset for th59 that is a copy of the
# "priced" to allow comparison with the extra data generated in
# our experiments.
th59 = let
    mock_purged = @where(th59, :model_variant .== "priced")
    mock_purged[!, :model_variant] .= "purged"
    vcat(th59, mock_purged; cols = :setequal)
end

# %%
# Just show how much data is missing for each variant.
qt_missing = combine(
    #groupby(th59, :model_variant),
    groupby(short_fr_df, :model_variant),
    :model_variant => length => :num_rows,
    :vars => (sum ∘ broadwrap(ismissing)) => :missing_vars,
    :plates => (sum ∘ broadwrap(ismissing) => :missing_plates)
)
showtable(qt_missing)

# %%
# Get full info for each run that has missing values:
missing_runs = filter([:vars, :plates] => ((v, p) -> ismissing(v) || ismissing(p)), short_fr_df)
#@show missing_runs
missing_ids = @with(missing_runs, tuple.(:model_variant, :instance_name))
#@show missing_ids
full_missing_runs = filter([:model_variant, :instance_name] => function (variant, name)
    return (variant, name) in missing_ids
end, fr_df)
showtable(full_missing_runs)

# %%
# This is a rare case. The LB and UB of the restricted model touch, and the
# priced restricted model is not solved.
showtable(@where(fr_df, (:instance_name .== "A1s") .& (:model_variant .== "priced")))

# %%
# Removes from both dataframes ('th59' and 'short_fr_df') all rows that have 'missing'
# in the 'vars' column of 'short_fr_df'
th59, short_fr_df = if any(ismissing, short_fr_df[!, :vars])
    let short_fr_df = deepcopy(short_fr_df), th59 = deepcopy(th59)
        expected_cols = sort!(["instance_name", "model_variant", "origin", "vars", "plates"])
        @assert sort!(names(th59)) == expected_cols
        @assert sort!(names(short_fr_df)) == expected_cols
        table_keys = ["instance_name", "model_variant", "origin"]
        sort!(th59, table_keys)
        sort!(short_fr_df, table_keys)
        to_delete = ismissing.(short_fr_df[!, "vars"])
        delete!(th59, to_delete)
        delete!(short_fr_df, to_delete)
        @assert th59[!, "instance_name"] == short_fr_df[!, "instance_name"]
        @assert th59[!, "model_variant"] == short_fr_df[!, "model_variant"]
        th59, short_fr_df
    end
else
    th59, short_fr_df
end

# %%
full_df = vcat(th59, short_fr_df; cols = :setequal)
nothing
#showtable(fulldf)

# %%
showtable(
    #@linq @where(full_df, (:model_variant .== "priced")) |>
    @linq @where(full_df, (:instance_name .== "gcut1")) |>
        sort!([:instance_name, :origin])
)

# %%
# Transforms the long format into a wide one and rename many columns and
# string cells to the final name they will have in the table.
function make_header(origin, variable)
    @assert origin ∈ ["reimplementation", "thomopulosthesis"]
    @assert variable ∈ ["vars", "plates"]
    prefix = origin == "reimplementation" ? "R. %" : "O. #"
    suffix = variable == "vars" ? "v" : "p"
    return prefix * suffix
end
make_headers(origins, variables) = make_header.(origins, variables)

#transform!(full_df, :origin => broadwrap(Symbol) => :origin)
#transform!(full_df, :model_variant => broadwrap(Symbol) => :model_variant)

fr_table = let full_df = deepcopy(full_df)
    fr_table = @linq full_df |>
        select(:model_variant, :origin, :vars, :plates) |>
        groupby([:model_variant, :origin]) |>
        based_on(; vars = sum(:vars), plates = sum(:plates)) |>
        stack([:vars, :plates]) |> # vars and plates become values in column 'variable', values in 'value'
        unstack([:origin, :variable], :model_variant, :value) |>
        select!([:origin, :variable] => make_headers => :header, Not([:origin, :variable])) |>
        stack(Not(:header)) |>
        unstack(:variable, :header, :value)

    fr_table[!, "R. %p"] = round.(100.0 .* fr_table[!, "R. %p"] ./ fr_table[!, "O. #p"]; digits = 2)
    fr_table[!, "R. %v"] = round.(100.0 .* fr_table[!, "R. %v"] ./ fr_table[!, "O. #v"]; digits = 2)
    #fr_table[:, ["variable", "O. #v", "R. %v", "O. #p", "R. %p"]]
    select!(fr_table, ["variable", "O. #v", "R. %v", "O. #p", "R. %p"])
    select!(fr_table, :variable => "Variant", Not(:variable))
    pretty_variant_names = Dict{String, Any}(
        "complete" => (longname = "Complete PP-G2KP", order = 1),
        "only_cut_position" => (longname = "Complete +Cut-Position", order = 2),
        "only_redundant_cut" => (longname = "Complete +Redundant-Cut", order = 3),
        "both_reductions" => (longname = "PP-G2KP (CP + RC)", order = 4),
                "restricted" => (longname = "Restricted PP-G2KP", order = 5),
        "after_restricted_pricing" => (longname = "Priced Restricted PP-G2KP", order = 6),
        "priced" => (longname = "(no purge) Priced PP-G2KP", order = 7),
        "purged" => (longname = "Priced PP-G2KP", order = 8),
    )
    sort!(fr_table, "Variant"; by = (name -> pretty_variant_names[name].order))
    fr_table[!, "Variant"] = getindex.(getindex.((pretty_variant_names,), fr_table[!, "Variant"]), :longname)

    fr_table
end

# %%
# Finally, the dataframe is transformed into a latex table.
let latex_df = deepcopy(fr_table)
    for col in ["O. #v", "R. %v", "O. #p", "R. %p"]
        latex_df[!, col] = number2latex.(latex_df[!, col])
    end
    # Needs to escape some symbols for latex.
    rename!(latex_df,
        ["Variant", "O. \\#v", "R. \\%v",  "O. \\#p", "R. \\%p"]
    )
    #=rename!(
        "Variant", "O. #v" => "O. \\#v", "R. %v" => "R. \\%v", "O. #p" => "O. \\#p",
        "R. %p" => "R. \\%p"
    )=#
    #=
    for int_col in ["O. #v", "O. #p"]
        latex_df[!, int_col] =
    end
    for float_col in ["R. %v", "R. %p"]
    end=#
    pretty_table(
        latex_df; backend = :latex, nosubheader = true, alignment = [:l, :r, :r, :r, :r]
    )
end


text:
#クロスバリデーションとグリッドチューンを同時実行する関数
function searchparm(Y,X;nfeature=0,ntree=[500,1000],nrate=[0.7],CV_n=5,rev=false)
  Ytraining=Y;
  Xtraining=X;
  CV_n=CV_n;

  #ランダムフォレストによる予測モデルの構築およびクロスバリデーションの実行
  estfun=Estfun(Ytraining,Xtraining,CV_n)

  #評価関数の定義
  evalfun=Evalfun(Ytraining,Xtraining,CV_n)

  #パラメータの自動算定
  if(nfeature==0)
    point1=Integer(round(size(Xtraining)[2]/3,0))
    point2=Integer(round(sqrt(size(Xtraining)[2]),0))
    nfeature=[1,2,Integer(round((point1+2)/2)),point1,Integer(round((point1+point2)/2)),point2,Integer(round((size(Xtraining)[2]+point2)/2)),size(Xtraining)[2]]
  end

     #特徴量の重複を除去
    nfeature=DataFrames.DataFrame(ind=nfeature)
    nfeature=DataFrames.by(nfeature, :ind, df -> size(df, 1))
    nfeature=convert(Array,nfeature[:,:ind])
    nfeature=sort(nfeature,rev=rev)

  # グリッドサーチの実行
  return  r=gridtune(estfun , evalfun,
                        ( "nfeature",nfeature ) ,
                        ( "ntree",ntree ) ,
                        (" nrate",nrate ) ;
                        verbose=true)

end

max_stars_repo_name: wakakusa/RFCVGT

max_issues_repo_name: wakakusa/RFCVGT

max_forks_repo_name: wakakusa/RFCVGT


text:
module WoodburyFactorizations
using LinearAlgebra
using LinearAlgebra: checksquare

include("util.jl")

# IDEA: can optimize threshold constant c which makes Woodbury Identity computationally advantageous over direct factorization
using LazyInverses

const AbstractMatOrFac{T} = Union{AbstractMatrix{T}, Factorization{T}}
export Woodbury, WoodburyFactorization

# represents A + αUCV
# things that prevent C from being a scalar: checkdims, factorize, ...
# the α is beneficial to preserve p.s.d.-ness during inversion (see inverse)
struct WoodburyFactorization{T, AT, UT, CT, VT, F, L, TT} <: Factorization{T}
    A::AT
    U::UT
    C::CT
    V::VT
	α::F # scalar, TODO: test if not ±1
	logabsdet::L
	temporaries::TT
end
const Woodbury = WoodburyFactorization
function Woodbury(A::AbstractMatOrFac, U::AbstractMatOrFac,
				  C::AbstractMatOrFac, V::AbstractMatOrFac,
				  α::Real = 1, logabsdet::Union{NTuple{2, Real}, Nothing} = nothing,
				  temporaries = allocate_temporaries(U, V);
				  check::Bool = true)
	(α == 1 || α == -1) || throw(DomainError("α ≠ ±1 not yet tested: α = $α"))
	check && checkdims(A, U, C, V)
	T = promote_type(typeof(α), eltype.((A, U, C, V))...)	# check promote_type
	F = typeof(α)
	AT, UT, CT, VT, LT = typeof.((A, U, C, V, logabsdet))
	TT = typeof(temporaries)
	Woodbury{T, AT, UT, CT, VT, F, LT, TT}(A, U, C, V, α, logabsdet, temporaries)
end

function Woodbury(A::AbstractMatOrFac, U::AbstractVector,
				  C::Real, V::Adjoint{<:Any, <:AbstractVector}, α::Real = 1,
				  logabsdet::Union{NTuple{2, Real}, Nothing} = nothing,
				  temporaries = allocate_temporaries(U, V);
				  check::Bool = true)
	(α == 1 || α == -1) || throw(DomainError("α ≠ ±1 not yet tested: α = $α"))
	check && checkdims(A, U, C, V)
	T = promote_type(typeof(α), eltype.((A, U, C, V))...)	# check promote_type
	F = typeof(α)
	AT, UT, CT, VT, LT = typeof.((A, U, C, V, logabsdet))
	TT = typeof(temporaries)
	Woodbury{T, AT, UT, CT, VT, F, LT, TT}(A, U, C, V, α, logabsdet, temporaries)
end

const RankOneCorrection = Woodbury{<:Any, <:Any, <:AbstractVector, <:Any, <:Adjoint{<:Any, <:AbstractVector}}

# casts all inputs to an equivalent Matrix
matrix(x::Real) = fill(x, (1, 1)) # could in principle extend Matrix
matrix(x::AbstractVector) = reshape(x, (:, 1))
matrix(x::Adjoint{<:Any, <:AbstractVector}) = reshape(x, (1, :))
matrix(x::AbstractMatOrFac) = x

function Woodbury(A, U, C, V, α::Real = 1, abslogdet = nothing)
	Woodbury(matrix.((A, U, C, V))..., α, abslogdet)
end

# low rank correction
# NOTE: cannot rid of type restriction on A without introducing ambiguities
function Woodbury(A::AbstractMatOrFac, U::AbstractVector, V::Adjoint = U',
	 			  α::Real = 1, logabsdet = nothing)
	T = eltype(A)
	Woodbury(A, U, one(T), V, α, logabsdet)
end
function Woodbury(A::AbstractMatOrFac, U::AbstractMatrix, V::AbstractMatrix = U',
	 			  α::Real = 1, logabsdet = nothing)
	T = eltype(A)
	Woodbury(A, U, (one(T)*I)(size(U, 2)), V, α, logabsdet)
end
# function Woodbury(A, L::LowRank, α::Real = 1, logabsdet = nothing)
#     Woodbury(A, L.U, 1.0I(size(L.U, 2)), L.V, α, logabsdet)
# end
function Woodbury(A, C::CholeskyPivoted, α::Real = 1, logabsdet = nothing)
	U = C.U[1:C.rank, invperm(C.p)]
	Woodbury(A, U', U, α, logabsdet)
end
LinearAlgebra.checksquare(::Number) = 1 # since size(::Number, ::Int) = 1
# checks if the dimensions of A, U, C, V are consistent to form a Woodbury factorization
function checkdims(A, U, C, V)
	n, m = A isa Number ? (1, 1) : size(A)
	k, l = C isa Number ? (1, 1) : size(C)
	s = "is inconsistent with A ($(size(A))) and C ($(size(C)))"
	!(size(U, 1) == n && size(U, 2) == k) && throw(DimensionMismatch("Size of U ($(size(U))) "*s))
	!(size(V, 1) == l && size(V, 2) == m) && throw(DimensionMismatch("Size of V ($(size(V))) "*s))
	return true
end

# pseudo-constructor?
# c is used to set threshold for conversion to Matrix
function woodbury(A, U, C, V, α = 1, c::Real = 1)
    W = Woodbury(A, U, C, V, α)
	#  only return Woodbury if it is efficient
	if size(W.C, 1) ≥ c * size(W.U, 1) || size(W.C, 2) ≥ c * size(W.V, 2)
		W = Matrix(W)
	end
	return W
end

################################## Base ########################################
Base.size(W::Woodbury) = size(W.A)
Base.size(W::Woodbury, d) = size(W.A, d)
Base.eltype(W::Woodbury{T}) where {T} = T
function Base.AbstractMatrix(W::Woodbury)
	Matrix(W.A) + W.α * *(W.U, W.C, W.V)
end
Base.Matrix(W::Woodbury) = AbstractMatrix(W)
function Base.copy(W::Woodbury)
	U = copy(W.U)
	V = W.U ≡ W.V' ? U' : copy(W.V)
	t = tuple(copy.(W.temporaries)...)
	Woodbury(copy(W.A), U, copy(W.C), V, W.α, W.logabsdet, t)
end
function Base.deepcopy(W::Woodbury)
	U = deepcopy(W.U)
	V = W.U ≡ W.V' ? U' : deepcopy(W.V)
	t = tuple(deepcopy.(W.temporaries)...)
	Woodbury(deepcopy(W.A), U, deepcopy(W.C), V, W.α, W.logabsdet, t)
end
Base.inv(W::Woodbury) = inv(factorize(W))
function Base.isequal(W1::Woodbury, W2::Woodbury)
	W1.A == W2.A && W1.U == W2.U && W1.C == W2.C && W1.V == W2.V && W1.α == W2.α
end

# sufficient but not necessary condition, useful for quick check
_ishermitian(A::Union{Real, AbstractMatOrFac}) = ishermitian(A)
function _ishermitian(W::Woodbury)
	try
		(W.U ≡ W.V' || W.U == W.V') && _ishermitian(W.A) && _ishermitian(W.C)
	catch e
		return false
	end
end
# if efficiently-verifiable sufficient condition fails, check matrix
LinearAlgebra.ishermitian(W::Woodbury) = _ishermitian(W) || ishermitian(Matrix(W))
LinearAlgebra.issymmetric(W::Woodbury) = eltype(W) <: Real && ishermitian(W)

function LinearAlgebra.isposdef(W::Woodbury)
	W.logabsdet isa Nothing ? isposdef(factorize(W)) : (W.logabsdet[2] > 0)
end
function LinearAlgebra.adjoint(W::Woodbury)
	ishermitian(W) ? W : Woodbury(W.A', W.V', W.C', W.U', W.α, W.logabsdet)
end
function LinearAlgebra.transpose(W::Woodbury)
	issymmetric(W) ? W : Woodbury(transpose.((W.A, W.V, W.C, W.U))..., W.α, W.logabsdet)
end

# indexed by two integers returns scalar
function Base.getindex(W::Woodbury, i, j)
	Ui = i isa Integer ? W.U[i, :]' : W.U[i, :]
	W.A[i, j] + W.α * *(Ui, W.C, W.V[:, j])
end
# indexed by two vectors other, returns woodbury
function Base.getindex(W::Woodbury, i::AbstractVector, j::AbstractVector)
	A = W.A[i, j]
	U = W.U[i, :]
	V = W.U ≡ W.V' && i == j ? U' : W.V[:, j]
	return Woodbury(A, U, W.C, V, W.α, nothing)
end

# WARNING: result of view references existing temporaries and is therefore not thread-safe
function Base.view(W::Woodbury, i, j)
	A = view(W.A, i, j)
	U = view(W.U, i, :)
	V = W.U ≡ W.V' && i == j ? U' : view(W.V, :, j)
	return Woodbury(A, U, W.C, V, W.α, nothing, W.temporaries)
end
function LinearAlgebra.tr(W::Woodbury)
	n = checksquare(W)
	trace_UCV = zero(eltype(W))
	for i in 1:n
		ui = @view W.U[i, :]
		vi = @view W.V[:, i]
		trace_UCV += dot(ui, W.C, vi)
	end
	return tr(W.A) + W.α * trace_UCV
end

######################## Linear Algebra primitives #############################
# IDEA: have efficient methods of hermitian WoodburyFactorization
Base.:*(a::Number, W::Woodbury) = Woodbury(a * W.A, W.U, a * W.C, W.V, W.α) # IDEA: could take over logabsdet efficiently
Base.:*(W::Woodbury, a::Number) = a * W

Base.:\(a::Number, W::Woodbury) = W / a # IDEA: could take over logabsdet efficiently
Base.:/(W::Woodbury, a::Number) = Woodbury(W.A / a, W.U, W.C / a, W.V, W.α)

Base.:*(W::Woodbury, x::AbstractVecOrMat) = mul!(zero(x), W, x)
Base.:*(B::AbstractMatrix, W::Woodbury) = adjoint(W'*B')
function LinearAlgebra.mul!(y::AbstractVector, W::Woodbury, x::AbstractVector, α::Real = 1, β::Real = 0)
	s, t = get_temporaries(W, x)
	mul!!(y, W, x, α, β, s, t)
end
function LinearAlgebra.mul!(y::AbstractMatrix, W::Woodbury, x::AbstractMatrix, α::Real = 1, β::Real = 0)
	s, t = get_temporaries(W, x)
	mul!!(y, W, x, α, β, s, t) # Pre-allocate!
end
# allocates s, t arrays for multiplication if not pre-allocated in W
allocate_temporaries(W::Woodbury, T::DataType = Float64) = allocate_temporaries(W.U, W.V, T)
function allocate_temporaries(W::Woodbury, n::Int, T::DataType = Float64)
	allocate_temporaries(W.U, W.V, n, T)
end
function allocate_temporaries(U, V, T::DataType = Float64)
	allocate_temporaries(size(U, 2), size(V, 1), T)
end
function allocate_temporaries(U, V, n::Int, T::DataType = Float64)
	allocate_temporaries(size(U, 2), size(V, 1), n, T)
end
function allocate_temporaries(nu::Int, nv::Int, T::DataType = Float64)
	s = zeros(T, nv)
	t = zeros(T, nu)
	return s, t
end
function allocate_temporaries(nu::Int, nv::Int, n::Int, T::DataType = Float64)
	s = zeros(T, nv, n)
	t = zeros(T, nv, n)
	return s, t
end
function get_temporaries(W::Woodbury, x::AbstractVector)
	T = promote_type(eltype(x), eltype(W))
	W.temporaries isa Tuple{<:AbstractVector{T}} ? W.temporaries : allocate_temporaries(W, eltype(x))
end
function get_temporaries(W::Woodbury, X::AbstractMatrix)
	T = promote_type(eltype(X), eltype(W))
	W.temporaries isa Tuple{<:AbstractMatrix{T}} ? W.temporaries : allocate_temporaries(W, size(X, 2), eltype(X))
end

function mul!!(y::AbstractVecOrMat, W::Woodbury, x::AbstractVecOrMat, α::Real, β::Real, s, t)
	mul!(s, W.V, x)
	mul!(t, W.C, s)
	mul!(y, W.U, t, α * W.α, β)
 	mul!(y, W.A, x, α, 1) # this allocates if D is diagonal due to broadcasting mechanism
end

# special case: rank one correction does not need additional temporaries for MVM
function LinearAlgebra.mul!(y::AbstractVector, W::RankOneCorrection, x::AbstractVector, α::Real = 1, β::Real = 0)
	s = dot(W.V', x)
	t = W.C * s
	@. y = (α * W.α) * W.U * t + β * y
 	mul!(y, W.A, x, α, 1)
end

# ternary dot
function LinearAlgebra.dot(x::AbstractVecOrMat, W::Woodbury, y::AbstractVector)
    Ux = W.U'x     # memory allocation can be avoided (with lazy arrays?)
    Vy = (x ≡ y && W.U ≡ W.V') ? Ux : W.V*y
    dot(x, W.A, y) + W.α * dot(Ux, W.C, Vy)
end

# ternary mul
function Base.:*(x::AbstractMatrix, W::Woodbury, y::AbstractVecOrMat)
    xU = x*W.U
    Vy = (x ≡ y' && W.U ≡ W.V') ? xU' : W.V*y
    *(x, W.A, y) + W.α * *(xU, W.C, Vy) # can avoid two temporaries
end

Base.:(\)(W::Woodbury, B::AbstractVector) = factorize(W)\B
Base.:(\)(W::Woodbury, B::AbstractMatrix) = factorize(W)\B
Base.:(/)(B::AbstractMatrix, W::Woodbury) = B/factorize(W)

function LinearAlgebra.diag(W::Woodbury)
    n = checksquare(W)
    d = zeros(eltype(W), n)
    for i in 1:n
        d[i] = W[i, i]
    end
    return d
end

########################## Matrix inversion lemma ##############################
# figure out constant c for which woodbury is most efficient
# could implement this in WoodburyMatrices and PR
function LinearAlgebra.factorize(W::Woodbury, c::Real = 1,
									compute_logdet::Val{T} = Val(true)) where T
	checksquare(W)
	if size(W.U, 1) < c*size(W.U, 2)
		return factorize(AbstractMatrix(W))
	else # only use Woodbury identity when it is beneficial to do so
		A = factorize(W.A)
		A⁻¹ = inverse(A)
		A⁻¹U = A⁻¹ * W.U
		VA⁻¹ = (W.U ≡ W.V' && ishermitian(A⁻¹)) ? A⁻¹U' : W.V * A⁻¹
		D = factorize_D(W, W.V*A⁻¹U)
		if T
			l, s = _logabsdet(A, W.C, D, W.α)
			W_logabsdet = (-l, s) # -l since the result is packaged in an inverse
		else
			W_logabsdet = nothing
		end
		α = -W.α # switch sign
		return inverse(Woodbury(A⁻¹, A⁻¹U, inverse(D), VA⁻¹, α, W_logabsdet))
	end
end

##################### conveniences for D = C⁻¹ ± V*A⁻¹*U ########################
compute_D(W::Woodbury) = compute_D!(W, *(W.V, inverse(W.A), W.U))
function compute_D!(W::Woodbury, VAU)
	invC = inv(W.C) # because we need the dense inverse matrix
	@. VAU = invC + W.α * VAU # could be made more efficient with 5 arg mul
end
compute_D!(W::Woodbury, VAU::Real) = inv(W.C) + W.α * VAU
factorize_D(W::Woodbury) = factorize_D(W, *(W.V, inverse(W.A), W.U))
function factorize_D(W::Woodbury, VAU)
	D = compute_D!(W, VAU)
	if size(D) == (1, 1)
		return D
	elseif _ishermitian(W)
		try return cholesky(Hermitian(D)) catch end
	end
	return factorize(D)
end

######################## Matrix determinant lemma ##############################
# if W.A = W.C = I, this is Sylvesters determinant theorem
# Determinant lemma for A + α*(UCV)
function LinearAlgebra.det(W::Woodbury)
	l, s = logabsdet(W)
	exp(l) * s
end
function LinearAlgebra.logdet(W::Woodbury)
	l, s = logabsdet(W)
	return s > 0 ? l : throw(DomainError("determinant negative: $(exp(l) * s)"))
end
function LinearAlgebra.logabsdet(W::Woodbury)
	W.logabsdet == nothing ? logabsdet(factorize(W)) : W.logabsdet
end

function _logabsdet(A, C, D, α::Real)
	n, m = checksquare(A), checksquare(D)
	la, sa = logabsdet(A)
	lc, sc = logabsdet(C)
	ld, sd = logabsdet(D)
	sα = (α == -1) && isodd(m) ? -1 : 1
	return +(la, lc, ld), *(sa, sc, sd, sα)
end

end # WoodburyFactorizations


text:
mutable struct GlobalApproximationValueIterationSolver{GFA <: GlobalFunctionApproximator, RNG <: AbstractRNG, F} <: Solver
    gfa::GFA
    num_samples::Int64
    belres::Float64
    num_iterations::Int64
    verbose::Bool
    rng::RNG
    is_mdp_generative::Bool
    n_generative_samples::Int64
    fv_type::F
end

function GlobalApproximationValueIterationSolver(gfa::GFA; num_samples::Int64, belres::Float64=1e-3,
                                                 max_iterations::Int64=1000, verbose::Bool=false,
                                                 rng::RNG=Random.GLOBAL_RNG, is_mdp_generative::Bool=false,
                                                 n_generative_samples::Int64=0, fv_type::F=Vector{Float64}) where {GFA <: GlobalFunctionApproximator, F, RNG <: AbstractRNG}
    return GlobalApproximationValueIterationSolver(gfa, num_samples, belres, max_iterations, verbose, rng, is_mdp_generative, n_generative_samples, fv_type)
end

function GlobalApproximationValueIterationSolver()
    throw(ArgumentError("GlobalApproximationValueIterationSolver needs a GlobalFunctionApproximator object for construction!"))
end

mutable struct GlobalApproximationValueIterationPolicy{GFA <: GlobalFunctionApproximator, RNG <: AbstractRNG, F} <: Policy
    gfa::GFA
    action_map::Vector
    mdp::Union{MDP,POMDP}
    is_mdp_generative::Bool
    n_generative_samples::Int64
    fv_type::F
    rng::RNG
end

function GlobalApproximationValueIterationPolicy(mdp::Union{MDP,POMDP},
                                                 solver::GlobalApproximationValueIterationSolver)
    return GlobalApproximationValueIterationPolicy(deepcopy(solver.gfa), ordered_actions(mdp), mdp,
                                                   solver.is_mdp_generative, solver.n_generative_samples, solver.fv_type, solver.rng)
end

function convert_featurevector(::Type{V}, s::S, mdp::Union{MDP,POMDP}) where {V <: AbstractArray, S}
    return convert_s(V, s, mdp)
end

# If global function approximator is non-linear, a default convert_s is required
@POMDP_require convert_featurevector(t::Type{V} where {V <: AbstractArray}, s::S where S, mdp::Union{MDP,POMDP}) begin
    @req convert_s(::Type{V} where {V <: AbstractArray}, ::S, ::typeof(mdp))
end

@POMDP_require solve(solver::GlobalApproximationValueIterationSolver, mdp::Union{MDP,POMDP}) begin

    P = typeof(mdp)
    S = statetype(P)
    A = actiontype(P)
    @req discount(::P)
    @subreq ordered_actions(mdp)

    @req actionindex(::P, ::A)
    @req actions(::P, ::S)
    as = actions(mdp)
    a = first(as)

    # Need to be able to sample states
    @req sample_state(::P, ::typeof(solver.rng))

    # Have different requirements depending on whether solver MDP is generative or explicit
    if solver.is_mdp_generative
        @req gen(::P, ::S, ::A, ::typeof(solver.rng))
    else
        @req transition(::P, ::S, ::A)
        ss = sample_state(mdp, solver.rng)
        dist = transition(mdp, ss, a)
        D = typeof(dist)
        @req support(::D)
    end

    # Feature vector conversion must be defined either directly or by default (through convert_s)
    @subreq convert_featurevector(solver.fv_type, sample_state(mdp, solver.rng), mdp)
end

function POMDPs.solve(solver::GlobalApproximationValueIterationSolver, mdp::Union{MDP,POMDP})

    @warn_requirements solve(solver,mdp)

    # Ensure that generative model has a non-zero number of samples
    if solver.is_mdp_generative
        @assert solver.n_generative_samples > 0
    end

    # Solver parameters
    belres = solver.belres
    num_iterations = solver.num_iterations
    num_samples = solver.num_samples
    discount_factor = discount(mdp)
    gfa_type = typeof(solver.gfa)

    # Initialize the policy
    policy = GlobalApproximationValueIterationPolicy(mdp, solver)

    total_time = 0.0
    iter_time = 0.0

    temp_s = sample_state(mdp, solver.rng)
    state_dim = length(convert_featurevector(solver.fv_type, temp_s, mdp))

    for iter = 1:num_iterations

        residual = 0.0

        # Setup input and outputs for fit functions
        state_matrix = zeros(eltype(solver.fv_type), num_samples, state_dim)
        val_vector = zeros(eltype(solver.fv_type), num_samples)

        iter_time = @elapsed begin

        # Loop over the chosen number of samples for approximation
        # and compute the current value function estimate at each sample
        for i = 1:num_samples

            s = sample_state(mdp, solver.rng)

            pt = convert_featurevector(solver.fv_type, s, mdp)
            state_matrix[i,:] = pt

            sub_aspace = actions(mdp,s)

            if isterminal(mdp, s)
                val_vector[i] = 0.0
            else
                old_util = value(policy, s)
                max_util = -Inf

                # Compute the approximate Q value for each action and choose the best
                for a in sub_aspace
                    iaction = actionindex(mdp,a)
                    u = 0.0

                    if solver.is_mdp_generative
                        for j in 1:solver.n_generative_samples
                            sp, r = @gen(:sp,:r)(mdp, s, a, solver.rng)
                            u += r

                            if !isterminal(mdp,sp)
                                sp_feature = convert_featurevector(solver.fv_type, sp, mdp)
                                u += discount_factor*compute_value(policy.gfa, sp_feature)
                            end
                        end
                        u = u / solver.n_generative_samples
                    else
                        dist = transition(mdp,s,a)
                        for (sp, p) in weighted_iterator(dist)
                            p == 0.0 ? continue : nothing
                            r = reward(mdp, s, a, sp)
                            u += p*r

                            # Only interpolate sp if it is non-terminal
                            if !isterminal(mdp,sp)
                                sp_feature = convert_featurevector(solver.fv_type, sp, mdp)
                                u += p * (discount_factor*compute_value(policy.gfa, sp_feature))
                            end
                        end
                    end

                    max_util = (u > max_util) ? u : max_util
                    util_diff = abs(max_util - old_util)
                    util_diff > residual ? (residual = util_diff) : nothing
                end #action

                val_vector[i] = max_util
            end
        end
        # Now fit!
        fit!(policy.gfa, state_matrix, val_vector)

        end # time

        total_time += iter_time
        solver.verbose ? @printf("[Iteration %-4d] residual: %10.3G | iteration runtime: %10.3f ms, (%10.3G s total)\n", iter, residual, iter_time*1000.0, total_time) : nothing
        residual < belres ? break : nothing
    end

    return policy
end

function POMDPs.value(policy::GlobalApproximationValueIterationPolicy, s::S) where S

    s_point = convert_featurevector(policy.fv_type, s, policy.mdp)
    val = compute_value(policy.gfa, s_point)
    return val
end

# Not explicitly stored in policy - extract from value function interpolation
function POMDPs.action(policy::GlobalApproximationValueIterationPolicy, s::S) where S

    mdp = policy.mdp
    best_a_idx = -1
    max_util = -Inf
    sub_aspace = actions(mdp,s)
    discount_factor = discount(mdp)

    for a in sub_aspace

        iaction = actionindex(mdp, a)
        u = value(policy,s,a)

        if u > max_util
            max_util = u
            best_a_idx = iaction
        end
    end

    return policy.action_map[best_a_idx]
end

function POMDPs.value(policy::GlobalApproximationValueIterationPolicy, s::S, a::A) where {S,A}

    mdp = policy.mdp
    discount_factor = discount(mdp)
    u = 0.0

    # As in solve(), do different things based on whether
    # mdp is generative or explicit
    if policy.is_mdp_generative
        for j in 1:policy.n_generative_samples
            sp, r = @gen(:sp,:r)(mdp, s, a, solver.rng)
            sp_point = convert_featurevector(policy.fv_type, sp, mdp)
            u += r + discount_factor*compute_value(policy.gfa, sp_point)
        end
        u = u / policy.n_generative_samples
    else
        dist = transition(mdp,s,a)
        for (sp, p) in weighted_iterator(dist)
            p == 0.0 ? continue : nothing
            r = reward(mdp, s, a, sp)
            u += p*r

            # Only interpolate sp if it is non-terminal
            if !isterminal(mdp,sp)
                sp_point = convert_featurevector(policy.fv_type, sp, mdp)
                u += p*(discount_factor*compute_value(policy.gfa, sp_point))
            end
        end
    end

    return u
end


text:
using ..Atmos
using ..Atmos: AbstractMoistureModel

# Helpers to gather the thermodynamic variables across the DG grid

function vars_thermo(atmos::AtmosModel, FT)
    @vars begin
        temp::FT
        pres::FT
        θ_dry::FT
        e_int::FT
        h_tot::FT
        h_int::FT

        moisture::vars_thermo(moisture_model(atmos), FT)
    end
end
vars_thermo(::AbstractMoistureModel, FT) = @vars()
function vars_thermo(m::Union{EquilMoist, NonEquilMoist}, FT)
    @vars begin
        q_liq::FT
        q_ice::FT
        q_vap::FT
        θ_vir::FT
        θ_liq_ice::FT
        has_condensate::Bool
    end
end
num_thermo(bl, FT) = varsize(vars_thermo(bl, FT))
thermo_vars(bl, array) = Vars{vars_thermo(bl, eltype(array))}(array)

# compute thermodynamic variables visitor function
function compute_thermo!(atmos::AtmosModel, state, aux, thermo)
    e_tot = state.energy.ρe / state.ρ
    ts = recover_thermo_state(atmos, state, aux)
    e_int = internal_energy(ts)

    thermo.temp = air_temperature(ts)
    thermo.pres = air_pressure(ts)
    thermo.θ_dry = dry_pottemp(ts)
    thermo.e_int = e_int

    thermo.h_tot = total_specific_enthalpy(ts, e_tot)
    thermo.h_int = specific_enthalpy(ts)

    compute_thermo!(moisture_model(atmos), state, aux, ts, thermo)

    return nothing
end
function compute_thermo!(::AbstractMoistureModel, state, aux, ts, thermo)
    return nothing
end
function compute_thermo!(
    moist::Union{EquilMoist, NonEquilMoist},
    state,
    aux,
    ts,
    thermo,
)
    thermo.moisture.has_condensate = has_condensate(ts)
    thermo.moisture.q_liq = liquid_specific_humidity(ts)
    thermo.moisture.q_ice = ice_specific_humidity(ts)
    thermo.moisture.q_vap = vapor_specific_humidity(ts)
    thermo.moisture.θ_vir = virtual_pottemp(ts)
    thermo.moisture.θ_liq_ice = liquid_ice_pottemp(ts)

    return nothing
end


text:
function test_yLP()

    model = Chemostat.Test.toy_model()
    glcidx = Chemostat.Utils.rxnindex(model, "gt")
    biomidx = Chemostat.Utils.rxnindex(model, "biom")

    d = zeros(length(model.c)); d[glcidx] = 1.0
    model.c[biomidx] = 1.0

    status, yflxm, yield = Chemostat.LP.yLP(model, d);

    @test all(isapprox.(model.S * yflxm, model.b; atol = 1e-4))
    @test yield > 0.0
    @test status == Chemostat.LP.JuMP.MOI.OPTIMAL
end
test_yLP()


text:
using JuLIP, PyPlot, PyCall

function myplot(at, ttl = nothing)
    x, y, _ = xyz(at)
    r = at.M
    N = length(at.Z)
    color = Array{String}(N)
    for i = 1:N
        z = at.Z[i];
        if z == 1
            color[i] = "blue";
        else
            color[i] = "red";
        end
    #    PyPlot.plot(mod(x[i],1), mod(y[i],1), color, markersize=20) #set to be mod(x, Lx), mod(y, Ly)
    end
    if ttl != nothing
        title(ttl)
    end

    cell_x = [0.0, 1.0, 1.0, 0.0, 0.0]
    cell_y = [0.0, 0.0, 1.0, 1.0, 0.0]

    shape = [circle(x0=mod(x[i],1)-r[i], y0=mod(y[i],1)-r[i],
                    x1=mod(x[i],1)+r[i], y1=mod(y[i],1)+r[i];
                    opacity = 0.8, fillcolor=color[i], line_color="transparent")
             for i in 1:N]

    PlotlyJS.plot(cell_x, cell_y, Layout(;height=600,width=600,showlegend=false,shapes=shape))

end


text:
"""
Variational partitioned Runge-Kutta integrator cache.

### Fields

* `n`: time step number
* `t`: time of current time step
* `t̅`: time of previous time step
* `q`: current solution of q
* `q̅`: previous solution of q
* `p`: current solution of p
* `p̅`: previous solution of p
* `v`: vector field of q
* `v̅`: vector field of q̅
* `f`: vector field of p
* `f̅`: vector field of p̅
* `q̃`: initial guess of q
* `p̃`: initial guess of p
* `ṽ`: initial guess of v
* `f̃`: initial guess of f
* `s̃`: holds shift due to periodicity of solution
* `Q`: internal stages of q
* `P`: internal stages of p
* `V`: internal stages of v
* `F`: internal stages of f
* `Y`: integral of vector field of internal stages of q
* `Z`: integral of vector field of internal stages of p
"""
mutable struct IntegratorCacheVPRK{ST,TT,D,S} <: IODEIntegratorCache{ST,D}
    n::Int
    t::TT
    t̅::TT

    q::Vector{ST}
    q̅::Vector{ST}
    p::Vector{ST}
    p̅::Vector{ST}
    λ::Vector{ST}
    λ̅::Vector{ST}

    qₑᵣᵣ::Vector{ST}
    pₑᵣᵣ::Vector{ST}

    q₋::Vector{ST}
    q̅₊::Vector{ST}
    p₋::Vector{ST}
    p̅₊::Vector{ST}

    v::Vector{ST}
    v̅::Vector{ST}
    f::Vector{ST}
    f̅::Vector{ST}

    u::Vector{ST}
    g::Vector{ST}

    q̃::Vector{ST}
    p̃::Vector{ST}
    ṽ::Vector{ST}
    f̃::Vector{ST}
    s̃::Vector{ST}

    Q::Vector{Vector{ST}}
    P::Vector{Vector{ST}}
    V::Vector{Vector{ST}}
    F::Vector{Vector{ST}}
    Λ::Vector{Vector{ST}}
    Φ::Vector{Vector{ST}}

    Y::Vector{Vector{ST}}
    Z::Vector{Vector{ST}}
    U::Vector{Vector{ST}}
    G::Vector{Vector{ST}}
    R::Vector{Vector{ST}}

    function IntegratorCacheVPRK{ST,TT,D,S}(projection::Bool=false) where {ST,TT,D,S}
        # create solution vectors
        q = zeros(ST,D)
        q̅ = zeros(ST,D)
        p = zeros(ST,D)
        p̅ = zeros(ST,D)

        # create error vectors
        qₑᵣᵣ = zeros(ST,D)
        pₑᵣᵣ = zeros(ST,D)

        # create update vectors
        v = zeros(ST,D)
        v̅ = zeros(ST,D)
        f = zeros(ST,D)
        f̅ = zeros(ST,D)

        # create temporary vectors
        q̃ = zeros(ST,D)
        p̃ = zeros(ST,D)
        ṽ = zeros(ST,D)
        f̃ = zeros(ST,D)
        s̃ = zeros(ST,D)

        # create internal stage vectors
        Q = create_internal_stage_vector(ST, D, S)
        P = create_internal_stage_vector(ST, D, S)
        V = create_internal_stage_vector(ST, D, S)
        F = create_internal_stage_vector(ST, D, S)
        Y = create_internal_stage_vector(ST, D, S)
        Z = create_internal_stage_vector(ST, D, S)

        # projection vectors
        if projection
            λ = zeros(ST,D)
            λ̅ = zeros(ST,D)

            q₋= zeros(ST,D)
            q̅₊= zeros(ST,D)
            p₋= zeros(ST,D)
            p̅₊= zeros(ST,D)

            u = zeros(ST,D)
            g = zeros(ST,D)

            Λ = create_internal_stage_vector(ST, D, S)
            Φ = create_internal_stage_vector(ST, D, S)

            U = create_internal_stage_vector(ST, D, 2)
            G = create_internal_stage_vector(ST, D, 2)
            R = create_internal_stage_vector(ST, D, S)
        else
            λ = Vector{ST}()
            λ̅ = Vector{ST}()

            q₋= Vector{ST}()
            q̅₊= Vector{ST}()
            p₋= Vector{ST}()
            p̅₊= Vector{ST}()

            u = Vector{ST}()
            g = Vector{ST}()

            Λ = create_internal_stage_vector(ST, 0, 0)
            Φ = create_internal_stage_vector(ST, 0, 0)

            U = create_internal_stage_vector(ST, 0, 0)
            G = create_internal_stage_vector(ST, 0, 0)
            R = create_internal_stage_vector(ST, 0, 0)
        end

        new(0, zero(TT), zero(TT),
            q, q̅, p, p̅, λ, λ̅,
            qₑᵣᵣ, pₑᵣᵣ,
            q₋, q̅₊, p₋, p̅₊,
            v, v̅, f, f̅, u, g, q̃, p̃, ṽ, f̃, s̃,
            Q, P, V, F, Λ, Φ, Y, Z, U, G, R)
    end
end

function IntegratorCacheVPRK(ST,TT,D,S)
    IntegratorCacheVPRK{ST,TT,D,S}(false)
end

function IntegratorCacheVPRKwProjection(ST,TT,D,S)
    IntegratorCacheVPRK{ST,TT,D,S}(true)
end

function update_params!(params::AbstractParametersVPRK, cache::IntegratorCacheVPRK)
    # set time for nonlinear solver and copy previous solution
    params.t̅  = cache.t
    params.q̅ .= cache.q
    params.p̅ .= cache.p
end

function update_solution!(int::AbstractIntegratorVPRK{DT,TT}, cache::IntegratorCacheVPRK{DT,TT}) where {DT,TT}
    update_solution!(cache.q, cache.qₑᵣᵣ, cache.V, tableau(int).q.b, tableau(int).q.b̂, timestep(int))
    update_solution!(cache.p, cache.pₑᵣᵣ, cache.F, tableau(int).p.b, tableau(int).p.b̂, timestep(int))
end

function project_solution!(int::AbstractIntegratorVPRK{DT,TT}, cache::IntegratorCacheVPRK, R::Vector{TT}) where {DT,TT}
    update_solution!(cache.q, cache.qₑᵣᵣ, cache.U, R, timestep(int))
    update_solution!(cache.p, cache.pₑᵣᵣ, cache.G, R, timestep(int))
end

function project_solution!(int::AbstractIntegratorVPRK{DT,TT}, cache::IntegratorCacheVPRK, RU::Vector{TT}, RG::Vector{TT}) where {DT,TT}
    update_solution!(cache.q, cache.qₑᵣᵣ, cache.U, RU, timestep(int))
    update_solution!(cache.p, cache.pₑᵣᵣ, cache.G, RG, timestep(int))
end

function CommonFunctions.set_solution!(cache::IntegratorCacheVPRK, sol, n=0)
    t, q, p = sol
    cache.n  = n
    cache.t  = t
    cache.q .= q
    cache.p .= p
    cache.v .= 0
    cache.f .= 0
end

function create_integrator_cache(int::AbstractIntegratorVPRK{DT,TT}) where {DT,TT}
    IntegratorCacheVPRK(DT, TT, ndims(int), nstages(int))
end

function create_integrator_cache(int::AbstractIntegratorVPRKwProjection{DT,TT}) where {DT,TT}
    IntegratorCacheVPRKwProjection(DT, TT, ndims(int), nstages(int))
end

function initialize!(int::AbstractIntegratorVPRK{DT}, cache::IntegratorCacheVPRK) where {DT,TT}
    cache.t̅ = cache.t - timestep(int)

    equation(int).v(cache.t, cache.q, cache.p, cache.v)
    equation(int).f(cache.t, cache.q, cache.p, cache.f)

    initialize!(int.iguess, cache.t, cache.q, cache.p, cache.v, cache.f,
                            cache.t̅, cache.q̅, cache.p̅, cache.v̅, cache.f̅)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

# Collection of sources required to build LASzipBuilder
sources = [
    "https://github.com/LASzip/LASzip.git" =>
    "002c42593a3a90609e12805e52350cfa4c6f6563",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
mkdir LASzip/build
cd LASzip/build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_VERBOSE_MAKEFILE=OFF ..
cmake --build . --target install --config Release

if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
mv $prefix/bin/liblaszip3.dll $prefix/bin/liblaszip.dll

fi
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    BinaryProvider.Linux(:i686, :glibc),
    BinaryProvider.Linux(:x86_64, :glibc),
    BinaryProvider.MacOS(),
    BinaryProvider.Windows(:i686),
    BinaryProvider.Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "liblaszip", :liblaszip),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "LASzipBuilder", sources, script, platforms, products, dependencies)

max_stars_repo_name: evetion/LASzipBuilder

max_issues_repo_name: evetion/LASzipBuilder

max_forks_repo_name: evetion/LASzipBuilder


text:
using PolygonOps
using Test
using StaticArrays

edge(p1, p2) = (SVector(p1),SVector(p2))

@testset "inpolygon" begin
    poly1 = SVector{2,Int}[(0,0),(0,10),(10,10),(10,0),(0,0)]

    @testset "Hao Sun" begin
        algo = HaoSun()
        for poly in (poly1, reverse(poly1))
            @test inpolygon(SVector(5,5),poly, algo) == 1 #in
            @test inpolygon(SVector(-5,-5),poly, algo) == 0 #out
            @test inpolygon(SVector(0,5),poly, algo) == -1 #on
            @test inpolygon(SVector(5,10),poly, algo) == -1 #on
            @test inpolygon(SVector(10,10),poly, algo) == -1 #on
            @test inpolygon(SVector(0,0),poly, algo) == -1 #on

            @test inpolygon(SVector(5,5),poly, algo, in=true,on=true,out=false) #in
            @test !inpolygon(SVector(-5,-5),poly, algo, in=true,on=true,out=false) #out
            @test inpolygon(SVector(0,5),poly, algo, in=true,on=true,out=false) #on
            @test inpolygon(SVector(5,10),poly, algo, in=true,on=true,out=false) #on
            @test inpolygon(SVector(10,10),poly, algo, in=true,on=true,out=false) #on
            @test inpolygon(SVector(0,0),poly, algo, in=true,on=true,out=false) #on
        end
    end
    @testset "Hormann Agathos" begin
        algo = HormannAgathos()

        poly = reverse(poly1)
        @test inpolygon(SVector(5,5),poly, algo) == 1 #in
        @test inpolygon(SVector(-5,-5),poly, algo) == 0 #out
        #@test_broken inpolygon(SVector(0,5),poly, algo) == -1 #on
        @test inpolygon(SVector(5,10),poly, algo) == -1 #on
        @test inpolygon(SVector(10,10),poly, algo) == -1 #on
        @test inpolygon(SVector(0,0),poly, algo) == -1 #on
    end
    @testset "API" begin
        poly = poly1
        @test inpolygon(SVector(5,5),poly, in=true,on=true,out=false) #in
        @test !inpolygon(SVector(-5,-5),poly, in=true,on=true,out=false) #out
        @test inpolygon(SVector(0,5),poly, in=true,on=true,out=false) #on
        @test inpolygon(SVector(5,10),poly, in=true,on=true,out=false) #on
        @test inpolygon(SVector(10,10),poly, in=true,on=true,out=false) #on
        @test inpolygon(SVector(0,0),poly, in=0x0,on=0x0,out=0x1) === 0x0 #on
    end
end

@testset "area/centroid" begin
    poly1 = SVector{2,Int}[(0,0),(0,10),(10,10),(10,0),(0,0)]
    poly2 = poly1.*5

    @testset "area" begin
        @test PolygonOps.area(poly1) == -100
        @test PolygonOps.area(reverse(poly1)) == 100
        @test PolygonOps.area(poly2) == -2500
        @test PolygonOps.area(reverse(poly2)) == 2500
    end

    @testset "centroid" begin
        @test PolygonOps.centroid(poly1) == SVector(5,5)
        @test PolygonOps.centroid(reverse(poly1)) == SVector(5,5)
        @test PolygonOps.centroid(poly2) == SVector(25,25)
        @test PolygonOps.centroid(reverse(poly2)) == SVector(25,25)
    end
end


text:
### VanillaMCTuner

# VanillaMCTuner is a dummy tuner type in the sense that it does not perform any tuning
# It is used only for determining whether the MCSampler will be verbose

immutable VanillaMCTuner <: MCTuner
  period::Integer # Tuning period over which acceptance rate is computed
  verbose::Bool # If verbose=false then the tuner is silent, else it is verbose

  function VanillaMCTuner(period::Integer, verbose::Bool)
    @assert period > 0 "Adaptation period should be positive"
    new(period, verbose)
  end
end

VanillaMCTuner(; period::Integer=100, verbose::Bool=false) = VanillaMCTuner(period, verbose)

Base.show(io::IO, tuner::VanillaMCTuner) = print(io, "VanillaMCTuner: period = $(tuner.period), verbose = $(tuner.verbose)")

Base.writemime(io::IO, ::MIME"text/plain", tuner::VanillaMCTuner) = show(io, tuner)


text:
export @MixedBag, eltypes, AbstractMixedBag

abstract AbstractMixedBag

Base.eachindex(iter::SimpleVector) = 1:length(iter)
function firstindex(condition, iter)
    for i in eachindex(iter)
        condition(iter[i]) && return i
    end
    throw(error("No index satisfying condition=$condition in iter=$iter"))
end

function fieldname(T, types)
    for i in eachindex(types)
        (types[i] == T) && return Symbol("_$i")
    end
    throw(error("No field T=$T in types=$types"))
end

function mixedbag_fields(types)
    fields = []
    for T in types
        fname = fieldname(T, types)
        field = :($(fname) :: Vector{$T})
        push!(fields, field)
    end
    fields
end

function typedef_mixedbag(name::Symbol, types)
    mutablility = false
    fields = mixedbag_fields(types)
    fieldblock = Expr(:block,fields...)
    Expr(:type, mutablility, :($name <: MixedBags.AbstractMixedBag), fieldblock)
end

eltypes{Bag <: AbstractMixedBag}(::Type{Bag}) = map(eltype, Bag.types)
eltypes(b::AbstractMixedBag) = eltypes(typeof(b))

# constructors
function empty_constructor_impl(name, types)
    args = [:($T[]) for T in types]
    rhs = Expr(:call, name, args...)
    lhs = Expr(:call, name, )
    Expr(Symbol("="), lhs, rhs)
end

# copy constructor
(::Type{Bag}){Bag <: AbstractMixedBag}(b::Bag) = b

function (::Type{Bag}){Bag <: AbstractMixedBag}(args...)
    b = Bag()
    for arg in args
        push!(b, arg)
    end
    b
end

macro MixedBag(name, types...)
    Expr(:block,
        typedef_mixedbag(name, types),
        empty_constructor_impl(name, types)
    ) |> esc
end


text:
"""
    checkfeasible(command::String, logicset::LogicalCombo; verbose=true, all=false, countany=false)

Is called when the user would like to check if a command produces a valid result,
possible result, or invalid result. The result is returned as a decimal from 0.0
to 1.0. With 0.0 being no matches and 1.0 being all matches.

### Arguments
* `verbose` : controls print
* `countall` : all sets have to be feasible or return 0
* `countany` : any set can be non-zero to return 1

### Examples
```julia
julia> myset = logicalparse("a, b, c ∈ red, blue, green")
a, b, c ∈ red, blue, green       feasible outcomes 27 ✓          :red blue blue

julia> myset = logicalparse("a != b,c; b = c ==> a = 'blue'", myset)
a != b,c                 feasible outcomes 12 ✓          :green blue blue
b = c ==> a = 'blue'     feasible outcomes 8 ✓           :blue green green

julia> checkfeasible("a = 'green' ==> b = 'red'", myset)
Check: a = 'green' ==> b = 'red' ... a = 'green' ==> b = 'red'   feasible outcomes 17 ✓          :blue red green
possible,  17 out of 21 possible combinations 'true'.
```
"""
function checkfeasible(command::String,
    logicset::LogicalCombo = LogicalCombo();
    verbose=true, countall=false, countany=false)

  clearcounter()

  countany && countall && throw("Both any and all can't be set to true")

  rowsin = sum(logicset.logical)

  if rowsin == 0
      println("No outcomes feasible with set!")
      return missing
  end

  verbose && !countall && !countany && print("check: $command ... ")
  verbose && countall && print("all: $command ... ")
  verbose && countany && print("any: $command ... ")

  logicsetout = logicalparse(command, logicset=logicset, verbose=verbose)

  rowsout = sum(logicsetout.logical)
  outcomeratio = rowsout/rowsin

  if countall
        outcomeratio != 1 && print("false")
        outcomeratio == 1 && print("true")
  elseif countany
        outcomeratio == 0 && print("false")
        outcomeratio != 0 && print("true")
  else
      outcomeratio == 1 && print("true")
      outcomeratio == 0 && print("false")
      (1 > outcomeratio > 0) && print("possible")
  end

  verbose && println(", $rowsout out of $rowsin possible combinations 'true'.")

  return [outcomeratio, logicsetout]
end


text:
# # In-place transforms

using FFTW
using MPI
using PencilFFTs

using Random: randn!

const INPUT_DIMS = (64, 32, 64)

MPI.Init()

dims = INPUT_DIMS

# Combine r2r and c2c in-place transforms.
transforms = (
    Transforms.R2R!(FFTW.REDFT01),
    Transforms.FFT!(),
    Transforms.R2R!(FFTW.DHT),
)

# MPI topology information
comm = MPI.COMM_WORLD
Nproc = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

# Let's do a 1D decomposition.
proc_dims = (Nproc, )

# Create in-place plan
plan = PencilFFTPlan(dims, transforms, proc_dims, comm)
rank == 0 && println(plan)
@assert Transforms.is_inplace(plan)

# Allocate data for the plan.
# This returns a `ManyPencilArray` container that holds multiple
# `PencilArray` views.
A = allocate_input(plan) :: PencilArrays.ManyPencilArray

# The input output `PencilArray`s are recovered using `first` and
# `last`.
u_in = first(A) :: PencilArray
u_out = last(A) :: PencilArray

# Initialise input data.
randn!(u_in)

# Apply in-place forward transform on the `ManyPencilArray` container.
plan * A

# After the transform, operations should be performed on the output view
# `u_out`. For instance, let's compute the global sum of the transformed data.
sum_global = sum(u_out)  # note that `sum` reduces over all processes

# Apply in-place backward transform.
plan \ A

# Now we can again perform operations on the input view `u_in`...


text:
cigar(x) = x[1]^2 + 1 * sum(abs2, view(x, 2:length(x)))   # 1e6

cigtab(x) = x[1]^2 + 1 * x[end]^2 + 1e4 * sum(abs2, view(x, 2:length(x)-1)) # 1e8

function griewank(x)
  n = length(x)
  1 + (1/4000)*sum(abs2, (x.+7)) - prod(cos.((x.+7) ./ sqrt.(1:n)))
end

function rosenbrock2d(x)
  return sum((x).*(x).*cos.((x)).*sin.((x.-5)))+200
end

function quartic(x)
  D = length(x)
  sum( (1:D) .* (x.-2).^4 )
end

function schwefel1_2(x)
  D = length(x)
  partsums = zeros(D)
  partsum = 0
  for i in 1:D
    partsum += x[i]-9
    partsums[i] = partsum
  end
  sum(abs2, partsums)
end

function ellipsoid(x)
    res = 0.0
    cumsum = 0.0
    for xx in x
        cumsum += xx-sqrt(2)
        res += cumsum^2
    end
    res
end

function elliptic(x)
  D = length(x)
  condition = 1e+3
  if D==1
      coefficients=1
  else
      coefficients = condition .^ range(0,stop=1,length=D)
  end
  sum(coefficients .* (x.+3/2).^2)
end

function rastrigin(x)
  D = length(x)
  10 * D + sum(abs2, (x.+0.7)) - 10 * sum(xx -> cos.(2π * xx), (x.+0.7))
end

function sphere(x)
  sum(abs2, (x.-1.3))
end

function sin_01_x_2(x)
  sum(sin.((x.+0.7))+0.01.*(x.+0.7).^2)
end

function x_i(x)
  ress=0
  for i=1:length(x)
    ress+=(x[i]-i-2.1)^2
  end
  return ress
end

function x_5_sq(x)
  sum((x.-5).^2).-5
end

function x_plus_y(x)
    ress=0
    for i=1:length(x)
      ress+=x[i]
    end
    return ress
end

#

function rand_poly(x)
    return sum(rand(length(x)).*x)
end

# rand_poly(rand(10))
# x=-10:0.001:10
# y=sin.(10rand()x).+(rand()/1000)x.^2
# plot(x,y)

# function limited_domain(x::Vector)::Float64
#   sum( sqrt.(1 ./ (x.+1) .+ (x.+1) ) )  ## min is at 1,1,1...
# end#function
# limited_domain(rand(9))

# x=-10:0.01:10

# y=rand()rastrigin.(x)-.010rand()elliptic.(x)
# scatter(x,y)

function alpine_1(x)
  return sum(abs.(x.*sin.(x)+0.1x))
end


text:
label(node::ArgumentNode) = String(node.name)
label(node::JuliaNode) = String(node.name)
label(node::RandomChoiceNode) = "$(node.dist) $(node.addr) $(node.name)"
label(node::GenerativeFunctionCallNode) = "$(node.addr) $(node.name)"

function draw_graph(ir::StaticIR, graphviz, fname)
    dot = graphviz.Digraph()
    nodes_to_name = Dict{StaticIRNode,String}()
    for node in ir.nodes
        nodes_to_name[node] = label(node)
    end
    for node in ir.nodes
        if isa(node, ArgumentNode)
            shape = "diamond"
            color = "white"
            parents = []
        elseif isa(node, RandomChoiceNode)
            shape = "ellipse"
            color = "white"
            parents = node.inputs
        elseif isa(node, GenerativeFunctionCallNode)
            shape = "star"
            color = "white"
            parents = node.inputs
        elseif isa(node, JuliaNode)
            shape = "box"
            color = "white"
            parents = values(node.inputs)
        end
        if node === ir.return_node
            color = "lightblue"
        end
        dot[:node](nodes_to_name[node], nodes_to_name[node], shape=shape, fillcolor=color, style="filled")
        for parent in parents
            dot[:edge](nodes_to_name[parent], nodes_to_name[node])
        end
    end
    dot[:render](fname, view=true)
end

function draw_graph(gen_fn::StaticIRGenerativeFunction, graphviz, fname)
    draw_graph(get_ir(typeof(gen_fn)), graphviz, fname)
end

export draw_graph


text:
function molecule(::Molecule"N₂")
  return """
  N	0.0000000	0.0000000	0.5502960
  N	0.0000000	0.0000000	-0.5502960
  """
end
molecule(m::Molecule"nitrogen") = molecule(alias(m))
alias(::Molecule"nitrogen") = Molecule"N₂"()


text:
export delete

function delete(out, mask::BitArray)
    return repack(out[:x][:, mask], out[:v][:, mask], out[:y][:, mask], out[:volume][mask], out[:type][mask])
end

"""
    delete!(obj::T, f::Function) where T

Delete mesh particle for object using boolean array from function f.
"""
function delete(obj::T, f::Function) where T
    # x, v, y, vol, type
    function func(out)
        mask = vec(f(out)) .== false
        out = delete(out, mask)
        return out
    end
    if isa(obj, Shape)
        return PostOpObj(obj, func)
    elseif isa(obj, PostOpObj)
        push!(obj.operations, func)
        return obj
    else
        error("Not allowed.")
    end
end

max_stars_repo_name: sajidmannan/PDMesh

max_issues_repo_name: sajidmannan/PDMesh

max_forks_repo_name: sajidmannan/PDMesh


text:
using Kafka
using Base.Test

io = IOBuffer()

writeobj(io, "test")
writeobj(io, 42)
writeobj(io, UInt8[1, 2, 3])
writeobj(io, ["hello", "world"])

seek(io, 0)

@test readobj(io, ASCIIString) == "test"
@test readobj(io, Int64) == 42
@test readobj(io, Vector{UInt8}) == UInt8[1, 2, 3]
@test readobj(io, Vector{ASCIIString}) == ["hello", "world"]

println("Ok.")


text:
#    Problem : GROUP C
#    *********
#    This is an error of a Power problem by Oren.
#
#    Origonal SIF Source:
#    S.S. Oren,
#    Self-scaling variable metric algorithms,
#    Part II: implementation and experiments"
#    Management Science 20(5):863-874, 1974.
#
#    The CUTE implementation specification of the source
#    formula is not correct.
#
#    POWER.SIF classification OUR2-AN-V-0
#
# Daniel Henderson, 08/2021

f = (x) -> begin
    fx = 0.0
    for i in firstindex(x):lastindex(x)
        fx += x[i]^2*i
    end
    return fx^2
end

g! = (g, x) -> begin
    α = 0.0
    for i in firstindex(x):lastindex(x)
        α += x[i]^2*i
    end
    for i in firstindex(x):lastindex(x)
        g[i] = 4α*i*x[i]
    end
    return g
end

fg! = (g, x) -> begin
    α = 0.0
    for i in firstindex(x):lastindex(x)
        α += x[i]^2*i
    end
    for i in firstindex(x):lastindex(x)
        g[i] = 4α*i*x[i]
    end
    return α^2, g
end

init = (n::Int=10000) -> begin
    return n, ones(n)
end

TestSet["POWER"] = UncProgram("POWER", f, g!, fg!, init)


text:
AbstractTrees.children(p::Problem) = (p.objective, p.constraints)

AbstractTrees.children(e::AbstractExpr) = e.children

AbstractTrees.children(v::AbstractVariable) = ()
AbstractTrees.children(c::Constant) = ()

AbstractTrees.children(C::Constraint) = (C.lhs, C.rhs)
AbstractTrees.children(C::SDPConstraint) = (C.child,)
AbstractTrees.children(C::SOCConstraint) = C.children
AbstractTrees.children(C::SOCElemConstraint) = C.children
AbstractTrees.children(C::ExpConstraint) = C.children

AbstractTrees.printnode(io::IO, node::AbstractExpr) = summary(io, node)
AbstractTrees.printnode(io::IO, node::Constraint) = summary(io, node)

function AbstractTrees.printnode(io::IO, node::Vector{Constraint})
    if length(node) == 0
        print(io, "no constraints")
    else
        print(io, "constraints")
    end
end

AbstractTrees.printnode(io::IO, node::Problem) = print(io, node.head)

function AbstractTrees.printnode(io::IO, node::Constant)
    if length(node.value) <= 3
        show(IOContext(io, :compact => true), node.value)
    else
        summary(io, node.value)
    end
end

AbstractTrees.printnode(io::IO, node::AbstractVariable) = summary(io, node)


text:
using TimerOutputs

function efficiency()

to = TimerOutput()

nₚ = 11
nₑ = nₚ-1

x = [1/nₑ*i for i in 0:nₑ]
k = zeros(nₚ,nₚ)
f = zeros(nₚ)

# poolx = DataPool([[1.0/nₑ*i for i in 0:10],zeros(nₚ),zeros(nₚ)],Dict(:x=>1,:y=>2,:z=>3))
# poolξ = DataPool([[0.0],[0.0],[0.0],[1.0]],Dict(:ξ=>1,:η=>2,:γ=>3,:w=>4))
b = (x,y,z)->x^2
data = Dict(:x=>x,:y=>zeros(nₚ),:z=>zeros(nₚ))
paradata = Dict(:ξ=>[0.0],:w=>[1.0],:σ=>[0.0],:b=>[0.0])
coefficients = Dict(:k=>1.0)
functions = Dict(:b=>b)
𝓒 = [Node(i,data) for i in 1:2]
𝓖 = [Point(1,1,paradata)]
@timeit to "node" begin
    ξ = 𝓖[1]
    @timeit to "construct node" node = Node(1,data)
    @timeit to "getindex id" node.I
    @timeit to "getindex x" node.x
    @timeit to "getindex y" node.y
    @timeit to "getindex z" node.z
    @timeit to "getindex ξ" ξ.ξ
    @timeit to "getindex σ" ξ.σ
end
@timeit to "cell" begin
    # @timeit to "AbstractPoi" begin
    #     @timeit to "construct Poi1" ap = Poi1(25,x)
    #     @timeit to "get_weight" get_weight(ap,ξ₁)
    #     @timeit to "get_coordinates" get_coordinates(ap,ξ₁)
    #     @timeit to "get_shape_functions Poi1" get_shape_functions(ap,ξ₁,Val(:∂1))
    # end
    @timeit to "AbstractSeg" begin
        ξ = 𝓖[1]
        @timeit to "construct Seg2" ap = Seg2(𝓒,𝓖)
        @timeit to "get Jacobe" ap.J(0.0)
        @timeit to "get coordinates" ap.coordinates(ξ)
        @timeit to "get shape functions Seg2 ∂1" ap.𝝭(ξ)
        @timeit to "get shape functions Seg2 ∂x" ap.∂𝝭∂x(ξ)
        @timeit to "get shape functions Seg2 ∂y" ap.∂𝝭∂y(ξ)
        @timeit to "get shape functions Seg2 ∂z" ap.∂𝝭∂z(ξ)
        @timeit to "get shape functions Seg2 ∇" ap.∇𝝭(ξ)
        # @timeit to "get_shape_functions Seg2 ∂1" get_shape_functions(ap,ξ,Val(:∂1))
    #     @timeit to "get_shape_functions Seg2 ∂x" get_shape_functions(ap,ξ,Val(:∂x))
    #     @timeit to "get_shape_functions Seg2 ∂y" get_shape_functions(ap,ξ,Val(:∂y))
        # @timeit to "get_shape_functions Seg2 ∂1 ∂x ∂y" get_shape_functions(ap,ξ₁,Val(:∂1),Val(:∂x),Val(:∂y))
        # @timeit to "get_shape_functions Seg2 ∂1 ∂x ∂y ∂z" get_shape_functions(ap,ξ₁,Val(:∂1),Val(:∂x),Val(:∂y),Val(:∂z))
    end
    # 𝒞 = [1,2,12]
    # @timeit to "AbstractTri" begin
    #     @timeit to "construct Tri3" ap = Tri3(𝒞,x)
    #     @timeit to "get_weight" get_weight(ap,ξ₂)
    #     @timeit to "get_coordinates" get_coordinates(ap,ξ₂)
    #     @timeit to "get_shape_functions Tri3 ∂1" get_shape_functions(ap,ξ₂,Val(:∂1))
    #     @timeit to "get_shape_functions Tri3 ∂x" get_shape_functions(ap,ξ₂,Val(:∂x))
    #     @timeit to "get_shape_functions Tri3 ∂y" get_shape_functions(ap,ξ₂,Val(:∂y))
    #     @timeit to "get_shape_functions Tri3 ∂1 ∂x ∂y" get_shape_functions(ap,ξ₂,Val(:∂1),Val(:∂x),Val(:∂y))
    #     @timeit to "get_shape_functions Tri3 ∂1 ∂x ∂y ∂z" get_shape_functions(ap,ξ₂,Val(:∂1),Val(:∂x),Val(:∂y),Val(:∂z))
    # end
    # 𝒞 = [1,2,13,12]
    # @timeit to "AbstractQuad" begin
    #     @timeit to "construct Quad" ap = Quad(𝒞,x)
    #     # @timeit to "get_weight" get_weight(ap,ξ₂)
    #     # @timeit to "get_coordinates" get_coordinates(ap,ξ₂)
    #     @timeit to "get_shape_functions Quad ∂1" get_shape_functions(ap,ξ₂,Val(:∂1))
    #     # @timeit to "get_shape_functions Quad ∂x ∂y" get_shape_functions(ap,ξ₂,Val(:∂x),Val(:∂y))
    #     # @timeit to "get_shape_functions Quad ∂x" get_shape_functions(ap,ξ₂,Val(:∂x))
    #     # @timeit to "get_shape_functions Quad ∂y" get_shape_functions(ap,ξ₂,Val(:∂y))
    #     # @timeit to "get_shape_functions Quad ∂1 ∂x ∂y" get_shape_functions(ap,ξ₂,Val(:∂1),Val(:∂x),Val(:∂y))
    #     # @timeit to "get_shape_functions Quad ∂1 ∂x ∂y ∂z" get_shape_functions(ap,ξ₂,Val(:∂1),Val(:∂x),Val(:∂y),Val(:∂z))
    # end
end
@timeit to "Operator" begin
    @timeit to "Construction" op = Operator(Val(:∫∇v∇udΩ),coefficients)
    # @timeit to "∇v∇u 1" op(ap,k,f,Val(:∇v∇u))
    @timeit to "∇v∇u 2" op(ap,k,f)
end
show(to)

end

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator