text:
"Holds the tableau of a stochastic explicit Runge-Kutta method."
struct TableauSERK{T} <: AbstractTableauERK{T}
    name::Symbol
    s::Int

    qdrift::CoefficientsRK{T}
    qdiff::CoefficientsRK{T}
    qdiff2::CoefficientsRK{T}

    # Order of the tableau is not included, because unlike in the deterministic
    # setting, it depends on the properties of the noise (e.g., the dimension of
    # the Wiener process and the commutativity properties of the diffusion matrix)
    #
    # Orders stored in qdrift, qdiff and qdiff2 are understood as the classical orders of these methods.

    function TableauSERK{T}(name, qdrift, qdiff, qdiff2) where {T}
        @assert qdrift.s == qdiff.s == qdiff2.s
        @assert qdrift.c[1] == qdiff.c[1] == qdiff2.c[1] == 0.
        @assert istrilstrict(qdrift.a)
        @assert istrilstrict(qdiff.a)
        @assert istrilstrict(qdiff2.a)
        @assert !(qdrift.s==1. && qdrift.a[1,1] ≠ 0.)
        @assert !(qdiff.s==1. && qdiff.a[1,1] ≠ 0.)
        @assert !(qdiff2.s==1. && qdiff2.a[1,1] ≠ 0.)
        new(name, qdrift.s, qdrift, qdiff, qdiff2)
    end
end

function TableauSERK(name, qdrift::CoefficientsRK{T}, qdiff::CoefficientsRK{T}, qdiff2::CoefficientsRK{T}) where {T}
    TableauSERK{T}(name, qdrift, qdiff, qdiff2)
end

function TableauSERK(name, qdrift::CoefficientsRK{T}, qdiff::CoefficientsRK{T}) where {T}
    TableauSERK{T}(name, qdrift, qdiff, CoefficientsRK(T, :NULL, 0, zero(qdrift.a), zero(qdrift.b), zero(qdrift.c)) )
end

function TableauSERK(name::Symbol, order_drift::Int, a_drift::Matrix{T}, b_drift::Vector{T}, c_drift::Vector{T},
                                    order_diff::Int, a_diff::Matrix{T}, b_diff::Vector{T}, c_diff::Vector{T},
                                    order_diff2::Int, a_diff2::Matrix{T}, b_diff2::Vector{T}, c_diff2::Vector{T}) where {T}
    TableauSERK{T}(name, CoefficientsRK(name, order_drift, a_drift, b_drift, c_drift), CoefficientsRK(name, order_diff, a_diff, b_diff, c_diff),
                    CoefficientsRK(name, order_diff2, a_diff2, b_diff2, c_diff2))
end

function TableauSERK(name::Symbol, order_drift::Int, a_drift::Matrix{T}, b_drift::Vector{T}, c_drift::Vector{T},
                                    order_diff::Int, a_diff::Matrix{T}, b_diff::Vector{T}, c_diff::Vector{T}) where {T}
    TableauSERK{T}(name, CoefficientsRK(name, order_drift, a_drift, b_drift, c_drift), CoefficientsRK(name, order_diff, a_diff, b_diff, c_diff),
                    CoefficientsRK(:NULL, 0, zero(a_drift), zero(b_drift), zero(c_drift)))
end

"Stochastic Explicit Runge-Kutta integrator."
struct IntegratorSERK{DT,TT,FT} <: StochasticIntegrator{DT,TT}
    equation::SDE{DT,TT,FT}
    tableau::TableauSERK{TT}
    Δt::TT
    ΔW::Vector{DT}
    ΔZ::Vector{DT}

    q ::Matrix{Vector{DT}}    # q[k,m]  - holds the previous time step solution (for k-th sample path and m-th initial condition)
    Q::Vector{Vector{DT}}     # Q[j][k] - the k-th component of the j-th internal stage
    V::Vector{Vector{DT}}     # V[j][k] - the k-th component of v(Q[j])
    B::Vector{Matrix{DT}}     # B[j]    - the diffusion matrix B(Q[j])

    ydiff::Vector{DT}

    function IntegratorSERK{DT,TT,FT}(equation, tableau, Δt) where {DT,TT,FT}
        D = equation.d
        M = equation.m
        NS= equation.ns
        NI= equation.n
        S = tableau.s

        # create solution vectors
        q = create_solution_vector(DT, D, NS, NI)

        # create internal stage vectors
        Q = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        B = create_internal_stage_vector(DT, D, M, S)

        new(equation, tableau, Δt, zeros(DT,M), zeros(DT,M), q, Q, V, B, zeros(DT,M))
    end
end

function IntegratorSERK(equation::SDE{DT,TT,FT}, tableau::TableauSERK{TT}, Δt::TT) where {DT,TT,FT}
    IntegratorSERK{DT,TT,FT}(equation, tableau, Δt)
end

function initialize!(int::IntegratorSERK, sol::SolutionSDE, k::Int, m::Int)
    check_solution_dimension_asserts(sol, k, m)

    # copy initial conditions from solution
    get_initial_conditions!(sol, int.q[k,m], k, m)
end

"""
Integrate SDE with explicit Runge-Kutta integrator.
  Calculating the n-th time step of the explicit integrator for the
  sample path r and the initial condition m
"""
function integrate_step!(int::IntegratorSERK{DT,TT,FT}, sol::SolutionSDE{DT,TT,NQ,NW}, r::Int, m::Int, n::Int) where {DT,TT,FT,NQ,NW}
    local tᵢ::TT
    local ydrift::DT

    # copy the increments of the Brownian Process
    if NW==1
        #1D Brownian motion, 1 sample path
        int.ΔW[1] = sol.W.ΔW[1,n-1]
        int.ΔZ[1] = sol.W.ΔZ[1,n-1]
    elseif NW==2
        #Multidimensional Brownian motion, 1 sample path
        for l = 1:sol.nm
            int.ΔW[l] = sol.W.ΔW[l,n-1]
            int.ΔZ[l] = sol.W.ΔZ[l,n-1]
        end
    elseif NW==3
        #1D or Multidimensional Brownian motion, r-th sample path
        for l = 1:sol.nm
            int.ΔW[l] = sol.W.ΔW[l,n-1,r]
            int.ΔZ[l] = sol.W.ΔZ[l,n-1,r]
        end
    end

    # calculates v(t,tQ) and assigns to the i-th column of V
    int.equation.v(sol.t[0] + (n-1)*int.Δt, int.q[r,m], int.V[1])
    # calculates B(t,tQ) and assigns to the matrix BQ[:,:,1]
    int.equation.B(sol.t[0] + (n-1)*int.Δt, int.q[r,m], int.B[1])

    @inbounds for i in 2:int.tableau.s
        for k in eachindex(int.Q[i])
            # contribution from the drift part
            ydrift = 0.
            for j = 1:i-1
                ydrift += int.tableau.qdrift.a[i,j] * int.V[j][k]
            end

            # ΔW contribution from the diffusion part
            int.ydiff .= 0.
            for j = 1:i-1
                for l = 1:sol.nm
                    int.ydiff[l] += int.tableau.qdiff.a[i,j] * int.B[j][k,l]
                end
            end

            int.Q[i][k] = int.q[r,m][k] + int.Δt * ydrift + dot(int.ydiff,int.ΔW)

            # ΔZ contribution from the diffusion part
            if int.tableau.qdiff2.name ≠ :NULL
                int.ydiff .= 0.
                for j = 1:i-1
                    for l = 1:sol.nm
                        int.ydiff[l] += int.tableau.qdiff2.a[i,j] * int.B[j][k,l]
                    end
                end

                int.Q[i][k] += dot(int.ydiff,int.ΔZ)/int.Δt
            end

        end
        tᵢ = sol.t[0] + (n-1)*int.Δt + int.Δt * int.tableau.qdrift.c[i]
        int.equation.v(tᵢ, int.Q[i], int.V[i])
        int.equation.B(tᵢ, int.Q[i], int.B[i])
    end

    # compute final update
    if int.tableau.qdiff2.name == :NULL
        update_solution!(int.q[r,m], int.V, int.B, int.tableau.qdrift.b, int.tableau.qdiff.b, int.Δt, int.ΔW)
    else
        update_solution!(int.q[r,m], int.V, int.B, int.tableau.qdrift.b, int.tableau.qdiff.b, int.tableau.qdiff2.b, int.Δt, int.ΔW, int.ΔZ)
    end

    # take care of periodic solutions
    cut_periodic_solution!(int.q[r,m], int.equation.periodicity)

    # copy to solution
    set_solution!(sol, int.q[r,m], n, r, m)
end


text:
# ---
# title: 72. Edit Distance
# id: problem72
# author: Qling
# date: 2021-03-17
# difficulty: Hard
# categories: String, Dynamic Programming
# link: <https://leetcode.com/problems/edit-distance/description/>
# hidden: true
# ---
#
# Given two strings `word1` and `word2`, return _the minimum number of
# operations required to convert`word1` to `word2`_.
#
# You have the following three operations permitted on a word:
#
#   * Insert a character
#   * Delete a character
#   * Replace a character
#
#
#
# **Example 1:**
#
#
#
#     Input: word1 = "horse", word2 = "ros"
#     Output: 3
#     Explanation:
#     horse -> rorse (replace 'h' with 'r')
#     rorse -> rose (remove 'r')
#     rose -> ros (remove 'e')
#
#
# **Example 2:**
#
#
#
#     Input: word1 = "intention", word2 = "execution"
#     Output: 5
#     Explanation:
#     intention -> inention (remove 't')
#     inention -> enention (replace 'i' with 'e')
#     enention -> exention (replace 'n' with 'x')
#     exention -> exection (replace 'n' with 'c')
#     exection -> execution (insert 'u')
#
#
#
#
# **Constraints:**
#
#   * `0 <= word1.length, word2.length <= 500`
#   * `word1` and `word2` consist of lowercase English letters.
#
#
## @lc code=start
using LeetCode

function min_distance(word1::String, word2::String)::Int32
    m, n = length(word1) + 1, length(word2) + 1
    dp = fill(0, m, n)

    for i = 1:m
        dp[i, 1] = i - 1
    end

    for j = 1:n
        dp[1, j] = j - 1
    end

    for i = 2:m, j = 2:n
        if word1[i-1] == word2[j-1]
            dp[i, j] = dp[i-1, j-1]
        else
            dp[i, j] = min(dp[i-1, j-1], dp[i, j-1], dp[i-1, j]) + 1
        end
    end

    return dp[m, n]
end

## @lc code=end

ext: jl


text:
# following https://support.hdfgroup.org/HDF5/doc/RM/RM_H5F.html#File-StartSwmrWrite
# and https://support.hdfgroup.org/HDF5/docNewFeatures/SWMR/HDF5_SWMR_Users_Guide.pdf
using HDF5
using Test
using Distributed

if nprocs() == 1
    procs = addprocs(1)
else
    procs = Int64[]
end
@everywhere using HDF5

@testset "swmr" begin
fname = tempname()

@testset "swmr modes" begin
    h5open(fname, "w", swmr=true) do h5
        h5["foo"] = collect(1:10)
    end
    h5open(fname, "r", swmr=true) do h5
        @test read(h5["foo"]) == collect(1:10)
    end
    h5open(fname, "r+", swmr=true) do h5
        @test read(h5["foo"]) == collect(1:10)
    end
end

@testset "h5d_oappend" begin
    h5open(fname, "w") do h5
        g = create_group(h5, "shoe")
        d = create_dataset(g, "bar", datatype(Float64), ((1,), (-1,)), chunk=(100,))
        dxpl_id = HDF5.get_create_properties(d)
        v = [1.0, 2.0]
        memtype = datatype(Float64)
        # @test HDF5.h5d_oappend(d, dxpl_id, 0, length(v), memtype, v)
    end
end

function dataset_write(d, ch_written, ch_read)
    for i = 1:10
        @assert take!(ch_read) == true
        HDF5.set_extent_dims(d, (i*10,))
        inds::UnitRange{Int} = (1:10) .+ (i - 1) * 10
        d[inds] = inds
        flush(d) # flush the dataset
        put!(ch_written,i)
    end
end

@everywhere function dataset_read(d, ch_written, ch_read)
    n = nlast = length(d)
    nbigger = 0
    i = 0
    put!(ch_read, true)
    while n < 100
        i = take!(ch_written)
        for j = 1:1000 # wait for new data to be available to avoid CI failures
            HDF5.refresh(d)
            nlast, n = n, length(d)
            n > nlast && break
            sleep(0.001)
        end
        vals = read(d)
        @assert vals == collect(1:n)
        n > nlast && (nbigger += 1)
        put!(ch_read,true)
    end
    return nbigger
end

@everywhere function swmr_reader(fname, ch_written, ch_read)
    h5open(fname, "r", swmr=true) do h5
        d = h5["foo"]
        dataset_read(d, ch_written, ch_read)
    end
end

# Spawn a reader function in a 2nd process, provide two channels for synchronization.
# Run a writing function in this process. The writing function writes,
# then notifies `ch_read`, then the reading function reads, and notifies `ch_read`. So read
# attempts should always follow writes, though there may be a delay before the data is available
# so there is a step that sleeps until data is available.
function remote_test(h5)
    ch_written, ch_read = RemoteChannel(1), RemoteChannel(1)
    a = @spawn(swmr_reader(fname, ch_written, ch_read))
    dataset_write(h5["foo"], ch_written, ch_read)
    nbigger = fetch(a)
    @test nbigger == 10
end

# create datasets and attributes before staring swmr writing
function prep_h5_file(h5)
    d = create_dataset(h5, "foo", datatype(Int), ((1,), (100,)), chunk=(1,))
    attributes(h5)["bar"] = "bar"
    g = create_group(h5, "group")
end

@testset "create by libver, then start_swmr_write" begin
    #test this h5open method with keyword arg
    h5open(fname, "w", libver_bounds=(:latest, :latest), swmr=false) do h5
        prep_h5_file(h5)
        HDF5.start_swmr_write(h5) # after creating datasets
        remote_test(h5)
    end
end

@testset "create by swmr mode, then close and open again" begin
    h5open(fname, "w", swmr=true) do h5
        prep_h5_file(h5)
    end
    # close the file after creating datasets, open again with swmr write access but not truncate
    h5open(fname, "r+", swmr=true) do h5
        remote_test(h5)
    end
