text:
using PyCall
using MAT
using PyPlot
using Statistics
using ADCME
using DelimitedFiles
mpl = pyimport("tikzplotlib")

n = 15
m = 2n
h = 0.01
NT = 20
Δt = 2/NT

function plot1(model_id, k)
    U0 = matread("data$model_id.mat")["U"][4]
    Sigma0 = matread("data$model_id.mat")["S"][4]

    d = matread("model$(model_id)nn$k.mat")
    Sigma_ = d["S"]
    U_ = d["U"]
    close("all")
    plot(LinRange(0, 2.0, NT+1), U0[:,1], "r--", label="\$u_x\$")  # reference
    plot(LinRange(0, 2.0, NT+1), U_[:,1], "ro")
    plot(LinRange(0, 2.0, NT+1), U0[:,1+(n+1)*(m+1)], "g--", label="\$u_y\$")
    plot(LinRange(0, 2.0, NT+1), U_[:,1+(n+1)*(m+1)], "go")
    xlabel("Time")
    ylabel("Displacement")
    legend()
    mpl.save("nlv$(model_id)_disp$k.tex")
    savefig("nlv$(model_id)_disp$k.png")

    close("all")
    plot(LinRange(0, 2.0, NT+1), mean(Sigma0[:,1:4,1], dims=2)[:],"r--", label="\$\\sigma_{xx}\$")
    plot(LinRange(0, 2.0, NT+1), mean(Sigma0[:,1:4,2], dims=2)[:],"b--", label="\$\\sigma_{yy}\$")
    plot(LinRange(0, 2.0, NT+1), mean(Sigma0[:,1:4,3], dims=2)[:],"g--", label="\$\\sigma_{xy}\$")
    plot(LinRange(0, 2.0, NT+1), mean(Sigma_[:,1:4,1], dims=2)[:],"ro")
    plot(LinRange(0, 2.0, NT+1), mean(Sigma_[:,1:4,2], dims=2)[:],"bo")
    plot(LinRange(0, 2.0, NT+1), mean(Sigma_[:,1:4,3], dims=2)[:],"go")
    legend()
    legend()
    xlabel("Time")
    ylabel("Stress")
    mpl.save("nlv$(model_id)_stress$k.tex")
    savefig("nlv$(model_id)_stress$k.png")
end

function plot2(model_id)
    m = 30
    n = 15
    θ = matread("nn$(model_id)_20.mat")["Variablecolon0"]
    U0 = matread("data$model_id.mat")["U"][4]
    Sigma0 = matread("data$model_id.mat")["S"][4]
    S = Sigma0[end,:,:]
    SS = zeros(1800÷4,3)
    for i = 1:450
        SS[i,:] = mean(S[4(i-1)+1:4i,:], dims=1)[:]
    end
    val = squeeze(ae(SS, [20,20,20,1], θ))
    σ = constant(SS)
    if model_id==1
        val2 = constant(10*ones(m*n)) + 5.0/(1+1000*sum(σ[:,1:2]^2, dims=2))
    else
        val2 = constant(10*ones(m*n)) + relu(50.0/(1+1000*sum(σ[:,1:2]^2, dims=2)) - 10.0)
    end
    sess = Session()
    V1, V2, S = run(sess, [val, val2, σ])

    close("all")
    plot((@. (S[:,1]^2 + S[:,2]^2)), V1, "o", label="Estimated")
    plot((@. (S[:,1]^2 + S[:,2]^2)), V2, "+", label="Reference")

    θ = matread("nn$(model_id)_0.mat")["Variablecolon0"]
    U0 = matread("data$model_id.mat")["U"][4]
    Sigma0 = matread("data$model_id.mat")["S"][4]
    S = Sigma0[end,:,:]
    SS = zeros(1800÷4,3)
    for i = 1:450
        SS[i,:] = mean(S[4(i-1)+1:4i,:], dims=1)[:]
    end
    val = squeeze(ae(SS, [20,20,20,1], θ))
    σ = constant(SS)
    val2 = constant(10*ones(m*n)) + 5.0/(1+1000*sum(σ[:,1:2]^2, dims=2))
    sess = Session()
    V1, V2, S = run(sess, [val, val2, σ])
    plot((@. (S[:,1]^2 + S[:,2]^2)), V1, ".", label = "Initial Guess")

    legend()
    xlabel("\$\\sqrt{\\sigma_{11}^2+\\sigma_{22}^2}\$")
    ylabel("\$\\eta\$")
    mpl.save("nlv$(model_id).tex")
    savefig("nlv$(model_id).png")
end

function normalize(x, n=201)
    T = zeros(n)
    T[1:length(x)] = log.(x)/log(10)
    for i = length(x)+1:n
        T[i] = T[length(x)]
    end
    return T
end

function get_data(model_id, width, depth, activation)
    loss = readdlm("data/loss_$(model_id)_$(width)_$(depth)_$(activation).txt")[:]
    return normalize(loss)
end

