text:
# TODO: Check math again!

"""
Marginal loglikelihood for sample i. Marginalizes over lambda and uses skew-t
likelihood. Used for updating beta.
"""
function marginal_loglike_beta(i::Char, state::State, data::Data)
  Nfinite = ref_Nfinite(data, 'T')
  Z = ref_Z(data, 'T')
  yfinite = ref_yfinite(data, 'T')

  loc = state.mu
  scale = state.sigma
  skew = state.phi
  df = state.nu

  eta = ref_eta(state, i)
  gamma = ref_gamma(state, i)

  f = logsumexp(Util.skewtlogpdf.(loc', scale', df', skew', yfinite) .+
                log.(eta)', dims=2)

  return Z * log(gamma) + Nfinite * log1p(-gamma) + sum(f)
end

"""
Marginal loglikelihood for sample i. Marginalizes over lambda and uses Normal
likelihood using latent variable representation of skew-t.  Used for updating
beta.
"""
function marginal_loglike_beta_latent_var(i::Char, state::State, data::Data)
  Nfinite = ref_Nfinite(data, 'T')
  yfinite = ref_yfinite(data, 'T')
  Z = ref_Z(data, 'T')

  K = length(state.mu)
  zeta = ref_zeta(state, 'T')
  v = ref_v(state, 'T')

  eta = ref_eta(state, i)
  gamma = ref_gamma(state, i)

  loc = state.mu' .+ state.psi' .* zeta
  scale = sqrt.(state.omega)' ./ sqrt.(v)
  g = logsumexp(normlogpdf.(loc, scale, yfinite) .+
                gammalogpdf.(state.nu' / 2, 2 ./ state.nu', v) .+
                log.(eta)', dims=2)

  return Z * log(gamma) + Nfinite * log1p(-gamma) + sum(g)
end

"""
Loglikelihood for sample i. Does not marginalize over lambda, and uses Normal
likelihood using latent variable representation of skew-t.
"""
function loglike_latent_var(i::Char, state::State, data::Data, use_star::Bool)
  Nfinite = ref_Nfinite(data, i)
  Z = ref_Z(data, i)
  yfinite = ref_yfinite(data, i)

  zeta = ref_zeta(state, i)
  v = ref_v(state, i)
  lambda = ref_lambda(state, i)

  j = (use_star && (i == 'T')) ? (state.beta ? 'T' : 'C') : i
  gamma = ref_gamma(state, j)

  loc = state.mu[lambda] + state.psi[lambda] .* zeta
  scale = sqrt.(state.omega)[lambda] ./ sqrt.(v)
  F = sum(normlogpdf.(loc, scale, yfinite))
  return Z * log(gamma) + Nfinite * log1p(-gamma) + sum(F)
end

"""
Loglikelihood for sample i. Does not marginalize over lambda, and uses Normal
likelihood using latent variable representation of skew-t. Arguments
are `state::State` and `data::Data` only.
"""
function loglike_latent_var(state::State, data::Data)
  llC = loglike_latent_var('C', state, data, false)
  llT = loglike_latent_var('T', state, data, true)
  return llC + llT
end

### TODO: Bayes Factors (Check all below here.) ###

# TODO: Marginal loglikelihood for computing Bayes factor
# and posterior odds. Check!
function marginal_loglike(data::Data, state::Dict{Symbol, Any})
  llC = marginal_loglike_C(data, state)
  llT = marginal_loglike_T(data, state)
  return llC + llT
end

# TODO: Check.
function marginal_loglike_C(data::Data, state::Dict{Symbol, Any})
  loc = state[:mu]
  scale = state[:sigma]
  skew = state[:phi]
  df = state[:nu]

  etaC = state[:etaC]
  gammaC = state[:gammaC]

  f = logsumexp(Util.skewtlogpdf.(loc', scale', df', skew', data.yC_finite) .+
                log.(etaC)', dims=2)

  return data.ZC * log(gammaC) + data.NC_finite * log1p(-gammaC) + sum(f)
end

# TODO: Check.
function marginal_loglike_T(data::Data, state::Dict{Symbol, Any})
  loc = state[:mu]
  scale = state[:sigma]
  skew = state[:phi]
  df = state[:nu]

  etaT_star = state[:beta] ? state[:etaT] : state[:etaC]
  gammaT_star = state[:beta] ? state[:gammaT] : state[:gammaC]

  f = logsumexp(Util.skewtlogpdf.(loc', scale', df', skew', data.yT_finite) .+
                log.(etaT_star)', dims=2)

  return data.ZT * log(gammaT_star) + data.NT_finite * log1p(-gammaT_star) + sum(f)
end

"""
Bayes factor in favor of model where β=1.
See: https://projecteuclid.org/download/pdfview_1/euclid.ba/1346158782
"""
function log_bayes_factor(data::Data,
                          chain0::Vector{Dict{Symbol, Any}},
                          chain1::Vector{Dict{Symbol, Any}})
  ll0 = [marginal_loglike(data, state) for state in chain0]
  ll1 = [marginal_loglike(data, state) for state in chain1]
  B0 = length(ll0)
  B1 = length(ll1)
  ll0_mean = logsumexp(-ll0) - log(B0)
  ll1_mean = logsumexp(-ll1) - log(B1)

  # Ratio of harmonic means.
  return ll0_mean - ll1_mean
end

function bayes_factor(data::Data,
                      chain0::Vector{Dict{Symbol, Any}},
                      chain1::Vector{Dict{Symbol, Any}})
  return exp(log_bayes_factor(data, chain0, chain1))
end

# TODO: Check.
"""
Posterior model probability of model where β=1.  To compute the posterior model
probability of the model where β=0, simply compute `1 - P(β=1 | y)`, where
`P(β=1 | y)` is the output to this function.
"""
function posterior_prob1(data::Data,
                         chain0::Vector{Dict{Symbol, Any}},
                         chain1::Vector{Dict{Symbol, Any}};
                         logpriorodds::Real=0)
  lbf = log_bayes_factor(data, chain0, chain1)
  return logistic(lbf + logpriorodds)
end

function posterior_prob1(lbf::Real; logpriorodds::Real=0)
  return logistic(lbf + logpriorodds)
end

# DIC.
function dic(chain, data)
  ll = [marginal_loglike(data, s) for s in chain]
  return MCMC.dic(ll)
end

function dic(chain0, chain1, pm1, data)
  @assert length(chain0) == length(chain1)
  B = length(chain0)
  chain = [pm1 > rand() ? chain1[b] : chain0[b] for b in 1:B]
  return MCMC.dic(chain, data)
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
module IndEco

# Write your package code here.

end

ext: jl

max_stars_repo_name: MathAero/IndEco

max_issues_repo_name: MathAero/IndEco

max_forks_repo_name: MathAero/IndEco


text:
using Graft
using StatsBase
import LightGraphs

# Fetch the dataset
# Uncompress the vertex metadata and convert to TSV
# Write the vertex metadata to vertex_data.txt
# Initialize the graph file, Graph.txt, with a header
include("build_dataset.jl")

# Remove duplicate vertex entries
;awk '!seen[$1]++' vertex_data.txt > vdata.txt

# Remove duplicate edge entries and convert to TSV
;awk '!seen[$0]++' gplus_combined.txt | tr ' ' '\t' > edata.txt

# Copy the metadata and edge data into the graph file
;cat vdata.txt edata.txt >> Graph.txt

# The graph dataset is now stored in Graph.txt
countlines("Graph.txt")

g = loadgraph("Graph.txt"; verbose=true)

size(g)

# Function to fetch the 5 most frequent entries
top5(x) = sort(collect(countmap(vcat(filter(y->length(y) > 0, collect(x))...))), by=x->x[2], rev=true)[1 : 5]

# Find the universities where alumni are well connected
@query(g |> filter(s.university == t.university) |> eachedge(s.university)) |> top5

# If you work for Google, which schools did people in your network go to?
network = hopgraph(g, @query(g |> filter("Google" in v.institution) |> eachvertex(v.label)), 1)
@query(network |> eachvertex(v.university)) |> top5

# Find the most popular schools in Los Angeles
@query(g |> filter("Los Angeles" in v.place) |> eachvertex(v.university)) |> top5

# Which cities are the most connect to New York?
@query(g |> filter("New York" in s.place) |> eachedge(t.place)) |> top5

# Run page rank, using LightGraphs, and set the result as a vertex property
M = export_adjacency(g)
setvprop!(g, :, LightGraphs.pagerank(LightGraphs.DiGraph(M)), :pagerank)

# Show a table containing a subset of vertex properties
VertexDescriptor(@query(g |> select(v.gender, v.last_name, v.pagerank)))

# Find the number of mutual friends between the source and target vertices for each edge
seteprop!(g, :, @query(g |> eachedge(e.mutualcount)), :mutual_friends);
EdgeDescriptor(g)


text:
using ControlSystems

abstract type AbstractController end

struct TVLQRController <: AbstractController
    nom_x
    nom_u
    K
    xf
    uf
    L
    N
    Δt
    LQR_eps
end

"""
Solves the discrete-time algebraic Riccati equation (DARE)

$(SIGNATURES)
"""
function TVDARE(model, nom_x, nom_u, Q, R, Qf, dt)
    P = Qf
    n = model.n
    m = model.m
    N = length(nom_x)
    K = [zeros(m,n) for i = 1:N-1]

    model_d = discretize_model(model, :rk4)

    J = TrajectoryOptimization.PartedMatrix(model_d)

    for i = 1:length(nom_u)
        t = N - i
        xt = nom_x[t]
        ut = nom_u[t]

        jacobian!(J, model_d, xt, ut, dt)

        A = J.xx
        B = J.xu

        APA = A' * P * A
        APB = A' * P * B
        BPA = B' * P * A
        BPB = B' * P * B

        Kt = -(R + BPB)\BPA
        K[t] .= Kt

        P = Q + Kt' * R * Kt + (A + B * Kt)' * P * (A + B * Kt)
    end

    return K
end

struct ALTRO
    x
    u
    K
end

"""
Performs direct-collocation using Augmented Lagrangian method

$(SIGNATURES)

**Arguments**:

- `model :: TrajectoryOptimization.Model`
- `x0 :: Vector{Float64}`
- `xf :: Vector{Float64}`
- `Q :: Matrix{Float64}`
- `R :: Matrix{Float64}`
- `Qf :: Matrix{Float64}`
- `N :: Int64`
- `dt :: Float64`

**Returns**:

`struct` containing the nominal tracking trajectory and the gains
"""
function ALTRO(model, x0, xf, Q, R, Qf, N, dt;
               ilqr_opts=(iterations=200,
                          iterations_linesearch=10,
                          verbose=true, cost_tolerance=1.0e-6),
               al_opts=(verbose=false,
                        iterations=30,
                        penalty_scaling=10.0,
                        cost_tolerance=1.0e-6,
                        cost_tolerance_intermediate=1.0e-5,
                        constraint_tolerance=1.0e-4),
               )
    model_d = discretize_model(model, :rk4)
    T = eltype(x0)
    n = model_d.n;
    m = model_d.m;

    goal = goal_constraint(xf);

    #N = 251;

    U = [ones(m) for k = 1:N-1];

    #dt = 0.01;

    obj = LQRObjective(Q,R,Qf,xf,N);
    con = Constraints(N);
    con[N] += goal;

    prob = Problem(model_d,obj,constraints=con,x0=x0,N=N,dt=dt,xf=xf);
    initial_controls!(prob, U);

    opts_ilqr = iLQRSolverOptions(;ilqr_opts...);

    opts_al = AugmentedLagrangianSolverOptions{T}(;opts_uncon=opts_ilqr, al_opts...);

    opts_altro = ALTROSolverOptions{T}(verbose=false,
                                       opts_al=opts_al);

    prob, solver = TrajectoryOptimization.solve(prob, opts_altro); # solve with ALTRO
    prob.X[end] = xf

    return ALTRO(prob.X, prob.U, solver.solver_al.solver_uncon.K)
end

"""
Generates the nominal trajectory to track and a TVLQR follwer for the same trajectory

$(SIGNATURES)
"""
function TVLQRController(model, x0, xf, Q, R, Qf, N, dt;
                         TVLQR_Q = nothing,
                         TVLQR_R = nothing,
                         TVLQR_Qf = nothing,
                         LQR_eps = 1e-1,
                         ilqr_opts=(iterations=200,
                                    iterations_linesearch=10,
                                    verbose=true, cost_tolerance=1.0e-6),
                         al_opts=(verbose=false,
                                  iterations=30,
                                  penalty_scaling=10.0,
                                  cost_tolerance=1.0e-6,
                                  cost_tolerance_intermediate=1.0e-5,
                                  constraint_tolerance=1.0e-4),
                         )

    if isnothing(TVLQR_Q)
        TVLQR_Q = Q
    end
    if isnothing(TVLQR_R)
        TVLQR_R = R
    end
    if isnothing(TVLQR_Qf)
        TVLQR_Qf = Qf
    end

    nom_traj = ALTRO(model, x0, xf, Q, R, Qf, N, dt, ilqr_opts=ilqr_opts, al_opts=al_opts)
    return TVLQRController(model, x0, xf, nom_traj, TVLQR_Q, TVLQR_R, TVLQR_Qf, N, dt, LQR_eps)
end

"""
Generates a TVLQR trajectory follower given the nominal trajectory `nom_traj`

$(SIGNATURES)

`Q`, `R`, `Qf` refer to the TVLQR costs
"""
function TVLQRController(model, x0, xf, nom_traj, Q, R, Qf, N, dt, LQR_eps)
    model_d = discretize_model(model, :rk4)

    X = nom_traj.x
    U = nom_traj.u
    K = TVDARE(model, nom_traj.x, nom_traj.u, Q, R, Qf, dt)
    # K = nom_traj.K

    J = TrajectoryOptimization.PartedMatrix(model_d)
    uf = zeros(model.m)

    jacobian!(J, model_d, xf, uf, dt)

    A = J.xx
    B = J.xu
    L = dlqr(A, B, Q, R)
    return TVLQRController(X, U, K, xf, uf, -L, N, dt, LQR_eps)
end

function (ctrl::TVLQRController)(t, x, thetamask)
    i = convert(Int, ceil(t/ctrl.Δt))
    i_ = min(i, ctrl.N-1)
    # i = min(t, ctrl.N-1)

    dx = x - ctrl.xf
    @views dx .= wraptopi(dx, thetamask)
    #@show norm(dx)

    if (i_ < i || norm(dx) < ctrl.LQR_eps)
        x_ = ctrl.xf
        u_ = ctrl.uf
        K = ctrl.L
        #@show dx
    else
        x_ = ctrl.nom_x[i_]
        u_ = ctrl.nom_u[i_]
        K  = ctrl.K[i_]
    end

    dx = x - x_
    @views dx .= wraptopi(dx, thetamask)

    du = K * dx
    u = u_ + du

    return u
end

struct RandomController <: AbstractController
end

"""
Generates a controller with random actions

$(SIGNATURES)
"""
function RandomController(model, std)
    return (t, x, thetamask) -> std * randn(model.m)
end

function addnoise(ctrl::AbstractController, noisestd)
    function noisy(t, x, thetamask)
        u = ctrl(t, x, thetamask)
        return u .+ randn(size(u)...) .* noisestd
    end
    return noisy
end

"""
Clip the output of a controller to a maximum force value `maxu`.

$(SIGNATURES)
"""
function clip(ctrl::AbstractController, maxu)
    function cu(t, x, thetamask)
        u = ctrl(t, x, thetamask)
        return clamp.(u, -maxu, maxu)
    end
    return cu
end

function wraptopi(radians)
     return radians .- 2*pi*floor.((radians .+ pi)/(2*pi))
end

function wraptopi(radians, mask)
    oldrad = copy(radians)
    return mask .* wraptopi(radians) .+ (1.0 .- mask) .* oldrad
end


text:
module Model
export make_model_data

#using Plots

function calculate_data(scale::Float64, r::Float64)

	q = 10 .^ (-1:0.01:1) # q-values [nm^-1]

	#scale = 0.02 # Number of spheres per nm^3 [nm^-3]
	#r = 0.5 # nm
	rho = 1 # nm^-2
	background = 0 # nm^-1

	qr = q .* r
	intensity = scale .* ((3*(4*pi*r^3/3)*rho).*(sin.(qr) - qr .* cos.(qr)) ./ qr.^3).^2 .+ background # nm^-1

	phi = scale * (4*pi*r^3/3)

	psi = 3*phi / (1 - phi)
	b = psi .* ((cos.(qr) + qr .* sin.(qr)) .* (sin.(qr) - qr .* cos.(qr))) ./ qr.^3
	c = psi .* (sin.(qr) - qr .* cos.(qr)).^2 ./ qr.^3
	d = 1 .+ psi .* (qr.^2 .* sin.(qr) .* cos.(qr)) ./ qr.^3
	e = psi .* (qr.^2 .* (sin.(qr).^2)) ./ qr.^3
	f = psi .* (qr .* sin.(qr) .* (sin.(qr) - qr .* cos.(qr))) ./ qr.^3
	g = - psi .* (qr .* cos.(qr) .* (sin.(qr) - qr .* cos.(qr))) ./ qr.^3

	X = 1 .+ b + (2 .* e .* f .* g + d .* (f.^2 - g.^2)) ./ (d.^2 + e.^2)
	Y = c + (2 .* d .* f .* g - e .* (f.^2 - g.^2)) ./ (d.^2 + e.^2)

	S = (Y ./ c) ./ (X.^2 + Y.^2)

	#gr()
	#global plt = plot(q, intensity; xaxis = :log, yaxis = :log)

	intensity = S .* intensity

	#plot!(q, intensity)

	return q::Array{Float64, 1}, intensity::Array{Float64, 1}

end

function save_data_to_file(data::Tuple{Array{Float64, 1}, Array{Float64, 1}}, outputFileName::String)

	io = open(outputFileName, "w")
	if isempty(read(outputFileName, String))
		println(io, "q-Value [nm^-1], Intensity [nm^-1]")
	end

	for i in 1:length(data[1])

		local qValue = data[1][i]
		local intensity = data[2][i]

		println(io, qValue, ", ", intensity)

	end

	close(io)

	return nothing

end

function make_model_data(scale::Float64, r::Float64)

	save_data_to_file(calculate_data(scale::Float64, r::Float64), "calculated.data")

	return nothing

end

end # Model

max_stars_repo_name: JonPallbo/Scatt

max_issues_repo_name: JonPallbo/Scatt

max_forks_repo_name: JonPallbo/Scatt


text:
@testset "Interactive mode" begin
    @testset "Template creation" begin
        write(stdin.buffer, "$me\n\n\r\n\n\n\n\nd")
        t = interactive_template()
        @test t.user == me
        @test t.host == "github.com"
        @test isempty(t.license)
        @test t.authors == LibGit2.getconfig("user.name", "")
        @test t.dir == default_dir
        @test t.julia_version == VERSION
        @test !t.ssh
        @test t.dev
        @test !t.manifest
        @test isempty(t.plugins)

        if isempty(LibGit2.getconfig("github.user", ""))
            write(stdin.buffer, "\n")
            @test_throws ArgumentError t = interactive_template()
        end

        down = '\x1b' * "[B"  # Down array key.
        write(stdin.buffer, "$me\ngitlab.com\n$down\r$me\n$test_file\n0.5\nyes\nno\nyes\n$down\r$down\rd\n\n")
        t = interactive_template()
        @test t.user == me
        @test t.host == "gitlab.com"
        # Not sure if the order the licenses are displayed in is consistent.
        @test !isempty(t.license)
        @test t.authors == me
        @test t.dir == abspath(test_file)
        @test t.julia_version == v"0.5.0"
        @test t.ssh
        @test !t.dev
        @test t.manifest
        # Like above, not sure which plugins this will generate.
        @test length(t.plugins) == 2

        write(stdin.buffer, "$me\nd")
        t = interactive_template(; fast=true)
        @test t.user == me
        @test t.host == "github.com"
        @test t.license == "MIT"
        @test t.authors == LibGit2.getconfig("user.name", "")
        @test t.dir == default_dir
        @test t.julia_version == VERSION
        @test !t.ssh
        @test !t.manifest
        @test isempty(t.plugins)
        println()

        # Host and SSH aren't prompted for when git is disabled.
        write(stdin.buffer, "$me\n\n\r\n\n\n\nd")
        t = interactive_template(; git=false)
        @test t.host == "github.com"
        @test !t.ssh
        println()
    end

    @testset "Package generation" begin
        write(stdin.buffer, "$me\n\n\r\n\n\n\n\n\n\n\nd")
        generate_interactive(test_pkg; gitconfig=gitconfig)
        @test isdir(joinpath(default_dir, test_pkg))
        rm(joinpath(default_dir, test_pkg); force=true, recursive=true)
    end

    @testset "Plugins" begin
        include("plugins.jl")
    end
end


text:
import Base: *, \, unsafe_getindex
using LinearAlgebra.BLAS: syr!, ger!, syrk!, syr2k!
using LinearAlgebra: BlasInt, BlasFloat, Cholesky, CholeskyPivoted
import LinearAlgebra: cholesky, cholesky!, ldlt, ldlt!

cholesky(::Type{Positive{T}}, A::AbstractMatrix, pivot=Val{false}; tol=default_tol(A), blocksize=default_blocksize(T)) where {T} = ldlt(Positive{T}, A, pivot; tol=tol, blocksize=blocksize)[1]
cholesky(::Type{Positive}, A::AbstractMatrix, pivot=Val{false}; tol=default_tol(A), blocksize=default_blocksize(floattype(eltype(A)))) = cholesky(Positive{floattype(eltype(A))}, A, pivot; tol=tol, blocksize=blocksize)

function ldlt(::Type{Positive{T}}, A::AbstractMatrix, pivot=Val{false}; tol=default_tol(A), blocksize=default_blocksize(T)) where {T}
    size(A, 1) == size(A, 2) || throw(DimensionMismatch("A must be square"))
    A0 = Array{floattype(T)}(undef, size(A))
    copyto!(A0, A)
    ldlt!(Positive{T}, A0, pivot; tol=tol, blocksize=blocksize)
end
ldlt(::Type{Positive}, A::AbstractMatrix, pivot=Val{false}; tol=default_tol(A), blocksize=default_blocksize(floattype(eltype(A)))) = ldlt(Positive{floattype(eltype(A))}, A, pivot; tol=tol, blocksize=blocksize)

cholesky!(::Type{Positive{T}}, A::AbstractMatrix{T}, pivot=Val{false}; tol=default_tol(A), blocksize=default_blocksize(T)) where {T<:AbstractFloat} = ldlt!(Positive{T}, A, pivot; tol=tol, blocksize=blocksize)[1]
cholesky!(::Type{Positive}, A::AbstractMatrix{T}, pivot=Val{false}; tol=default_tol(A), blocksize=default_blocksize(T)) where {T<:AbstractFloat} = cholesky!(Positive{T}, A; tol=tol, blocksize=blocksize)

# Blocked, cache-friendly algorithm (unpivoted)
function ldlt!(::Type{Positive{T}}, A::AbstractMatrix{T}, pivot::Type{Val{false}}=Val{false}; tol=default_tol(A), blocksize=default_blocksize(T)) where {T}
    size(A,1) == size(A,2) || error("A must be square")
    eltype(A)<:Real || error("element type $(eltype(A)) not yet supported")
    K = size(A, 1)
    d = Array{signtype(T)}(undef, K)
    for j = 1:blocksize:K
        # Split A into
        #            |
        #       B11  |
        #            |
        # A = ----------------
        #            |
        #       B21  |   B22
        #            |
        jend = min(K, j+blocksize-1)
        B11 = view(A, j:jend, j:jend)
        d1 = view(d, j:jend)
        solve_diagonal!(B11, d1, tol)
        if jend < K
            B21 = view(A, jend+1:K, j:jend)
            solve_columns!(B21, d1, B11)
            B22 = view(A, jend+1:K, jend+1:K)
            update_columns!(B22, d1, B21)
        end
    end
    @static if VERSION >= v"0.7.0-DEV.393"
        return Cholesky(A, :L, BLAS.BlasInt(0)), d
    else
        return Cholesky(A, :L), d
    end
end

# Version with pivoting
function ldlt!(::Type{Positive{T}}, A::AbstractMatrix{T}, pivot::Type{Val{true}}; tol=default_tol(A), blocksize=default_blocksize(T)) where {T}
    size(A,1) == size(A,2) || error("A must be square")
    eltype(A)<:Real || error("element type $(eltype(A)) not yet supported")
    K = size(A, 1)
    d = Array{signtype(T)}(undef, K)
    piv = convert(Vector{BlasInt}, 1:K)
    Ad = diag(A)
    for j = 1:blocksize:K
        jend = min(K, j+blocksize-1)
        solve_columns_pivot!(A, d, piv, Ad, tol, j:jend)
        if jend < K
            d1 = view(d, j:jend)
            B21 = view(A, jend+1:K, j:jend)
            B22 = view(A, jend+1:K, jend+1:K)
            update_columns!(B22, d1, B21)
        end
    end
    CholeskyPivoted(A, 'L', piv, BLAS.BlasInt(K), tol, BLAS.BlasInt(0)), d
end

ldlt!(::Type{Positive}, A::AbstractMatrix{T}, pivot=Val{false}; tol=default_tol(A), blocksize=default_blocksize(T)) where {T} = ldlt!(Positive{T}, A; tol=tol, blocksize=blocksize)

function solve_diagonal!(B, d, tol)
    K = size(B, 1)
    for j = 1:K
        Bjj = B[j,j]
        if abs(Bjj) > tol
            # compute ℓ (as the jth column of B)
            d[j] = sign(Bjj)
            s = sqrt(abs(Bjj))
            B[j,j] = s
            f = d[j]/s
            for i = j+1:K
                B[i,j] *= f
            end
            # subtract ℓ[j+1:end]⊗ℓ[j+1:end] from the lower right quadrant
            update_columns!(view(B, j+1:K, j+1:K), d[j], view(B, j+1:K, j))
        else
            # For the zero diagonals, replace them with 1. In a Newton step,
            # this corresponds to following the gradient (i.e., H = eye).
            d[j] = 0
            B[j,j] = 1
            for i = j+1:K
                B[i,j] = 0
            end
        end
    end
    B
end

function solve_columns!(B21, d, B11)
    I, J = size(B21)
    for j = 1:J
        dj = d[j]
        dj == 0 && continue
        s = B11[j,j]
        f = dj/s
        for i = 1:I
            B21[i,j] *= f
        end
        update_columns!(view(B21, :, j+1:J), dj, view(B21, :, j), view(B11, j+1:J, j))
    end
    B21
end

# Here, pivoting applies to the whole matrix, so we don't pass in a view.
# The jrange input describes the columns we're supposed to handle now.
function solve_columns_pivot!(A, d, piv, Ad, tol, jrange)
    K, KA = last(jrange), size(A, 1)
    jmin = first(jrange)
    for j in jrange
        # Find the remaining diagonal with largest magnitude
        Amax = zero(eltype(A))
        jmax = j-1
        for jj = j:KA
            tmp = abs(Ad[jj])
            if tmp > Amax
                Amax = tmp
                jmax = jj
            end
        end
        if jmax > j
            pivot!(A, j, jmax)
            Ad[j], Ad[jmax] = Ad[jmax], Ad[j]
            piv[j], piv[jmax] = piv[jmax], piv[j]
        end
        Ajj = A[j,j]
        for k = jmin:j-1
            tmp = A[j,k]
            Ajj -= d[k]*tmp*tmp
        end
        if abs(Ajj) > tol
            # compute ℓ (as the jth column of A)
            d[j] = sign(Ajj)
            s = sqrt(abs(Ajj))
            A[j,j] = s
            f = d[j]/s
        else
            d[j] = 0
            A[j,j] = 1
            f = zero(eltype(A))
        end
        for k = jmin:j-1
            @inbounds ck = d[k]*A[j,k]
            @simd for i = j+1:KA
                @inbounds A[i,j] -= ck*A[i,k]
            end
        end
        dj = d[j]
        @simd for i = j+1:KA
            @inbounds tmp = A[i,j]
            tmp *= f
            @inbounds A[i,j] = tmp
            @inbounds Ad[i] -= dj*tmp*tmp
        end
    end
    A
end

# Computes dest -= d*c*c', in the lower diagonal
@inline function update_columns!(dest::StridedMatrix{T}, d::Number, c::StridedVector{T}) where {T<:BlasFloat}
    syr!('L', convert(T, -d), c, dest)
    dest
end

# Computes dest -= d*x*y'
@inline function update_columns!(dest::StridedMatrix{T}, d::Number, x::StridedVector{T}, y::StridedVector{T}) where {T<:BlasFloat}
    ger!(convert(T, -d), x, y, dest)
    dest
end

# Computes dest -= C*diagm(d)*C', in the lower diagonal
function update_columns!(dest::StridedMatrix{T}, d::AbstractVector, C::StridedMatrix{T}) where {T<:BlasFloat}
    isempty(d) && return dest
    # If d is homogeneous, we can use syr rather than syr2
    allsame = true
    d1 = d[1]
    for i = 2:length(d)
        allsame &= (d[i] == d1)
    end
    allsame && d1 == 0 && return dest
    if allsame
        syrk!('L', 'N', convert(T, -d1), C, one(T), dest)
    else
        Cd = C*Diagonal(d)
        syr2k!('L', 'N', -one(T)/2, C, Cd, one(T), dest)
    end
    dest
end

# Pure-julia fallbacks for the above routines
# Computes dest -= d*c*c', in the lower diagonal
function update_columns!(dest, d::Number, c::AbstractVector)
    K = length(c)
    for j = 1:K
        dcj = d*c[j]
        @simd for i = j:K
            @inbounds dest[i,j] -= dcj*c[i]
        end
    end
    dest
end

# Computes dest -= d*x*y'
function update_columns!(dest, d::Number, x::AbstractVector, y::AbstractVector)
    I, J = size(dest)
    for j = 1:J
        dyj = d*y[j]
        @simd for i = 1:I
            @inbounds dest[i,j] -= dyj*x[i]
        end
    end
    dest
end

# Computes dest -= C*diagm(d)*C', in the lower diagonal
function update_columns!(dest, d::AbstractVector, C::AbstractMatrix)
    Ct = C'
    Cdt = Diagonal(d) * Ct
    K = size(dest, 1)
    nc = size(C, 2)
    for j = 1:K
        for i = j:K
            tmp = zero(eltype(dest))
            @simd for k = 1:nc
                @inbounds tmp += Ct[k,i]*Cdt[k,j]
            end
            @inbounds dest[i,j] -= tmp
        end
    end
    dest
end

# Diagonal pivoting (row&column swap) for a lower triangular matrix
function pivot!(A, i::Integer, j::Integer)
    i, j = min(i,j), max(i,j)
    for k = 1:i-1
        A[i,k], A[j,k] = A[j,k], A[i,k]
    end
    A[i,i], A[j,j] = A[j,j], A[i,i]
    for k = i+1:j-1
        A[k,i], A[j,k] = A[j,k], A[k,i]
    end
    for k = j+1:size(A,1)
        A[k,i], A[k,j] = A[k,j], A[k,i]
    end
    A
end

floattype(::Type{T}) where {T<:AbstractFloat} = T
floattype(::Type{T}) where {T<:Integer} = Float64
floattype(::Type{T}) where {T<:Real} = T

signtype(::Type{T}) where {T<:AbstractFloat} = Int8
signtype(::Type{T}) where {T<:Integer} = Int8
signtype(::Type{T}) where {T<:Real} = T

const cachesize = 2^15

default_δ(A) = 10 * size(A, 1) * eps(floattype(real(eltype(A))))
default_tol(A) = default_δ(A) * maximum(abs,A)
default_blocksize(::Type{T}) where {T} = max(4, floor(Int, sqrt(cachesize/sizeof(T)/4)))


text:
function gradient!{T <: AbstractFloat}(u::AbstractVector{T}, dist::Divergence, a::AbstractVector{T}, b::AbstractVector{T})
    @inbounds for i = eachindex(a, b)
        u[i] = gradient(dist, a[i], b[i])
    end
    u
end

function gradient!{T <: AbstractFloat}(u::AbstractVector{T}, dist::Divergence, a::AbstractVector{T})
    @inbounds for i = eachindex(a)
        u[i] = gradient(dist, a[i])
    end
    u
end

function gradient{T <: AbstractFloat}(dist::Divergence, a::AbstractVector{T}, b::AbstractVector{T})
    if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    gradient!(Array(T, length(a)), dist, a, b)
end

function gradient{T <: AbstractFloat}(dist::Divergence, a::AbstractVector{T})
    gradient!(Array(T, length(a)), dist, a)
end

function hessian!{T <: AbstractFloat}(u::Vector{T}, dist::Divergence, a::AbstractVector{T}, b::AbstractVector{T})
    @inbounds for i = eachindex(a, b)
        u[i] = hessian(dist, a[i], b[i])
    end
    u
end

function hessian!{T <: AbstractFloat}(u::Vector{T}, dist::Divergence, a::AbstractVector{T})
    @inbounds for i = eachindex(a)
        u[i] = hessian(dist, a[i])
    end
    u
end

function hessian{T <: AbstractFloat}(dist::Divergence, a::AbstractVector{T}, b::AbstractVector{T})
    if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    hessian!(Array(T, length(a)), dist, a, b)
end

function hessian{T <: AbstractFloat}(dist::Divergence, a::AbstractVector{T})
    hessian!(Array(T, length(a)), dist, a)
end


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
using ArgParse
s = ArgParseSettings()
@add_arg_table s begin
    "--gpu"
        help = "Use GPU (cuda)"
        action = :store_true
    "--llvm"
        help = "Use LLVM to gernerate GPU code"
        action = :store_true
    "--profile"
        help = "Use profiler"
        action = :store_true
    "--minibatch"
        help = "Minibatch size"
        arg_type = Int
        default = 100
    "--epochs"
        help = "Epochs to train"
        arg_type = Int
        default = 10
    "--dropout"
        help = "Use dropout"
        action = :store_true
    "--once"
        help = "Don't compile the evaluator"
        action = :store_true
    "--devectorize"
        help = "Devectorize the graph before evaluation"
        action = :store_true
    "--test"
        help = "Evaluate test set"
        action = :store_true
    "--seed"
        help = "Random seed"
        arg_type = Int
    "--graphviz"
        help = "Dump the network structure to mnist.gv"
        action = :store_true
    "--load"
        help = "Load a snapshot before training"
        #arg_type = Int
    "--store"
        help = "Store a snapshot after training"
        #arg_type = Int
    "--64"
        help = "Use 64 bits"
        action = :store_true
    "--solver"
        help = "Select solver (sgd, momentum, adadelta, adam)"
        default = "momentum"
end

const args = parse_args(ARGS, s)
const gpu = args["gpu"]
const llvm = args["llvm"]
const profile = args["profile"]
const minibatch_size = args["minibatch"]
const epochs = args["epochs"]
const dropout = args["dropout"]
const once = args["once"]
const solver = args["solver"]

if !gpu
end

const T = args["64"] ? Float64 : Float32

if gpu
    using CUDArt
    CUDArt.init(0)
    device(0)

    using CUBLAS, MNIST, Seep, JLD

    make_node(n, a::Array) = ANode(n, CudaArray(map(T, a)))
    if once
      make_node(n, x::Int...) = ANode(n, CudaArray(T, x...))
    else
      make_node(n, x::Int...) = ANode(n, x...)
    end
    get_data(x) = to_host(x)
    const pool = BuddyPool(CudaArray(T, 1<<23))
else
    ENV["SEEP_NO_GPU"]=1
    using MNIST, Seep, JLD

    make_node(n, a::Array) = ANode(n, map(T, a))
    if once
      make_node(n, x::Int...) = ANode(n, Array(T, x...))
    else
      make_node(n, x::Int...) = ANode(n, x...)
    end
    get_data(x) = x
    const pool = BuddyPool(Array(T, 1<<23))
end

if args["seed"] != nothing
  srand(args["seed"])
end

const ftrain = map(T, traindata()[1]/255)
const ltrain = map(T, traindata()[2])
const ytrain = zeros(T, 10, length(ltrain))
for i in 1:size(ytrain, 2)
  ytrain[1+round(Int, ltrain[i]),i] = 1
end

const n0 = 28*28
const x = ANode("x", n0, minibatch_size)
a = x

reg_terms = ANode[]
params = ANode[]

if true
  const n1 = 25
  const w1 = make_node("w1", randn(n1, n0)/sqrt(n0))
  const b1 = make_node("b1", zeros(n1, 1))
  a = tanh(w1*a.+b1)

  if dropout
    const d = make_node("d1", zeros(n1, minibatch_size))
    a = a.*d
  end
  push!(params, w1)
  push!(params, b1)
  push!(reg_terms, w1)
else
  const n1 = n0
end

if false
  const w2 = make_node("w2", randn(n2, n1)/sqrt(n1))
  const b2 = make_node("b2", zeros(n2, 1))
  a = tanh(w2*a.+b2)
  push!(params, w2)
  push!(params, b2)
  push!(reg_terms, w2)
else
  const n2 = n1
end

const n3 = 10
const w3 = make_node("w3", randn(n3, n2)/sqrt(n2))
const b3 = make_node("b3", zeros(n3, 1))
a = softmax(w3*a.+b3, 1)
push!(params, w3)
push!(params, b3)
push!(reg_terms, w3)

const y = ANode("y", 10, minibatch_size)
const e = -log(a).*y

const obj= sum(e) + 1e-3*reduce(+, map(w->dot(w,w), reg_terms))

const grad = gradients(obj)

if solver == "sgd"
  const sol = gradient_descent
elseif solver == "momentum"
  const sol = momentum
elseif solver == "adadelta"
  const sol = Seep.adadelta
elseif solver == "adam"
  const sol = adam
elseif solver == "adamfast"
  const sol = Seep.adam_fast
else
  error("Unknown solver: $solver")
end

upd = map(x->sol(x, grad[x]/sqrt(minibatch_size)), params)

const E = Float64[]
const times = Float64[]

if args["graphviz"]
  open(f->graphviz(f,  e, upd...), "mnist.gv", "w")
end

if args["devectorize"]
  Seep.devectorize(e, upd...)
end

if once
  run() = global dict = evaluate(pool, e, a, upd...)
else
  const run,dict = (gpu && llvm ? Seep.llvm_instance : instance)(pool, ANode[e; a; upd])
end

Seep.print_stats(STDOUT, pool)

if args["load"] !== nothing
  load_snapshot(args["load"], e)
end

n = 0

if once
  const xa,ya = x.arg,y.arg
else
  const xa,ya = dict[x],dict[y]
end

if dropout && once
  const da = d.arg
elseif dropout
  const da = dict[d]
end

function minibatch(ii)
    if gpu
       copy!(xa, ftrain[:,ii])
       copy!(ya, ytrain[:,ii])
       if dropout
         copy!(da, 2. * bitrand(n1,length(ii)))
       end
    else
        for i in 1:length(ii)
          for j in 1:(28*28)
            xa[j,i] = ftrain[j,ii[i]]
          end

          for j in 1:10
            ya[j,i] = ytrain[j,ii[i]]
          end

          if dropout
            for j in 1:n1
              da[j,i] = 2. * rand(0:1)
            end
          end
        end
    end

    run()

    push!(E, sum(get_data(dict[e]))/length(ii))
    push!(times, time())
end

function epoch()
  perm = randperm(60000)
  for i in 1:minibatch_size:60000
    minibatch(perm[i+(1:minibatch_size)-1])
  end
end

epoch()

Profile.clear_malloc_data()

if profile
    Profile.init(delay=0.01)
    for i in 1:epochs
        @profile epoch()
    end

    data,lidict = Profile.retrieve()
    @save "mnist-profile.jld" data lidict
else
    for i in 1:epochs
        @time epoch()
    end
end

if args["test"]
  ftest = map(T, testdata()[1]/255)
  ltest = map(T, testdata()[2])

  pp = zeros(10, 10000)
  if dropout
    da[:] = ones(n1, minibatch_size)
  end

  for i in 1:div(10000,minibatch_size)
    n += 1
    ix = (i-1)*minibatch_size+(1:minibatch_size)
    copy!(xa, vec(ftest[:,ix]))
    run()
    pp[:,ix] = get_data(dict[a])
  end

  correct = countnz([indmax(pp[:,i]) == (1+ltest[i]) for i in 1:10000])/10000
  @show correct

  gmean = exp(mean([log(pp[round(Int, 1+ltest[i]), i]) for i in 1:10000]))
  @show gmean
end

if args["store"] !== nothing
  save_snapshot(args["store"], e)
end

close(pool)

if gpu
  CUDArt.close(0)
end

@show E[end]


text:
using MarketData

facts("Volume") do

    context("obv") do
        @fact obv(ohlcv).values[1]     --> 4783900 # TTR value is 4783900
        @fact obv(ohlcv).values[12]    --> 9390200 # TTR value is 2000-01-19   9390200
        @fact obv(ohlcv).timestamp[12] --> Date(2000,1,19)
    end

    context("vwap") do
       @fact vwap(ohlcv).values[1]    --> roughly(97.9215, atol=.01)  # TTR value 2000-01-14 97.92154
       @fact vwap(ohlcv).values[end]  --> roughly(21.4446, atol=.01)  # TTR value  21.44458
       @fact vwap(ohlcv).timestamp[1] --> Date(2000,1,14)
    end
end


text:
using Test

include("EndmemberLibrary.jl")

datafile = "data/basic_endmember_library.csv"
classname = "Class"

@info "Basic Loading"
lib = load_data!(SpectralLibrary(datafile, classname, 18))
@test isnothing(lib) == false

@time @test isnothing(lib.spectra) == false
@time @test length(size(lib.spectra)) == 2

@info "Filter By Class"
tmplib = deepcopy(lib)
filter_by_class!(tmplib)
@time @test length(lib.classes) == length(tmplib.classes)

tmplib =  load_data!(SpectralLibrary(datafile, classname, 18, [unique(lib.classes)[1]]))
filter_by_class!(tmplib)
@time @test length(lib.classes) != length(tmplib.classes)
@time @test length(unique(tmplib.classes)) == 1
tmplib=nothing

@info "Remove wavelength region inplace"

tmplib =  load_data!(SpectralLibrary(datafile, classname, 18))

remove_wavelength_region_inplace!(tmplib)
@time @test size(lib.spectra) != size(tmplib.spectra)

@info "Wavelength interpolation"
tmplib = deepcopy(lib)
new_wl = 400 .+ Vector{Float64}(1:200)*10.
interpolate_library_to_new_wavelengths!(tmplib, new_wl)
@time @test size(tmplib.spectra)[2] == length(new_wl)
@time @test size(lib.spectra) != size(tmplib.spectra)

@info "Scale Library"
tmplib = deepcopy(lib)
scale_library!(tmplib)
@time @test lib.spectra == tmplib.spectra

scale_library!(tmplib, 2)
@time @test lib.spectra != tmplib.spectra

tmplib = deepcopy(lib)
brightness_normalize!(tmplib)
@time @test lib.spectra != tmplib.spectra


text:
module GW
    include("../../examples/gridworld.jl")
end

@testset "gridworld" begin
    env = GW.GridWorld()
    reset!(env)
    while !terminated(env)
        act!(env, rand(actions(env)))
        render(env)
    end
    @test haskey(env.rewards, state(env))

    env2 = GW.GridWorld()
    @test state(clone(env2)) == state(env2)
    setstate!(env, state(env2))
    @test state(env) == state(env2)
    @test observe(env) == observe(env2)
    @test observations(env) == observations(env2)
end


text:
mutable struct QuadraticModel <: TransistorModel
    mos::MOSStructure
    Vth
    Cox
    W
    L
    μ
    T0
end

setW!(Wn,m::QuadraticModel) = begin m.W = Wn end
setL!(Ln,m::QuadraticModel) = begin m.L = Ln end
function setT!(T,m::QuadraticModel)
    Sdist = m.mos.semiconductor.dist
    Ssemi = m.mos.semiconductor.material
    nMos = MOSStructure(m.mos.metal,m.mos.oxide,SemiconductorData(T,Sdist,Ssemi),m.mos.tox,m.mos.tsi)
    nModel = QuadraticModel(m.W,m.L,nMos;μ=m.μ)
    m.mos = nMos
    m.Vth = nModel.Vth
    m.Cox = nModel.Cox
    m.W = nModel.W
    m.L = nModel.L
    m.μ = nModel.μ
end
function Vth_Quadratic(Vfb,Na,tox,T,phif,sigma=0.0)
    phit = kb*T
    G = γ(Na,tox)
    k = 0#log(n/(n-1))
    #uprintln(n,",",phit,",",G,",",k,",",Vfb,",",Na,",",tox,",",T,",",phif)
    return Vfb+2*phif+phit*k+G*sqrt(2*phif)
end

function Vth_Quadratic(mos::MOSStructure)
    return Vth_Quadratic(VFB(mos),abs(C(mos)),mos.tox,Temperature(mos),ϕb(mos))
end

QuadraticModel(W,L,mos::MOSStructure;μ=520.0) = QuadraticModel(mos,Vth_Quadratic(mos),Cox(mos),W,L,μ,Temperature(mos))

function Id(Vg,Vd,Vs,Vb,mdl::QuadraticModel)
    Vgs = Vg-Vs
    Vds = Vd-Vs
    ϕt = Temperature(mdl.mos)*kb
    @unpack mos,Vth,Cox,W,L,μ,T0 = mdl
    Vov = Vgs-Vth
    S = W/L
    if Vov < 0.0
        return S*μ*(T0/300.0)^(-3/2.0)*Cox*ϕt^2*(exp(Vov/ϕt))*(1-exp(-Vds/ϕt))
    elseif Vds < Vov
        return S*μ*(T0/300.0)^(-3/2.0)*Cox*(Vov*Vds-Vds^2*0.5)
    else
        return S*μ*(T0/300.0)^(-3/2.0)*Cox*0.5*(Vov)^2
    end

end

function id(Vg,Vd,Vs,Vb,mdl::QuadraticModel)
    Vgs = Vg-Vs
    Vds = Vd-Vs
    ϕt = Temperature(mdl.mos)*kb
    @unpack mos,Vth,Cox,W,L,μ = mdl
    Vov = Vgs-Vth
    S = W/L
    if Vov < 0.0
        return ϕt^2*(exp(Vov/ϕt))*(1-exp(-Vds/ϕt))
    elseif Vds < Vov
        return (Vov*Vds-Vds^2*0.5)
    else
        return 0.5*(Vov)^2
    end

end

Id(Vg,Vd,Vs,m::QuadraticModel) = Id(Vg,Vd,Vs,Vs,m::QuadraticModel)
# Symbolics.@register Id(Vg,Vd,Vs,m::QuadraticModel)
Id(Vg,Vd,m::QuadraticModel) = Id(Vg,Vd,0.0,m::QuadraticModel)
gm(Vg,Vd,Vs,Vb,m::QuadraticModel) = Zygote.ForwardDiff.derivative(x->Id(x,Vd,Vs,Vb,m),Vg)
IdNR(vg,vd,vs,m::QuadraticModel) = Id(vg,vd,vs,m)-gm(vg,vd,vs,m)*(vg-vs)-gds(vg,vd,vs,m)*(vd-vs)
# Symbolics.@register IdNR(vg,vd,vs,m::QuadraticModel)
# Symbolics.@register gm(Vg,Vd,Vs,Vb,m::QuadraticModel)
gm(Vg,Vd,Vs,m::QuadraticModel) = Zygote.ForwardDiff.derivative(x->Id(x,Vd,Vs,m),Vg)
gm(Vg,Vd,m::QuadraticModel) = Zygote.ForwardDiff.derivative(x->Id(x,Vd,0.0,m),Vg)
gmId(Vg,Vd,m::QuadraticModel) = gm(Vg,Vd,m)/Id(Vd,Vg,m)
gds(Vg,Vd,Vs,m::QuadraticModel) = Zygote.ForwardDiff.derivative(x->Id(Vg,x,Vs,m),Vd)
# ModelingToolkit.@register gds(Vg,Vd,Vs,Vb,m::QuadraticModel)

gdsId(Vg,Vd,m::QuadraticModel;order=5) = gds(Vg,Vd,m)/Id(Vd,Vg,m)


text:
module TestSet

using Folds
using Test
include("utils.jl")

@testset begin
    @test Folds.unique(x for x in 1:100 if x < 3) ==ₜ [1, 2]
    @test Folds.unique(x for x in 1:100 if x < 3; basesize = 1) ==ₜ [1, 2]
end

end  # module


text: #using .DualNumbers


text:
"""
$(TYPEDEF)

Two sample superiority test for proportion

Constructors
------------
* `TwoSamplePropSuperior(p1::Real, p2::Real, k::Real, delta::Real)`

Arguments
---------
* `p1`: Proportion of group 1

* `p2`: Proportion of group 2

* `k`: Allocation ratio of the groups, k = n(group1) / n(group 2)

* `delta`: Non-inferiority Margin

Fields
------
$(FIELDS)
"""
type TwoSamplePropSuperior <: TrialTest
    p1::Real
    p2::Real
    k::Real
    delta::Real

    # Validator
    function TwoSamplePropSuperior(p1, p2, k, delta)

        if !((0 < p1 < 1) & (0 < p2 < 1))

            error("Proportion values must be in (0, 1)")

        end # end if

        if !(0 < k < Inf)

            error("Allocation ratio must be in (0, Inf)")

        end # end if

        if !(0 <= delta < 1)

            error("The superiority margin δ must be in [0, 1)")

        end # end if

        new(p1, p2, k, delta)

    end # end function

end # function

# Two sample superiority test for proportion
function hypotheses{T <: TwoSamplePropSuperior}(test::T, n::Real, std::Union{Real, Tuple{Real, Real}, Void}, alpha::Real, side::String)

    diff = test.p1 - test.p2 - test.delta
    se = sqrt(1 / (test.k * n) * test.p1 * (1 - test.p1) + 1 / n * test.p2 * (1 - test.p2))
    z = diff / se
    p = cdf(ZDIST, z - quantile(ZDIST, 1 - alpha)) + cdf(ZDIST, -z - quantile(ZDIST, 1 - alpha))
    return p

end # end function


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct DeploymentWhatIfSettings <: SwaggerModel
    resultFormat::Any # spec type: Union{ Nothing, String } # spec name: resultFormat

    function DeploymentWhatIfSettings(;resultFormat=nothing)
        o = new()
        validate_property(DeploymentWhatIfSettings, Symbol("resultFormat"), resultFormat)
        setfield!(o, Symbol("resultFormat"), resultFormat)
        o
    end
end # type DeploymentWhatIfSettings

const _property_map_DeploymentWhatIfSettings = Dict{Symbol,Symbol}(Symbol("resultFormat")=>Symbol("resultFormat"))
const _property_types_DeploymentWhatIfSettings = Dict{Symbol,String}(Symbol("resultFormat")=>"String")
Base.propertynames(::Type{ DeploymentWhatIfSettings }) = collect(keys(_property_map_DeploymentWhatIfSettings))
Swagger.property_type(::Type{ DeploymentWhatIfSettings }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DeploymentWhatIfSettings[name]))}
Swagger.field_name(::Type{ DeploymentWhatIfSettings }, property_name::Symbol) =  _property_map_DeploymentWhatIfSettings[property_name]

const _allowed_DeploymentWhatIfSettings_resultFormat = ["ResourceIdOnly", "FullResourcePayloads"]

function check_required(o::DeploymentWhatIfSettings)
    true
end

function validate_property(::Type{ DeploymentWhatIfSettings }, name::Symbol, val)
    if name === Symbol("resultFormat")
        Swagger.validate_param(name, "DeploymentWhatIfSettings", :enum, val, _allowed_DeploymentWhatIfSettings_resultFormat)
    end
end


text:
function Grid=DGNodes(Grid,Trans,DG)
for iE=1:Grid.NumEdges
  Grid.Edges(iE).DGPoints=Trans(DG.xw,P1,P2,Grid);
end
end


text:
using BinDeps
using Compat

@BinDeps.setup

world = library_dependency("libworld", aliases=["libworld", "world-0"])

const version = "0.2.1-4"

github_root = "https://github.com/r9y9/WORLD"
arch = WORD_SIZE == 64 ? "x86_64" : "i686"
major = version[1]
provides(Binaries,
         URI("$(github_root)/releases/download/v$(version)/world-$(major)_mingw$(WORD_SIZE)_$(arch).zip"),
         world, unpacked_dir = "usr/lib", os = :Windows)

provides(Sources,
         URI("$(github_root)/archive/v$(version).tar.gz"),
         world,
         unpacked_dir="WORLD-$(version)")

prefix = joinpath(BinDeps.depsdir(world), "usr")
srcdir = joinpath(BinDeps.depsdir(world), "src", "WORLD-$(version)")

provides(SimpleBuild,
          (@build_steps begin
              GetSources(world)
              @build_steps begin
                  ChangeDirectory(srcdir)
                  `./waf configure --prefix=$prefix`
                  `./waf build`
                  `./waf install`
              end
           end), world, os = :Unix)

@BinDeps.install @compat Dict(:libworld => :libworld)


text:
mutable struct Solve{V} <: State{V}
    value::V
end

Solve(; unit, _type, _...) = begin
    U = value(unit)
    V = valuetype(_type, U)
    v = zero(V)
    Solve{V}(v)
end

constructortags(::Val{:Solve}) = (:unit,)

# \Equal opeartor for both solve/bisect
⩵(x, y) = x - y
export ⩵

using SymPy: SymPy, sympy
genpolynomial(v::VarInfo; scope) = begin
    x = v.name
    V = extractfuncargpair.(v.args) .|> first
    push!(V, x)
    #HACK: use scope of caller, avoiding eval into precompiled Cropbox module
    p = scope.eval(@q let $(V...)
        Cropbox.SymPy.@vars $(V...)
        Cropbox.sympy.Poly(begin
            $(v.body)
        end, $x)
    end)
    Q = p.coeffs()
    #HACK: normalize coefficient to avoid runtime unit generation
    Q = Q / Q[1] |> reverse .|> SymPy.simplify
    Q .|> repr .|> Meta.parse
end
genpolynomialunits(U, n) = [@q($U^$(i-1)) for i in n:-1:1]

import PolynomialRoots
solvepolynomial(p, pu) = begin
    sp = [deunitfy(q, qu) for (q, qu) in zip(p, pu)]
    r = PolynomialRoots.roots(sp)
    real.(filter(isreal, r))
end
solvequadratic(a, b, c) = begin
    (a == 0) && return (-c/b,)
    Δ = b^2 - 4a*c
    if Δ > 0
        s = √Δ
        ((-b - s) / 2a, (-b + s) / 2a)
    elseif Δ == 0
        (-b/2a,)
    else
        #HACK: make up a non-complex solution
        x = -b/2a
        @warn "ignore complex roots for quadratic equation: $a*x^2 + $b*x + $c = 0" x
        (x,)
    end
end

updatetags!(d, ::Val{:Solve}; _...) = begin
    !haskey(d, :lower) && (d[:lower] = -Inf)
    !haskey(d, :upper) && (d[:upper] = Inf)
    !haskey(d, :pick) && (d[:pick] = QuoteNode(:maximum))
end

genvartype(v::VarInfo, ::Val{:Solve}; V, _...) = @q Solve{$V}

gendefault(v::VarInfo, ::Val{:Solve}) = nothing

genupdate(v::VarInfo, ::Val{:Solve}, ::MainStep; scope, kw...) = begin
    U = gettag(v, :unit)
    isnothing(U) && (U = @q(u"NoUnits"))
    P = genpolynomial(v; scope)
    n = length(P)
    PU = genpolynomialunits(U, n)
    lower = gettag(v, :lower)
    upper = gettag(v, :upper)
    pick = gettag(v, :pick).value
    body = if n == 2 # linear
        @gensym a b xl xu
        @q let $a = $(esc(P[2])),
               $b = $(esc(P[1])),
               $xl = $C.unitfy($C.value($lower), $U),
               $xu = $C.unitfy($C.value($upper), $U)
            clamp(-$b / $a, $xl, $xu)
        end
    elseif n == 3 # quadratic
        @gensym a b c X xl xu l
        @q let $a = $C.deunitfy($(esc(P[3])), $(PU[3])),
               $b = $C.deunitfy($(esc(P[2])), $(PU[2])),
               $c = $C.deunitfy($(esc(P[1])), $(PU[1])),
               $X = $C.unitfy($C.solvequadratic($a, $b, $c), $U),
               $xl = $C.unitfy($C.value($lower), $U),
               $xu = $C.unitfy($C.value($upper), $U)
            #TODO: remove duplication
            $l = filter(x -> $xl <= x <= $xu, collect($X))
            isempty($l) && ($l = clamp.($X, $xl, $xu))
            $l |> $pick
        end
    else # generic polynomials (slow!)
        @gensym X xl xu l
        @q let $X = $C.unitfy($C.solvepolynomial([$(esc.(P)...)], [$(PU...)]), $U),
               $xl = $C.unitfy($C.value($lower), $U),
               $xu = $C.unitfy($C.value($upper), $U)
            $l = filter(x -> $xl <= x <= $xu, $X)
            #TODO: better report error instead of silent clamp?
            isempty($l) && ($l = clamp.($X, $xl, $xu))
            $l |> $pick
        end
    end
    val = genbody(v, body)
    genstore(v, val)
end


text:
using JSON

const LERNA_CONFIG = joinpath(@__DIR__, "..", "packages", "lerna.json")
const VERSION = open(LERNA_CONFIG, "r") do io
    VersionNumber(JSON.parse(io)["version"])
end


text:
function get_ax2o_a(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getAx2oA", RealFieldElement, ())
end

function get_bo_a(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getBoA", RealFieldElement, ())
end

function get_bo_a_bpo(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getBoABpo", RealFieldElement, ())
end

function get_chi(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getChi", RealFieldElement, ())
end

function get_chi2(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getChi2", RealFieldElement, ())
end

function get_co2ab(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getCo2AB", RealFieldElement, ())
end

function get_e2(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getE2", RealFieldElement, ())
end

function get_max_ecc_pow(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getMaxEccPow", jint, ())
end

function get_mean_motion(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getMeanMotion", RealFieldElement, ())
end

function get_moa(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getMoa", RealFieldElement, ())
end

function get_oo_ab(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getOoAB", RealFieldElement, ())
end

function get_orbit_period(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getOrbitPeriod", RealFieldElement, ())
end

function get_ratio(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getRatio", RealFieldElement, ())
end

function get_res_orders(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getResOrders", List, ())
end

function get_roa(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getRoa", RealFieldElement, ())
end

function get_theta(obj::FieldDSSTTesseralContext)
    return jcall(obj, "getTheta", RealFieldElement, ())
end


text:
import Base: reducedim, sum, maximum
export argmax, argmin

function reshape3d(x::CuArray, dim::Int)
    # dim == 0 && return (1, length(x), 1)
    dim1, dim2, dim3 = 1, size(x,dim), 1
    for i = 1:dim-1
        dim1 *= size(x, i)
    end
    for i = dim+1:ndims(x)
        dim3 *= size(x, i)
    end
    reshape(x, dim1, dim2, dim3)
end

function reducedim_f(T, v0, op)
    CuFunction("""
    template<typename T>
    __inline__ __device__ void warpReduce(T &a) {
        for (int delta = warpSize/2; delta > 0; delta /= 2) {
            T b = __shfl_down(a, delta);
            a = $op;
        }
    }

    template<typename T>
    __inline__ __device__ void blockReduce(T &val) {
        static __shared__ T shared[32]; // 32*32=1024
        int laneId = threadIdx.x % warpSize;
        int warpId = threadIdx.x / warpSize;
        warpReduce<T>(val);
        if (laneId == 0) shared[warpId] = val;
        __syncthreads();

        val = shared[laneId];
        if (warpId == 0) warpReduce<T>(val);
    }

    __global__ void reduce($T *x, $T *y, int size1, int size2) {
        int idxX = blockIdx.x * blockDim.x + threadIdx.x;
        int idxY = blockIdx.y * blockDim.y + threadIdx.y;
        if (idxX >= size1) return;
        if (idxY >= size2) return;

        $T val = x[idxX + idxY * size2];
        blockReduce<$T>(val);
        if (idxX == 0) y[idxY] = val;
        //if (threadIdx.x == 0) y[blockIdx.x] = val;
    }
    """)
end

@generated function sum2(x::CuArray{T}, dim::Int) where T
    f = reducedim_f(cstring(T), "0", "a + b")
    quote
        #x3d = reshape3d(x, dim)
        #xsize = [size(x)...]
        #xsize[dim] = 1
        #y = similar(x, xsize...)

        y = similar(x, 1, size(x,2))
        #n = size(x,2)
        bx = 1024
        #gx = n <= bx ? 1 : ceil(Int, n/bx)
        gdims, bdims = (1,size(x,2),1), (bx,1,1)
        # gdims, bdims = cudims(length(x))
        culaunch($f, gdims, bdims, x.ptr, y.ptr, size(x,1), size(x,2))
        y
    end
end

function find_f{T}(::Type{T}, v0::String, op::String)
    CuFunction("""
    template<typename T>
    __inline__ __device__ void warpReduce(T &a, int &i) {
        for (int delta = warpSize/2; delta > 0; delta /= 2) {
            T b = __shfl_down(a, delta);
            int j = __shfl_down(i, delta);
            if ($op) {
                a = b;
                i = j;
            }
        }
    }
    template<typename T>
    __inline__ __device__ void blockReduce(T &value, int &index) {
        static __shared__ T values[32]; // 32*32=1024
        static __shared__ int indexs[32];
        int laneId = threadIdx.y % warpSize;
        int warpId = threadIdx.y / warpSize;
        warpReduce<T>(value, index);
        if (laneId == 0) {
            values[warpId] = value;
            indexs[warpId] = index;
        }
        __syncthreads();
        value = (laneId < 32) ? values[laneId] : $v0;
        index = (laneId < 32) ? indexs[laneId] : -1;
        warpReduce<T>(value, index);
    }
    __global__ void f(Array<$T,3> x, Array<int,3> y) {
        int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
        int idx_y = blockIdx.y * blockDim.y + threadIdx.y;
        int idx_z = blockIdx.z * blockDim.z + threadIdx.z;
        $T value = (idx_y < x.dims[1]) ? x(idx_x,idx_y,idx_z) : $v0;
        int index = (idx_y < x.dims[1]) ? blockIdx.y * blockDim.y + threadIdx.y : -1;
        //for (int i = idx_y+blockDim.y*gridDim.y; i < x.dims[1]; i += blockDim.y*gridDim.y) {
        //    $T b = x(idx_x, i, idx_z);
        //    a = $op;
        //}
        blockReduce<$T>(value, index);
        if (threadIdx.y == 0) y(blockIdx.x, blockIdx.y, blockIdx.z) = index + 1;
    }
    """)
end

function reducedim_f2(T, op::String)
    blocksize = 1024
    CuFunction("""
    __global__ void reduce($T *x, $T *y) {
        static __shared__ $T temp[1024];
        int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
        //int idx_y = blockIdx.y * blockDim.y + threadIdx.y;
        int idx_z = blockIdx.z * blockDim.z + threadIdx.z;
        //if (idx_y >= x.dims[1]) return;
        int tid = threadIdx.y;
        int idx_y = blockIdx.y * $blocksize * 2 + threadIdx.y;
        $T a = (idx_y < x.dims[1]) ? x(idx_x, idx_y, idx_z) : 0;

        if (idx_y+$blocksize < x.dims[1]) {
            a += x[idx_y+$blocksize];
        }
        temp[tid] = a;
        $T b;
        __syncthreads();

        if (($blocksize >= 512) && (tid < 256)) {
            b = temp[tid+256];
            a = $op;
            temp[tid] = a;
            __syncthreads();
        }
        if (($blocksize >= 256) &&(tid < 128)) {
            b = temp[tid+128];
            a = $op;
            temp[tid] = a;
            __syncthreads();
        }
        if (($blocksize >= 128) && (tid < 64)) {
            b = temp[tid+64];
            a = $op;
            temp[tid] = a;
            __syncthreads();
        }
        if (tid < 32) {
            if ($blocksize >= 64) a += temp[tid+32];
            for (int offset = warpSize/2; offset > 0; offset /= 2) {
                a += __shfl_down(a, offset);
            }
        }
        if (threadIdx.y == 0) y(blockIdx.x, blockIdx.y, blockIdx.z) = a;
    }""")
end


text:
using Mea
using Base.Test

tests = [
        "green",
        "periodize",
        #"sigmadc",
        "transport"
        ]

if length(ARGS) > 0
    tests = ARGS
end

for t in tests
    test_file = "test_$t.jl"
    print_with_color(:green, "* $test_file\n")
    include(test_file)
end


text:
@testset "create_struct_AST" begin
    # Mutable
    eval(GraphQLClient.create_struct_AST(
        :NewType,
        Number,
        Dict(:field1 => Int)
    ))
    @test fieldnames(NewType) == (:field1,)
    new_type = NewType(1)
    @test ismutable(new_type) # Check mutability
    @test new_type isa Number # Check supertype
    new_type.field1 = nothing # Check Union
    @test isnothing(new_type.field1)

    # Immutable
    eval(GraphQLClient.create_struct_AST(
        :NewTypeImmutable,
        Number,
        Dict(:field1 => Int),
        mutable=false
    ))
    @test fieldnames(NewTypeImmutable) == (:field1,)
    @test !ismutable(NewTypeImmutable(1)) # Check mutability
    @test NewTypeImmutable(1) isa Number # Check supertype
end

@testset "build_name_to_type introspection" begin
    # Simple scalar types
    client = Client("url", "ws"; introspect=false)
    client.type_to_fields_map = Dict(
        "MyObject" => Dict(
            "my_scalar" => build_type(nothing, "my_scalar", type=
                build_type("SCALAR", "String")),
            "my_nonnull_scalar" => build_type("", "my_nonnull_scalar", type=
                build_type("NON_NULL", nothing, ofType=
                    build_type("SCALAR", "Boolean")))
        )
    )
    name_to_type = GraphQLClient.build_name_to_type(client, "MyObject")
    @test name_to_type[:my_nonnull_scalar] == Bool
    @test name_to_type[:my_scalar] == String

    # Custom scalar
    push!(client.type_to_fields_map, "MyObjectWithCustomScalar" => Dict(
        "my_custom_scalar_field" => build_type(nothing, "my_custom_scalar_field", type=
                build_type("SCALAR", "CustomScalarType"))
    ))
    @test_throws ArgumentError GraphQLClient.build_name_to_type(client, "MyObjectWithCustomScalar")
    name_to_type = GraphQLClient.build_name_to_type(client, "MyObjectWithCustomScalar", scalar_types=Dict("CustomScalarType" => Int8))
    @test name_to_type[:my_custom_scalar_field] == Int8

    # Object - stopped by recursion
    push!(client.type_to_fields_map, "MyTopLevelObject" => Dict(
        "object_field" => build_type(nothing, "object_field", type=
            build_type("OBJECT", "MyObject"))
    ))

    @test_logs (:warn, r"Cannot introspect field") GraphQLClient.build_name_to_type(
        client,
        "MyTopLevelObject",
        objects_being_introspected=["MyObject"],
    )

    @test_logs (:warn, r"Cannot introspect field") GraphQLClient.build_name_to_type(
        client,
        "MyTopLevelObject",
        allowed_level=1,
    )

    # Object - introspected
    name_to_type = GraphQLClient.build_name_to_type(client, "MyTopLevelObject")
    @test name_to_type[:object_field] == GraphQLClient.get_introspected_type(client, "MyObject")

    # Already introspected vector
    push!(client.type_to_fields_map, "MyTopLevelVectorObject" => Dict(
        "vector_object_field" => build_type(nothing, "vector_object_field", type=
            build_type("LIST", nothing, ofType=
                build_type("OBJECT", "MyObject")))))

    name_to_type = GraphQLClient.build_name_to_type(client, "MyTopLevelVectorObject")
    @test name_to_type[:vector_object_field] == Vector{GraphQLClient.get_introspected_type(client, "MyObject")}
end

@testset "_instrospect_object" begin
    # Object with scalar types
    client = Client("url", "ws"; introspect=false)
    client.type_to_fields_map = Dict(
        "MyObject" => Dict(
            "my_scalar" => build_type(nothing, "my_scalar", type=
                build_type("SCALAR", "String")),
            "my_nonnull_scalar" => build_type("", "my_nonnull_scalar", type=
                build_type("NON_NULL", nothing, ofType=
                    build_type("SCALAR", "Boolean")))
        )
    )

    T = GraphQLClient._instrospect_object(client, "MyObject")
    @test hasfield(T, :my_scalar)
    @test hasfield(T, :my_nonnull_scalar)
    @test T <: GraphQLClient.AbstractIntrospectedStruct

    # check doesn't get introspected again
    T2 = GraphQLClient._instrospect_object(client, "MyObject")
    @test T===T2
end

@testset "introspect_object" begin
    client = Client("url", "ws"; introspect=false)
    client.type_to_fields_map = Dict(
        "MyObject" => Dict(
            "my_scalar" => build_type(nothing, "my_scalar", type=
                build_type("SCALAR", "String")),
            "my_sub_object" => build_type("", "my_sub_object", type=
                build_type("OBJECT", "MySubObject"))
        ),
        "MySubObject" => Dict(
            "my_scalar_2" => build_type(nothing, "my_scalar_2", type=
                build_type("SCALAR", "Float"))
        )
    )
    # Test objects are not re-introspected
    T1 = GraphQLClient.introspect_object(client, "MyObject")
    T1_sub = GraphQLClient.get_introspected_type(client, "MySubObject")
    T2 = GraphQLClient.introspect_object(client, "MyObject")
    T2_sub = GraphQLClient.get_introspected_type(client, "MySubObject")
    @test T1 === T2
    @test T1_sub === T2_sub
    # Test global client methods
    global_graphql_client(client)
    T2_global = GraphQLClient.introspect_object("MyObject")
    @test T1 === T2_global
    @test GraphQLClient.get_introspected_type("MySubObject") === T1_sub
    # Test that force re-introspects all objects that are requested and are subobjects and no others
    T3 = GraphQLClient.introspect_object(client, "MyObject", force=true)
    T3_sub = GraphQLClient.get_introspected_type(client, "MySubObject")
    @test T3 !== T2
    @test T3_sub !== T2_sub
    T4_sub = GraphQLClient.introspect_object(client, "MySubObject", force=true)
    T4 = GraphQLClient.get_introspected_type(client, "MyObject")
    @test T4 == T3
    @test T4_sub !== T3_sub
    # Test reset all resets all
    T5_sub = GraphQLClient.introspect_object(client, "MySubObject", reset_all=true)
    @test T5_sub !== T4_sub
    @test !haskey(client.introspected_types, "MyObject") # deleted

    # Test that already introspected subobject is used
    GraphQLClient.introspect_object(client, "MyObject")
    @test GraphQLClient.get_introspected_type(client, "MySubObject") === T5_sub

    # Test that mutable cascades to all sub types
    @test ismutable(GraphQLClient.initialise_introspected_struct(client, "MyObject"))
    @test ismutable(GraphQLClient.initialise_introspected_struct(client, "MySubObject"))
    GraphQLClient.introspect_object(client, "MyObject", force=true, mutable=false)
    @test !ismutable(GraphQLClient.initialise_introspected_struct(client, "MyObject"))
    @test !ismutable(GraphQLClient.initialise_introspected_struct(client, "MySubObject"))

    # Test parent types
    T4 = GraphQLClient.introspect_object(client, "MyObject", force=true, parent_type=Number)
    T4_sub = GraphQLClient.get_introspected_type(client, "MySubObject")
    @test T4 <: Number
    @test T4_sub <: GraphQLClient.AbstractIntrospectedStruct # only top level set
    T5 = GraphQLClient.introspect_object(client, "MyObject", force=true, parent_map=Dict("MyObject" => Real, "MySubObject" => AbstractString))
    T5_sub = GraphQLClient.get_introspected_type(client, "MySubObject")
    @test T5 <: Real
    @test T5_sub <: AbstractString
    @test_logs (:warn, r"Parent type for MyObject supplied in both parent_map and parent_type kwarg") GraphQLClient.introspect_object(client, "MyObject", force=true, parent_type=Number, parent_map=Dict("MyObject" => AbstractString, "MySubObject" => AbstractString))
    T6 = GraphQLClient.introspect_object(client, "MyObject", force=true, parent_type=Number, parent_map=Dict("MyObject" => AbstractString, "MySubObject" => AbstractString))
    T6_sub = GraphQLClient.get_introspected_type(client, "MySubObject")
    @test T6 <: Number
    @test T6_sub <: AbstractString

    # Test force only changes structs that are re-introspected and not others
    client.type_to_fields_map["MyNewObject"] =  Dict(
            "my_scalar" => build_type(nothing, "my_scalar", type=
                build_type("SCALAR", "String")))
    GraphQLClient.introspect_object(client, "MyNewObject", force=true)
    T7_sub = GraphQLClient.get_introspected_type(client, "MySubObject")
    @test T7_sub === T6_sub

    # Test recursion
    client = Client("url", "ws"; introspect=false)
    client.type_to_fields_map = Dict(
        "MyObject" => Dict(
            "my_sub_object" => build_type("", "my_sub_object", type=
                build_type("OBJECT", "MySubObject"))
        ),
        "MySubObject" => Dict(
            "my_object" => build_type(nothing, "my_object", type=
                build_type("OBJECT", "MyObject")),
            "my_scalar" => build_type(nothing, "my_scalar", type=
                build_type("SCALAR", "String"))
        )
    )
    GraphQLClient.introspect_object(client, "MyObject", force=true)
    T_subobject = GraphQLClient.get_introspected_type(client, "MySubObject")
    @test !hasfield(T_subobject, :my_object)

    # Test scalar type
    client = Client("url", "ws"; introspect=false)
    push!(client.type_to_fields_map, "MyObjectWithCustomScalar" => Dict(
        "my_custom_scalar_field" => build_type(nothing, "my_custom_scalar_field", type=
                build_type("SCALAR", "CustomScalarType")),
        "my_scalar_field" => build_type(nothing, "my_scalar_field", type=
                build_type("SCALAR", "String"))
    ))
    T = GraphQLClient.introspect_object(client, "MyObjectWithCustomScalar", custom_scalar_types=Dict("CustomScalarType" => Int), force=true)
    @test hasfield(T, :my_custom_scalar_field)
    @test hasfield(T, :my_scalar_field)
    @test fieldtypes(T) == (Union{Nothing, Int64}, Union{Nothing, String})
end

@testset "Initialising" begin
    client = Client("url", "ws"; introspect=false)
    client.type_to_fields_map = Dict(
        "MyObject" => Dict(
            "my_scalar" => build_type(nothing, "my_scalar", type=
                build_type("SCALAR", "String")),
            "my_sub_object" => build_type("", "my_sub_object", type=
                build_type("OBJECT", "MySubObject"))
        ),
        "MySubObject" => Dict(
            "my_scalar_2" => build_type(nothing, "my_scalar_2", type=
                build_type("SCALAR", "Float"))
        )
    )
    GraphQLClient.introspect_object(client, "MyObject")

    # Test Initialise
    my_object = GraphQLClient.initialise_introspected_struct(client, SubString("MyObject"))
    @test all(isnothing(getproperty(my_object, property)) for property in propertynames(my_object))
    my_object = GraphQLClient.initialise_introspected_struct(SubString("MyObject"))
    @test all(isnothing(getproperty(my_object, property)) for property in propertynames(my_object))

    # Test create
    my_sub_object = GraphQLClient.create_introspected_struct(
        client,
        "MySubObject",
        Dict("my_scalar_2" => 3.0)
    )

    my_object = GraphQLClient.create_introspected_struct(
        client,
        "MyObject",
        Dict("my_scalar" => "string", "my_sub_object" => my_sub_object)
    )

    @test my_object.my_sub_object.my_scalar_2 == 3.0

    # create with global client
    global_graphql_client(client)
    my_sub_object = GraphQLClient.create_introspected_struct(
        "MySubObject",
        Dict("my_scalar_2" => 3.0)
    )
    my_object = GraphQLClient.create_introspected_struct(
        "MyObject",
        Dict("my_scalar" => "string", "my_sub_object" => my_sub_object)
    )
    @test my_object.my_sub_object.my_scalar_2 == 3.0

    # Only can create for mutable
    GraphQLClient.introspect_object(client, "MyObject"; force=true, mutable=false)

    @test_throws GraphQLClient.GraphQLClientException GraphQLClient.create_introspected_struct(
        client,
        "MyObject",
        Dict("my_scalar" => "string", "my_sub_object" => my_sub_object)
    )
end


text:
export Layer, LayerState, AbstractParameter
export copy_solver_state!

@compat abstract type Layer end
@compat abstract type LayerState end

# forward declaration
@compat abstract type AbstractParameter end


text:
export build_SDP_Instance_from_SDPPrimal

function build_SDP_Instance_from_SDPPrimal(sdpprimal::SDPPrimal)
    sdp_pb = SDP_Problem()

    # Build numerical descirption of problem, indexed on strings
    set_matrices!(sdp_pb, sdpprimal)
    set_linear!(sdp_pb, sdpprimal)
    set_const!(sdp_pb, sdpprimal)

    # Build structural information and string to id maps
    set_blockvartypes!(sdp_pb, sdpprimal)

    set_constraints!(sdp_pb)
    set_blocks!(sdp_pb)
    set_scalvars!(sdp_pb)

    return sdp_pb
end

function set_matrices!(sdp_pb::SDP_Problem, sdpprimal::SDPPrimal)
    for ((moment, block, γ, δ), coeff) in sdpprimal.blocks
        γ_str, δ_str = format_string(γ), format_string(δ)

        ctr_name = (format_string(moment.conj_part), format_string(moment.expl_part), moment.clique)
        block_name = block
        row_name = max(γ_str, δ_str)
        col_name = min(γ_str, δ_str)

        # sanity check
        @assert !haskey(sdp_pb.matrices, (ctr_name, block_name, row_name, col_name))
        @assert imag(coeff) == 0

        sdp_pb.matrices[(ctr_name, block_name, row_name, col_name)] = coeff
    end
end

function set_linear!(sdp_pb::SDP_Problem, sdpprimal::SDPPrimal)
    # Building symmetric terms into linear terms
    for ((moment, block, var), coeff) in sdpprimal.linsym
        ctr_name = (format_string(moment.conj_part), format_string(moment.expl_part), moment.clique)
        var_name = format_string(var, block)

        # sanity check
        @assert !haskey(sdp_pb.linear, (ctr_name, var_name))
        @assert imag(coeff) == 0

        sdp_pb.linear[(ctr_name, var_name)] = coeff
    end

    # Building scalar linear terms
    for ((moment, var), coeff) in sdpprimal.lin
        ctr_name = (format_string(moment.conj_part), format_string(moment.expl_part), moment.clique)
        var_name = format_string(var)

        # sanity check
        @assert !haskey(sdp_pb.linear, (ctr_name, var_name))
        @assert imag(coeff) == 0

        sdp_pb.linear[(ctr_name, var_name)] = coeff
    end
end

function set_const!(sdp_pb::SDP_Problem, sdpprimal::SDPPrimal)
    for (moment, coeff) in sdpprimal.cst
        ctr_name = (format_string(moment.conj_part), format_string(moment.expl_part), moment.clique)

        # sanity check
        @assert !haskey(sdp_pb.cst_ctr, ctr_name)
        @assert imag(coeff) == 0

        sdp_pb.cst_ctr[ctr_name] = coeff
    end
end

"""
set_blockvartypes!(sdp_pb::SDP_Problem, sdpprimal::SDPPrimal)

Set attributes `name_to_sdpblock` and `id_to_sdpblock` given a primal sdp.
"""
function set_blockvartypes!(sdp_pb::SDP_Problem, sdpprimal::SDPPrimal)
    n_sdp = 0

    for (blockname, blocktype) in sdpprimal.block_to_vartype
        if blocktype == :SDP
            n_sdp += 1

            block = SDP_Block(n_sdp, blockname)
            sdp_pb.name_to_sdpblock[blockname] = block
            sdp_pb.id_to_sdpblock[n_sdp] = block
        elseif blocktype ∉ Set([:Sym])
            error("set_blockvartypes!(): Unknown block variable type '$blocktype' for variable $blockname")
        end
    end
end