end

rm(fname) # cleanup file created by swmr tests

if nprocs() > 1
    rmprocs(procs)
end

end # testset swmr


text:
push!(LOAD_PATH, pwd()*"/../../src/")
using Distributions
using PyPlot
using BVB

lower = -100.
upper = 100.

max_iter = 50

# example 0
true_d = MixtureModel(Normal, [(5.0, 2), (10.0, 2), (20.0, 1), (-2, 1)], [1.0/3, 1.0/4, 1.0/4, 1.0/6])
ll(x) = logpdf(true_d, x)

init_approx_d =  Normal(0, 40)
x = linspace(-10,30,500)

include("run_laplacian_gradient_boosting_exclusive.jl")
show()

# example 1
true_d = MixtureModel(Normal, [(-2.0, 1), (2.0, 1)], [0.5, 0.5])
ll(x) = logpdf(true_d, x)

init_approx_d = Normal(0, 20)
x = linspace(-10,10,500)

include("run_laplacian_gradient_boosting_exclusive.jl")
show()

println("Done.")

# example 2
true_d = MixtureModel(Normal, [(-3.5, 1), (-1.0, 0.4), (1.0, 0.5), (6.0, 1)], [0.2, 0.3, 0.3, 0.2])
ll(x) = logpdf(true_d, x)
x = linspace(-10,10,500)

include("run_laplacian_gradient_boosting_exclusive.jl")
show()

println("Done.")

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
"""
    ijth_coordinate_update!(i::Int64, j::Int64, T_minus_p::Int64, estim_settings::EstimSettings, Y_i::SubArray{Float64}, Y_lagged::FloatArray, Ψ::FloatArray)

Update the (i,j)-th element of Ψ the coordinate descent algorithm (Friedman et al., 2010) as in Pellegrino (2020).

# References
Friedman et al. (2010) and Pellegrino (2020)
"""
function ijth_coordinate_update!(i::Int64, j::Int64, T_minus_p::Int64, estim_settings::EstimSettings, Y_i::SubArray{Float64}, Y_lagged::FloatArray, Ψ::FloatArray)

    # Views
    Y_lagged_j = @view Y_lagged[j,:];

    # Partial residuals
    Ψ_i_ex_j = Ψ[i,:]; # Julia creates a copy of Ψ[i,:]
    Ψ_i_ex_j[j] = 0.0;
    V_ij = Y_i - Y_lagged'*Ψ_i_ex_j;

    # Scalar product between Y_lagged_j and V_ij divided by (T-p)
    scalar_Y_lagged_j = sum(V_ij .* Y_lagged_j)/T_minus_p;

    # Norm squared of Y_lagged_j divided by (T-p)
    normsq_Y_lagged_j = sum(Y_lagged_j .^ 2)/T_minus_p;

    # Soft-thresholding operator: update for the coefficients
    hyper_prod = estim_settings.λ*estim_settings.β^fld(j-1, estim_settings.n);
    Ψ[i, j] = soft_thresholding(scalar_Y_lagged_j, estim_settings.α*hyper_prod)/(normsq_Y_lagged_j + (1-estim_settings.α)*hyper_prod);

    # Return update for objfun_new
    update_objfun_new = hyper_prod*(0.5*(1-estim_settings.α)*Ψ[i,j]^2 + estim_settings.α*abs(Ψ[i,j]));
    return update_objfun_new;
end