function plot3(model_id, depth)
    U = zeros(20, 201)
    close("all")
    k = 0
    for activation in ["tanh", "relu", "selu", "elu"]
        for width in [1 5 10 20 40]
            k += 1
            U[k,:] = get_data(model_id, width, depth, activation)
            @info activation, width
        end
    end
    # subplot(120+v)
    close("all")
    pcolormesh(U', rasterized=true, cmap="hot")
    xticks([collect(0.5:1.0:19.5);2.5;7.5;12.5;17.5],
        [repeat([1 5 10 20 40], 1,4)[:];"tanh"; "relu"; "selu"; "elu"])
    c = colorbar()
    c.set_ticks([-2,-1,0,1,2,3])
    c.set_ticklabels(["\$10^{-2}\$", "\$10^{-1}\$",
            "\$10^0\$", "\$10^{1}\$", "\$10^{2}\$", "\$10^{3}\$"])
    va = -0.05*ones(5)
    ts = gca().get_xticklabels( )
    for i = 21:24
        ts[i].set_y( va[i-20] )
    end
    z = 0:200
    for i = 0:20
        plot(ones(201)*i, z, "k", alpha=0.2)
    end
    for i in [5, 10, 15]
        plot(ones(201)*i, z, "k")
    end
    # gca().tick_params(axis="both", which="major", labelsize=20)
    gca().tick_params(axis="both", which="major", labelsize=10)
    ylabel("Iterations")

    savefig("nn$(model_id)$depth.png")
    savefig("nn$(model_id)$depth.pdf")
end

plot1(1, 0)
plot1(1, 20)
plot1(2, 0)
plot1(2, 20)

plot2(1)
plot2(2)

for i in [1,3,5,10]
    plot3(1, i)
    plot3(2, i)
end

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
module DataInterpolations

### Interface Functionality

abstract type AbstractInterpolation{FT,T} <: AbstractVector{T} end

Base.size(A::AbstractInterpolation) = size(A.u)
Base.size(A::AbstractInterpolation{true}) = length(A.u) .+ size(A.t)
Base.getindex(A::AbstractInterpolation,i) = A.u[i]
Base.getindex(A::AbstractInterpolation{true},i) =
    i<=length(A.u) ? A.u[i] : A.t[i-length(A.u)]
Base.setindex!(A::AbstractInterpolation,x,i) = A.u[i] = x
Base.setindex!(A::AbstractInterpolation{true},x,i) =
    i <= length(A.u) ? (A.u[i] = x) : (A.t[i-length(A.u)] = x)

using ChainRulesCore, LinearAlgebra, RecursiveArrayTools, RecipesBase, Reexport
@reexport using Optim

include("interpolation_caches.jl")
include("interpolation_utils.jl")
include("interpolation_methods.jl")
include("plot_rec.jl")
include("derivatives.jl")
include("integrals.jl")
include("online.jl")

function ChainRulesCore.rrule(::typeof(_interpolate),
                              A::Union{LagrangeInterpolation,AkimaInterpolation,
                                       BSplineInterpolation,BSplineApprox}, t::Number)
    interpolate_pullback(Δ) = (NoTangent(), NoTangent(), derivative(A, t) * Δ)
    return _interpolate(A, t), interpolate_pullback
end

ChainRulesCore.frule((_, _, Δt), ::typeof(_interpolate), A::AbstractInterpolation,
                     t::Number) = _interpolate(A, t), derivative(A, t) * Δt

(interp::AbstractInterpolation)(t::Number) = _interpolate(interp, t)

using Symbolics: Num, unwrap, SymbolicUtils
(interp::AbstractInterpolation)(t::Num) = SymbolicUtils.term(interp, unwrap(t))
SymbolicUtils.promote_symtype(t::AbstractInterpolation, _...) = Real

export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation,
    AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline,
    BSplineInterpolation, BSplineApprox, Curvefit

# added for RegularizationSmooth, JJS 11/27/21
import RegularizationTools
const RT = RegularizationTools
include("regularization_smooth.jl")
export RegularizationSmooth

# Deprecated April 2020
export ZeroSpline

end # module

ext: jl


text:
@gen function f_dynamic(x, y)
    z ~ normal(0, 1)
    return x + y + z
end
@gen (static) function f_static(x, y)
    z ~ normal(0, 1)
    return x + y + z
end
@load_generated_functions()

for (lang, f) in [:dynamic => f_dynamic,
                  :static  => f_static]
    @testset "update(...) shorthand assuming unchanged args ($lang modeling lang)" begin
        trace0 = simulate(f, (5, 6))

        constraints = choicemap((:z, 0))
        trace1, _, _, discard = update(trace0, constraints)
        # The main test is that the shorthand version runs without crashing,
        # which is already shown by the time we get here.  Beyond that, let's
        # sanity-check that `update` did what it's supposed to.
        @test get_args(trace1) == (5, 6)
        @test trace1[:z] == 0
        @test :z in keys(get_values_shallow(discard))
    end

    @testset "regenerate(...) shorthand assuming unchanged args ($lang modeling lang)" begin
        trace0 = simulate(f, (5, 6))
        trace1, _, _ = regenerate(trace0, select(:z))
        # The main test is that the shorthand version runs without crashing,
        # which is already shown by the time we get here.  Beyond that, let's
        # sanity-check that `regenerate` did what it's supposed to.
        @test get_args(trace1) == (5, 6)
        @test trace1[:z] != trace0[:z]
    end
end

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
@testset "Poisson 1D" begin

    alpha = 0.1
    kx = 0.5
    xmin, xmax = 0.0, 2π / kx
    nx = 64

    mesh = OneDGrid( xmin, xmax, nx )

    rho = zeros(nx)
    expected_ex = zeros(nx)

    for i in 1:nx
        rho[i] = alpha * cos(kx * mesh.x[i])
        expected_ex[i] = alpha / kx * sin(kx * mesh.x[i])
    end

    poisson = OneDPoisson( mesh )

    computed_ex = zero(expected_ex)

    ParticleInCell.compute_e_from_rho!( computed_ex, poisson, rho )

    @show maximum(abs.( computed_ex .- expected_ex ))

    @test computed_ex ≈ expected_ex

end


text:
module AugmentQTests

include(joinpath(dirname(@__FILE__), "..", "src", "AndersonMoore.jl"))
# test augmentQ
using .AndersonMoore

# test augmentQ! firmvalue example
function firmvalue()::Bool

neq=4::Int64;nlag=1::Int64;nlead=1::Int64
qRows=(neq*nlead)::Int64;qCols=(neq*(nlag+nlead))::Int64

qFinal=[0.  -0.4  0.  0.  0.  1.  -1.  0.;
0.  0.  0.  0.  0.  0.  1.  0.;
0.  0.  0.  -1.  0.  0.  0.  1.;
0.  0.  0.  0.  0.868243  -0.496139  0.  0.]::Array{Float64,2}

qNew=[0.  -0.4  0.  0.  0.  1.  -1.  0.;
0.  0.  0.  0.  0.  0.  1.  0.;
0.  0.  0.  -1.  0.  0.  0.  1.;
0.  0.  0.  0.  0.  0.  0.  0.]::Array{Float64,2}

ww=[0.868243  0.  0.;
-0.496139  0.  1.;
0.  1.  0.]::Array{Float64,2}

iqNew=3::Int64

js=[5  6  8]::Array{Int64,2}

(qqJulia)=augmentQ!(qNew,ww,js,iqNew,qRows)
isapprox(qqJulia,qFinal,rtol=0.::Float64,atol=1e-16::Float64)
end;

# test augmentQ! firmvalue3Leads2Lags example
function firmvalue3Leads2Lags()::Bool

neq=4::Int64;nlag=2::Int64;nlead=3::Int64
qRows=(neq*nlead)::Int64;qCols=(neq*(nlag+nlead))::Int64

qFinal=[0.  -0.4  0.  0.  0.  0.  0.  0.  0.  1.  -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  -0.4  0.  0.  0.  0.  0.  0.  0.  1.  -1.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  1.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  -0.4  0.  0.  0.  0.  0.  0.  0.  1.  -1.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  1.;
0.  0.  0.  0.  0.  0.  0.  0.  0.556804  0.  0.  0.  -0.691462  -0.42174  0.  0.  0.168433  0.392762  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.556804  0.  0.  0.  0.185277  -0.255383  0.  0.  -0.628602  -0.21595  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.455955  0.  0.  0.  0.414505  -0.511734  0.  0.  0.376823  -0.465213  0.  0.]::Array{Float64,2}

qNew=[0.  -0.4  0.  0.  0.  0.  0.  0.  0.  1.  -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  -0.4  0.  0.  0.  0.  0.  0.  0.  1.  -1.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  1.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  -0.4  0.  0.  0.  0.  0.  0.  0.  1.  -1.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  1.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]::Array{Float64,2}

ww=[0.556804  0.556804  0.455955  0.  -0.00000000000000000884102  0.00000000000000000251725;
-0.691462  0.185277  0.414505  0.  0.0000000000000000128181  0.0000000000000000214521;
-0.42174  -0.255383  -0.511734  0.  0.534522  0.534522;
0.168433  -0.628602  0.376823  0.  -0.00000000000000000533916  -0.0000000000000000148456;
0.392762  -0.21595  -0.465213  0.  0.845154  -0.845154;
0.  0.  0.  1.  0.  0.]::Array{Float64,2}

iqNew=9::Int64

js=[9  13  14  17  18  20]::Array{Int64,2}

(qJulia)=augmentQ!(qNew,ww,js,iqNew, qRows)
isapprox(qJulia,qFinal,rtol=0.::Float64,atol=1e-16::Float64)
end;

# test augmentQ! example7 example
function example7()::Bool

neq=4::Int64;nlag=1::Int64;nlead=1::Int64
qRows=(neq*nlead)::Int64;qCols=(neq*(nlag+nlead))::Int64

qFinal=[0.  0.  0.  1.  0.  0.  1.  -1.;
0.  0.  0.  -0.66  0.  -1.1  0.  1.;
0.  0.  0.  0.  -0.0975567  0.663306  0.  0.391168;
0.  0.  0.  0.  -0.665933  0.663306  0.  -0.716975]::Array{Float64,2}

qNew=[0.  0.  0.  1.  0.  0.  1.  -1.;
0.  0.  0.  -0.66  0.  -1.1  0.  1.;
0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.]::Array{Float64,2}

ww=[-0.0975567  -0.665933  0.136197;
0.663306  0.663306  0.897059;
0.391168  -0.716975  -0.420399]::Array{Float64,2}

iqNew=2::Int64

js=[5  6  8]::Array{Int64,2}

(qJulia)=augmentQ!(qNew,ww,js,iqNew,qRows)
isapprox(qJulia,qFinal,rtol=0.::Float64,atol=1e-16::Float64)
end;

# test augmentQ! oneEquationNoLead example
function oneEquationNoLead()::Bool

neq=1::Int64;nlag=1::Int64;nlead=1::Int64
qRows=(neq*nlead)::Int64;qCols=(neq*(nlag+nlead))::Int64

qFinal=[2.  3.]::Array{Float64,2}

qNew=[2.  3.]::Array{Float64,2}

ww = hcat([1.])::Array{Float64,2}

iqNew=1::Int64

js = hcat([2])::Array{Int64,2}

(qJulia)=augmentQ!(qNew,ww,js,iqNew,qRows)
isapprox(qJulia,qFinal,rtol=0.::Float64,atol=1e-16::Float64)
end;

# test augmentQ! reliablePaperExmpl example
function reliablePaperExmpl()::Bool

neq=5::Int64;nlag=1::Int64;nlead=1::Int64
qRows=(neq*nlead)::Int64;qCols=(neq*(nlag+nlead))::Int64

qFinal=[0.  0.  0.  -0.5  0.  0.  -1.  0.7  -0.5  1.;
0.  0.  0.  0.  0.  1.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  1.  0.  0.  0.;
0.  0.  0.  0.  0.  -0.548572  0.  0.548572  0.475037  -0.544611;
0.  0.  0.  0.  0.  -0.548572  0.  0.548572  -0.508387  0.124726]::Array{Float64,2}

qNew=[0.  0.  0.  -0.5  0.  0.  -1.  0.7  -0.5  1.;
0.  0.  0.  0.  0.  1.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  1.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]::Array{Float64,2}

ww=[-0.548572  -0.548572  -0.604781  1.;
0.548572  0.548572  0.604781  0.;
0.475037  -0.508387  -0.291542  0.;
-0.544611  0.124726  0.428349  0.]::Array{Float64,2}

iqNew=3::Int64

js=[6  8  9  10]::Array{Int64,2}

(qJulia)=augmentQ!(qNew,ww,js,iqNew,qRows)
isapprox(qJulia,qFinal,rtol=0.::Float64,atol=1e-16::Float64)
end;

# test augmentQ! athan example
function athan()::Bool

neq=9::Int64;nlag=1::Int64;nlead=1::Int64
qRows=(neq*nlead)::Int64;qCols=(neq*(nlag+nlead))::Int64

qFinal=[-0.293  -0.764  0.293  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  -1.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  -2.  -1.  1.  1.  0.  0.  2.  0.  0.;
0.  0.  0.  -0.911  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  -1.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  0.  0.  0.  0.  0.  1.;
-0.381771  0.  0.  0.  0.  0.  0.  0.  0.  0.447207  -0.250613  0.0434036  -0.081062  -0.763541  0.  0.  0.  0.]::Array{Float64,2}

qNew=[-0.293  -0.764  0.293  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  -1.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  -2.  -1.  1.  1.  0.  0.  2.  0.  0.;
0.  0.  0.  -0.911  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  -1.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  0.  0.  0.  0.  0.  1.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]::Array{Float64,2}

ww=[-0.381771  0.  0.  0.330041  0.330041  0.0000000000000000360124  0.;
0.447207  0.  0.  -0.454296  0.269733  -0.120158  0.;
-0.250613  0.  0.  -0.147243  0.605615  -0.0600791  0.;
0.0434036  0.  0.  0.102035  -0.186692  0.0600791  0.;
-0.081062  0.  1.  -0.015374  -0.444561  0.0600791  0.;
-0.763541  0.  0.  0.660083  0.660083  0.987286  1.;
0.  1.  0.  0.  0.  0.  0.]::Array{Float64,2}

iqNew=8::Int64

js=[1  10  11  12  13  14  18]::Array{Int64,2}

(qJulia)=augmentQ!(qNew,ww,js,iqNew,qRows)
isapprox(qJulia,qFinal,rtol=0.::Float64,atol=1e-16::Float64)
end;

# test augmentQ! habitmod example
function habitmod()::Bool

neq=12::Int64;nlag=4::Int64;nlead=1::Int64
qRows=(neq*nlead)::Int64;qCols=(neq*(nlag+nlead))::Int64

qFinal=[0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  -1.  -1.  0.  0.  3.45419  0.  0.  -0.74118  0.00208517  0.494116  -0.1149;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -0.88879  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.52802  0.  0.11279  0.031833  0.024734  0.  0.  0.  0.  0.  0.  0.  -0.71526  0.  0.042076  -0.22829  -0.012267  0.  0.  0.  0.  0.  0.  0.  0.22215  0.  -0.14172  0.26679  -0.11959  0.  0.  0.  0.  0.  0.  0.  -0.42308  0.  -0.61076  -0.10946  0.17172  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  -0.00011277;
0.46436  0.  0.16131  0.02744  0.17087  0.  0.  0.  0.  0.  0.  0.  -1.3425  0.  0.16794  -0.35818  -0.09837  0.  0.  0.  0.  0.  0.  0.  1.06382  0.  -0.02255  0.51594  0.08949  0.  0.  0.  0.  0.  0.  0.  -0.46672  0.  -0.11369  -1.10657  -0.26021  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  -0.00141;
-0.23292  0.  0.08512  -0.04213  0.14149  0.  0.  0.  0.  0.  0.  0.  0.03313  0.  0.15223  0.01572  -0.34762  0.  0.  0.  0.  0.  0.  0.  0.70552  0.  -0.24697  0.21888  -0.07805  0.  0.  0.  0.  0.  0.  0.  -0.56316  0.  -0.08268  -0.16054  -0.68637  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  -0.00348;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -0.99853  0.  0.  0.  0.  0.  -0.00147  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.00000847733  0.  0.00000117376  0.00000064385  -0.000000322128  0.  0.  0.  0.  0.  0.  0.  -0.0000101418  0.  -0.000000129057  -0.00000326377  0.00000147074  0.  0.  0.  0.  0.  0.  0.  -0.000000239717  0.  -0.000000811654  0.00000268973  -0.00000129793  0.  0.  0.  0.  0.  0.  0.  -0.0370804  -0.0000485593  -0.0000488121  -0.00000549337  0.00000563828  0.000193477  0.998571  0.0384857  -0.0000406731  0.000000114426  0.0000271152  -0.00000626373;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -0.0489971  0.  -0.0112006  -0.00419088  -0.0119854  0.  0.  0.  0.  0.  0.  0.  0.0895547  0.  -0.0207267  0.0311734  -0.0160863  0.  0.  0.  0.  0.  0.  0.  -0.000792326  0.  -0.0311743  -0.01287  -0.0312023  0.  0.  0.  0.  0.  0.  0.  0.0246579  0.0958625  0.127132  0.12544  -0.0259403  -0.944771  -0.109495  -0.000509914  0.162178  -0.000456259  -0.108118  0.0240777;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.00289063  0.  0.000619471  0.000197516  0.000288522  0.  0.  0.  0.  0.  0.  0.  -0.0021565  0.  0.000949653  -0.0012804  0.000358271  0.  0.  0.  0.  0.  0.  0.  0.000033234  0.  0.000632742  0.000431941  0.000187315  0.  0.  0.  0.  0.  0.  0.  -0.0000000000000109556  -0.000000000000104638  0.429473  -0.00191937  0.000869072  0.000000000000557064  -0.561381  -0.00525455  0.431073  -0.00523235  -0.560807  -0.0000161223;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.00337162  0.  0.000722547  0.000230381  0.00033653  0.  0.  0.  0.  0.  0.  0.  -0.00251533  0.  0.00110767  -0.00149345  0.000417886  0.  0.  0.  0.  0.  0.  0.  0.0000387639  0.  0.000738027  0.000503813  0.000218483  0.  0.  0.  0.  0.  0.  0.  -0.000000000000000619476  -0.000000000000000412551  0.500935  -0.00223874  0.00101368  -0.000000000000159624  -0.4981  0.000143356  0.502801  0.00014275  -0.498119  -0.000018805;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -0.00296384  0.  0.0000660393  -0.00028174  0.000922354  0.  0.  0.  0.  0.  0.  0.  -0.0024055  0.  0.000894541  0.000153904  -0.000584756  0.  0.  0.  0.  0.  0.  0.  0.0022348  0.  0.00092268  0.00083439  0.000286668  0.  0.  0.  0.  0.  0.  0.  0.00016314  0.0012582  -0.5851  -0.00344982  -0.0042903  0.0275468  0.393258  0.00166552  -0.589634  0.00165883  0.393087  0.000519754]::Array{Float64,2}

qNew=[0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  -1.  -1.  0.  0.  3.45419  0.  0.  -0.74118  0.00208517  0.494116  -0.1149;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -0.88879  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.52802  0.  0.11279  0.031833  0.024734  0.  0.  0.  0.  0.  0.  0.  -0.71526  0.  0.042076  -0.22829  -0.012267  0.  0.  0.  0.  0.  0.  0.  0.22215  0.  -0.14172  0.26679  -0.11959  0.  0.  0.  0.  0.  0.  0.  -0.42308  0.  -0.61076  -0.10946  0.17172  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  -0.00011277;
0.46436  0.  0.16131  0.02744  0.17087  0.  0.  0.  0.  0.  0.  0.  -1.3425  0.  0.16794  -0.35818  -0.09837  0.  0.  0.  0.  0.  0.  0.  1.06382  0.  -0.02255  0.51594  0.08949  0.  0.  0.  0.  0.  0.  0.  -0.46672  0.  -0.11369  -1.10657  -0.26021  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  -0.00141;
-0.23292  0.  0.08512  -0.04213  0.14149  0.  0.  0.  0.  0.  0.  0.  0.03313  0.  0.15223  0.01572  -0.34762  0.  0.  0.  0.  0.  0.  0.  0.70552  0.  -0.24697  0.21888  -0.07805  0.  0.  0.  0.  0.  0.  0.  -0.56316  0.  -0.08268  -0.16054  -0.68637  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  -0.00348;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -0.99853  0.  0.  0.  0.  0.  -0.00147  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]::Array{Float64,2}

ww=[0.00000847733  -0.0489971  0.00289063  0.00337162  -0.00296384  0.  0.00289826  0.108025  0.  -0.0438951  -0.017043  -0.0201472  -0.0791448  -0.0676478  0.204738  0.00741613  0.0457941  -0.0564212  -0.00691027  -0.00120837  0.000000000269016  -0.00000000000411353  -0.00000000000450183  -0.000000000277616;
0.00000117376  -0.0112006  0.000619471  0.000722547  0.0000660393  0.  -0.0176289  -0.0074902  0.  -0.0033955  -0.00895106  -0.00876825  -0.0129104  -0.0103006  0.028083  0.0304513  0.00587952  -0.0187752  -0.000579761  0.00136933  -0.000000000000250339  -0.0000000000000146502  -0.0000000000000142445  0.000000000000221426;
0.00000064385  -0.00419088  0.000197516  0.000230381  -0.00028174  0.  0.000265449  0.0158931  0.  -0.00228812  -0.000766883  -0.00000997697  -0.00614689  -0.00623485  0.0148889  -0.00363062  0.00467283  0.00161117  -0.00049169  -0.000532766  0.0000000000000819604  0.00000000000000523411  0.00000000000000510892  -0.0000000000000716094;
-0.000000322128  -0.0119854  0.000288522  0.00033653  0.000922354  0.  -0.0441305  -0.010256  0.  0.0149254  -0.0121877  -0.00423626  -0.00265616  -0.00671859  -0.0126441  0.0452943  0.00474333  0.0127896  0.0014004  0.00101101  0.0000000000000103699  0.0000000000000104046  0.0000000000000103302  0.0000000000000103699;
-0.0000101418  0.0895547  -0.0021565  -0.00251533  -0.0024055  0.  0.200931  -0.107545  0.  -0.0900449  0.0628586  -0.0156808  0.254502  0.0432761  -0.397022  0.133855  -0.0655314  0.115807  -0.0154165  -0.0167395  -0.00000132535  0.00000141058  0.00000141025  -0.00000150715;
-0.000000129057  -0.0207267  0.000949653  0.00110767  0.000894541  0.  -0.0567803  -0.0272877  0.  0.00609649  -0.0186839  -0.0242131  0.0165386  -0.0222346  -0.0754314  0.0606591  0.0329998  0.047426  -0.00405635  0.00410942  0.00000000130129  -0.00000000121979  -0.00000000121974  0.00000000113833;
-0.00000326377  0.0311734  -0.0012804  -0.00149345  0.000153904  0.  0.0443083  -0.039071  0.  -0.00331164  0.0177129  0.00839306  0.0485082  0.0314834  -0.0792245  -0.0148154  -0.0394007  0.0012426  0.00120493  -0.00215158  -0.0000000004235  0.000000000397062  0.000000000397044  -0.000000000370604;
0.00000147074  -0.0160863  0.000358271  0.000417886  -0.000584756  0.  -0.0213281  0.0800915  0.  0.0159632  -0.00532987  0.0000349571  -0.000670496  -0.019063  -0.0240676  -0.121402  -0.0196196  0.0286035  0.000173991  0.00151099  -0.0000000000000629115  0.0000000000000138601  0.0000000000000138823  0.000000000000090653;
-0.000000239717  -0.000792326  0.000033234  0.0000387639  0.0022348  0.  -0.0287681  -0.110363  0.  -0.0442625  0.0259932  -0.00592813  -0.180197  0.323735  0.106488  -0.446833  -0.091663  -0.259238  0.00224181  -0.0837174  -0.000657423  -0.0139054  0.0153282  -0.000773483;
-0.000000811654  -0.0311743  0.000632742  0.000738027  0.00092268  0.  -0.0906195  0.0653455  0.  0.0382729  -0.0108365  -0.0250779  0.0303865  0.0158253  -0.0580298  -0.193109  -0.139833  0.0894769  -0.0202529  -0.00355066  -0.000000207769  0.0000128265  -0.0000124039  -0.00000021471;
0.00000268973  -0.01287  0.000431941  0.000503813  0.00083439  0.  -0.0300244  -0.0284717  0.  0.00447302  -0.0108379  0.00467599  -0.0998469  0.0278434  0.146725  -0.00860273  0.0905782  -0.0943029  0.00879296  -0.00123068  0.0000000672228  -0.00000417477  0.00000403798  0.0000000695332;
-0.00000129793  -0.0312023  0.000187315  0.000218483  0.000286668  0.  -0.0907005  0.151454  0.  0.0594932  -0.0159445  0.0111903  0.0325489  -0.01021  0.153668  0.0785177  -0.0105616  0.00218828  -0.000218761  0.00114412  0.000000000423522  -0.00000000039704  -0.000000000397022  0.000000000370627;
-0.0370804  0.0246579  -0.0000000000000109556  -0.000000000000000619476  0.00016314  0.  0.00390901  0.042236  0.  -0.0708795  0.0430473  0.17663  -0.121307  -0.402926  0.168196  0.378419  0.222077  0.350449  0.162999  0.244079  0.252045  0.252094  0.251997  0.252043;
-0.0000485593  0.0958625  -0.000000000000104638  -0.000000000000000412551  0.0012582  0.  -0.369591  0.161292  1.  0.361963  -0.329714  -0.229698  0.23356  0.209559  -0.209672  -0.194319  -0.231622  -0.219961  -0.26965  -0.250136  -0.252046  -0.252038  -0.252051  -0.252044;
-0.0000488121  0.127132  0.429473  0.500935  -0.5851  0.  -0.434604  -0.00160089  0.  0.346801  -0.203295  -0.183594  0.193904  0.218294  -0.0525509  -0.0682873  -0.110412  -0.370963  -0.229611  -0.256183  -0.252045  -0.252057  -0.252032  -0.252043;
-0.00000549337  0.12544  -0.00191937  -0.00223874  -0.00344982  0.  0.253458  -0.250185  0.  -0.0925811  0.0341928  -0.0312099  0.05905  -0.0406851  -0.096241  0.0167473  -0.061666  0.0923319  -0.000184476  0.00342656  0.0000000000000517042  0.0000000000000518429  0.0000000000000515574  0.0000000000000517101;
0.00000563828  -0.0259403  0.000869072  0.00101368  -0.0042903  0.  0.00882998  0.295823  0.  0.0250455  -0.00230957  0.0163765  -0.0328191  0.0263861  -0.10227  0.0107738  0.0347595  -0.0520895  -0.010955  -0.0036157  -0.0000000672229  0.00000417477  -0.00000403798  -0.0000000695333;
0.000193477  -0.944771  0.000000000000557064  -0.000000000000159624  0.0275468  0.  0.764362  0.764362  0.  -0.808862  0.888581  0.888581  -0.804108  -0.804108  0.733356  0.733356  0.807099  0.807099  0.872336  0.872336  0.870612  0.870612  0.870605  0.870605;
0.998571  -0.109495  -0.561381  -0.4981  0.393258  0.  0.132405  0.104521  0.  -0.117758  0.129437  0.116763  -0.10817  -0.0916235  0.0952083  0.0831025  0.102589  0.0954434  0.117021  0.112195  0.000659587  0.0138518  -0.0152792  0.000775859;
0.0384857  -0.000509914  -0.00525455  0.000143356  0.00166552  0.  0.000553558  0.000363793  0.  -0.000366298  0.000468001  0.000164372  -0.000233277  0.000362695  0.0000913327  -0.000352556  0.0000236021  -0.000248118  0.000192018  0.0000171507  0.00000000166111  -0.000000103175  0.0000000997946  0.00000000171824;
-0.0000406731  0.162178  0.431073  0.502801  -0.589634  0.  -0.198838  -0.160111  0.  0.182241  -0.197634  -0.189195  0.172734  0.170126  -0.157146  -0.155475  -0.173005  -0.171768  -0.188855  -0.187001  -0.186811  -0.18681  -0.18681  -0.18681;
0.000000114426  -0.000456259  -0.00523235  0.00014275  0.00165883  0.  0.000559392  0.000450443  0.  -0.000512702  0.000556006  0.000532264  -0.000485953  -0.000478617  0.0004421  0.000437399  0.000486716  0.000483238  0.000531309  0.000526092  0.000525558  0.000525556  0.000525556  0.000525554;
0.0000271152  -0.108118  -0.560807  -0.498119  0.393087  0.  0.132557  0.10674  0.  -0.121493  0.131755  0.126129  -0.115155  -0.113416  0.104763  0.103649  0.115336  0.114511  0.125903  0.124666  0.12454  0.124539  0.124539  0.124539;
-0.00000626373  0.0240777  -0.0000161223  -0.000018805  0.000519754  1.  -0.0390233  -0.026122  0.  0.0288168  -0.031094  -0.0292324  0.0267982  0.0262171  -0.0240533  -0.0241842  -0.0268358  -0.0264974  -0.0293165  -0.028979  -0.02896  -0.0289599  -0.0289599  -0.0289598]::Array{Float64,2}

iqNew=7::Int64

js=[13  15  16  17  25  27  28  29  37  39  40  41  49  50  51  52  53  54  55  56  57  58  59  60]::Array{Int64,2}

(qJulia)=augmentQ!(qNew,ww,js,iqNew,qRows)
isapprox(qJulia,qFinal,rtol=0.::Float64,atol=1e-16::Float64)
end;

end


text:
# src path
src_path="C:\\Users\\Koike\\Desktop\\git\\Qtree3d\\src\\"

# main (変更しないこと)
src_read="read_uns.jl"
include(src_path*src_read)
src_read="read_para.jl"
include(src_path*src_read)
src_read="misc.jl"
include(src_path*src_read)
src_read="pre.jl"
include(src_path*src_read)

src_read="main_Qtree.jl"
include(src_path*src_read)


text:
function getSubpopulationAtTime( mpSim::MPsim, timePoint::Float64,
    subpopulations::Vector{Subpopulation} )

    result = fill( Vector{String}(), length( subpopulations ) )
    hasAttributeConds = map( subpopulation -> !isempty(
        subpopulation.attributeConds ), subpopulations )

    # Find all personnel members that satisfy the source node, the history
    #   conditions, and the time conditions.
    for ii in eachindex( subpopulations )
        subpopulation = subpopulations[ii]
        result[ii] = subpopulation.sourceNode == "active" ?
            getActiveAtTime( mpSim, timePoint, subpopulation ) :
            getActiveAtTime( mpSim, timePoint, subpopulation, false )  # ! To change
    end  # for ii in eachindex( subpopulations )

    if !any( hasAttributeConds )
        return result
    end  # if !any( hasAttributeConds )

    # Get the personnel database snapshot of the personnel members for whom a
    #   database state reconstruction is required.
    idsToReconstruct = vcat( result[hasAttributeConds]... )
    queryCmd = string( "SELECT * FROM `", mpSim.persDBname, "` WHERE",
        "\n    `", mpSim.idKey, "` IN ( '" , join( idsToReconstruct, "', '" ),
        "' )" )
    activeAtTime = DataFrame( DBInterface.execute( mpSim.simDB, queryCmd ) )

    if isempty( activeAtTime )
        return result
    end  # if isempty( activeAtTime )

    getSubpopulationStateAtTime!( activeAtTime, mpSim, timePoint )
    idSymbol = Symbol( mpSim.idKey )

    # For each request, filter out personnel members satisfying the conditions
    #   of the request.
    for ii in filter( ii -> hasAttributeConds[ii], eachindex( result ) )
        subpopulation = subpopulations[ii]
        satisfyConds = map( id -> id ∈ result[ii],
            activeAtTime[:, idSymbol] )

        for condition in subpopulation.attributeConds
            satisfyConds .&= map(  activeAtTime[:,
                Symbol( condition.attribute )] ) do val
                return condition.operator( val, condition.value )
            end  # map( ... )
        end  # for cond in attribConds

        result[ii] = activeAtTime[satisfyConds, idSymbol]
    end  # for ii in filter( ..., eachindex( result ) )

    return result

end  # getSubpopulationAtTime( mpSim, timePoint, subpopulations )

function getActiveAtTime( mpSim::MPsim, timePoint::Float64,
    subpopulation::Subpopulation )

    # Generate the SQLite query.
    queryCmd = string( "SELECT *, ",
        timePoint, " - timeEntered tenure, ",
        timePoint, " - timeEntered + ageAtRecruitment age FROM `",
        mpSim.persDBname, "` WHERE",
        "\n    timeEntered <= ", timePoint, " AND ( timeExited > ", timePoint,
        " OR timeExited IS NULL )" )

    if !isempty( subpopulation.timeConds )
        condSQLite = conditionToSQLite.( subpopulation.timeConds, Ref( mpSim ) )

        # Replace "time in node" by "tenure"
        for ii in eachindex( condSQLite )
            if startswith( condSQLite[ii], "`time in node`" )
                condSQLite[ii] = replace( condSQLite[ii],
                    "`time in node`" => "tenure" )
            end  # if startswith( condSQLite[ii], "`time in node`" )
        end  # for ii in eachindex( condSQLite )

        queryCmd = string( queryCmd, " AND ", join( condSQLite, " AND " ) )
    end  # if !isempty( timeConds )

    if !isempty( subpopulation.historyConds )
        # Identify which history conditions are negatives and adjust for them.
        histConds = deepcopy( subpopulation.historyConds )
        isNegHistCond = map( cond -> cond.operator ∈ [!=, ∉], histConds )

        for ii in eachindex( histConds )
            if isNegHistCond[ii]
                condition = histConds[ii]
                histConds[ii] = MPcondition( condition.attribute,
                    condition.operator == Base.:∈ ? Base.:∈ : Base.:(==),
                    condition.value )
            end  # if isNegHistCond[ii]
        end  # for ii in eachindex( histConds )

        # Generate history condition queries.
        condSQLite = conditionToSQLite.( histConds, Ref( mpSim ) )
        condSQLite = string.( "`", mpSim.idKey, "`",
            map( bVal -> bVal ? " NOT" : "", isNegHistCond ), " IN ( SELECT `",
            mpSim.idKey, "` FROM `" , mpSim.transDBname, "`
                WHERE timeIndex <= ", timePoint, condSQLite,
            " )" )
        queryCmd = string( queryCmd, " AND ", join( condSQLite, " AND " ) )
    end  # if !isempty( subpopulation.historyConds )

    queryCmd = string( queryCmd, "\n    ORDER BY `", mpSim.idKey, "`" )
    result = DataFrame( DBInterface.execute( mpSim.simDB, queryCmd ) )
    return isempty( result ) ? Vector{String}() :
        string.( result[:, Symbol( mpSim.idKey )] )

end  # getActiveAtTime( MPsim, timePoint, subpopulation )

function getActiveAtTime( mpSim::MPsim, timePoint::Float64,
    subpopulation::Subpopulation, fullPop::Bool )

    # Only for safety.
    if fullPop
        return getActiveAtTime( mpSim, timePoint, subpopulation )
    end  # if fullPop

    # Identify which history conditions are negatives and adjust for them.
    histConds = deepcopy( subpopulation.historyConds )
    isNegHistCond = map( cond -> cond.operator ∈ [!=, ∉], histConds )

    for ii in eachindex( histConds )
        if isNegHistCond[ii]
            condition = histConds[ii]
            histConds[ii] = MPcondition( condition.attribute,
                condition.operator == Base.:∈ ? Base.:∈ : Base.:(==),
                condition.value )
        end  # if isNegHistCond[ii]
    end  # for ii in eachindex( histConds )

    # Generate history condition queries.
    queryCmd = conditionToSQLite.( histConds, Ref( mpSim ) )
    queryCmd = string.( "`", mpSim.idKey, "`",
        map( bVal -> bVal ? " NOT" : "", isNegHistCond ), " IN ( SELECT `",
        mpSim.idKey, "` FROM `" , mpSim.transDBname, "` WHERE",
            "\n    timeIndex <= ", timePoint, queryCmd, " )" )

    # Generate subquery.
    nodeName = subpopulation.sourceNode
    nodeCond = haskey( mpSim.baseNodeList, nodeName ) ? nodeName :
        join( mpSim.compoundNodeList[nodeName].baseNodeList, "', '" )
    queryCmd = string( "SELECT `", mpSim.idKey, "` tmpID, ",
        "max(timeIndex) lastTransTime FROM `", mpSim.transDBname, "` WHERE",
        "\n    ", mpSim.tNode, " IS NOT 'active' AND",
        "\n    timeIndex <= ", timePoint, isempty( queryCmd ) ? "" : " AND ", join( queryCmd, " AND " ), "\n    GROUP BY `", mpSim.idKey, "`",
        "\n    HAVING ", mpSim.tNode, " IN ( '", nodeCond, "' )" )

    # Generate full query.
    queryCmd = string( "SELECT *, ",
        timePoint, " - timeEntered tenure, ",
        timePoint, " - timeEntered + ageAtRecruitment age, ",
        timePoint, " - lastTransTime `time in node` FROM `",
        mpSim.persDBname, "`",
        "\n    INNER JOIN ( ", queryCmd, " ) tmpList ON `", mpSim.idKey,
        "` IS tmpID" )

    if !isempty( subpopulation.timeConds )
        condSQLite = conditionToSQLite.( subpopulation.timeConds, Ref( mpSim ) )
        queryCmd = string( queryCmd, "\n    WHERE ",
            join( condSQLite, " AND " ) )
    end  # if !isempty( timeConds )

    queryCmd = string( queryCmd, "\n    ORDER BY `", mpSim.persDBname, "`.`",
        mpSim.idKey, "`" )
    result = DataFrame( DBInterface.execute( mpSim.simDB, queryCmd ) )
    return isempty( result ) ? Vector{String}() :
        string.( result[:, Symbol( mpSim.idKey )] )

end  # getActiveAtTime( mpSim, timePoint, subpopulation, fullPop )

function getSubpopulationStateAtTime!( activeAtTime::DataFrame, mpSim::MPsim,
    timePoint::Float64 )

    # Get the value of the attributes at time tPoint.
    activeIDs = activeAtTime[:, Symbol( mpSim.idKey )]
    queryCmd = string( "SELECT `", mpSim.idKey, "`, attribute, ", mpSim.valName,
        ", max(timeIndex) lastChangeTime FROM `", mpSim.histDBname, "` WHERE",
        "\n    timeIndex <= ", timePoint, " AND ",
        "\n    attribute IS NOT 'status' AND ",
        "\n   `", mpSim.idKey, "` IN ( '",
        join( activeIDs, "', '" ) , "')",
        "\n    GROUP BY `", mpSim.idKey, "`, attribute",
        "\n    ORDER BY attribute, `", mpSim.idKey, "`" )
    currentAttributeVals = DataFrame( DBInterface.execute( mpSim.simDB,
        queryCmd ) )

    # Reconstruct the state of the personnel members satisfying the node and
    #   time conditions.
    attributes = unique( currentAttributeVals[:, :attribute] )

    for attribute in attributes
        attributeInds = currentAttributeVals[:, :attribute] .== attribute
        activeAtTime[:, Symbol( attribute )] =
            currentAttributeVals[attributeInds, Symbol( mpSim.valName )]
    end  # for attrib in attribshelp

end  # getSubpopulationStateAtTime!( activeAtTime, mpSim, timePoint )


text:
__precompile__(true)
module GLVisualize

using Compat
using GLFW
using GLWindow
using GLAbstraction
using ModernGL
using FixedSizeArrays
using GeometryTypes
using ColorTypes
using Colors
using Reactive
using Quaternions
using FixedPointNumbers
using FileIO
using Packing
using SignedDistanceFields
using FreeType
using Iterators
import ColorVectorSpace

import Images
using Base.Markdown
using Compat
import Compat.unsafe_wrap
import Compat.String
import Compat.unsafe_string

if VERSION < v"0.5.0-dev+4612"
	function Base.checkbounds(::Type{Bool}, array::AbstractArray, indexes...)
		checkbounds(Bool, size(array), indexes...)
	end
end

typealias GLBoundingBox AABB{Float32}

import Base: merge, convert, show

if VERSION < v"0.5.0-dev+4612"
	function Base.checkbounds(::Type{Bool}, array::AbstractArray, indexes...)
		checkbounds(Bool, size(array), indexes...)
	end
else
	import Base: view
end

export renderloop

function assetpath(folders...)
    path = joinpath(dirname(@__FILE__), "..", "assets", folders...)
    isfile(path) || isdir(path) || error("Could not locate file at $path")
    path
end
loadasset(folders...; kw_args...) = load(assetpath(folders...); kw_args...)
export assetpath, loadasset

include("FreeTypeAbstraction.jl")
using .FreeTypeAbstraction

include("StructsOfArrays.jl")
using .StructsOfArrays

include("types.jl")
export CIRCLE, RECTANGLE, ROUNDED_RECTANGLE, DISTANCEFIELD, TRIANGLE

include("config.jl")
using .Config
import .Config: default

include("boundingbox.jl")

include("visualize_interface.jl")
export _view #push renderobject into renderlist of the default screen, or supplied screen
export visualize    # Visualize an object
export visualize_default # get the default parameter for a visualization

include("utils.jl")
export y_partition
export x_partition
export loop, bounce
export clicked, dragged_on, is_hovering
export OR, AND, isnotempty
export color_lookup

include("renderloop.jl")

include("texture_atlas.jl")

include(joinpath("gui", "color_chooser.jl"))
include(joinpath("gui", "numbers.jl"))
include(joinpath("gui", "line_edit.jl"))
include(joinpath("gui", "buttons.jl"))
include(joinpath("gui", "options.jl"))

export widget # edits some value, name should be changed in the future!

include(joinpath("visualize", "lines.jl"))
include(joinpath("visualize", "containers.jl"))
include(joinpath("visualize", "image_like.jl"))
include(joinpath("visualize", "mesh.jl"))
include(joinpath("visualize", "particles.jl"))
include(joinpath("visualize", "surface.jl"))
include(joinpath("visualize", "text.jl"))

include("camera.jl")
export cubecamera

# Compose/Gadfly only work for 0.4 right now
if VERSION.minor == 4
include("compose_backend.jl")
end

include("videotool.jl")
export create_video

include("documentation.jl")
export get_docs

#include("glv_userimg.jl")

end # module


text:
function bipartite_matching_primal_dual(rp::Vector{Int64}, ci::Vector{Int64},
                    ai::Vector{T}, m::Int64, n::Int64) where T

    # variables used for the primal-dual algorithm
    # normalize ai values # updated on 2-19-2019
    #ai ./= maximum(abs.(ai))
    alpha=zeros(Float64,m)
    bt=zeros(Float64,m+n)#beta
    queue=zeros(Int64,m)
    t=zeros(Int64,m+n)
    match1=zeros(Int64,m)
    match2=zeros(Int64,m+n)
    tmod = zeros(Int64,m+n)
    ntmod=0

    # initialize the primal and dual variables
    for i=1:m
        for rpi=rp[i]:rp[i+1]-1
            if ai[rpi] > alpha[i]
               alpha[i]=ai[rpi]
            end
        end
    end

    # dual variables (bt) are initialized to 0 already
    # match1 and match2 are both 0, which indicates no matches

    i=1
    while i<=m
        for j=1:ntmod
            t[tmod[j]]=0
        end
        ntmod=0
        # add i to the stack
        head=1
        tail=1
        queue[head]=i
        while head <= tail && match1[i]==0
            k=queue[head]
            for rpi=rp[k]:rp[k+1]-1
                j = ci[rpi]
                if ai[rpi] < alpha[k] + bt[j] - 1e-8
                    continue
                end # skip if tight
                if t[j]==0
                    tail=tail+1
                    if tail <= m
                        queue[tail]=match2[j]
                    end
                    t[j]=k
                    ntmod=ntmod+1
                    tmod[ntmod]=j
                    if match2[j]<1
                        while j>0
                            match2[j]=t[j]
                            k=t[j]
                            temp=match1[k]
                            match1[k]=j
                            j=temp
                        end
                        break
                    end
                end
            end
            head=head+1
        end
        if match1[i] < 1
            theta=Inf
            for j=1:head-1
                t1=queue[j]
                for rpi=rp[t1]:rp[t1+1]-1
                    t2=ci[rpi]
                    if t[t2] == 0 && alpha[t1] + bt[t2] - ai[rpi] < theta
                        theta = alpha[t1] + bt[t2] - ai[rpi]
                    end
                end
            end
            for j=1:head-1
                alpha[queue[j]] -= theta
            end
            for j=1:ntmod
                bt[tmod[j]] += theta
            end
            continue
        end
        i=i+1
    end
    val=0
    for i=1:m
        for rpi=rp[i]:rp[i+1]-1
            if ci[rpi]==match1[i]
                val=val+ai[rpi]
            end
        end
    end
    noute = 0
    for i=1:m
        if match1[i]<=n
            noute=noute+1
        end
    end

    M_output = MatrixNetworks.Matching_output(m,n,val,noute,match1)
    return M_output
end


text:
# Implementation taken from p316 RL Textbook and incorporating TB
Base.@kwdef mutable struct ETB{O, T<:AbstractTraceUpdate} <: LearningUpdate
    lambda::Float64
    opt::O
    hdpi::HordeDPI
    clip_threshold::Float64
    normalize_interest::Bool
    trace::T = AccumulatingTraces()
    e::IdDict = IdDict()
    followon::IdDict = IdDict()
    prev_discounts::IdDict = IdDict()
    emphasis_logging::Any = zeros(4) #I'm lazy and this is used for tracking hard to recalculate values
    rho_logging::Any = zeros(4) #I'm lazy and this is used for tracking hard to recalculate values
end

function get_demon_parameters(lu::ETB, learner, demons, obs, state, action, next_obs, next_state, next_action, env_reward)
    # C, next_discounts, _ = get(demons, obs, action, next_obs, next_action)
    C, next_discounts, _ = get(demons; state_t = obs, action_t = action, state_tp1 = next_obs, action_tp1 = next_action, reward = env_reward)
    target_pis = get_demon_pis(demons, learner.num_actions, state, obs)
    next_target_pis = get_demon_pis(demons, learner.num_actions, next_state, next_obs)
    C, next_discounts, target_pis, next_target_pis
end

function update!(lu::ETB,
                 learner::QLearner{M, LU},
                 demons,
                 obs,
                 next_obs,
                 state,
                 action,
                 next_state,
                 next_action,
                 is_terminal,
                 behaviour_pi_func,
                 env_reward) where {M<:AbstractMatrix, LU<:ETB}

    weights = learner.model
    λ = lu.lambda

    behaviour_pis = behaviour_pi_func(state, obs)
    C, next_discounts, target_pis, next_target_pis =
        get_demon_parameters(lu,
                             learner,
                             demons,
                             obs,
                             state,
                             action,
                             next_obs,
                             next_state,
                             next_action,
                             env_reward)

    discounts = get!(()->zero(next_discounts), lu.prev_discounts, learner)::typeof(next_discounts)
    e = get!(()->zero(weights), lu.e, weights)::typeof(weights)
    followon = get!(()->zeros(learner.num_actions), lu.followon, weights)::typeof(zeros(learner.num_actions))

    #Broadcast the policy and pseudotermination of each demon across the actions
    inds = get_action_inds(action, learner.num_actions, learner.num_demons)
    state_action_row_ind = inds

    # interest = get_interest(learner, lu, obs, action)

    # using dpi for interest instead
    interest = lu.hdpi(obs, action)
    if (lu.normalize_interest)
        interest = [ceil(i) for i in interest]
    end

    # getting IS ratio
    if (behaviour_pis[action] == 0)
        ρ = zeros(size(target_pis[:, action]))
    else
        ρ = target_pis[:, action] / behaviour_pis[action]
    end

    # accumulating followon
    followon[:] .+= interest

    # Get Emphasis vector - size is # demons
    # Correcting with ρ in the beginning since it is the action-value emphasis rather than the state-value emphasis
    emphasis = ρ .* (λ * behaviour_pis[action] * interest + (1 - λ * behaviour_pis[action]) * followon)

    # I'm lazy, here we're just logging the emphasis and rho as variables.
    lu.emphasis_logging = emphasis
    lu.rho_logging = ρ

    clamp!(emphasis, 0, lu.clip_threshold)

    # Update eligibility trace
    update_trace!(lu.trace,
                  e,
                  state,
                  λ,
                  repeat(discounts, inner = learner.num_actions),
                  repeat(target_pis[:,action], inner = learner.num_actions),
                  inds;
                  emphasis=emphasis)

    # decaying followon
    followon[:] .*= ρ .* next_discounts

    pred = learner(next_state)
    Qs = reshape(pred, (learner.num_actions, learner.num_demons))'

    # Target Pi is num_demons  x num_actions
    backup_est_per_demon = vec(sum((Qs .* next_target_pis ), dims = 2))
    target = C + next_discounts .* backup_est_per_demon

     # TD error per demon is the td error on Q
    td_err = target - (weights * state)[inds]
    td_err_across_demons = repeat(td_err, inner=learner.num_actions)
    # update discounts
    discounts .= next_discounts

    if lu.opt isa Auto
        next_state_action_row_ind = get_action_inds(next_action, learner.num_actions, learner.num_demons)
        state_discount = zero(e)
        state_discount[state_action_row_ind,:] .+= state'
        state_discount[next_state_action_row_ind,:] .-= next_discounts * next_state'
        abs_phi = abs.(e)
        update!(lu.opt, weights, e, td_err_across_demons, abs_phi .* max.(state_discount, abs_phi), learner.num_demons, learner.num_actions)
    else
        Flux.Optimise.update!(lu.opt, weights,  -(e .* td_err_across_demons))
    end
end

function update!(lu::ETB,
                 learner::Union{SRLearner,GPI},
                 demons,
                 obs,
                 next_obs,
                 state::SparseVector,
                 action,
                 next_state,
                 next_action,
                 is_terminal,
                 behaviour_pi_func,
                 env_reward)

    ψ = learner.ψ
    w = learner.r_w

    C, next_discounts, target_pis, next_target_pis =
        get_demon_parameters(lu,
                             learner,
                             demons,
                             obs,
                             state,
                             action,
                             next_obs,
                             next_state,
                             next_action,
                             env_reward)

    discounts = get!(()->zero(next_discounts), lu.prev_discounts, learner)::typeof(next_discounts)
    e_nz = get!(()->Int[], lu.e, learner)::Array{Int, 1}
    e_w_nz = get!(()->Int[], lu.e, learner)::Array{Int, 1}
    e_ψ = get!(()->zero(ψ), lu.e, ψ)::typeof(ψ)
    e_w = get!(()->zero(w), lu.e, w)::typeof(w)
    λ = lu.lambda

    next_active_state_action = get_active_action_state_vector(next_state, next_action,length(next_state), learner.num_actions)
    active_state_action = get_active_action_state_vector(state, action, length(state), learner.num_actions)

    projected_state = learner.feature_projector(obs, action, next_obs)

    (reward_C, SF_C) = C[1:learner.num_tasks] , C[learner.num_tasks + 1:end]
    (reward_discounts, SF_discounts) = discounts[1:learner.num_tasks], discounts[learner.num_tasks+1:end]
    (reward_next_discounts, SF_next_discounts) = next_discounts[1:learner.num_tasks], next_discounts[learner.num_tasks+1:end]
    (reward_target_pis, SF_target_pis) = target_pis[1:learner.num_tasks,:], target_pis[learner.num_tasks+1:end,:]
    (reward_next_target_pis, SF_next_target_pis) = next_target_pis[1:learner.num_tasks,:], next_target_pis[learner.num_tasks+1:end, :]

    # Followon
    followon = get!(()->zeros(learner.num_demons), lu.followon, weights)::typeof(zeros(learner.num_demons))

    # Setting interest to constant ones for now for each demon
    # interest = get_interest(learner, lu, obs, action)

    demon_interest = lu.hdpi(obs, action)

    if (lu.normalize_interest)
        demon_interest = [ceil(i) for i in demon_interest]
    end

    num_repeat = convert(Integer, (learner.num_demons - learner.num_tasks) / learner.num_tasks)
    if (rem(learner.num_demons - learner.num_tasks, learner.num_tasks) != 0)
        println("The number to repeat for SF demons aren't even with the number of tasks for SR. Something might be wrong :(")
    end

    SF_interest = repeat(demon_interest, inner=num_repeat)
    reward_interest = demon_interest

    interest = vcat(reward_interest, SF_interest)

    # getting IS ratio
    behaviour_pis = behaviour_pi_func(state, obs)
    if (behaviour_pis[action] == 0)
        ρ = zeros(size(target_pis[:, action]))
    else
        ρ = target_pis[:, action] / behaviour_pis[action]
    end
    reward_ρ, SF_ρ = ρ[1:learner.num_tasks], ρ[learner.num_tasks+1:end]

    # accumulating followon
    followon[:] .+= interest

    # Get Emphasis vector - size is # demons
    # Correcting with ρ in the beginning since it is the action-value emphasis rather than the state-value emphasis
    emphasis = ρ .* (λ * behaviour_pis[action] * interest + (1 - λ * behaviour_pis[action]) * followon)

    # I'm lazy, here we're just logging the emphasis and rho as variables.
    lu.emphasis_logging = emphasis
    lu.rho_logging = ρ

    clamp!(emphasis, 0, lu.clip_threshold)
    reward_emphasis, SF_emphasis = emphasis[1:learner.num_tasks], emphasis[learner.num_tasks+1:end]

    # Update Traces: See update_utils.jl
    update_trace!(lu.trace, e_ψ, active_state_action, λ, SF_discounts, SF_target_pis[:, action]; emphasis=SF_emphasis)
    update_trace!(lu.trace, e_w, projected_state, λ, reward_discounts, reward_target_pis[:, action]; emphasis=reward_emphasis)
    # e_nz = e_nz ∪ active_state_action.nzind
    # e_w_nz = e_w_nz ∪ projected_state_action.nzind
    if λ == 0.0
        e_nz = active_state_action.nzind
        e_w_nz = projected_state.nzind
    else
        e_nz = e_nz ∪ active_state_action.nzind
        e_w_nz = e_w_nz ∪ projected_state.nzind
    end

    # decaying followon
    followon[:] .*= ρ .* next_discounts

    pred = ψ * next_active_state_action
    reward_feature_backup = zeros(length(SF_C))
    for a in 1:learner.num_actions
        next_possible_active_state_action = get_active_action_state_vector(next_state, a, length(next_state), learner.num_actions)
        reward_feature_backup += SF_next_target_pis[:,a] .* (ψ * next_possible_active_state_action)
    end

    target = SF_C + SF_next_discounts .* reward_feature_backup
    td_err = target - ψ * active_state_action

    pred_err = reward_C - w * projected_state

    # This should always be true as this is immediate next step prediction which is equivalent to having discounts of 0 for all states
    @assert sum(reward_discounts) == 0
    # TD err is applied across rows

    if lu.opt isa Auto
        # next_state_action_row_ind = get_action_inds(next_action, learner.num_actions, learner.num_demons)
        state_discount = -SF_next_discounts * next_active_state_action'
        state_discount .+= active_state_action'
        abs_ϕ_ψ = if λ == 0.0
            abs.(repeat(active_state_action, outer=(1, length(td_err)))')
        else
            abs.(e_ψ)
        end
        z = abs_ϕ_ψ .* max.(abs_ϕ_ψ, state_discount)
        update!(lu.opt, ψ, e_ψ, td_err, z,  learner.num_demons - learner.num_tasks, 1)

        # reward state discount is always  0 so no need for next_state_action.
        # state_discount_r = -reward_next_discounts * projected_next_state_action'
        state_discount_r = projected_state'
        abs_ϕ_w = if λ == 0.0
            abs.(repeat(projected_state, outer=(1, length(pred_err)))')
        else
            abs.(e_w)
        end
        z_r = abs_ϕ_w .* max.(abs_ϕ_w, state_discount_r)
        update!(lu.opt, w, e_w, pred_err, z_r, learner.num_tasks, 1)
    elseif lu.opt isa Flux.Descent
        α = lu.opt.eta
        if λ == 0.0
            ψ[:, active_state_action.nzind] .+= (α  * td_err) * active_state_action.nzval'
            w .= w .+ α * pred_err * projected_state'
        else
            ψ[:, e_nz] .+= (α  * td_err) .* e_ψ[:, e_nz]
            w[:, e_w_nz] .+= (α  * pred_err) .* e_w[:, e_w_nz]
        end
    else
        update!(lu.opt, ψ, - e_ψ .* td_err)
        update!(lu.opt, w, - e_w .* pred_err)
    end
    discounts .= next_discounts
end

function zero_eligibility_traces!(lu::ETB)
    for (k, v) ∈ lu.e
        if eltype(v) <: Integer
            lu.e[k] = Int[]
        elseif eltype(v) <: Vector #Used for LSTD
            lu.e[k] = zero.(v)
        else
            v .= 0
        end
    end

end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
struct Element
    name::AbstractString
    refstate::AbstractString # reference state
    mass::Float64
    H298::Float64
    S298::Float64
end

abstract type AbstractGFunction end

mutable struct GFunction <: AbstractGFunction
    name::AbstractString
    temp::Tuple{<:Real,<:Real}
    temps::Vector{<:Real}
    exprs::Vector{AbstractString}
    funcname::AbstractString
    funcstr::AbstractString
end

function GFunction(name::AbstractString,
                   temp::Tuple{<:Real,<:Real},
                   funcstr::AbstractString)
    return GFunction(name, temp, Real[], AbstractString[], getfuncname(name), funcstr)
end

function GFunction(name::AbstractString,
                   temps::Vector{<:Real},
                   exprs::Vector{<:AbstractString})
    temp = (minimum(temps), maximum(temps))
    return GFunction(name, temp, temps, exprs, getfuncname(name), "")
end

const Sublattice = Vector{AbstractString}
const Constitution = Vector{Sublattice}

mutable struct Parameter <: AbstractGFunction
    name::AbstractString
    symbol::Char
    comb::Constitution
    order::Int
    temp::Tuple{<:Real,<:Real}
    temps::Vector{<:Real}
    exprs::Vector{AbstractString}
    funcname::AbstractString
    funcstr::AbstractString
end

function Parameter(name::AbstractString,
                   symbol::Char,
                   comb::Vector{<:Vector},
                   order::Int,
                   temps::Vector{<:Real},
                   exprs::Vector{<:AbstractString},
                   funcname::AbstractString)
    temp = (minimum(temps), maximum(temps))
    return Parameter(name, symbol, comb, order, temp, temps, exprs, funcname, "")
end

function Parameter(name::AbstractString,
                   symbol::Char,
                   comb::Vector{<:Vector},
                   order::Int,
                   temp::Tuple{<:Real,<:Real},
                   funcname::AbstractString,
                   funcstr::AbstractString)
    return Parameter(name, symbol, comb, order, temp,
                     Real[], AbstractString[], funcname, funcstr)
end

mutable struct Phase
    name::AbstractString
    state::Char
    type::AbstractString
    sites::Vector{<:Real}
    cons::Constitution
    params::Vector{Parameter}
end

function Phase(name::AbstractString,
               state::Char,
               sites::Vector{<:Real})
    return Phase(name, state, "%", sites, Sublattice[], Parameter[])
end

function Phase(name::AbstractString,
               state::Char,
               sites::Vector{<:Real},
               cons::Vector{<:Vector},
               params::Vector{Parameter})
    return Phase(name, state, "%", sites, cons, params)
end

function Phase(name::AbstractString,
               state::Char,
               type::AbstractString,
               sites::Vector{<:Real})
    phas = Phase(name, state, type, sites, Sublattice[], Parameter[])
end

struct TypeDefinition
    code::Char
    key::AbstractString # SEQ, GES, etc.
    args::AbstractString
end

struct Database
    elems::Vector{Element}
    funcs::Vector{GFunction}
    phass::Vector{Phase}
    types::Vector{TypeDefinition}
end

function Database(elems::Vector{Element},
                  funcs::Vector{GFunction},
                  phass::Vector{Phase})
    return Database(elems, funcs, phass, TypeDefinition[])
end

function Base.display(db::Database)
    println("Database Summary:")
    nelem = length(db.elems)
    println("\tElements: $nelem")
    nfunc = length(db.funcs)
    println("\tFunctions: $nfunc")
    nphas = length(db.phass)
    println("\tPhases: $nphas")
end

function Base.print(db::Database)
    println("Database:")

    nelem = length(db.elems)
    println("\tElements: $nelem")
    for elem in db.elems
        println("\t\t$(elem.name)")
    end

    nfunc = length(db.funcs)
    println("\tFunctions: $nfunc")
    for func in db.funcs
        println("\t\t$(func.name)")
    end

    nphas = length(db.phass)
    println("\tPhases: $nphas")
    for k in 1:nphas
        phas = db.phass[k]
        print("\t\t$k: $(phas.name); ")
        println(constitutionstring(phas))
        for param in phas.params
            println("\t\t\t" * param.symbol * '(' * localname(param) * ')')
        end
    end
end


text:
module FETermsWithAutodiff

using LinearAlgebra
using Gridap.FESpaces
using Gridap.Arrays
using Gridap.Fields
using Gridap.Geometry
using Gridap.CellData

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

V = FESpace(model=model,valuetype=Float64,reffe=:Lagrangian,order=1,conformity=:H1)
U = TrialFESpace(V)

dv = get_cell_basis(V)
du = get_cell_basis(U)
uh = FEFunction(U,rand(num_free_dofs(U)))

trian = Triangulation(model)
quad = CellQuadrature(trian,2)

res(uh,v) = ∇(uh)⋅∇(v)
jac(uh,u,v) = ∇(u)⋅∇(v)

t_Ω = FETerm(res,jac,trian,quad)
t_auto_Ω = FETerm(res,trian,quad)

cell_r = get_cell_residual(t_Ω,uh,dv)
cell_j = get_cell_jacobian(t_Ω,uh,du,dv)

cell_r_auto = get_cell_residual(t_auto_Ω,uh,dv)
cell_j_auto = get_cell_jacobian(t_auto_Ω,uh,du,dv)

test_array(cell_r_auto,cell_r)
test_array(cell_j_auto,cell_j)

trian_Γ = BoundaryTriangulation(model)
quad_Γ = CellQuadrature(trian_Γ,2)

t_Γ = FETerm(res,jac,trian_Γ,quad_Γ)
t_auto_Γ = FETerm(res,trian_Γ,quad_Γ)

cell_r_Γ = get_cell_residual(t_Γ,uh,dv)
cell_j_Γ = get_cell_jacobian(t_Γ,uh,du,dv)

cell_r_Γ_auto = get_cell_residual(t_auto_Γ,uh,dv)
cell_j_Γ_auto = get_cell_jacobian(t_auto_Γ,uh,du,dv)

test_array(cell_r_Γ_auto,cell_r_Γ)
test_array(cell_j_Γ_auto,cell_j_Γ)

ener(uh) = 0.5*∇(uh)⋅∇(uh)
t_auto_Ω = FEEnergy(ener,trian,quad)

cell_r_auto = get_cell_residual(t_auto_Ω,uh,dv)
cell_j_auto = get_cell_jacobian(t_auto_Ω,uh,du,dv)

test_array(cell_r_auto,cell_r)
test_array(cell_j_auto,cell_j)

t_auto_Γ = FEEnergy(ener,trian_Γ,quad_Γ)

cell_r_Γ_auto = get_cell_residual(t_auto_Γ,uh,dv)
cell_j_Γ_auto = get_cell_jacobian(t_auto_Γ,uh,du,dv)

test_array(cell_r_Γ_auto,cell_r_Γ)
test_array(cell_j_Γ_auto,cell_j_Γ)

end # module


text:
module TestNChisq

using Base.Test
using Distributions
using TimeData
using NLopt
using JuMP

import Econometrics

## reload(joinpath(Pkg.dir(), "Econometrics", "src", "Econometrics.jl"))

## constructors / display
##-----------------------

ncs = Econometrics.NChiSq([3.2, 0.2])
ncs = Econometrics.NChiSq(3.2, 0.2)

println("\nTesting display output:")
display(ncs)

@test_throws Exception Econometrics.NChiSq(-0.2, 0.2)
@test_throws Exception Econometrics.NChiSq(0.2, -1.2)

## dof
##----

ncs = Econometrics.NChiSq([3.2, 0.2])
@test Econometrics.dof(ncs) == 3.2

## getParams
##----------

@test Econometrics.getParams(ncs) == (3.2, 0.2)

## rand
##-----

nSim = 5000
ncs = Econometrics.NChiSq([3.2, 0.2])
simVals = rand(ncs, nSim)

## pdf
##----

d = Econometrics.NChiSq([3.2, 0.2])
pdfVal = Econometrics.pdf(d, 0.2)

@test isa(pdfVal, Float64)

pdfVals = Econometrics.pdf(d, simVals)

@test isa(pdfVals, Array{Float64, 1})
@test all(pdfVals .>= 0)

## cdf
##----

d = Econometrics.NChiSq([3.2, 0.2])
cdfVal = Econometrics.cdf(d, 0.2)

@test isa(cdfVal, Float64)

cdfVals = Econometrics.cdf(d, simVals)

@test isa(cdfVals, Array{Float64, 1})
@test all(cdfVals .>= 0)
@test all(cdfVals .<= 1)

## quantile
##---------

d = Econometrics.NChiSq([3.2, 0.2])
pVal = Econometrics.quantile(d, 0.2)

@test isa(pVal, Float64)

## test cdf(quantile(x)) and quantile(cdf(x))
##-------------------------------------------

nSim = 5000
d = Econometrics.NChiSq([3.2, 0.2])
simVals = rand(d, nSim)
cdfVals = Econometrics.cdf(d, simVals)
pVals = Econometrics.quantile(d, cdfVals)
@test_approx_eq pVals simVals

grid = [0.01:0.01:0.99]
pVals = Econometrics.quantile(d, grid)
cdfVals = Econometrics.cdf(d, pVals)
@test_approx_eq cdfVals grid

## nllh
##-----

nSim = 5000
d = Econometrics.NChiSq([3.2, 0.2])
simVals = rand(d, nSim)

nllhVal1 = Econometrics.nllh(d, simVals)
nllhVal2 = Econometrics.nllh_nchis([Econometrics.getParams(d)...], simVals)

@test_approx_eq nllhVal1 nllhVal2

@test_throws Exception Econometrics.nllh_nchis([3.2, 0.2, 1.2, 4.2], simVals)

## fit
##----

nSim = 1000
params = [3.2, 0.2]
d = Econometrics.NChiSq(params)

@time dHat = Econometrics.fit(Econometrics.NChiSq, simVals)

@test isa(dHat, Econometrics.NChiSq)

## compare parameters
[params'
[Econometrics.getParams(dHat)...]']

end


text:
module Rank

using DataStructures
using Base: tail

export extrema_filter

# Max-min filter

# This is an implementation of the Lemire max-min filter
# http://arxiv.org/abs/cs.DS/0610046

# Monotonic wedge
immutable Wedge{T}
    L::CircularDeque{T}
    U::CircularDeque{T}
end
(::Type{Wedge{T}}){T}(n::Integer) = Wedge(CircularDeque{T}(n), CircularDeque{T}(n))

function Base.push!(W::Wedge, i::Integer)
    push!(W.L, i)
    push!(W.U, i)
    W
end

function addtoback!(W::Wedge, A, i, J)
    mn, mx = A[i, J]
    @inbounds while !isempty(W.L) && mn < A[back(W.L), J][1]
        pop!(W.L)
    end
    @inbounds while !isempty(W.U) && mx > A[back(W.U), J][2]
        pop!(W.U)
    end
    push!(W.L, i)
    push!(W.U, i)
    W
end

function Base.empty!(W::Wedge)
    empty!(W.L)
    empty!(W.U)
    W
end

@inline function getextrema(A, W::Wedge, J)
    (A[front(W.L), J][1], A[front(W.U), J][2])
end

"""
    extrema_filter(A, window) --> Array{(min,max)}

Calculate the running min/max over a window of width `window[d]` along
dimension `d`, centered on the current point. The returned array has
the same indices as the input `A`.
"""
function extrema_filter{T,N}(A::AbstractArray{T,N}, window::NTuple{N,Integer})
    _extrema_filter!([(a,a) for a in A], window...)
end
extrema_filter(A::AbstractArray, window::AbstractArray) = extrema_filter(A, (window...,))
extrema_filter(A::AbstractArray, window) = error("`window` must have the same number of entries as dimensions of `A`")

extrema_filter{T,N}(A::AbstractArray{T,N}, window::Integer) = extrema_filter(A, ntuple(d->window, Val{N}))

function _extrema_filter!(A::Array, w1, w...)
    if w1 > 1
        a = first(A)
        cache = ntuple(i->a, w1>>1)
        _extrema_filter1!(A, w1, cache)
    end
    _extrema_filter!(permutedims(A, [2:ndims(A);1]), w...)
end
_extrema_filter!(A::Array) = A

# Extrema-filtering along "columns" (dimension 1). This implements Lemire
# Algorithm 1, with the following modifications:
#   - multidimensional array support by looping over trailing dimensions
#   - working with min/max pairs rather than plain values, to
#     facilitate multidimensional processing
#   - output for all points of the array, handling the edges as max-min
#     over halfwindow on either side
function _extrema_filter1!{T}(A::AbstractArray{Tuple{T,T}}, window::Int, cache)
    # Initialise the internal wedges
    # U[1], L[1] are the location of the global (within the window) maximum and minimum
    # U[2], L[2] are the maximum and minimum over (U1, end] and (L1, end], respectively
    W = Wedge{Int}(window+1)
    tmp = Array{Tuple{T,T}}(window)
    c = z = first(cache)

    inds = indices(A)
    inds1 = inds[1]
    halfwindow = window>>1
    iw = min(last(inds1), first(inds1)+window-1)
    for J in CartesianRange(tail(inds))
        empty!(W)
        # Leading edge. We can't overwrite any values yet in A because
        # we'll need them again in later computations.
        for i = first(inds1):iw
            addtoback!(W, A, i, J)
            c, cache = cyclecache(cache, getextrema(A, W, J))
        end
        # Process the rest of the "column"
        for i = iw+1:last(inds1)
            A[i-window, J] = c
            if i == window+front(W.U)
                shift!(W.U)
            end
            if i == window+front(W.L)
                shift!(W.L)
            end
            addtoback!(W, A, i, J)
            c, cache = cyclecache(cache, getextrema(A, W, J))
        end
        for i = last(inds1)-window+1:last(inds1)-1
            if i >= first(inds1)
                A[i, J] = c
            end
            if i == front(W.U)
                shift!(W.U)
            end
            if i == front(W.L)
                shift!(W.L)
            end
            c, cache = cyclecache(cache, getextrema(A, W, J))
        end
        A[last(inds1), J] = c
    end
    A
end

# This is slightly faster than a circular buffer
@inline cyclecache(b, x) = b[1], (Base.tail(b)..., x)

end


text:
# Functions for discretizing continuous values using various algorithms

using Discretizers

export get_bin_ids!, get_frequencies, get_frequencies_from_bin_ids

# Parameters:
# 	- values_x, arrays of floats (multiple arrays supported)
function get_root_n(values_x...)
	return round(Int, sqrt(size(values_x[1])[1]))
end

# Parameters:
# 	- values_x, arrays of floats (multiple arrays supported)
# 	- mode, string
#	- number_of_bins, number
function get_frequencies(mode, number_of_bins, values_x)
	bin_ids = Array{Int}(size(values_x))
	number_of_bins = get_bin_ids!(values_x, mode, number_of_bins, bin_ids)
	frequencies = zeros(Int, (number_of_bins, 1))
	for id in bin_ids
		frequencies[id] += 1
	end
	return frequencies
end
function get_frequencies(mode, number_of_bins, values_x, values_y)
	# n is the number of data points
	n = size(values_x)[1]
	bin_ids_x = Array{Int}(size(values_x))
	bin_ids_y = Array{Int}(size(values_y))
	number_of_bins_x = get_bin_ids!(values_x, mode, number_of_bins, bin_ids_x)
	number_of_bins_y = get_bin_ids!(values_y, mode, number_of_bins, bin_ids_y)
	frequencies = zeros(Int, (number_of_bins_x, number_of_bins_y))
	for i in 1:n
		frequencies[bin_ids_x[i], bin_ids_y[i]] += 1
	end
	return frequencies
end
function get_frequencies(mode, number_of_bins, values_x, values_y, values_z)
	# n is the number of data points
	n = size(values_x)[1]
	bin_ids_x = Array{Int}(size(values_x))
	bin_ids_y = Array{Int}(size(values_y))
	bin_ids_z = Array{Int}(size(values_z))
	number_of_bins_x = get_bin_ids!(values_x, mode, number_of_bins, bin_ids_x)
	number_of_bins_y = get_bin_ids!(values_y, mode, number_of_bins, bin_ids_y)
	number_of_bins_z = get_bin_ids!(values_z, mode, number_of_bins, bin_ids_z)
	frequencies = zeros(Int, (number_of_bins_x, number_of_bins_y, number_of_bins_z))
	for i in 1:n
		frequencies[bin_ids_x[i], bin_ids_y[i], bin_ids_z[i]] += 1
	end
	return frequencies
end
function get_frequencies(mode, number_of_bins, values_x...)
	# n is the number of data points
	n = size(values_x[1])[1]
	# d is the number of dimensions
	d = length(values_x)
	bin_ids = Array{Array{Int}}(d)
	all_number_of_bins = Array{Int}(d)
	for i in 1:d
		current_values = values_x[i]
		bin_ids[i] = Array{Int}(size(current_values))
		all_number_of_bins[i] = get_bin_ids!(current_values, mode, number_of_bins, bin_ids[i])
	end
	frequencies = zeros(Int, tuple(all_number_of_bins...))
	for i in 1:n
		frequencies[ntuple(x -> bin_ids[x][i], d)...] += 1
	end
	return frequencies
end

# Parameters:
# 	- bin_ids_x, array of ints
# 	- number_of_bins_x, number
function get_frequencies_from_bin_ids(bin_ids_x, number_of_bins_x)
	n = size(bin_ids_x)[1]
	frequencies = zeros(Int, number_of_bins_x)
	for i in 1:n
		frequencies[bin_ids_x[i]] += 1
	end
	return frequencies
end
function get_frequencies_from_bin_ids(bin_ids_x, bin_ids_y, number_of_bins_x, number_of_bins_y)
	n = size(bin_ids_x)[1]
	frequencies = zeros(Int, (number_of_bins_x, number_of_bins_y))
	for i in 1:n
		frequencies[bin_ids_x[i], bin_ids_y[i]] += 1
	end
	return frequencies
end
function get_frequencies_from_bin_ids(bin_ids_x, bin_ids_y, bin_ids_z, number_of_bins_x, number_of_bins_y, number_of_bins_z)
	n = size(bin_ids_x)[1]
	frequencies = zeros(Int, (number_of_bins_x, number_of_bins_y, number_of_bins_z))
	for i in 1:n
		frequencies[bin_ids_x[i], bin_ids_y[i], bin_ids_z[i]] += 1
	end
	return frequencies
end

# Parameters:
# 	- values_x, array of floats
# 	- mode, string
#	- number_of_bins, integer
#	- bin_ids, 1-dimensional array of bin ids for each value
function get_bin_ids!(values_x, mode, number_of_bins, bin_ids)
	min, max = extrema(values_x)
	if min == max
		# If values are all the same, assign them all to bin 1
		number_of_bins = 1
		bin_ids[1:end] = 1
	elseif mode == "uniform_width"
		bin_ids[1:end] = encode(LinearDiscretizer(binedges(DiscretizeUniformWidth(number_of_bins), values_x)), values_x)
	elseif mode == "binarize"
		number_of_bins = 2
		bin_ids[1:end] = map(v -> v == 0 ? 1 : 2, values_x)
	elseif mode == "uniform_count"
		try
			bin_ids[1:end] = encode(LinearDiscretizer(binedges(DiscretizeUniformCount(number_of_bins), reshape(values_x, length(values_x)))), values_x)
		catch
			bin_ids = encode(LinearDiscretizer(binedges(DiscretizeUniformWidth(number_of_bins), values_x)), values_x)
			println("Uniform count failed, fell back to uniform width")
		end
	elseif mode == "bayesian_blocks"
		try
			edges = binedges(DiscretizeBayesianBlocks(), values_x)
			bin_ids[1:end] = encode(LinearDiscretizer(edges), values_x)
			number_of_bins = length(edges) - 1
		catch
			bin_ids[1:end] = encode(LinearDiscretizer(binedges(DiscretizeUniformWidth(number_of_bins), values_x)), values_x)
			println("Bayesian blocks failed, fell back to uniform width")
		end
	else
		println("Mode doesn't exist, fell back to uniform width")
		bin_ids[1:end] = encode(LinearDiscretizer(binedges(DiscretizeUniformWidth(number_of_bins), values_x)), values_x)
	end

	return number_of_bins
end


text:
"""
test 1
"""
@recipe function plot(layer::T) where {T <: SimpleSDMLayer}
    eltype(layer) <: Number || throw(ArgumentError("Plotting is only supported for layers with number values ($(eltype(layer)))"))
    seriestype --> :heatmap
    if get(plotattributes, :seriestype, :heatmap) in [:heatmap, :contour]
        aspect_ratio --> 1
        xlims --> extrema(longitudes(layer))
        ylims --> extrema(latitudes(layer))
        lg = copy(layer.grid)
        replace!(lg, nothing => NaN)
        lg = convert(Matrix{Float64}, lg)
        longitudes(layer), latitudes(layer), lg
    elseif get(plotattributes, :seriestype, :histogram) in [:histogram, :density]
        collect(layer)
    elseif get(plotattributes, :seriestype, :surface) in [:surface, :wireframe]
        aspect_ratio --> 1
        xlims --> extrema(longitudes(layer))
        ylims --> extrema(latitudes(layer))
        lg = copy(layer.grid)
        replace!(lg, nothing => minimum(layer))
        lg = convert(Matrix{Float64}, lg)
        longitudes(layer), latitudes(layer), lg
    end
end

"""
test 2
"""
@recipe function plot(l1::FT, l2::ST) where {FT <: SimpleSDMLayer, ST <: SimpleSDMLayer}
    eltype(l1) <: Number || throw(ArgumentError("Plotting is only supported for layers with number values ($(eltype(l1)))"))
    eltype(l2) <: Number || throw(ArgumentError("Plotting is only supported for layers with number values ($(eltype(l2)))"))
    seriestype --> :scatter
    if get(plotattributes, :seriestype, :scatter) in [:scatter, :histogram2d]
        SimpleSDMLayers._layers_are_compatible(l1, l2)
        valid_i = findall(.!isnothing.(l1.grid) .& .!isnothing.(l2.grid))
        l1.grid[valid_i], l2.grid[valid_i]
    end
end


text:
"""
    pytype(x)

The Python `type` of `x`.
"""
pytype(x) = pynew(errcheck(@autopy x C.PyObject_Type(getptr(x_))))
pytype(name, bases, dict) = pybuiltins.type(name, ispy(bases) && pyistype(bases) ? pytuple((bases,)) : pytuple(bases), pydict(dict))
export pytype

"""
    pyclass(name, bases=(); members...)

Construct a new Python type with the given name, bases and members.

Equivalent to `pytype(name, bases, members)`.
"""
pyclass(name, bases=(); dict...) = pytype(name, bases, pystrdict_fromiter(dict))
export pyclass

pyistype(x) = pytypecheckfast(x, C.Py_TPFLAGS_TYPE_SUBCLASS)

pytypecheck(x, t) = (@autopy x t C.Py_TypeCheck(getptr(x_), getptr(t_))) == 1
pytypecheckfast(x, f) = (@autopy x C.Py_TypeCheckFast(getptr(x_), f)) == 1


text:
typ(n) = Symbol(string("Series", n))
elt(n) = Symbol(string("c", n))
mem(s,n) = Expr(:., s, Expr(:quote, elt(n)))

function generate_type(n::Integer)
  @assert n > 0
  LastType = typ(n-1)
  Typ = typ(n)

  defn = :(immutable $Typ{T} <: AbstractSeries{T, $n} c0::T end)

  for i = 1:n
    push!(defn.args[3].args, :($(elt(i))::T))
  end

  eval(defn)

  if (n > 1)
    @eval restrict(x::$Typ) = $LastType($([mem(:x, i) for i = 0:n-1]...))
    @eval polyder(x::$Typ) = $LastType($([:($i*$(mem(:x, i))) for i = 1:n]...))
    @eval polyint{T}(x::$LastType{T}) = $Typ(zero(T),
      $([:($(mem(:x, i-1))/$i) for i = 1:n]...)
    )
  else
    @eval restrict(x::$Typ) = x.c0
    @eval polyder(x::$Typ) = x.c1
    @eval polyint(x::Real) = Series1(zero(typeof(x)), x)
  end

  @eval $Typ(x::Real) = $Typ(x, $([:(zero(x)) for i = 1:n]...))

  syms = [gensym() for i=0:n]

  @eval series($(syms...)) = $Typ(promote($(syms...))...)

  @eval convert{T<:Real}(::Type{$Typ{T}}, x::Real) = $Typ(convert(T, x))

  #Horner's method
  polyval_body = mem(:x, n)
  for m = n-1:-1:0
    polyval_body = :($(mem(:x, m)) + eps*$polyval_body)
  end
  @eval polyval(x::$Typ, eps::Real) = $polyval_body

  @eval constant(x::$Typ) = x.c0

  @eval +(x::$Typ, y::$Typ) = $Typ(
    $([:($(mem(:x, i)) + $(mem(:y, i))) for i = 0:n]...)
  )
  @eval +{T<:Real, S}(c::T, x::$Typ{S}) = $Typ{promote_type(T, S)}(
    c + $(mem(:x, 0)),
    $([mem(:x, i) for i = 1:n]...)
  )
  @eval +{T<:Real, S}(x::$Typ{S}, c::T) = $Typ{promote_type(T, S)}(
    $(mem(:x, 0)) + c,
    $([mem(:x, i) for i = 1:n]...)
  )

  @eval -(x::$Typ) = $Typ($([:(-$(mem(:x, i))) for i = 0:n]...))

  @eval -(x::$Typ, y::$Typ) = $Typ(
    $([:($(mem(:x, i)) - $(mem(:y, i))) for i = 0:n]...)
  )
  @eval -{T<:Real, S}(c::T, x::$Typ{S}) = $Typ{promote_type(T, S)}(
    c - $(mem(:x, 0)),
    $([:(-$(mem(:x, i))) for i = 1:n]...)
  )
  @eval -{T<:Real, S}(x::$Typ{S}, c::T) = $Typ{promote_type(T, S)}(
    $(mem(:x, 0)) - c,
    $([:($(mem(:x, i))) for i = 1:n]...)
  )

  times_term(i) = :(+($([:($(mem(:x, j))*$(mem(:y, i-j))) for j = 0:i]...)))
  @eval *(x::$Typ, y::$Typ) = $Typ($([times_term(i) for i = 0:n]...))

  # Disambiguate bool multiplication. Begrudgingly...
  # See https://github.com/JuliaLang/julia/issues/5611
  @eval *(x::Bool, z::$Typ) = ifelse(x,z,zero(z))
  @eval *(z::$Typ, x::Bool) = x * z

  @eval *(c::Real, x::$Typ) = $Typ($([:(c*$(mem(:x, i))) for i = 0:n]...))
  @eval *(x::$Typ, c::Real) = $Typ($([:($(mem(:x, i))*c) for i = 0:n]...))

  @eval /(x::$Typ, c::Real) = $Typ($([:($(mem(:x, i))/c) for i = 0:n]...))

  # TODO probably a better way than eval to get a type from a string
  return eval(Typ)
end


text:
using Weave

markdown(file) = begin
  try
    weave(file, out_path = "$(pwd())/public")
  catch
    println("wait to be modified")
  end
end

watch(file) = begin
  watch_file(file)
  markdown(file)
  watch(file)
end

max_stars_repo_name: moldray/jQuant

max_issues_repo_name: moldray/jQuant

max_forks_repo_name: moldray/jQuant