"""
    coordinate_descent(estim_settings::EstimSettings)

Estimate an elastic-net VAR(p) with the coordinate descent algorithm (Friedman et al., 2010) as in Pellegrino (2020).

# Arguments
- `estim_settings`: settings used for the estimation

# References
Friedman et al. (2010) and Pellegrino (2020)
"""
function coordinate_descent(estim_settings::EstimSettings)

    # Check inputs
    check_bounds(estim_settings.p, 1);
    check_bounds(estim_settings.λ, 0);
    check_bounds(estim_settings.α, 0, 1);
    check_bounds(estim_settings.β, 1);
    check_bounds(estim_settings.max_iter, 3);
    check_bounds(estim_settings.max_iter, estim_settings.prerun);
    check_bounds(estim_settings.n, 2); # It supports only multivariate models (for now ...)

    # Interpolate data
    data = interpolate(estim_settings.Y, estim_settings.n, estim_settings.T);

    # VAR(p) data
    Y, Y_lagged = lag(data, estim_settings.p);
    T_minus_p = estim_settings.T - estim_settings.p;

    # Memory pre-allocation
    Ψ = zeros(estim_settings.n, estim_settings.np);
    objfun_old = -Inf;
    objfun_new = -Inf;

    # Run coordinate_descent, in turn, for each target variable
    for i=1:estim_settings.n

        # Views
        Y_i = @view Y[i,:];
        Ψ_i = @view Ψ[i,:];

        # Loop until max_iter is reached
        for iter=1:estim_settings.max_iter

            # Reset objfun_new
            objfun_new = 0.0;

            # Loop over the predictors
            for j=1:estim_settings.np
                objfun_new += ijth_coordinate_update!(i, j, T_minus_p, estim_settings, Y_i, Y_lagged, Ψ);
            end

            # Update objfun_new
            objfun_new += (1/(2*T_minus_p))*sum((Y_i - Y_lagged'*Ψ_i).^2);
            verb_message(estim_settings.verb, "coordinate_descent ($i-th row) > iter=$(iter), objfun=$(round(objfun_new, digits=5))");

            # Stop when the algorithm converges
            if iter > 1
                if isconverged(objfun_new, objfun_old, estim_settings.tol, estim_settings.ε, false)
                    verb_message(estim_settings.verb, "coordinate_descent ($i-th row) > converged!\n");
                    break;
                end
            end

            # Store current run information
            objfun_old = copy(objfun_new);
        end
    end

    # Estimate var-cov matrix of the residuals
    V = Y - Ψ*Y_lagged;
    Σ = Symmetric((V*V')./T_minus_p)::SymMatrix;

    # Return output
    return Ψ, Σ;
end


text:
"""
    quadgk(f,s::AbstractReferenceShape;kwargs...)

Integrate `f` over the reference shape `s`.
"""
quadgk(f,::Type{ReferenceLine};kwargs...)     = quadgk(f,0,1;kwargs...)

quadgk(f,l::AbstractReferenceShape;kwargs...) = quadgk(f,typeof(l);kwargs...)

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
module CustomSpinnerFrostHelper

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/IceSpinner" IceSpinner(x::Integer, y::Integer, attachToSolid::Bool=false, directory::String="danger/FrostHelper/icecrystal", spritePathSuffix::String="", destroyColor::String="b0eaff", tint::String="ffffff", borderColor::String="000000", moveWithWind::Bool=false, dashThrough::Bool=false, bloomAlpha::Number=0.0, bloomRadius::Number=0.0, collidable::Bool=true, rainbow::Bool=false, drawOutline::Bool=true, attachGroup::Integer=-1)

const placements = Ahorn.PlacementDict(
   "Custom Spinner (Rainbow Spinner Sprite, Frost Helper)" => Ahorn.EntityPlacement(
        IceSpinner,
		"point",
		Dict{String, Any}(),
        function(entity)
			entity.data["directory"] = "danger/crystal"
			entity.data["spritePathSuffix"] = "_white"
        end
    ),
	"Custom Spinner (Frost Helper)" => Ahorn.EntityPlacement(
		IceSpinner
	)
)

function Ahorn.selection(entity::IceSpinner)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 8, y - 8, 16, 16)
end

# TODO - Add support for background
function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::IceSpinner, room::Maple.Room)
	texture = get(entity.data, "directory", "danger/FrostHelper/icecrystal")
	spritePathSuffix = get(entity.data, "spritePathSuffix", "")

	realTexture = "$texture/fg$(spritePathSuffix)03.png"

	rawTint = Ahorn.argb32ToRGBATuple(parse(Int, lstrip(get(entity.data, "tint", "ffffff"), [ '#' ]), base=16))[1:3] ./ 255
	realTint = (rawTint..., 1.0)
    Ahorn.drawSprite(ctx, realTexture, 0, 0, tint = realTint)
end

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
@testset "UMLS" begin
    user = ""
    psswd = ""
    apikey = ""

    try
        user = ENV["UMLS_USER"]
        psswd = ENV["UMLS_PSSWD"]
        apikey = ENV["UMLS_APIKEY"]
    catch
        @warn "UMLS tests could not run: Set up credentials as enviroment variables: UMLS_USER, UMLS_PSSWD and UMLS_APIKEY"
    end

    if (user != "") && (psswd != "") && (apikey != "")
        term = "obesity"
        query = Dict("string"=>term, "searchType"=>"exact" )

        tgt = get_tgt(force_new=true, username=user, password=psswd)
        #do it again to test reading from file
        tgt = get_tgt(username=user, password=psswd)
        #test apikey
        tgt = get_tgt(force_new=true, apikey=apikey)

        @testset "Testing Search/Content/CUI" begin
            all_results= search_umls(tgt, query)
            @test length(all_results[1]["result"]["results"]) == 2
            cui = best_match_cui(all_results)
            @test cui == "C0028754"
            sm = get_semantic_types(tgt, cui)
            @test sm[1] == "Disease or Syndrome"
        end
    end
end


text:
module PosteriorFactorizationTest

using Test
using ForneyLab

using ForneyLab: Cluster, deterministicEdges, deterministicEdgeSet, isDeterministic

@testset "PosteriorFactorization" begin
    pfz = PosteriorFactorization()
    @test pfz.posterior_factors == Dict{Symbol, PosteriorFactor}()
    @test pfz.edge_to_posterior_factor == Dict{Edge, PosteriorFactor}()
    @test pfz.node_edge_to_cluster == Dict{Tuple{FactorNode, Edge}, Cluster}()
    @test currentPosteriorFactorization() === pfz
end

@testset "isDeterministic" begin
    # Clamp
    fg = FactorGraph()
    @RV x = constant(0.0)
    @test isDeterministic(fg.nodes[:clamp_1].i[:out], Dict{Interface, Bool}())

    # Stochastic
    fg = FactorGraph()
    @RV x ~ GaussianMeanVariance(0.0, 1.0)
    @test !isDeterministic(fg.nodes[:gaussianmeanvariance_1].i[:out], Dict{Interface, Bool}())

    # Equality
    fg = FactorGraph()
    @RV x = equal(constant(0.0), constant(0.0))
    @test isDeterministic(fg.nodes[:equality_1].i[1],
                          Dict{Interface, Bool}(fg.nodes[:clamp_1].i[:out] => true,
                                                fg.nodes[:clamp_2].i[:out] => false))
    @test !isDeterministic(fg.nodes[:equality_1].i[1],
                           Dict{Interface, Bool}(fg.nodes[:clamp_1].i[:out] => false,
                                                 fg.nodes[:clamp_2].i[:out] => false))
    # Deterministic
    fg = FactorGraph()
    @RV x ~ Addition(0.0, 0.0)
    @test isDeterministic(fg.nodes[:addition_1].i[:out],
                          Dict{Interface, Bool}(fg.nodes[:clamp_1].i[:out] => true,
                                                fg.nodes[:clamp_2].i[:out] => true))
    @test !isDeterministic(fg.nodes[:addition_1].i[:out],
                           Dict{Interface, Bool}(fg.nodes[:clamp_1].i[:out] => true,
                                                 fg.nodes[:clamp_2].i[:out] => false))
end

@testset "deterministicEdgeSet and deterministicEdges" begin
    fg = FactorGraph()

    @RV x ~ GaussianMeanVariance(0.0, 1.0)
    @RV y ~ GaussianMeanVariance(0.0, 1.0)
    @RV w = x + y
    @RV a = constant(1.0)
    @RV z = a*w
    placeholder(z, :z)

    e_x = x.edges[1]
    e_y = y.edges[1]
    e_w = w.edges[1]
    e_a = a.edges[1]
    e_z = z.edges[1]
    e_1 = fg.nodes[:clamp_1].i[:out].edge
    e_2 = fg.nodes[:clamp_2].i[:out].edge
    e_3 = fg.nodes[:clamp_3].i[:out].edge
    e_4 = fg.nodes[:clamp_4].i[:out].edge

    @test deterministicEdgeSet(e_x) == Set{Edge}([e_w, e_a, e_z])
    @test deterministicEdges(fg) == Set{Edge}([e_w, e_a, e_z, e_1, e_2, e_3, e_4])
end

@testset "deterministicEdges for graph with dangling edge" begin
    fg = FactorGraph()

    @RV x ~ GaussianMeanVariance(1.0, 1.0)
    @RV y = 4.0 * x

    e_1 = fg.nodes[:clamp_1].i[:out].edge
    e_2 = fg.nodes[:clamp_2].i[:out].edge
    e_3 = fg.nodes[:clamp_3].i[:out].edge

    @test deterministicEdges(fg) == Set{Edge}([e_1, e_2, e_3])
end

end # module


text:
# let

	mat = [0.0 0.1 0.2 0.3]'
	dset = Dataset()
	indA = 0

	add_float_var(dset, "A")
	@test get_number_of_variables(dset) == 1
	@test !is_discrete(dset, indA)
	@test get_id(dset, indA) == "A"
	set_id(dset, indA, "Bubba")
	@test get_id(dset, indA) == "Bubba"
	set_id(dset, indA, "A")

	set_number_of_records(dset, length(mat))
	@test get_number_of_records(dset) == length(mat)

	for (i,v) in enumerate(mat)
		set_float(dset, indA, i-1, v)
	end
	for (i,v) in enumerate(mat)
		@test isapprox(get_float(dset, indA, i-1), v, atol=1e-8)
	end

	test_vec_approx_equal( dset[indA, 0:2], [0.0, 0.1, 0.2], 1e-7 )
	test_vec_approx_equal( dset[indA], [0.0, 0.1, 0.2, 0.3], 1e-7 )
	@test has_missing_data( dset, indA ) == false
	@test is_constant( dset, indA ) == false

	# add and remove a variable
	add_float_var(dset, "B")
	@test get_number_of_variables(dset) == 2
	remove_var(dset, 1)
	@test get_number_of_variables(dset) == 1
	for (i,v) in enumerate(mat)
		@test isapprox(get_float(dset, indA, i-1), v, atol=1e-8)
	end

	# create a dset copy
	dset2 = Dataset(dset)
	@test get_number_of_variables(dset2) == 1
	@test !is_discrete(dset2, indA)
	for (i,v) in enumerate(mat)
		@test isapprox(get_float(dset2, indA, i-1), v, atol=1e-8)
	end

	# test DatasetParseParams
	pp = DatasetParseParams()
	@test_throws ErrorException pp[:sadasdf]
	@test getindex(pp, :columnIdsPresent) == true
	@test pp[:missingValueToken] == "*"
	@test pp[:missingInt] == -1
	@test isnan(pp[:missingFloat])
	@test pp[:columnIdsPresent] == true
	setindex!(pp, "h", :missingValueToken)
	@test pp[:missingValueToken] == "h"
	pp[:missingInt] = 5
	@test pp[:missingInt] == 5
	pp[:missingFloat] = 1.2
	@test isapprox(pp[:missingFloat], 1.2, atol=1e-7)
	pp[:columnIdsPresent] = false
	@test pp[:columnIdsPresent] == false
	@test_throws ErrorException pp[:sadasdf]=5
	@test_throws AssertionError pp[:missingValueToken]=5

	# test DatasetWriteParams
	wp = DatasetWriteParams()
	@test_throws ErrorException pp[:sadasdf]
	@test getindex(wp, :columnIdsPresent) == true
	@test wp[:useStateIndices] == false
	@test wp[:separator] == '\t'
	@test wp[:missingValueToken] == "*"
	@test wp[:floatFormat] == "%g"
	wp[:columnIdsPresent] = false
	@test wp[:columnIdsPresent] == false
	wp[:useStateIndices] = true
	@test wp[:useStateIndices] == true
	wp[:separator] = ' '
	@test wp[:separator] == ' '
	wp[:missingValueToken] = "howaboutthat"
	@test wp[:missingValueToken] == "howaboutthat"
	wp[:floatFormat] = "%.3f"
	@test wp[:floatFormat] == "%.3f"

	dset3 = Dataset()
	write_file( dset, "out.dat" )
	@test isfile("out.dat")
	read_file( dset3, "out.dat" )
	rm("out.dat")
	@test get_number_of_variables(dset3) == 1
	@test !is_discrete(dset3, indA)
	for (i,v) in enumerate(mat)
		@test isapprox(get_float(dset3, indA, i-1), v, atol=1e-8)
	end

	dset4 = Dataset()
	write_file( dset, "out.dat", wp )
	read_file( dset4, "out.dat", pp )
	@test get_number_of_variables(dset4) == 1
	@test !is_discrete(dset4, indA)
	for (i,v) in enumerate(mat)
		@test isapprox(get_float(dset4, indA, i-1), v, atol=1e-8)
	end
	rm("out.dat")

	# ------------

	@test find_variable(dset, "A") == 0

	indB = 1
	add_int_var(dset, "B")
	for (i,v) in enumerate([0,1,0,2])
		set_int(dset, indB, i-1, v)
	end
	@test get_int_data(dset, indB) == [0, 1, 0, 2]
	test_vec_approx_equal( get_float_data(dset, indA), [0.0, 0.1, 0.2, 0.3], 1e-7 )

	# discretize
	discretize( dset, indA, 2, statePrefix="state", algorithm=DSL_DISCRETIZE_UNIFORMWIDTH )
	test_vec_approx_equal( dset[indA], [0,0,1,1], 1e-7 )

	indC = 2
	add_float_var( dset, "C" )
	dset[indC,0:3] = [0.0, 0.3, 0.2, 0.1]
	test_vec_approx_equal( dset[indC], [0.0, 0.3, 0.2, 0.1], 1e-7 )
	discretize( dset, indC, 2, algorithm=DSL_DISCRETIZE_UNIFORMCOUNT, edges=[0.05, 5.05, 5.15] )
	test_vec_approx_equal( dset[indC], [0, 1, 1, 1], 1e-7 )

	indD = 3
	add_float_var( dset, "D" )
	dset[indD,0:3] = [0.0, 0.3, 0.2, 0.1]
	test_vec_approx_equal( dset[indD], [0.0, 0.3, 0.2, 0.1], 1e-7 )
	edges = discretize_with_info( dset, indD, 2 )
	test_vec_approx_equal( dset[indD], [0, 1, 1, 0], 1e-7 )

# end


text:
using RDatasets, Gadfly

set_default_plot_size(6inch, 3inch)

plot(dataset("car", "SLID"), color="Language", y="Wages",
     Geom.histogram(position=:stack, orientation=:horizontal))


text:
module Database

using MySQL

const HOST = "localhost"
const USER = "root"
const PASS = "root"
const DB = "six_degrees"

const CONN = MySQL.connect(HOST, USER, PASS, db = DB)

export CONN

disconnect() = MySQL.disconnect(CONN)

atexit(disconnect)

end

max_stars_repo_name: PacktPublishing/Julia-Projects

max_issues_repo_name: PacktPublishing/Julia-Projects

max_forks_repo_name: PacktPublishing/Julia-Projects


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 81cdc1aa-04fd-11eb-3235-bff0e3f52535
using Pkg, DrWatson

# ╔═╡ 81cdee8c-04fd-11eb-29f5-85927396b6dc
begin
	@quickactivate "StatisticalRethinkingTuring"
	using Turing
	using StatisticalRethinking
end

# ╔═╡ 5e06a11c-04f6-11eb-0116-1fd51d614f99
md"## clip-04-69-71t"

# ╔═╡ b8b8de5e-0512-11eb-0932-830c3055e00a
#include(srcdir("tools.jl"))

md"### snippet 4.64"

# ╔═╡ 81ce6d32-04fd-11eb-2c4f-77bdf6a3f9e2
begin
	df = CSV.read(sr_datadir("Howell1.csv"), DataFrame)
	scale!(df, [:weight])
end;

# ╔═╡ 82407c7c-04fd-11eb-2bba-3783e79b0671
md"### snippet 4.69"

# ╔═╡ da70abc2-0527-11eb-31fa-d5b029193ba5
f_cube(weight_s, a, b1, b2, b3) = a + b1 * weight_s + b2 * weight_s^2 + b3 * weight_s^3

# ╔═╡ 82483e30-04fd-11eb-1290-b98ae5d06a45
@model function cube(weight_s, heights)
    a ~ Normal(178, 20)
    b1 ~ LogNormal(0, 1)
    b2 ~ Normal(0, 10)
    b3 ~ Normal(0, 10)
    σ ~ Uniform(0, 50)
    μ = f_cube.(weight_s, a, b1, b2, b3)
    heights ~ MvNormal(μ, σ)
end

# ╔═╡ 8250f05c-04fd-11eb-0370-899604fffdfd
begin
	quap4_6t = quap(cube(df.weight_s, df.height), NelderMead())
	quap4_6t.coef
end

# ╔═╡ 825e8410-04fd-11eb-3e36-6d3bac027215
md"### snippets 4.70, 4.71"

# ╔═╡ 423cd608-0529-11eb-126e-7fde3de2165b
begin
	weight_seq_2 = range(-2.2, 2, length = 30)
	post_2 = DataFrame(rand(quap4_6t.distr, 1_000)', quap4_6t.params)
	mu_2 = f_cube.(weight_seq_2, post_2.a', post_2.b1', post_2.b2', post_2.b3') |> meanlowerupper
	sim_2 = rand.(Normal.(mu_2.raw, post_2.σ')) |> meanlowerupper
	weight_seq_rescaled = weight_seq_2 .* std(df.weight) .+ mean(df.weight)
	scatter(df.weight, df.height, ms = 3, alpha = 0.7, legend = false)
	plot!(weight_seq_rescaled, mu_2.mean, ribbon = (mu_2.mean .- mu_2.lower, mu_2.upper .- mu_2.mean))
	plot!(weight_seq_rescaled, sim_2.lower, fillrange = sim_2.upper, alpha = 0.3, la = 0.0, c = 2)
end

# ╔═╡ 4f386d5c-0573-11eb-084c-1758217d06f9
md"## End of clip-04-69-71t.jl"

# ╔═╡ Cell order:
# ╟─5e06a11c-04f6-11eb-0116-1fd51d614f99
# ╠═81cdc1aa-04fd-11eb-3235-bff0e3f52535
# ╠═81cdee8c-04fd-11eb-29f5-85927396b6dc
# ╟─b8b8de5e-0512-11eb-0932-830c3055e00a
# ╠═81ce6d32-04fd-11eb-2c4f-77bdf6a3f9e2
# ╟─82407c7c-04fd-11eb-2bba-3783e79b0671
# ╠═da70abc2-0527-11eb-31fa-d5b029193ba5
# ╠═82483e30-04fd-11eb-1290-b98ae5d06a45
# ╠═8250f05c-04fd-11eb-0370-899604fffdfd
# ╟─825e8410-04fd-11eb-3e36-6d3bac027215
# ╠═423cd608-0529-11eb-126e-7fde3de2165b
# ╟─4f386d5c-0573-11eb-084c-1758217d06f9


text:
"""
Abstract type for all Components
"""
abstract type Component end

struct ComponentPoolIndex
    type::Type{<:Component}
    index::Int
end
ComponentPoolIndex(t::Component, i::Int) = ComponentPoolIndex(typeof(t), i)

struct ComponentSet{T<:Tuple} <: Component
    _raw::Vector{ComponentPoolIndex}
end

struct ComponentPool
    _raw::Dict{Type{<:Component}, Vector{Component}}
    _sets::Vector{Type{<:ComponentSet}}
end
ComponentPool()  = ComponentPool(Dict(), Vector())
Base.getindex(c::ComponentPool, t::Type{<:Component}, i::Int) = c._raw[t][i]
Base.getindex(c::ComponentPool, t::Type{T}) where T<:Component = if t in keys(c._raw) c._raw[t] else Vector() end
function Base.getindex(c::ComponentPool, i::ComponentPoolIndex)
    res = getindex(c, i.type, i.index)
    if typeof(res) <: ComponentSet
        res = Tuple(c[v] for v in res._raw)
    end
    res
end

ComponentSet(T::Vararg{Type{<:Component}}) = ComponentSet{T}
ComponentSet(v::Vararg{ComponentPoolIndex}) = ComponentSet(collect(v))
ComponentSet{T}(args...) where T<:Tuple = ComponentSet(args...)
function ComponentSet(c::ComponentPool, v::Vararg{Tuple{<:Component}})
    set = ComponentSet(v)
    if !(v in c._sets) push!(c._sets, set) end
end

function indexify(v::Vector{Component})
    i = 1
    a = []
    for x in v
        push!(a, ComponentPoolIndex(x, i))
        i+=1
    end
    a
end

tupify(v) = (v,)
veccify(v) = [v]
#+------------------------------------------------------+
#| This is the main idea                                |
#| Combinations are only generated if needed            |
#+------------------------------------------------------+
function create_combinations(scheme, c::ComponentPool)
    arrs = []
    # Get a ComponentPoolIndex version of the array in the ComponentPool
    # for the Type of each element in the scheme.
    # since we're iterating the scheme in order,
    # there is no need to worry about passing the
    # args to the system function in the correct order.
    for x in scheme
        # If one of these has no components defined for its type,
        # then the overall combination will just be empty,
        # so we might as well just return []
        if c[x] == [] return c[x] end
        push!(arrs, indexify(c[x]))
    end
    # if arrs is length 1 then a will be returned straight away,
    # so we need each of its elements to be a tuple for compat.
    a = veccify.(first(arrs))
    k = 2
    # Create pairs until you reach the end of the scheme.
    # This will be ordered since arrs is based on scheme.
    while k<=length(arrs)
        a = create_pair_combinations(a, arrs[k])
        k+=1
    end
    return a
end
create_pair_combinations(a, b) = reduce(vcat, [[ae..., be...] for ae in a for be in b])


text:
module MinkowskiReduction

using LinearAlgebra, Random
export GaussReduce, RandUnimodMat, RandLowerTri, minkReduce, DeviousMat, isMinkReduced, orthogonalityDefect
"""
    minkReduce(U, V, W, debug=false)

Find the shortest equivalent basis of that lattice formed by {`U`, `V`, `W`}

```jldoctest
julia> U = [1, 2, 3]; V = [-1, 2, 3]; W = [3, 0, 4]; minkReduce(U,V,W)
([-2.0, 0.0, 0.0], [0.0, -2.0, 1.0], [-1.0, 2.0, 3.0])
```
"""
function minkReduce(U, V, W, debug=false)
    i = 0
    while true
        i +=1
        norms = [norm(U), norm(V), norm(W)]
        p = sortperm(norms)
        U,V,W = (U,V,W)[p] # sort into ascending order
        U,V,W = shortenW_in_UVW(U, V, W)
        i > 15 && error("minkReduce: Too many iterations")
#        println(U,V,W,i,"det",det(hcat(U,V,W)))
        norm(W) ≥ norm(V) && break
    end
    if !debug return U, V, W end
    return U, V, W, i
end

"""
    shortenW_in_UVW

Reduce vector W so that it is as close to the origin as possible.

Subtract multiples of U and V from W. W will remain in an affine plane, which
is parallel to the U-V plane but which passes through the end of the W vector.
(See Lecture notes in computer science, ISSN 0302-974, ANTS - VI : algorithmic
number theory, 2004, vol. 3076, pp. 338-357 ISBN 3-540-22156-5)
"""
function shortenW_in_UVW(U,V,W)
    # If U, V are themselves mink reduced, then the projection of W (shifted by multiples of U,V)
    # that is the closest to the origin, will be contained in the parallelogram
    # formed by the new U, V (GLWH Aug 2020. Not sure exactly what I meant here...)
    U, V = GaussReduce(U,V)
    # Project W into the U-V plane, call it T
    # n̂ is a unit vector ⟂ to U-V plane. Subtract multiples of thi
    # from W to move it into the U-V plane. Dot product gives number of multiples.
    n̂ = (U×V)/norm(U×V)
    T = W - W⋅n̂ * n̂
    # find multiples of U, V that move T inside the parallelogram formed by U, V
    # (See notes from Rod Forcade in the support folder.)
    denom = (U⋅U)*(V⋅V)-(U⋅V)^2
    a = floor(Int, ((U⋅W)*(V⋅V)-(V⋅W)*(U⋅V))/denom)
    b = floor(Int, ((V⋅W)*(U⋅U)-(U⋅W)*(U⋅V))/denom)
    # Find the corner of the U-V cell closest to shifted T and shift T to be closest to origin
    W = W - a*U - b*V # Try the corner in "first quadrant"

    # Now try the other three corners, keep the shortest that is found
    norm(W) > norm(W - U) && (W = W - U)
    norm(W) > norm(W - V) && (W = W - V)
    norm(W) > norm(W - U - V) && (W = W - U - V)
    return U, V, W
end

"""
    GaussReduce(U, V)

Reduce the basis vectors {`U`, `V`} to the shortest possible basis.

# Examples
```jldoctest
julia> GaussReduce([5 8], [8 13])
([-1.0 0.0], [0.0 -1.0])
```
"""
function GaussReduce(U, V)
    maxval = max(abs.(U)...,abs.(V)...)
    if norm(U) > norm(V) U, V = V, U end
    i = 0
    while true
        V, U = U, V - round((U⋅V)/(U⋅U))*U
        i += 1
        if norm(U) > norm(V) || norm(U)≈norm(V) break; end
        isnan(norm(U)) && error("GaussReduce: input vectors were linearly independent")
        i > 50 && error("GaussReduce: Too many iterations") # failsafe to break out if not converging
    end
    return V, U
end

"""
    orthogonalityDefect(a,b,c)

Compute the orthogonality defect of three basis vectors.

# Examples
```jldoctest
julia> orthogonalityDefect([1,1,0],[1,0,1],[0,1,1])
1.4142135623730954
```
"""
function orthogonalityDefect(a,b,c)
    return prod(norm.([a,b,c]))/abs((a×b)⋅c)
end

"""
    RandUnimodMat(n)

Generate a random unimodular 2x2 matrix. `n` is a small integer (number of row and column operations).

See also: `RandLowerTri(n)`, `FibonacciMat(n)`, `DeviousMat(n)`
"""
function RandUnimodMat(n)
    mat = RandLowerTri(1)
    for i ∈ 1:n
        mat = mat*RandLowerTri(1)
        mat = mat*transpose(RandLowerTri(1))
    end
    return mat
end

"""
    RandLowerTri(n)
Generate a random 2x2 matrix of the form [1 0; 0 ±n].

See also: `RandUnimodMat(n)`, `FibonacciMat(n)`, `DeviousMat(n)`
"""
function RandLowerTri(n)
    return [1 0; rand(-n:n) 1]
end

"""
    FibonacciMat(k)
Generate a 2x2 matrix of the form [f2 f3; f1 f2] where f1, f2, f3 are consecutive Fibonacci-like numbers

See also: `RandUnimodMat(n)`, `FibonacciMat(n)`, `DeviousMat(n)`

"""
function FibonacciMat(k)
    f1 = round(Int64,1.61803398875^k/sqrt(5))
    f2 = round(Int64,1.61803398875^(k+1)/sqrt(5))
    f3 = f1 + f2
    any(i -> i < 1, [f1 f2 f3]) && error("Overflow in FibonacciMat function")
    return [f2 f3; f1 f2]
end

"""
    DeviousMat(n)

Make a unimodular 3x3 matrix that requires a large number of steps to reduce
(See email from Rod Feb 1 2020)

`n` dictates the size of the entries
"""
function DeviousMat(n)
    n < 3 && error("for DeviousMat, n > 2")
    u,v = round(Int64,(2+√3)^n/(2*√3)), round(Int64,(2+√3)^n/2)
    a,b = convert(Int64,(u+v+1)/2), -u
    c,d = a-1, v-u
    return [a b c; b d b; c b a]
end

"""
    isMinkReduced(U,V,W)

Check if the basis {`U`,`V`,`W`} is Minkoswki reduced.

"""
function isMinkReduced(U,V,W)
    if norm(U) > norm(V)+eps() println("Condition 1 failed"); return false end
    if norm(V) > norm(W)+eps() println("Condition 2 failed"); return false end
    if norm(V) > norm(U+V)+eps() println("Condition 3 failed"); return false end
    if norm(V) > norm(U-V)+eps() println("Condition 4 failed"); return false end
    if norm(W) > norm(U+W)+eps() println("Condition 5 failed"); return false end
    if norm(W) > norm(U-W)+eps() println("Condition 6 failed"); return false end
    if norm(W) > norm(V+W)+eps() println("Condition 7 failed"); return false end
    if norm(W) > norm(V-W)+eps() println("Condition 8 failed"); return false end
    if norm(W) > norm(U+V+W)+eps() println("Condition 9 failed"); return false end
    if norm(W) > norm(U-V+W)+eps() println("Condition 10 failed"); return false end
    if norm(W) > norm(U+V-W)+eps() println("Condition 11 failed"); return false end
    if norm(W) > norm(U-V-W)+eps() println("Condition 12 failed"); return false end
    return true
end

 end


text:
# NOTE: Carefully review all the NOTES here.

include(joinpath(@__DIR__, "../../imports.jl"))
include(joinpath(@__DIR__, "../../simulatedata.jl"))

function simfn(settings::Dict{Symbol, Any})
  println("pid: $(getpid())")
  println("Threads: $(Threads.nthreads())")
  println("pwd: $(pwd())")
  flush(stdout)

  println("settings:")
  for (k, v) in settings
    println("$k => $v")
  end

  # Results dir
  results_dir = settings[:results_dir]
  mkpath(results_dir)

  # Repulsive FAM penalty scale
  phi = settings[:repfam_dist_scale]
  sim_z = CytofRepFAM.Model.sim_fn_exp_decay_generator(phi)
  use_repulsive = phi > 0

  function init_state_const_data(simdat; K, L)
    deltaz_prior = TruncatedNormal(1.0, 0.1, 0.0, Inf)
    d = CytofRepFAM.Model.Data(simdat[:y_complete])
    c = CytofRepFAM.Model.defaultConstants(d, K, L,
                                      sig2_prior=InverseGamma(3, 1),
                                      delta0_prior=deltaz_prior,
                                      delta1_prior=deltaz_prior,
                                      alpha_prior=Gamma(0.1, 10.0),
                                      # NOTE:  These 3 items are important
                                      yQuantiles=[.0, .25, .5],
                                      pBounds=[.05, .8, .05],
                                      y_grid=collect(range(-10,
                                                           stop=4,
                                                           length=100)),
                                      probFlip_Z=1.0,
                                      similarity_Z=sim_z)

    # NOTE: This is new. We want to spread weights evenly.
    c.eta_prior = Dict(z => Dirichlet(L[z], 1) for z in 0:1)

    # Initialize state
    # s = CytofRepFAM.Model.smartInit(c, d)  # mclust init
    s = CytofRepFAM.Model.smartInit(c, d, modelNames="kmeans",
                                    seed=settings[:mcmcseed],
                                    iterMax=30)  # kmeans init
    # s = CytofRepFAM.Model.genInitialState(c, d,  # NOTE: random inits
    #                                       sb_ibp=false,
    #                                       allow_repeated_Z_columns=false)

    t = CytofRepFAM.Model.Tuners(d.y, c.K)
    X = CytofRepFAM.Model.eye(Float64, d.I)

    cfs = CytofRepFAM.Model.ConstantsFS(c)

    # NOTE: These priors are important
    # similar to p ~ Beta(1, 9), weakly informative
    cfs.omega_prior = Normal(-3, 1.3)
    cfs.W_star_prior = Gamma(1.0, 2) # shape, scale

    dfs = CytofRepFAM.Model.DataFS(d, X)
    sfs = CytofRepFAM.Model.StateFS{Float64}(s, dfs)
    tfs = CytofRepFAM.Model.TunersFS(t, s, X)

    return Dict(:dfs => dfs, :cfs => cfs, :sfs => sfs, :tfs => tfs,
                :simdat => simdat)
  end

  @time simdat = simulatedata1(Z=settings[:Z],
                               N=settings[:N],
                               W=settings[:W],  # Categorical.
                               sig2=[.5, .5],
                               seed=settings[:dataseed],
                               propmissingscale=.6,
                               sortLambda=true,
                               mus=Dict(0=>[-1.0], 1=>[1.0]),
                               eps_mus_dist=Uniform(-.3, .3));

  # Parameters to monitor
  monitor1 = [:theta__Z, :theta__v, :theta__alpha,
              :omega, :r, :theta__lam, :W_star, :theta__eta,
              :theta__W, :theta__delta, :theta__sig2]
  monitor2 = [:theta__y_imputed, :theta__gam]

  # MCMC Specs
  nsamps_to_thin(nsamps::Int, nmcmc::Int) = max(1, div(nmcmc, nsamps))
  nsamps = settings[:nsamps]  # Number of samples
  thin_samps = settings[:thin_samps]  # Factor to thin the primary parameters
  mcmc_iter = nsamps * thin_samps  # Number of MCMC iterations

  # LPML / DIC are computed based on `mcmc_iter` samples
  nburn = settings[:nburn]  # burn-in time

  # Configurations: priors, initial state, data, etc.
  config = init_state_const_data(simdat,
                                 K=settings[:Kmcmc],
                                 L=settings[:Lmcmc])

  # Print constants
  println("N: $(config[:dfs].data.N)")
  println("J: $(config[:dfs].data.J)")
  CytofRepFAM.Model.printConstants(config[:cfs])
  flush(stdout)

  # Fit model
  @time out = CytofRepFAM.Model.fit_fs_tp!(
    config[:sfs], config[:cfs], config[:dfs],
    nmcmc=mcmc_iter, nburn=nburn,
    batchprop=settings[:batchprop],
    prior_thin=settings[:pthin],
    temper=1.0, anneal=true, mb_update_burn_prop=0.7,
    Z_marg_lamgam=1.0,
    Z_marg_lamgam_decay_rate=100.0,
    Z_marg_lamgam_min=1.0,
    printFreq=1,
    seed=settings[:mcmcseed],
    computedden=true,
    computeDIC=true,
    computeLPML=true,
    thins=[settings[:thin_samps],
           nsamps_to_thin(10, mcmc_iter)],
    thin_dden=nsamps_to_thin(200, mcmc_iter),
    time_updates=false,
    verbose=3)

  # Dump output
  BSON.bson("$(results_dir)/output.bson", out)
  BSON.bson("$(results_dir)/simdat.bson", Dict(:simdat => config[:simdat]))

  println("Completed!")
end  # simfn


text:
#=
This is a script that renumbers the IDs for reactions in a chemical network
The script takes 2 arguments, the first of which is the file to be renumbered,
and the second is the number to start to count at (inclusive).
=#

infile = open(ARGS[1])
outfile = open("renumbered_$(ARGS[1])","w")
start  = int(ARGS[2])

println("The infile is $(ARGS[1])")
println("The outfile is renumbered_$(ARGS[1])")
println("The starting number is $(ARGS[2])")

num = start
for ln in eachline(infile)
#  println(ln)
  newln = ln
  numstr = @sprintf("%5d",num)
  outline = "$(newln[1:166])$numstr$(newln[172:end])"
  write(outfile,outline)
  num = num + 1
end

close(infile)
close(outfile)


text:
# MadNLPHSL.jl
# Created by Sungho Shin (sungho.shin@wisc.edu)

using Pkg.Artifacts, BinaryProvider

 if haskey(ENV,"MADNLP_HSL_BLAS")
    blasvendor = (ENV["MADNLP_HSL_BLAS"]=="openblas") ? :openblas : :mkl
else
    blasvendor =  artifact_hash("MKL",joinpath(@__DIR__, "..", "Artifacts.toml")) != nothing ? :mkl : :openblas
 end

const verbose = "--verbose" in ARGS
const prefix = Prefix(@__DIR__)
const so = BinaryProvider.platform_dlext()
const rpath = `-Wl,-rpath,`
const whole_archive= Sys.isapple() ? `-Wl,-all_load` : `-Wl,--whole-archive`
const no_whole_archive = Sys.isapple() ? `-Wl,-noall_load` : `-Wl,--no-whole-archive`
const libdir     = mkpath(joinpath(@__DIR__, "lib"))
const hsl_library_path = haskey(ENV,"MADNLP_HSL_LIBRARY_PATH") ? ENV["MADNLP_HSL_LIBRARY_PATH"] : ""
const hsl_source_path = haskey(ENV,"MADNLP_HSL_SOURCE_PATH") ? ENV["MADNLP_HSL_SOURCE_PATH"] : ""
const hsl_version = haskey(ENV,"MADNLP_HSL_VERSION_NUMBER") ? ENV["MADNLP_HSL_VERSION_NUMBER"] : "2015.06.23"
const FC = haskey(ENV,"MADNLP_FC") ? ENV["MADNLP_FC"] : `gfortran`
const libmetis_dir = joinpath(artifact"METIS", "lib")
const with_metis = `-L$libmetis_dir $rpath$libmetis_dir -lmetis`
const libblas_dir = joinpath(blasvendor == :mkl ? artifact"MKL" : artifact"OpenBLAS32","lib")
const with_blas = blasvendor == :mkl ? `-L$libblas_dir $rpath$libblas_dir -lmkl_rt` :
    `-L$libblas_dir $rpath$libblas_dir -lopenblas`
const installer = Sys.isapple() ? "brew install" : Sys.iswindows() ? "pacman -S" : "sudo apt install"

rm(libdir;recursive=true,force=true)
mkpath(libdir)
isvalid(cmd::Cmd)=(try run(cmd) catch e return false end; return true)

# HSL
if hsl_source_path != ""
    if isvalid(`$FC --version`)
        OC = OutputCollector[]
        cd(hsl_source_path)
        push!(OC,OutputCollector(`$FC -fopenmp -fPIC -c -O3 -o common/deps.o common/deps.f`,verbose=verbose))
        push!(OC,OutputCollector(`$FC -fopenmp -fPIC -c -O3 -o common/deps90.o common/deps90.f90`,verbose=verbose))
        wait.(OC); empty!(OC)
        push!(OC,OutputCollector(`$FC -fPIC -c -O3 -o mc19/mc19d.o mc19/mc19d.f`,verbose=verbose))
        push!(OC,OutputCollector(`$FC -fPIC -c -O3 -o ma27/ma27d.o ma27/ma27d.f`,verbose=verbose))
        push!(OC,OutputCollector(`$FC -fopenmp -fPIC -c -O3 -o ma57/ma57d.o ma57/ma57d.f`,verbose=verbose))
        wait.(OC); empty!(OC)
        push!(OC,OutputCollector(`$FC -fopenmp -fPIC -c -O3 -o hsl_ma77/hsl_ma77d.o hsl_ma77/hsl_ma77d.f90`,verbose=verbose))
        push!(OC,OutputCollector(`$FC -fopenmp -fPIC -c -O3 -o hsl_ma86/hsl_ma86d.o hsl_ma86/hsl_ma86d.f90`,verbose=verbose))
        push!(OC,OutputCollector(`$FC -fopenmp -fPIC -c -O3 -o hsl_ma97/hsl_ma97d.o hsl_ma97/hsl_ma97d.f90`,verbose=verbose))
        wait.(OC); empty!(OC)
        push!(OC,OutputCollector(`$FC -fopenmp -fPIC -c -O3 -o hsl_ma77/C/hsl_ma77d_ciface.o hsl_ma77/C/hsl_ma77d_ciface.f90`,verbose=verbose))
        push!(OC,OutputCollector(`$FC -fopenmp -fPIC -c -O3 -o hsl_ma86/C/hsl_ma86d_ciface.o hsl_ma86/C/hsl_ma86d_ciface.f90`,verbose=verbose))
        push!(OC,OutputCollector(`$FC -fopenmp -fPIC -c -O3 -o hsl_mc68/C/hsl_mc68i_ciface.o hsl_mc68/C/hsl_mc68i_ciface.f90`,verbose=verbose))
        push!(OC,OutputCollector(`$FC -fopenmp -fPIC -c -O3 -o hsl_ma97/C/hsl_ma97d_ciface.o hsl_ma97/C/hsl_ma97d_ciface.f90`,verbose=verbose))
        wait.(OC); empty!(OC)
        wait(OutputCollector(`$FC -o$(libdir)/libhsl.$so -shared -fPIC -O3 common/deps90.o common/deps.o mc19/mc19d.o ma27/ma27d.o ma57/ma57d.o hsl_ma77/hsl_ma77d.o hsl_ma77/C/hsl_ma77d_ciface.o hsl_ma86/hsl_ma86d.o hsl_ma86/C/hsl_ma86d_ciface.o hsl_mc68/C/hsl_mc68i_ciface.o hsl_ma97/hsl_ma97d.o hsl_ma97/C/hsl_ma97d_ciface.o -fopenmp $with_metis $with_blas`,verbose=verbose))
        cd("$(@__DIR__)")
        product = FileProduct(prefix,joinpath(libdir,"libhsl.$so"), :libhsl)
    end
else
    product = FileProduct(hsl_library_path, :libhsl)
end

# write deps.jl
if satisfied(product)
    @info "Building HSL (Ma27, Ma57, Ma77, Ma86, Ma97) succeeded."
    write_deps_file(joinpath(@__DIR__, "deps.jl"),Product[product], verbose=verbose)
else
    @error "Building HSL (Ma27, Ma57, Ma77, Ma86, Ma97) failed."
    write_deps_file(joinpath(@__DIR__, "deps.jl"),Product[], verbose=verbose)
end


text:
@info "Running functionality tests..."
@testset "Stochastic Programs: Functionality" begin
    @testset "SP Constructs: $name" for (sp,res,name) in problems
        tol = 1e-2
        @test optimize!(sp) == :Optimal
        @test isapprox(optimal_decision(sp), res.x̄, rtol = tol)
        @test isapprox(optimal_value(sp), res.VRP, rtol = tol)
        @test isapprox(EWS(sp), res.EWS, rtol = tol)
        @test isapprox(EVPI(sp), res.EVPI, rtol = tol)
        @test isapprox(VSS(sp), res.VSS, rtol = tol)
        @test isapprox(EV(sp), res.EV, rtol = tol)
        @test isapprox(EEV(sp), res.EEV, rtol = tol)
    end
    @testset "Inequalities: $name" for (sp,res,name) in problems
        @test EWS(sp) <= VRP(sp)
        @test VRP(sp) <= EEV(sp)
        @test VSS(sp) >= 0
        @test EVPI(sp) >= 0
        @test VSS(sp) <= EEV(sp)-EV(sp)
        @test EVPI(sp) <= EEV(sp)-EV(sp)
    end
    @testset "Deferred model creation" begin
        @test decision_length(deferred) == 0
        @test nscenarios(deferred) == 2
        @test nsubproblems(deferred) == 0
        @test optimize!(deferred) == :Optimal
        @test decision_length(deferred) == 2
        @test nscenarios(deferred) == 2
        @test nsubproblems(deferred) == 2
        @test isapprox(optimal_value(deferred), -855.83, rtol = 1e-2)
    end
    @testset "Copying: $name" for (sp,res,name) in problems
        tol = 1e-2
        sp_copy = copy(sp)
        add_scenarios!(sp_copy, scenarios(sp))
        @test nscenarios(sp_copy) == nscenarios(sp)
        generate!(sp_copy)
        @test nsubproblems(sp_copy) == nsubproblems(sp)
        @test optimize!(sp_copy) == :Optimal
        optimize!(sp)
        @test isapprox(optimal_decision(sp_copy), optimal_decision(sp), rtol = tol)
        @test isapprox(optimal_value(sp_copy), optimal_value(sp), rtol = tol)
        @test isapprox(EWS(sp_copy), EWS(sp), rtol = tol)
        @test isapprox(EVPI(sp_copy), EVPI(sp), rtol = tol)
        @test isapprox(VSS(sp_copy), VSS(sp), rtol = tol)
        @test isapprox(EV(sp_copy), EV(sp), rtol = tol)
        @test isapprox(EEV(sp_copy), EEV(sp), rtol = tol)
    end
    @testset "Sampling" begin
        sampled_sp = sample(simple_model, sampler, 100, solver=GLPKSolverLP())
        @test nscenarios(sampled_sp) == 100
        @test nsubproblems(sampled_sp) == 100
        @test isapprox(stage_probability(sampled_sp), 1.0)
        sample!(sampled_sp, sampler, 100)
        @test nscenarios(sampled_sp) == 200
        @test nsubproblems(sampled_sp) == 200
        @test isapprox(stage_probability(sampled_sp), 1.0)
    end
    @testset "Confidence intervals" begin
        glpk = GLPKSolverLP()
        try
            CI = confidence_interval(simple_model, sampler, solver = glpk, N = 200, log = false)
            @test lower(CI) <= upper(CI)
        catch end
        sol = optimize!(simple_model, sampler, solver = glpk, confidence = 0.95, tol = 1e-1, log = false)
        @test lower(confidence_interval(sol)) <= upper(confidence_interval(sol))
    end
end


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################

# Overloads
#
# Different objects that must all interact:
# 1. Number
# 2. Variable
# 3. [Generic]Norm
# 4. [Generic]AffExpr
# 5. QuadExpr
# 6. [Generic]NormExpr

# Number
# Number--Number obviously already taken care of!
# Number--Variable
(+)(lhs::Number, rhs::Variable) = AffExpr([rhs],[+1.],convert(Float64,lhs))
(-)(lhs::Number, rhs::Variable) = AffExpr([rhs],[-1.],convert(Float64,lhs))
(*)(lhs::Number, rhs::Variable) = AffExpr([rhs],[convert(Float64,lhs)], 0.)
# Number--GenericNorm
(+){P,C,V}(lhs::Number, rhs::GenericNorm{P,C,V}) = GenericNormExpr{P,C,V}(copy(rhs),  one(C), convert(C,lhs))
(-){P,C,V}(lhs::Number, rhs::GenericNorm{P,C,V}) = GenericNormExpr{P,C,V}(copy(rhs), -one(C), convert(C,lhs))
(*){P,C,V}(lhs::Number, rhs::GenericNorm{P,C,V}) = GenericNormExpr{P,C,V}(copy(rhs),     lhs, zero(GenericAffExpr{C,V}))
# Number--GenericAffExpr
(+)(lhs::Number, rhs::GenericAffExpr) = GenericAffExpr(copy(rhs.vars),copy(rhs.coeffs),lhs+rhs.constant)
(-)(lhs::Number, rhs::GenericAffExpr) = GenericAffExpr(copy(rhs.vars),    -rhs.coeffs ,lhs-rhs.constant)
(*)(lhs::Number, rhs::GenericAffExpr) = GenericAffExpr(copy(rhs.vars),[lhs*rhs.coeffs[i] for i=1:length(rhs.coeffs)],lhs*rhs.constant)
# Number--QuadExpr
(+)(lhs::Number, rhs::QuadExpr) = QuadExpr(copy(rhs.qvars1),copy(rhs.qvars2),copy(rhs.qcoeffs),lhs+rhs.aff)
(-)(lhs::Number, rhs::QuadExpr) = QuadExpr(copy(rhs.qvars1),copy(rhs.qvars2),    -rhs.qcoeffs ,lhs-rhs.aff)
(*)(lhs::Number, rhs::QuadExpr) = QuadExpr(copy(rhs.qvars1),copy(rhs.qvars2), lhs*rhs.qcoeffs ,lhs*rhs.aff)
# Number--GenericNormExpr
(+){T<:GenericNormExpr}(lhs::Number, rhs::T) = T(copy(rhs.norm),     rhs.coeff, lhs+rhs.aff)
(-){T<:GenericNormExpr}(lhs::Number, rhs::T) = T(copy(rhs.norm),    -rhs.coeff, lhs-rhs.aff)
(*){T<:GenericNormExpr}(lhs::Number, rhs::T) = T(copy(rhs.norm), lhs*rhs.coeff, lhs*rhs.aff)

# Variable (or, AbstractJuMPScalar)
(+)(lhs::AbstractJuMPScalar) = lhs
(-)(lhs::Variable) = AffExpr([lhs],[-1.0],0.0)
(*)(lhs::AbstractJuMPScalar) = lhs # make this more generic so extensions don't have to define unary multiplication for our macros
# Variable--Number
(+)(lhs::Variable, rhs::Number) = (+)( rhs,lhs)
(-)(lhs::Variable, rhs::Number) = (+)(-rhs,lhs)
(*)(lhs::Variable, rhs::Number) = (*)(rhs,lhs)
(/)(lhs::Variable, rhs::Number) = (*)(1./rhs,lhs)
# Variable--Variable
(+)(lhs::Variable, rhs::Variable) = AffExpr([lhs,rhs], [1.,+1.], 0.)
(-)(lhs::Variable, rhs::Variable) = AffExpr([lhs,rhs], [1.,-1.], 0.)
(*)(lhs::Variable, rhs::Variable) = QuadExpr([lhs],[rhs],[1.],AffExpr(Variable[],Float64[],0.))
# Variable--Norm
(+){C,V<:Variable}(lhs::Variable, rhs::GenericNorm{2,C,V}) = GenericSOCExpr{C,V}(copy(rhs),  one(C), GenericAffExpr{C,V}(lhs))
(-){C,V<:Variable}(lhs::Variable, rhs::GenericNorm{2,C,V}) = GenericSOCExpr{C,V}(copy(rhs), -one(C), GenericAffExpr{C,V}(lhs))
# Variable--AffExpr
(+){C,V<:JuMPTypes}(lhs::V, rhs::GenericAffExpr{C,V}) =
    GenericAffExpr{C,V}(vcat(rhs.vars,lhs),vcat(rhs.coeffs,one(C)), rhs.constant)
(-){C,V<:JuMPTypes}(lhs::V, rhs::GenericAffExpr{C,V}) =
    GenericAffExpr{C,V}(vcat(rhs.vars,lhs),vcat(-rhs.coeffs,one(C)),-rhs.constant)
function (*)(lhs::Variable, rhs::AffExpr)
    n = length(rhs.vars)
    if rhs.constant != 0.
        ret = QuadExpr([lhs for i=1:n],copy(rhs.vars),copy(rhs.coeffs),AffExpr([lhs], [rhs.constant], 0.))
    else
        ret = QuadExpr([lhs for i=1:n],copy(rhs.vars),copy(rhs.coeffs),zero(AffExpr))
    end
end
(/)(lhs::Variable, rhs::AffExpr) = error("Cannot divide a variable by an affine expression")
# Variable--QuadExpr
(+)(v::Variable, q::QuadExpr) = QuadExpr(copy(q.qvars1),copy(q.qvars2),copy(q.qcoeffs),v+q.aff)
(-)(v::Variable, q::QuadExpr) = QuadExpr(copy(q.qvars1),copy(q.qvars2),    -q.qcoeffs ,v-q.aff)
# Variable--GenericSOCExpr
(+){C,V<:Variable}(lhs::Variable, rhs::GenericSOCExpr{C,V}) = GenericSOCExpr{C,V}(copy(rhs.norm),  rhs.coeff, lhs+rhs.aff)
(-){C,V<:Variable}(lhs::Variable, rhs::GenericSOCExpr{C,V}) = GenericSOCExpr{C,V}(copy(rhs.norm), -rhs.coeff, lhs-rhs.aff)

# GenericNorm
(+)(lhs::GenericNorm) = lhs
(-){P,C,V}(lhs::GenericNorm{P,C,V}) = GenericNormExpr{P,C,V}(copy(lhs), -one(C), zero(GenericAffExpr{C,V}))
(*)(lhs::GenericNorm) = lhs
# GenericNorm--Number
(+){P,C,V}(lhs::GenericNorm{P,C,V},rhs::Number) =
    GenericNormExpr{P,C,V}(copy(lhs), one(C), GenericAffExpr{C,V}( rhs))
(-){P,C,V}(lhs::GenericNorm{P,C,V},rhs::Number) =
    GenericNormExpr{P,C,V}(copy(lhs), one(C), GenericAffExpr{C,V}(-rhs))
(*){P,C,V}(lhs::GenericNorm{P,C,V},rhs::Number) =
    GenericNormExpr{P,C,V}(copy(lhs), rhs, zero(GenericAffExpr{C,V}))
(/){P,C,V}(lhs::GenericNorm{P,C,V},rhs::Number) =
    GenericNormExpr{P,C,V}(copy(lhs), 1/rhs, zero(GenericAffExpr{C,V}))
# Norm--Variable
(+)(lhs::Norm,rhs::Variable) = SOCExpr(copy(lhs), 1.0, AffExpr( rhs))
(-)(lhs::Norm,rhs::Variable) = SOCExpr(copy(lhs), 1.0, AffExpr(-rhs))
# GenericNorm--GenericNorm
# GenericNorm--GenericAffExpr
(+){P,C,V}(lhs::GenericNorm{P,C,V},rhs::GenericAffExpr{C,V}) =
    GenericNormExpr{P,C,V}(copy(lhs), one(C), copy(rhs))
(-){P,C,V}(lhs::GenericNorm{P,C,V},rhs::GenericAffExpr{C,V}) =
    GenericNormExpr{P,C,V}(copy(lhs), one(C),     -rhs)
# Norm--QuadExpr
# Norm--SOCExpr

# GenericAffExpr
(+)(lhs::GenericAffExpr) = lhs
(-)(lhs::GenericAffExpr) = GenericAffExpr(lhs.vars, -lhs.coeffs, -lhs.constant)
(*)(lhs::GenericAffExpr) = lhs
# AffExpr--Number
(+)(lhs::GenericAffExpr, rhs::Number) = (+)(+rhs,lhs)
(-)(lhs::GenericAffExpr, rhs::Number) = (+)(-rhs,lhs)
(*)(lhs::GenericAffExpr, rhs::Number) = (*)(rhs,lhs)
(/)(lhs::GenericAffExpr, rhs::Number) = (*)(1.0/rhs,lhs)
function (^)(lhs::Union{Variable,AffExpr}, rhs::Integer)
    if rhs == 2
        return lhs*lhs
    elseif rhs == 1
        return QuadExpr(lhs)
    elseif rhs == 0
        return QuadExpr(1)
    else
        error("Only exponents of 0, 1, or 2 are currently supported. Are you trying to build a nonlinear problem? Make sure you use @NLconstraint/@NLobjective.")
    end
end
(^)(lhs::Union{Variable,AffExpr}, rhs::Number) = error("Only exponents of 0, 1, or 2 are currently supported. Are you trying to build a nonlinear problem? Make sure you use @NLconstraint/@NLobjective.")
# AffExpr--Variable
(+){C,V<:JuMPTypes}(lhs::GenericAffExpr{C,V}, rhs::V) = (+)(rhs,lhs)
(-){C,V<:JuMPTypes}(lhs::GenericAffExpr{C,V}, rhs::V) = GenericAffExpr{C,V}(vcat(lhs.vars,rhs),vcat(lhs.coeffs,-one(C)),lhs.constant)
(*)(lhs::AffExpr, rhs::Variable) = (*)(rhs,lhs)
(/)(lhs::AffExpr, rhs::Variable) = error("Cannot divide affine expression by a variable")
# AffExpr--Norm
(+){P,C,V}(lhs::GenericAffExpr{C,V}, rhs::GenericNorm{P,C,V}) = GenericNormExpr{P,C,V}(copy(rhs),  one(C), copy(lhs))
(-){P,C,V}(lhs::GenericAffExpr{C,V}, rhs::GenericNorm{P,C,V}) = GenericNormExpr{P,C,V}(copy(rhs), -one(C), copy(lhs))
# AffExpr--AffExpr
(+){C,V<:JuMPTypes}(lhs::GenericAffExpr{C,V}, rhs::GenericAffExpr{C,V}) = (operator_warn(lhs,rhs); GenericAffExpr(vcat(lhs.vars,rhs.vars),vcat(lhs.coeffs, rhs.coeffs),lhs.constant+rhs.constant))
(-){C,V<:JuMPTypes}(lhs::GenericAffExpr{C,V}, rhs::GenericAffExpr{C,V}) = GenericAffExpr(vcat(lhs.vars,rhs.vars),vcat(lhs.coeffs,-rhs.coeffs),lhs.constant-rhs.constant)
function (*)(lhs::AffExpr, rhs::AffExpr)
    ret = QuadExpr(Variable[],Variable[],Float64[],AffExpr(Variable[],Float64[],0.))

    # Quadratic terms
    n = length(lhs.coeffs)
    m = length(rhs.coeffs)
    sizehint!(ret.qvars1, n*m)
    sizehint!(ret.qvars2, n*m)
    sizehint!(ret.qcoeffs, n*m)
    for i = 1:n
        for j = 1:m
            push!(ret.qvars1,  lhs.vars[i])
            push!(ret.qvars2,  rhs.vars[j])
            push!(ret.qcoeffs, lhs.coeffs[i]*rhs.coeffs[j])
        end
    end

    # Try to preallocate space for aff
    if lhs.constant != 0 && rhs.constant != 0
        sizehint!(ret.aff.vars, n+m)
        sizehint!(ret.aff.coeffs, n+m)
    elseif lhs.constant != 0
        sizehint!(ret.aff.vars, n)
        sizehint!(ret.aff.coeffs, n)
    elseif rhs.constant != 0
        sizehint!(ret.aff.vars, m)
        sizehint!(ret.aff.coeffs, m)
    end

    # [LHS constant] * RHS
    if lhs.constant != 0
        c = lhs.constant
        for j = 1:m
            push!(ret.aff.vars,   rhs.vars[j])
            push!(ret.aff.coeffs, rhs.coeffs[j] * c)
        end
        ret.aff.constant += c * rhs.constant
    end

    # [RHS constant] * LHS
    if rhs.constant != 0
        c = rhs.constant
        for i = 1:n
            push!(ret.aff.vars,   lhs.vars[i])
            push!(ret.aff.coeffs, lhs.coeffs[i] * c)
        end
        # Don't do the following line
        #ret.aff.constant += c * lhs.constant
        # If lhs.constant is 0, its a waste of time
        # If lhs.constant is non-zero, its already done
    end

    return ret
end
# AffExpr--QuadExpr
(+)(a::AffExpr, q::QuadExpr) = QuadExpr(copy(q.qvars1),copy(q.qvars2),copy(q.qcoeffs),a+q.aff)
(-)(a::AffExpr, q::QuadExpr) = QuadExpr(copy(q.qvars1),copy(q.qvars2),    -q.qcoeffs ,a-q.aff)
# AffExpr--SOCExpr
(+){C,V}(lhs::GenericAffExpr{C,V}, rhs::GenericSOCExpr{C,V}) = GenericSOCExpr{C,V}(copy(rhs.norm),  copy(rhs.coeff), lhs+rhs.aff)
(-){C,V}(lhs::GenericAffExpr{C,V}, rhs::GenericSOCExpr{C,V}) = GenericSOCExpr{C,V}(copy(rhs.norm), -copy(rhs.coeff), lhs-rhs.aff)

# QuadExpr
(+)(lhs::QuadExpr) = lhs
(-)(lhs::QuadExpr) = 0.0-lhs
(*)(lhs::QuadExpr) = lhs
# QuadExpr--Number
(+)(lhs::QuadExpr, rhs::Number) = (+)(+rhs,lhs)
(-)(lhs::QuadExpr, rhs::Number) = (+)(-rhs,lhs)
(*)(lhs::QuadExpr, rhs::Number) = (*)(rhs,lhs)
(/)(lhs::QuadExpr, rhs::Number) = (*)(1.0/rhs,lhs)
# QuadExpr--Variable
(+)(q::QuadExpr, v::Variable) = (+)(v,q)
(-)(q::QuadExpr, v::Variable) = QuadExpr(copy(q.qvars1),copy(q.qvars2),copy(q.qcoeffs),q.aff-v)
(*)(q::QuadExpr, v::Variable) = error("Cannot multiply a quadratic expression by a variable")
(/)(q::QuadExpr, v::Variable) = error("Cannot divide a quadratic expression by a variable")
# QuadExpr--AffExpr
(+)(q::QuadExpr, a::AffExpr) = QuadExpr(copy(q.qvars1),copy(q.qvars2),copy(q.qcoeffs),q.aff+a)
(-)(q::QuadExpr, a::AffExpr) = QuadExpr(copy(q.qvars1),copy(q.qvars2),copy(q.qcoeffs),q.aff-a)
(*)(q::QuadExpr, a::AffExpr) = error("Cannot multiply a quadratic expression by an aff. expression")
(/)(q::QuadExpr, a::AffExpr) = error("Cannot divide a quadratic expression by an aff. expression")
# QuadExpr--QuadExpr
(+)(q1::QuadExpr, q2::QuadExpr) = QuadExpr( vcat(q1.qvars1, q2.qvars1),     vcat(q1.qvars2, q2.qvars2),
                                            vcat(q1.qcoeffs, q2.qcoeffs),   q1.aff + q2.aff)
(-)(q1::QuadExpr, q2::QuadExpr) = QuadExpr( vcat(q1.qvars1, q2.qvars1),     vcat(q1.qvars2, q2.qvars2),
                                            vcat(q1.qcoeffs, -q2.qcoeffs),  q1.aff - q2.aff)

# GenericNormExpr
(+)(lhs::GenericNormExpr) = lhs
(-){T<:GenericNormExpr}(lhs::T) = T(copy(lhs.norm), -lhs.coeff, -lhs.aff)
(*)(lhs::GenericNormExpr) = lhs
# GenericNormExpr--Number
(+){T<:GenericNormExpr}(lhs::T,rhs::Number) = T(copy(lhs.norm), lhs.coeff, lhs.aff+rhs)
(-){T<:GenericNormExpr}(lhs::T,rhs::Number) = T(copy(lhs.norm), lhs.coeff, lhs.aff-rhs)
(*){T<:GenericNormExpr}(lhs::T,rhs::Number) = T(copy(lhs.norm), lhs.coeff*rhs, lhs.aff*rhs)
(/){T<:GenericNormExpr}(lhs::T,rhs::Number) = T(copy(lhs.norm), lhs.coeff/rhs, lhs.aff/rhs)
# SOCExpr--Variable
(+)(lhs::SOCExpr,rhs::Variable) = SOCExpr(copy(lhs.norm), lhs.coeff, lhs.aff+rhs)
(-)(lhs::SOCExpr,rhs::Variable) = SOCExpr(copy(lhs.norm), lhs.coeff, lhs.aff-rhs)
# GenericNormExpr--GenericNorm
# GenericNormExpr--GenericAffExpr
(+){P,C,V}(lhs::GenericNormExpr{P,C,V},rhs::GenericAffExpr{C,V}) =
    GenericNormExpr{P,C,V}(copy(lhs.norm), lhs.coeff, lhs.aff+rhs)
(-){P,C,V}(lhs::GenericNormExpr{P,C,V},rhs::GenericAffExpr{C,V}) =
    GenericNormExpr{P,C,V}(copy(lhs.norm), lhs.coeff, lhs.aff-rhs)
# SOCExpr--QuadExpr
# SOCExpr--SOCExpr

(==)(lhs::AffExpr,rhs::AffExpr) = (lhs.vars == rhs.vars) && (lhs.coeffs == rhs.coeffs) && (lhs.constant == rhs.constant)
(==)(lhs::QuadExpr,rhs::QuadExpr) = (lhs.qvars1 == rhs.qvars1) && (lhs.qvars2 == rhs.qvars2) && (lhs.qcoeffs == rhs.qcoeffs) && (lhs.aff == rhs.aff)

#############################################################################
# Helpers to initialize memory for AffExpr/QuadExpr
#############################################################################

_sizehint_expr!(q::GenericAffExpr, n::Int) = begin
        sizehint!(q.vars,   length(q.vars)   + n)
        sizehint!(q.coeffs, length(q.coeffs) + n)
        nothing
end

_sizehint_expr!(q::GenericQuadExpr, n::Int) = begin
        sizehint!(q.qvars1,  length(q.qvars1)  + n)
        sizehint!(q.qvars2,  length(q.qvars2)  + n)
        sizehint!(q.qcoeffs, length(q.qcoeffs) + n)
        _sizehint_expr!(q.aff, n)
        nothing
end
_sizehint_expr!(q, n) = nothing

#############################################################################
# High-level operators
# Currently supported
#  - sum
#  - dot
#############################################################################

Base.sum(j::JuMPArray) = sum(j.innerArray)
Base.sum(j::JuMPDict)  = sum(values(j.tupledict))
Base.sum(j::JuMPArray{Variable}) = AffExpr(vec(j.innerArray), ones(length(j.innerArray)), 0.0)
Base.sum(j::JuMPDict{Variable})  = AffExpr(collect(values(j.tupledict)), ones(length(j.tupledict)), 0.0)
Base.sum(j::Array{Variable}) = AffExpr(vec(j), ones(length(j)), 0.0)
function Base.sum{T<:GenericAffExpr}(affs::Array{T})
    new_aff = zero(T)
    for aff in affs
        append!(new_aff, aff)
    end
    return new_aff
end

import Base.vecdot

_dot_depr() = warn("dot is deprecated for multidimensional arrays. Use vecdot instead.")

# Base Julia's generic fallback vecdot requires that dot be defined
# for scalars, so instead of defining them one-by-one, we will
# fallback to the multiplication operator
Base.dot(lhs::JuMPTypes, rhs::JuMPTypes) = lhs*rhs
Base.dot(lhs::JuMPTypes, rhs::Number)    = lhs*rhs
Base.dot(lhs::Number,    rhs::JuMPTypes) = lhs*rhs

Base.dot{T,S,N}(lhs::Array{T,N}, rhs::JuMPArray{S,N})    = begin _dot_depr(); vecdot(lhs,rhs); end
Base.dot{T,S,N}(lhs::JuMPArray{T,N},rhs::Array{S,N})     = begin _dot_depr(); vecdot(lhs,rhs); end
Base.dot{T,S,N}(lhs::JuMPArray{T,N},rhs::JuMPArray{S,N}) = begin _dot_depr(); vecdot(lhs,rhs); end
Base.dot{T<:JuMPTypes,S,N}(lhs::Array{T,N}, rhs::Array{S,N}) = begin _dot_depr(); vecdot(lhs,rhs); end
Base.dot{T<:JuMPTypes,S<:JuMPTypes,N}(lhs::Array{T,N}, rhs::Array{S,N}) = begin _dot_depr(); vecdot(lhs,rhs); end
Base.dot{T,S<:JuMPTypes,N}(lhs::Array{T,N}, rhs::Array{S,N}) = begin _dot_depr(); vecdot(lhs,rhs); end

Base.dot{T<:JuMPTypes,S<:JuMPTypes}(lhs::Vector{T},rhs::Vector{S}) = _dot(lhs,rhs)
Base.dot{T<:JuMPTypes,S}(lhs::Vector{T},rhs::Vector{S}) = _dot(lhs,rhs)
Base.dot{T,S<:JuMPTypes}(lhs::Vector{T},rhs::Vector{S}) = _dot(lhs,rhs)

# TODO: qualify Base.vecdot once v0.3 support is dropped
vecdot{T<:JuMPTypes,S,N}(lhs::Array{T,N},rhs::Array{S,N}) = _dot(lhs,rhs)
vecdot{T<:JuMPTypes,S<:JuMPTypes,N}(lhs::Array{T,N},rhs::Array{S,N}) = _dot(lhs,rhs)
vecdot{T,S<:JuMPTypes,N}(lhs::Array{T,N},rhs::Array{S,N}) = _dot(lhs,rhs)

function _dot{T,S}(lhs::Array{T}, rhs::Array{S})
    size(lhs) == size(rhs) || error("Incompatible dimensions")
    ret = zero(one(T)*one(S))
    for (x,y) in zip(lhs,rhs)
        ret = addtoexpr(ret, x, y)
    end
    ret
end

###############################################################################
# A bunch of operator junk to make matrix multiplication and friends act
# reasonably sane with JuMP types

Base.promote_rule{R<:Real}(::Type{Variable},::Type{R}       ) = AffExpr
Base.promote_rule(         ::Type{Variable},::Type{AffExpr} ) = AffExpr
Base.promote_rule(         ::Type{Variable},::Type{QuadExpr}) = QuadExpr
Base.promote_rule{R<:Real}(::Type{AffExpr}, ::Type{R}       ) = AffExpr
Base.promote_rule(         ::Type{AffExpr}, ::Type{QuadExpr}) = QuadExpr
Base.promote_rule{R<:Real}(::Type{QuadExpr},::Type{R}       ) = QuadExpr

_throw_transpose_error() = error("Transpose not currently implemented for JuMPArrays with arbitrary index sets.")
Base.transpose(x::AbstractJuMPScalar) = x
Base.transpose( x::JuMPArray) = _throw_transpose_error()
Base.ctranspose(x::JuMPArray) = _throw_transpose_error()

# Can remove the following code once == overloading is removed

function Compat.issymmetric{T<:JuMPTypes}(x::Matrix{T})
    (n = size(x,1)) == size(x,2) || return false
    for i in 1:n, j in (i+1):n
        isequal(x[i,j], x[j,i]) || return false
    end
    true
end

# Special-case because the the base version wants to do fill!(::Array{Variable}, zero(AffExpr))
Base.diagm(x::Vector{Variable}) = diagm(convert(Vector{AffExpr}, x))

###############
# The _multiply!(buf,y,z) adds the results of y*z into the buffer buf. No bounds/size
# checks are performed; it is expected that the caller has done this, has ensured
# that the eltype of buf is appropriate, and has zeroed the elements of buf (if desired).

function _multiply!{T<:JuMPTypes}(ret::Array{T}, lhs::Array, rhs::Array)
    m, n = size(lhs,1), size(lhs,2)
    r, s = size(rhs,1), size(rhs,2)
    for i in 1:m, j in 1:s
        q = ret[i,j]
        _sizehint_expr!(q, n)
        for k in 1:n
            tmp = convert(T, lhs[i,k]*rhs[k,j])
            append!(q, tmp)
        end
    end
    ret
end

function _multiply!{T<:Union{GenericAffExpr,GenericQuadExpr}}(ret::Array{T}, lhs::SparseMatrixCSC, rhs::Array)
    nzv = lhs.nzval
    rv  = lhs.rowval
    for col in 1:lhs.n
        for k in 1:size(ret, 2)
            for j in lhs.colptr[col]:(lhs.colptr[col+1]-1)
                append!(ret[rv[j], k], nzv[j] * rhs[col,k])
            end
        end
    end
    ret
end

function _multiply!{T<:Union{GenericAffExpr,GenericQuadExpr}}(ret::Array{T}, lhs::Matrix, rhs::SparseMatrixCSC)
    rowval = rhs.rowval
    nzval  = rhs.nzval
    for multivec_row in 1:size(lhs,1)
        for col in 1:rhs.n
            idxset = rhs.colptr[col]:(rhs.colptr[col+1]-1)
            q = ret[multivec_row, col]
            _sizehint_expr!(q, length(idxset))
            for k in idxset
                append!(q, lhs[multivec_row, rowval[k]] * nzval[k])
            end
        end
    end
    ret
end

# TODO: implement sparse * sparse code as in base/sparse/linalg.jl (spmatmul)
_multiply!{T<:JuMPTypes}(ret::AbstractArray{T}, lhs::SparseMatrixCSC, rhs::SparseMatrixCSC) = _multiply!(ret, lhs, full(rhs))

_multiply!(ret, lhs, rhs) = A_mul_B!(ret, lhs, ret)

(*){T<:JuMPTypes}(             A::Union{Matrix{T},SparseMatrixCSC{T}}, x::Union{Matrix,   Vector,   SparseMatrixCSC})    = _matmul(A, x)
(*){T<:JuMPTypes,R<:JuMPTypes}(A::Union{Matrix{T},SparseMatrixCSC{T}}, x::Union{Matrix{R},Vector{R},SparseMatrixCSC{R}}) = _matmul(A, x)
(*){T<:JuMPTypes}(             A::Union{Matrix,   SparseMatrixCSC},    x::Union{Matrix{T},Vector{T},SparseMatrixCSC{T}}) = _matmul(A, x)

function _matmul(A, x)
    m, n = size(A,1), size(A,2)
    r, s = size(x,1), size(x,2)
    n == r || error("Incompatible sizes")
    ret = _return_array(A, x)
    _multiply!(ret, A, x)
    ret
end

_multiply_type(R,S) = typeof(one(R) * one(S))

# Don't do size checks here in _return_array, defer that to (*)
_return_array{R,S}(A::AbstractMatrix{R}, x::AbstractVector{S}) = _fillwithzeros(Array(_multiply_type(R,S), size(A,1)))
_return_array{R,S}(A::AbstractMatrix{R}, x::AbstractMatrix{S}) = _fillwithzeros(Array(_multiply_type(R,S), size(A,1), size(x,2)))

# helper so we don't fill the buffer array with the same object
function _fillwithzeros{T}(arr::Array{T})
    for I in eachindex(arr)
        arr[I] = zero(T)
    end
    arr
end

# Let's be conservative and only define arithmetic for the basic types
typealias ArrayOrSparseMat{T} Union{Array{T}, SparseMatrixCSC{T}}

for op in [:+, :-]; @eval begin
    function $op{T<:JuMPTypes}(lhs::Number,rhs::ArrayOrSparseMat{T})
        ret = Array(typeof($op(lhs, zero(T))), size(rhs))
        for I in eachindex(ret)
            ret[I] = $op(lhs, rhs[I])
        end
        ret
    end
    function $op{T<:JuMPTypes}(lhs::ArrayOrSparseMat{T},rhs::Number)
        ret = Array(typeof($op(zero(T), rhs)), size(lhs))
        for I in eachindex(ret)
            ret[I] = $op(lhs[I], rhs)
        end
        ret
    end
    function $op{T<:JuMPTypes,S}(lhs::T,rhs::ArrayOrSparseMat{S})
        ret = Array(typeof($op(lhs, zero(S))), size(rhs))
        for I in eachindex(ret)
            ret[I] = $op(lhs, rhs[I])
        end
        ret
    end
    function $op{T<:JuMPTypes,S}(lhs::ArrayOrSparseMat{S},rhs::T)
        ret = Array(typeof($op(zero(S), rhs)), size(lhs))
        for I in eachindex(ret)
            ret[I] = $op(lhs[I], rhs)
        end
        ret
    end
end; end

for op in [:*, :/]; @eval begin
    function $op{T<:JuMPTypes}(lhs::Number,rhs::Array{T})
        ret = Array(typeof($op(lhs, zero(T))), size(rhs))
        for I in eachindex(ret)
            ret[I] = $op(lhs, rhs[I])
        end
        ret
    end
    function $op{T<:JuMPTypes}(lhs::Array{T},rhs::Number)
        ret = Array(typeof($op(zero(T), rhs)), size(lhs))
        for I in eachindex(ret)
            ret[I] = $op(lhs[I], rhs)
        end
        ret
    end
    function $op{T<:JuMPTypes,S}(lhs::T,rhs::Array{S})
        ret = Array(typeof($op(lhs, zero(S))), size(rhs))
        for I in eachindex(ret)
            ret[I] = $op(lhs, rhs[I])
        end
        ret
    end
    function $op{T<:JuMPTypes,S}(lhs::Array{S},rhs::T)
        ret = Array(typeof($op(zero(S), rhs)), size(lhs))
        for I in eachindex(ret)
            ret[I] = $op(lhs[I], rhs)
        end
        ret
    end
end; end

# Special-case sparse matrix scalar multiplication/division
(*){T<:JuMPTypes}(lhs::Number, rhs::SparseMatrixCSC{T}) =
    SparseMatrixCSC(rhs.m, rhs.n, copy(rhs.colptr), copy(rhs.rowval), lhs .* rhs.nzval)
(*)(lhs::JuMPTypes, rhs::SparseMatrixCSC) =
    SparseMatrixCSC(rhs.m, rhs.n, copy(rhs.colptr), copy(rhs.rowval), lhs .* rhs.nzval)
(*){T<:JuMPTypes}(lhs::SparseMatrixCSC{T}, rhs::Number) =
    SparseMatrixCSC(lhs.m, lhs.n, copy(lhs.colptr), copy(lhs.rowval), lhs.nzval .* rhs)
(*)(lhs::SparseMatrixCSC, rhs::JuMPTypes) =
    SparseMatrixCSC(lhs.m, lhs.n, copy(lhs.colptr), copy(lhs.rowval), lhs.nzval .* rhs)
(/){T<:JuMPTypes}(lhs::SparseMatrixCSC{T}, rhs::Number) =
    SparseMatrixCSC(lhs.m, lhs.n, copy(lhs.colptr), copy(lhs.rowval), lhs.nzval ./ rhs)

# The following are primarily there for internal use in the macro code for @constraint
for op in [:(+), :(-)]; @eval begin
    function $op(lhs::Array{Variable},rhs::Array{Variable})
        (sz = size(lhs)) == size(rhs) || error("Incompatible sizes for $op: $sz $op $(size(rhs))")
        ret = Array(AffExpr, sz)
        for I in eachindex(ret)
            ret[I] = $op(lhs[I], rhs[I])
        end
        ret
    end
end; end

for (dotop,op) in [(:.+,:+), (:.-,:-), (:.*,:*), (:./,:/)]
    @eval begin
        $dotop(lhs::Number,rhs::JuMPTypes) = $op(lhs,rhs)
        $dotop(lhs::JuMPTypes,rhs::Number) = $op(lhs,rhs)
    end
    for (T1,T2) in [(:JuMPTypes,:Number),(:JuMPTypes,:JuMPTypes),(:Number,:JuMPTypes)]
        # Need these looks over S1,S2 for v0.3 because Union{Array,SparseMatrix}
        # gives ambiguity warnings
        for S1 in (:Array,:SparseMatrixCSC)
            @eval $dotop{S<:$T1}(lhs::$S1{S},rhs::$T2) = $op(lhs,rhs)
        end
        for (S1,S2) in [(:Array,:Array),(:Array,:SparseMatrixCSC),(:SparseMatrixCSC,:Array)]
            @eval begin
                function $dotop{S<:$T1,T<:$T2}(lhs::$S1{S},rhs::$S2{T})
                    size(lhs) == size(rhs) || error("Incompatible dimensions")
                    arr = Array(typeof($op(zero(S),zero(T))), size(rhs))
                    @inbounds for i in eachindex(lhs)
                        arr[i] = $op(lhs[i], rhs[i])
                    end
                    arr
                end
            end
        end
        for S2 in (:Array,:SparseMatrixCSC)
            @eval $dotop{T<:$T2}(lhs::$T1,rhs::$S2{T}) = $op(lhs,rhs)
        end
    end
end

(+){T<:JuMPTypes}(x::Array{T}) = x
function (-){T<:JuMPTypes}(x::Array{T})
    ret = similar(x, typeof(-one(T)))
    for I in eachindex(ret)
        ret[I] = -x[I]
    end
    ret
end
(*){T<:JuMPTypes}(x::Array{T}) = x

###############################################################################
# Add nonlinear function fallbacks for JuMP built-in types
const op_hint = "Are you trying to build a nonlinear problem? Make sure you use @NLconstraint/@NLobjective."
for (func,_) in Calculus.symbolic_derivatives_1arg(), typ in [:Variable,:AffExpr,:QuadExpr]
    errstr = "$func is not defined for type $typ. $op_hint"
    @eval Base.$(func)(::$typ) = error($errstr)
end

*{T<:QuadExpr,S<:Union{Variable,AffExpr,QuadExpr}}(::T,::S) =
    error( "*(::$T,::$S) is not defined. $op_hint")
(*)(lhs::QuadExpr, rhs::QuadExpr) =
    error( "*(::QuadExpr,::QuadExpr) is not defined. $op_hint")
*{T<:QuadExpr,S<:Union{Variable,AffExpr,QuadExpr}}(::S,::T) =
    error( "*(::$S,::$T) is not defined. $op_hint")
/{S<:Union{Number,Variable,AffExpr,QuadExpr},T<:Union{Variable,AffExpr,QuadExpr}}(::S,::T) =
    error( "/(::$S,::$T) is not defined. $op_hint")


text:
module ImpvolOutput
	using JLD2
	using FileIO
	using CSV
	using Missings
	using Plots
	using HypothesisTests
	using Distributions
	using DataFrames
	using Logging
	include("calibration_utils.jl")
	Logging.configure(level=DEBUG)

	function read_results(path = "experiments/baseline/actual/results.jld2")
		file = jldopen(path, "r")
		return file["results"]
	end

	function sort_results(results)
		return sort(collect(results), by = x -> x[1])
	end

	function list_keys(results)
		for (key, value) in results[1][2]
			println(key)
		end
	end

	function make_series(results, key = :real_GDP)
		series = zeros(size(results[1][2][key])[1], size(results[1][2][key])[2], size(results[1][2][key])[3], length(results))
		for t in 1:length(results)
			# Only the first element in the shock dimension is interesting, the rest are there only for optimizaton purposes used in the algorithm
			series[:,:,:,t] = results[t][2][key][:,:,:,1]
		end
		return series
	end

	function calculate_volatilities(x, parameters, bool_detrend::Bool, range=:)
		if bool_detrend
			x_c, x_t = DetrendUtilities.detrend(log.(x),parameters[:bp_weights])
		else
			x_c = log.(x)
		end

		# time is the last dimension
		return var(x_c[:,:,:,range],ndims(x_c))
	end

	function plot_model_vs_data(plot_data::Tuple, title::String)
		# data and model are expected to be of dimension N x T
		# length of label should be equal N
		data = plot_data[1]
		model = plot_data[2]
		label = plot_data[3]

		data = permutedims(data,(2,1))
		model = permutedims(model,(2,1))

		size(data,2) == size(model,2) || error("You can only compare matching time series between the model outcome and data")

		colors = distinguishable_colors(size(data,2))

		ENV["GKSwstype"] = "gksqt"
		fig = Plots.plot()
		for i in 1:size(data,2)
			Plots.plot!([data[:,i] model[:,i]], color = colors[i,1], ls = [:solid :dash], label = [label[i,1] ""], title = title)
		end
		return fig
	end

	function plot_data(path = "experiments/baseline/actual/results.jld2", key = :real_GDP, country_range = ":")
		# output: "data, model, label", as an input for the function 'plot_model_vs_data'
		data = load("data/impvol_data.jld2")
		gdp_d = squeeze(sum(data["va"],3), (1,3))
		gdp_d = Float64.(collect(Missings.replace(gdp_d, NaN)))
		gdp_d = gdp_d[eval(parse(country_range)),:]

		pwt = squeeze(data["pwt"], (1,3))
		pwt = Float64.(collect(Missings.replace(pwt, NaN)))
		pwt = pwt[eval(parse(country_range)),:]

		cpi = CSV.read("data/cpi.csv", header = true)
		cpi = permutedims(convert(Array, cpi)[:,2:end], (2,1))
		cpi = cpi ./ cat(2, cpi[:,24]) # 24 = 1995-base
		cpi = Float64.(collect(Missings.replace(cpi, NaN)))
		cpi = cpi[eval(parse(country_range)),:]

		xr = CSV.read("data/exchange_rates.csv", header = true)
		xr = permutedims(convert(Array, xr)[:,2:end], (2,1))
		xr = Float64.(collect(Missings.replace(xr, NaN)))
		xr = xr[eval(parse(country_range)),:]

		country_names = CSV.read("data/country_name.txt", header = false, types = [String])
		country_names = country_names[:1][eval(parse(country_range))]

		results = sort_results(read_results(path))
		gdp_m = make_series(results, key)
		gdp_m = sum(gdp_m[1,eval(parse(country_range)),:,:],2)
		gdp_m = squeeze(gdp_m, 2)

		gdp_d = gdp_d ./ cpi .* xr
		#gdp_d = gdp_d ./ pwt

		# normalization: model(1972) = data(1972)
		gdp_m = gdp_m .* (gdp_d[:,1] ./ gdp_m[:,1])

		return log.(gdp_d), log.(gdp_m), country_names
	end

	################ Running the 'plot_data' function ################
	# include("output.jl")
	# plt_dta = ImpvolOutput.plot_data() # with the desired arguments if needed
	# Base.print_matrix(IOContext(STDOUT, :limit => false), plt_dta[3]) # to check the countries' position
	# idx = zeros(length(plt_dta[3]))
	# idx[6], idx[10], idx[11], idx[13], idx[15], idx[16], idx[24], idx[25] = 1, 1, 1, 1, 1, 1, 1, 1
	# idx = idx .== 1
	# ImpvolOutput.plot_model_vs_data((plt_dta[1][idx,:], plt_dta[2][idx,:], plt_dta[3][idx]), "GDP")
	# # solid line = data, dashed line = model
	##################################################################

	function write_results(parameters, rootpath = "experiments/baseline/", key = :real_GDP, bool_detrend = true, dirsdown = 1, pattern = r"jld2$")
		# Volatility of the desired variable found in the last folders of depth 'dirsdown' from the 'rootpath' is calculated by this function

		# Create 'stats' array to store volatilities
		stats = DataFrame([String, Real, Real, Real, Real, Real, Real, Real, Real, Real], [:country_names, :actual, :kappa1972, :nosectoral, :nosectoral_kappa1972, :onlyglobal, :onlyglobal_kappa1972, :trade_barriers, :diversification, :specialization], 25)

		stats[:country_names] = CSV.read("data/country_name.txt", header = false, types = [String], nullable=false)[1]

		# Fill 'stats' array
		for (root, dirs, files) in walkdir(dirname(rootpath))
			dir_depth = length(matchall(r"/", root)) - 1

			if dir_depth == dirsdown

				for file in files
					if ismatch(pattern, file)
						sectoral_series = make_series(sort_results(read_results(joinpath(root, file))), key)

						# Summation over the sectors
						series = cat(ndims(sectoral_series), sum(sectoral_series, 3))
						stats[eval(parse(":" * SubString(root,findin(root,'/')[end] + 1, length(root))))] = squeeze(calculate_volatilities(series, parameters, bool_detrend), (1,3,4))
					end
				end
			end
		end

		Logging.debug(stats)

		# Trade barriers
		stats[:trade_barriers] = 100 * (stats[:actual] - stats[:kappa1972]) ./ stats[:kappa1972]

		# Diversification
		stats[:diversification] = 100 * (stats[:nosectoral] - stats[:nosectoral_kappa1972]) ./ stats[:kappa1972]

		# Specialization
		stats[:specialization] = 100 * (stats[:actual] - stats[:kappa1972] - stats[:nosectoral] + stats[:nosectoral_kappa1972]) ./ stats[:kappa1972]

		CSV.write(rootpath * "/output_table.csv", stats)
	end

	################ Running the 'write_results' function ################
	# include("output.jl")
	# parameters = Dict{Symbol, Any}()
	# parameters[:N] = 25
	# parameters[:bp_weights] = [0.774074394803123, -0.201004684236153, -0.135080548288772,-0.050951964876636]
	# stats = ImpvolOutput.write_results(parameters)
	######################################################################

end


text:
using LatPhysPlottingSVG
using Test

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


text:
# experiment_utils
#
# A set of utility functions for running experiments.
# These are useful for removing ugly/repetitive code.

using DataFrames
using DataStructures
using JSON

# This is a utility function for parsing commandline arguments to a script.
# The doc uses the following inputs for this explanation: -s7 or --seed=7
# @arguments:
# arg:      the entire commandline argument to be parsed (e.g. -s7 or --seed=7)
# shortopt: the root for short options (in this example its 's')
# longopt:  the root for long option ('seed' in this case)
# @returns: (?string) The matching parsed argument as a string (e.g. '7') else Nothing
function matchcommandlinearg(arg::AbstractString, shortopt::AbstractString, longopt::AbstractString)
    shortmatch = match(Regex("^-$(shortopt)(.*)"), arg)
    if shortmatch != nothing
        return shortmatch.captures[1]
    end
    longmatch = match(Regex("^--$(longopt)=(.*)"), arg)
    if longmatch != nothing
        return longmatch.captures[1]
    end
end

# This is a utility function for saving a Dictionary of mixed types (instancedata)
# as a json in the given project. It constructs the pathname for the saved
# data as follows:
#
# ./results/<dir>/data/<prefix>_kwarg1=kwval1_..._kwargn=kwvaln.json
#
# where kwarg1=kwval1, ..., kwargn=kwvaln are all the keyword arguments passed in
# If any of the directories doesn't exists, it just creates that directory.
#
# CAREFUL! This will always construct the directories for the current dir.
function save_json(instancedata::Dict; prefix="julia", dir="misc", kwargs...)
    # first make sure the directory exists
    relpath_dir = joinpath(".", "results", dir, "data")
    if !isdir(relpath_dir)
        mkpath(relpath_dir)
    end
    # now build up the filename
    running_filename = prefix
    for (chunkkey, chunkvalue) in kwargs
        chunk = string(chunkkey, "=", chunkvalue)
        running_filename = string(running_filename, "_", chunk)
    end
    json_filename = string(running_filename, ".json")
    # and lets write it
    relpath_json_filename = joinpath(relpath_dir, json_filename);
    fh = open(relpath_json_filename, "w");
    write(fh, JSON.json(instancedata));
    close(fh);
end

# This is a utility function for laoding json data saved at filename
function load_json(filename)
    fh = open(filename);
    content_string = readlines(fh)[1];
    result = JSON.parse(content_string);
    close(fh);
    return result;
end

# This is a very specific macro for running try/catch blocks
# and printing a set of parameters to STDERR when something
# goes wrong.
#
# If it hits the catch portion, it will always continue (i.e.
# we assume this is always called within a loop).
macro trycatchcontinue(trycode, catchcode)
    quote
        try
            $(esc(trycode))
        catch err
            $(esc(catchcode))
            showerror(STDERR,err); Base.show_backtrace(STDERR,catch_backtrace())
            println()
            continue
        end
    end
end

# This is a macro that times the execution of a block of code.
# It returns a tuple with the amount of time as the first argument
# and the value returned by the code (if applicable) as the second.
# @arguments
# ex: Any block of code (an AST)
# @returns: (<time_delta>, <code_return_value>)
macro gettime(ex)
    quote
        local tstart = time();
        local val = $(esc(ex));
        local tend = time();
        (tend - tstart, val)
    end
end

# This is a macro that can be used to set the seed before executing code.
# It can be used in the following ways:
#
# @setseed <some_nondeterministic_code>
# @setseed <seednumber> <some_nondeterministic_code>
#
# In the former case, it searches the above scope to look for a value
# bound to the variable 'seed.'
macro setseed(args...)
    if length(args) == 1
        quote
            srand($(esc(seed)))
            $(esc(args[1]))
        end
    elseif length(args) == 2
        quote
            local seed = $(esc(args[1]))
            srand(seed)
            $(esc(args[2]))
        end
    else
        error("@setseed only accepts 1 or 2 arguments")
    end
end

# This is a private helper function for parsestringcli and parseintcli.
# It basically cycles through all the commandline args in ARGS and parses
# them according to the acceptable short and long options given in 'opts'
function _parsecli(opts...)
    for arg in ARGS
        climatch = matchcommandlinearg(arg, opts...)
        if climatch != nothing
           return climatch
        end
    end
    return nothing
end

# This is a utility function for parsing a commandline arg and binding
# that _string_ to the first argument here called 'variable.'
#
# @arguments
# variable: This is the symbol that the parsed cli chunk will be bound to.
# args:     The first two elements of this list must be the short and long
#           options for the commandline argument.
#           If there is a third element, it is the default in case no match is found
#
# Example usage: @parsestringcli foo "n" "num" "bar"
macro parsestringcli(variable, args...)
    if length(args) < 2 || length(args) > 3
        error("@parsecli only takes 2 possible flags and optional default")
    end
    local climatch = _parsecli(args[1], args[2])
    if climatch != nothing
        return :($(esc(variable)) = $(esc(climatch)))
    elseif length(args) == 3
        return :($(esc(variable)) = $(esc(args[3])))
    end
end

# This is a utility function for parsing a commandline arg and binding
# that _int_ to the first argument here called 'variable.'
#
# @arguments
# variable: This is the symbol that the parsed cli chunk will be bound to.
# args:     The first two elements of this list must be the short and long
#           options for the commandline argument.
#           If there is a third element, it is the default in case no match is found
#
# Example usage: @parseintcli seed "s" "seed" 8
macro parseintcli(variable, args...)
    if length(args) < 2 || length(args) > 3
        error("@parsecli only takes 2 possible flags and optional default")
    end
    local climatch = _parsecli(args[1], args[2])
    if climatch != nothing
        return :($(esc(variable)) = parse(Int, ($(esc(climatch)))))
    elseif length(args) == 3
        return :($(esc(variable)) = $(esc(args[3])))
    end
end


text:
using AutoExpe
using Test

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

max_stars_repo_name: ZacharieALES/AutoExpe

max_issues_repo_name: ZacharieALES/AutoExpe

max_forks_repo_name: ZacharieALES/AutoExpe


text:
module Comlineoption
    include("comlineoption_module.jl")
    using ArgParse
    using .Comlineoption_module

    const DEFINPNAME = "input.inp"

    function construct(args)
        parsed_args = parse_commandline(args)
        opt = Comlineoption_module.comlineoption_param(parsed_args["inputfile"], parsed_args["usethread"])

        return opt
    end

    function parse_commandline(args)
        s = ArgParseSettings()

        @add_arg_table! s begin
            "--inputfile", "-I"
                help = "インプットファイル名を指定します"
                arg_type = String
                default = DEFINPNAME

            "--usethread", "-T"
                help = "スレッド並列化するかどうかを指定します"
                arg_type = Bool
                default = true
        end

        return parse_args(args, s)
    end
end


text:
a = 3
b = 2
c = 1

fs = 100
T = 1/fs
f1 = 5
f2 = 10
f3 = 30

t = T:T:fs

x1 = a*sin.(2*pi*f1*T*t)+b*sin.(2*pi*f2*T*t)+c*sin.(2*pi*f3*T*t);

@testset "Minima/maxima" begin
    @test length(maxima(x1)) == 30
    @test length(minima(x1)) == 30

    @test length(maxima(x1, 1, true)) == 31
    @test length(minima(x1, 1, true)) == 31

    @test length(maxima(x1, 1000)) == 4
    @test length(minima(x1, 1000)) == 4

    @test length(maxima(x1, 1000, true)) == 5
    @test length(minima(x1, 1000, true)) == 5
end


text:
# Problem 6 - Project Euler
# http://projecteuler.net/index.php?section=problems&id=6

function p006()
    sqsum = 0
    sum = 0

    for i in 1:100
        sqsum += i * i
        sum += i
    end

    (sum * sum) - sqsum
end

println(p006())

max_stars_repo_name: emergent/ProjectEuler

max_issues_repo_name: emergent/ProjectEuler

max_forks_repo_name: emergent/ProjectEuler