text:
# The functions below check if a matrix is approximately Hermitian
using Compat.LinearAlgebra: RealHermSymComplexHerm

const ERR = 1E-13

isapproxhermitian(A::RealHermSymComplexHerm) = true

function isapproxhermitian(A::AbstractMatrix)
    m, n = size(A)
    m == n || return false
    @inbounds for j = 1:m
        for i = 1:j-1
            Aij, Aji = A[i,j], A[j,i]
            if abs(real(Aij)-real(Aji)) > ERR || abs(imag(Aij)+imag(Aji)) > ERR
                return false
            end
        end
        abs(imag(A[j,j])) > ERR && return false
    end
    return true
end

function isapproxhermitian(A::SparseMatrixCSC)
    # Modified from base
    m, n = size(A)
    m == n || return false
    colptr = A.colptr
    rowval = A.rowval
    nzval = A.nzval
    @inbounds for col = 1:m
        for p = colptr[col]:colptr[col+1]-1
            val = nzval[p]
            row = rowval[p]
            if val == 0 # ignore stored zeros
                continue
            end
            if row == col # diagonal element
                abs(imag(val)) > ERR && return false
            else # off-diagonal element
                val2 = A[col,row]
                if abs(real(val)-real(val2)) > ERR || abs(imag(val)+imag(val2)) > ERR
                    return false
                end
            end
        end
    end
    return true
end

function hermitianize!(A::AbstractMatrix)
    m = checksquare(A)
    @inbounds for j = 1:m
         for i = 1:j-1
            A[i,j] = conj(A[j,i])
        end
        A[j,j] = real(A[j,j])
    end
    return A
end

function isunitary(A::AbstractMatrix{T}) where {T}
    # check efficiently if A'*A = I
    m, n = size(A)
    m == n || return false
    @inbounds for i = 1:m, j = 1:m
        a = zero(T)
        for k = 1:m
            a += conj(A[k,i]) * A[k,j]
        end
        if i == j
            a == 1 || return false
        else
            a == 0 || return false
        end
    end
    return true
end

function isunitary(A::SparseMatrixCSC{T,S}) where {T,S}
    # check efficiently if A'*A = I
    m, n = size(A)
    m == n || return false
    colptr = A.colptr
    rowval = A.rowval
    nzval = A.nzval
    @inbounds for i = 1:m, j = 1:m
        a = zero(T)
        for p1 = colptr[i]:colptr[i+1]-1, p2 = colptr[j]:colptr[j+1]-1
            if rowval[p1] < rowval[p2]
                continue
            elseif rowval[p1]==rowval[p2]
                a += conj(nzval[p1]) * nzval[p2]
                continue
            end
        end
        if i == j
            a == 1 || return false
        else
            a == 0 || return false
        end
    end
    return true
end

function isapproxunitary(A::AbstractMatrix{T}) where {T}
    # check efficiently if A'*A ≈ I
    m, n = size(A)
    m == n || return false
    @inbounds for i = 1:m, j = 1:m
        a = zero(T)
        for k = 1:m
            a += conj(A[k,i]) * A[k,j]
        end
        if i == j
            if abs(1-real(a)) > ERR || abs(imag(a)) > ERR
                return false
            end
        else
            if abs(real(a)) > ERR || abs(imag(a)) > ERR
                return false
            end
        end
    end
    return true
end

function isapproxunitary(A::SparseMatrixCSC{T,S}) where {T,S}
    # check efficiently if A'*A ≈ I
    m, n = size(A)
    m == n || return false
    colptr = A.colptr
    rowval = A.rowval
    nzval = A.nzval
    @inbounds for i = 1:m, j = 1:m
        a = zero(T)
        for p1 = colptr[i]:colptr[i+1]-1, p2 = colptr[j]:colptr[j+1]-1
            if rowval[p1] < rowval[p2]
                continue
            elseif rowval[p1]==rowval[p2]
                a += conj(nzval[p1]) * nzval[p2]
                continue
            end
        end
        if i == j
            if abs(1-real(a)) > ERR || abs(imag(a)) > ERR
                return false
            end
        else
            if abs(real(a)) > ERR || abs(imag(a)) > ERR
                return false
            end
        end
    end
    return true
end


text:
include("deepmap.jl")
include("fileio.jl")
include("preprocess.jl")
include("resnet.jl"); export ResNet
include("mobilenet.jl"); export MobileNet
include("efficientnet.jl"); export EfficientNet

ext: jl


text:
Base.@kwdef struct DGSBPNormal <: AbstractModel
    # Data
    y0::Vector{Float64}
    X0::Matrix{Float64}
    y1::Vector{Float64}
    X1::Matrix{Float64}
    mapping::Vector{Vector{Int}} = [[i] for i in 1:size(X0, 2)]
    update_g::Vector{Bool} = ones(Bool, length(mapping))
    # Transformed data
    D0::Int = size(X0, 2)
    # Hyperparameters
    m_μ0::Float64 = 0.0
    c_μ0::Float64 = 1.0
    a_τ0::Float64 = 0.1
    b_τ0::Float64 = 0.1
    m0::Vector{Float64} = zeros(D0)
    B0::Symmetric{Float64, Matrix{Float64}} = Symmetric(9.0 * I(D0))
    # Parameters
    τ::Vector{Float64} = ones(1)
    μ::Vector{Float64} = zeros(1)
    b::Vector{Vector{Float64}} = [zeros(D0)]
    # Transformed parameters
    ȳ::Vector{Float64} = zeros(1)
    v::Vector{Float64} = zeros(1)
    # Skeleton
    skl::Skeleton = Skeleton(; y0, y1, X0, X1, mapping, update_g)
end

function skeleton(m::DGSBPNormal)
    return m.skl
end

function kernel_pdf(m::DGSBPNormal, yi::Float64, xi::Vector{Float64}, j::Int)
    kernel = Normal(m.μ[j], 1 / √m.τ[j])
    return pdf(kernel, yi)
end

function update_atoms!(m::DGSBPNormal)
    (; m_μ0, c_μ0, a_τ0, b_τ0, μ, τ, ȳ, v, skl) = m
    update_suffstats!(m)
    n = cluster_sizes(skl)
    while length(μ) < rmax(skl)
        new_τ = rand(Gamma(a_τ0, 1 / b_τ0))
        new_μ = rand(Normal(m_μ0, √(c_μ0 / new_τ)))
        push!(τ, new_τ)
        push!(μ, new_μ)
    end
    for j in 1:rmax(skl)
        m_μ1 = (c_μ0 * n[j] * ȳ[j] + m_μ0) / (c_μ0 * n[j] + 1)
        c_μ1 = c_μ0 / (c_μ0 * n[j] + 1)
        a_τ1 = a_τ0 + n[j] / 2
        b_τ1 = b_τ0 + n[j] * (ȳ[j] - m_μ0)^2 / 2 / (c_μ0 * n[j] + 1) + v[j] / 2
        τ[j] = rand(Gamma(a_τ1, 1 / b_τ1))
        μ[j] = rand(Normal(m_μ1, √(c_μ1 / τ[j])))
    end
    return nothing
end

function update_suffstats!(m::DGSBPNormal)
    (; y0, ȳ, v, skl) = m
    d = cluster_labels(skl)
    n = cluster_sizes(skl)
    y = y0
    if length(ȳ) < rmax(skl)
        resize!(ȳ, rmax(skl))
        resize!(v, rmax(skl))
    end
    for j in 1:rmax(skl)
        ȳ[j] = 0.0
        v[j] = 0.0
    end
    for i in 1:length(y)
        ȳ[d[i]] += y[i]
    end
    for j in 1:rmax(skl)
        iszero(n[j]) && continue
        ȳ[j] /= n[j]
    end
    for i in 1:length(y)
        v[d[i]] += (y[i] - ȳ[d[i]])^2
    end
end


text:
### Start 31.07.2017

# "Check if the condition holds"
# function condForReduce2(subsig, delta)
#     n = length(subsig)
#     @assert n > 1
#     for i = 1 : n-1
#         s0 = subsig[1:i]
#         s1 = subsig[(i+1):n]
#         ecut = calcEcut(length(s0), length(s1), delta)
#         #ecut = calcEcut(subsig[1:n], length(s0), length(s1), delta)
#         m0 = mean(s0)
#         m1 = mean(s1)
#         if abs(m1-m0) >= ecut
#             return true
#         end
#     end
#     return false
# end
#
# function condForReduce(subsig, ecut_custom)
#     n = length(subsig)
#     @assert n > 1
#     for i = 1 : n-1
#         s0 = subsig[1:i]
#         s1 = subsig[(i+1):n]
#         m0 = mean(s0)
#         m1 = mean(s1)
#         if abs(m1-m0) >= ecut_custom
#             return true
#         end
#     end
#     return false
# end

# function test_cond()
#     println( cond([1, 2;], 1))
#     println( cond([1, 2;], 0.9))
# end
# test_cond()
#
# function recReduce!(w, ecut_custom)
#     """OLD"""
#     if length(w) > 1
#         if condForReduce(w, ecut_custom)
#             recReduce!(w[1:(end-1)], ecut_custom)
#         end
#     else
#         return w
#     end
# end
#
# function test_recReduce()
#     function loctest(_v; e = 1.0)
#         recReduce!(_v, e)
#         println(_v)
#     end
#     loctest([1.0])
#     loctest([1.0, 3.0, 3.1])
#     loctest([1.0, 3.0, 3.1, 3.1, 3.1])
#     loctest([1.0, 1.0, 1.0, 1.0, 1.0])
#     loctest([1.0, 1.0, 1.0, 4.0, 4.0])
# end
# #test_recReduce()

### End 31.07.2017

include("./ListItem.jl")
include("./List.jl")

mutable struct AdwinObject
    DELTA                   :: Float64
    mintMinimLongitudWindow :: Int
    mdbldelta               :: Float64
    mintTime                :: Int
    mintClock               :: Int
    mdblWidth               :: Float64  # Mean of Width = mdblWidth/Number of items
    #
    #BUCKET
    MAXBUCKETS              :: Int
    lastBucketRow           :: Int
    TOTAL                   :: Float64
    VARIANCE                :: Float64
    WIDTH                   :: Int
    BucketNumber            :: Int
    #
    Detect                  :: Int
    numberDetections        :: Int
    DetectTwice             :: Int
    blnBucketDeleted        :: Bool
    BucketNumberMAX         :: Int
    mintMinWinLength        :: Int

    listRowBuckets          :: List
end

function initAdwinObject()
    AO = AdwinObject()

end

function adwin(d::Float64)
    DELTA::Float64 = 0.002
    mintMinimLongitudWindow::Int = 10
    mdbldelta::Float64 = 0.002
    mintTime::Int = 0
    mintClock::Int = 32
    mdblWidth::Float64 = 0 # Mean of Width = mdblWidth/Number of items

    #BUCKET
    MAXBUCKETS::Int = 5
    lastBucketRow::Int = 0
    TOTAL::Float64 = 0
    VARIANCE::Float64 = 0
    WIDTH::Int = 0
    BucketNumber::Int = 0

    Detect::Int = 0
    numberDetections::Int = 0
    DetectTwice::Int  = 0
    blnBucketDeleted::Bool = false
    BucketNumberMAX::Int = 0
    mintMinWinLength::Int = 5

    listRowBuckets::List = List()

    function insertElement!(value::Float64)
        WIDTH +=1
    end

    #function setInput!(intEntrada::Float64)
    #    setInput!(intEntrada, mdbldelta)
    #end

    function setInput!(intEntrada::Float64, delta::Float64 = mdbldelta)
        blnChange::Bool = false
        blnExit::Bool = false
        cursor::ListIt
        #mintTime += 1

        # 1,2)Increment window in one element
        #insertElement!(intEntrada)
    end

    for p in 0:999
        x = p < 500 ? 1000.0 : 500.0
        setInput!(x)
    end

end

adwin(0.01)

function a()
    a = 0
    function setInput(b)
        a = a+1
        c = b
    end
    println(a)
    setInput(5)
    println(a)
end
a()

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
function num_var(model::CentroidalQuadrupedIC)
    nq = model.nq
    nc = model.nc
    nb = 2

    nq + nc + nc + nc + nc * nb + nc + nc * nb
end

function indices_z(model::CentroidalQuadrupedIC)
    nq = model.nq
    nc = model.nc
    nb = 2
    q = collect(1:nq)
    γ = collect(nq .+ (1:nc))
    sγ = collect(nq + nc .+ (1:nc))
    ψ = collect(nq + nc + nc .+ (1:nc))
    b = collect(nq + nc + nc + nc .+ (1:(nc * nb)))
    sψ = collect(nq + nc + nc + nc + nc * nb .+ (1:nc))
    sb = collect(nq + nc + nc + nc + nc * nb + nc .+ (1:(nc * nb)))
    IndicesZ(q, γ, sγ, ψ, b, sψ, sb)
end

function indices_optimization(model::CentroidalQuadrupedIC)
    nq = model.nq
    nc = model.nc
    nb = 2
    nz = num_var(model)

    IndicesOptimization(
        nz,
        nz,
        [collect(nq .+ (1:nc)), collect(nq + nc .+ (1:nc))],
        [collect(nq .+ (1:nc)), collect(nq + nc .+ (1:nc))],
        [
            [collect([nq + 2nc + 1, (nq + 2nc + nc + 0 .+ (1:2))...]), collect([nq + 2nc + nc + nc * nb + 1, (nq + 2nc + nc + nc * nb + nc + 0 .+ (1:2))...])],
            [collect([nq + 2nc + 2, (nq + 2nc + nc + 2 .+ (1:2))...]), collect([nq + 2nc + nc + nc * nb + 2, (nq + 2nc + nc + nc * nb + nc + 2 .+ (1:2))...])],
            [collect([nq + 2nc + 3, (nq + 2nc + nc + 4 .+ (1:2))...]), collect([nq + 2nc + nc + nc * nb + 3, (nq + 2nc + nc + nc * nb + nc + 4 .+ (1:2))...])],
            [collect([nq + 2nc + 4, (nq + 2nc + nc + 6 .+ (1:2))...]), collect([nq + 2nc + nc + nc * nb + 4, (nq + 2nc + nc + nc * nb + nc + 6 .+ (1:2))...])],
            [collect([nq + 2nc + 5, (nq + 2nc + nc + 8 .+ (1:2))...]), collect([nq + 2nc + nc + nc * nb + 5, (nq + 2nc + nc + nc * nb + nc + 8 .+ (1:2))...])],
        ],
        [
            [collect([nq + 2nc + 1, (nq + 2nc + nc + 0 .+ (1:2))...]), collect([nq + 2nc + nc + nc * nb + 1, (nq + 2nc + nc + nc * nb + nc + 0 .+ (1:2))...])],
            [collect([nq + 2nc + 2, (nq + 2nc + nc + 2 .+ (1:2))...]), collect([nq + 2nc + nc + nc * nb + 2, (nq + 2nc + nc + nc * nb + nc + 2 .+ (1:2))...])],
            [collect([nq + 2nc + 3, (nq + 2nc + nc + 4 .+ (1:2))...]), collect([nq + 2nc + nc + nc * nb + 3, (nq + 2nc + nc + nc * nb + nc + 4 .+ (1:2))...])],
            [collect([nq + 2nc + 4, (nq + 2nc + nc + 6 .+ (1:2))...]), collect([nq + 2nc + nc + nc * nb + 4, (nq + 2nc + nc + nc * nb + nc + 6 .+ (1:2))...])],
            [collect([nq + 2nc + 5, (nq + 2nc + nc + 8 .+ (1:2))...]), collect([nq + 2nc + nc + nc * nb + 5, (nq + 2nc + nc + nc * nb + nc + 8 .+ (1:2))...])],
        ],
        collect(1:(nq + nc + nc * nb + nc)),
        collect(nq + nc + nc * nb + nc .+ (1:nc)),
        collect(nq + nc + nc * nb + nc + nc .+ (1:(nc * 3))),
        [collect(nq + nc + nc * nb + nc + nc + (i - 1) * 3 .+ (1:3)) for i = 1:nc],
        collect(nq + nc + nc * nb + nc .+ (1:(nc + nc * 3))))
end

Trajectory(model::CentroidalQuadrupedIC, T) = Trajectory(model, T, nc=model.nc, nb=(2 * model.nc))
GradientTrajectory(model::CentroidalQuadrupedIC, T) = GradientTrajectory(model, T, nc=model.nc, nb=(2 * model.nc))

function residual(model::CentroidalQuadrupedIC, mass_matrix, dynamics_bias, kinematics, kinematics_jacobians, z, θ, μ)

    # dimensions
    nq = model.nq
    nu = model.nu
    nw = model.nw
    nc = model.nc
    nb = 2
    nf = length(friction_coefficients(model))

    # unpack data
    q0 = θ[1:nq]
    q1 = θ[nq .+ (1:nq)]
    u1 = θ[2nq .+ (1:nu)]
    w1 = θ[2nq + nu .+ (1:nw)]
    friction = θ[2nq + nu + nw .+ (1:nf)]
    h = θ[2nq + nu + nw + nf .+ (1:1)]

    # unpack variables
    q2 = z[1:nq]
    γ1 = z[nq .+ (1:nc)]
    sγ1 = z[nq + nc .+ (1:nc)]
    ψ1 = z[nq + nc + nc .+ (1:nc)]
    b1 = z[nq + nc + nc + nc .+ (1:(nc * nb))]
    sψ1 = z[nq + nc + nc + nc + (nc * nb) .+ (1:nc)]
    sb1 = z[nq + nc + nc + nc + (nc * nb) + nc .+ (1:(nc * nb))]

    # rotation matrix at surface
    rotation = [I(3) for i = 1:nc] #TODO non-flat surfaces

    # tangential velocity
    v1 = (q2 - q1) / h[1]
    vT = [(rotation[i] * kinematics_jacobians[i](q2) * v1)[1:2] for i = 1:nc]

    # contact forces
    J = contact_jacobian(model, q1)
    λ1 = transpose(J) * vcat([transpose(rotation[i]) * [b1[(i - 1) * nb .+ (1:nb)]; γ1[i]] for i = 1:nc]...)

    # residual
    [
    dynamics(model, mass_matrix, dynamics_bias, h, q0, q1, u1, w1, λ1, q2);
    signed_distance(model, q2) - sγ1;
    vcat([vT[i] - sb1[(i - 1) * nb .+ (1:nb)] for i = 1:nc]...);
    vcat([friction[i] * γ1[i] - ψ1[i] for i = 1:nc]...);
    γ1 .* sγ1 .- μ[1];
    vcat([cone_product([ψ1[i]; b1[(i - 1) * nb .+ (1:nb)]], [sψ1[i]; sb1[(i - 1) * nb .+ (1:nb)]]) - [μ[1]; 0.0; 0.0] for i = 1:nc]...);
    ]
end

function codegen_dynamics(model::CentroidalQuadrupedIC)
    mm = q -> mass_matrix(model, q)
    db = (q, q̇) -> dynamics_bias(model, q, q̇)
    return mm, db
end


text:
# This is test of the code for computing a transit model integrated over a time step.
# Computes derivatives over the timestep.
using PyPlot

#include("../src/integrate_transit_cubature.jl")
#include("../src/integrate_transit_simpson_vec.jl")
include("integrate_transit_simpson.jl")

# Test it out:

include("../../test/loglinspace.jl")

t1 = -1.5; t2 = 1.5; nt = 1000; dt = 0.3
t = zeros(nt)
t .= linearspace(t1,t2,nt)
#t = [t[1235]] ; nt =1
# The following compares different tolerances and maxdepths:
r = 0.1; b0 = 0.5; u_n = [0.3,0.3]; nu = length(u_n)
favg0 = zeros(nt)
favg1 = zeros(nt)
neval1= zeros(Int64,nt)
depthmax1= zeros(Int64,nt)
favg2 = zeros(nt)
neval2= zeros(Int64,nt)
depthmax2= zeros(Int64,nt)
favg3 = zeros(nt)
neval3= zeros(Int64,nt)
depthmax3= zeros(Int64,nt)
trans = transit_init(r,b0,u_n,false)
param = [0.0,1.0,b0]   # [t_0,v,b_0]

# Now, carry out speed test:

function compute_lightcurve!(trans::Transit_Struct{T},param::Array{T,1},t::Array{T,1},favg0::Array{T,1},nt::Int64) where {T <: Real}
@inbounds for i=1:nt
  trans.b = sqrt(param[3]^2+(param[2]*(t[i]-param[1]))^2)
  favg0[i] =  transit_poly!(trans)
end
return
end

compute_lightcurve!(trans,param,t,favg0,nt)
@time compute_lightcurve!(trans,param,t,favg0,nt)

function integrate_lightcurve!(trans::Transit_Struct{T},param::Array{T,1},t::Array{T,1},dt::T,favg1::Array{T,1},nt::Int64,tol::T,maxdepth::Int64,neval_t::Array{Int64,1},depthmax::Array{Int64,1}) where {T <: Real}
# First, find the points of contact:
b0 = abs(param[3]); r = trans.r; t0 = param[1]; v = param[2]
if b0 > (1.0+r)
  println("No transit")
  favg1 .= one(T)
  return
elseif (b0+r) > 1.0
  # Two points of contact:
  tc = [t0 - sqrt((1.0+r)^2-b0^2)/v,t0 + sqrt((1.0+r)^2-b0^2)/v]
  nc = 2
else
  # Four points of contact:
  tc = [t0 - sqrt((1.0+r)^2-b0^2)/v,t0 - sqrt((1.0-r)^2-b0^2)/v,t0 + sqrt((1.0-r)^2-b0^2)/v,t0 + sqrt((1.0+r)^2-b0^2)/v]
  nc = 4
end
dtinv = inv(dt)
fint  = zeros(T,6+trans.n)
@inbounds for i=1:nt
  neval_t[i] = 0
  depthmax[i] = 0
  t1 = t[i]-0.5*dt ; t2 = t[i]+0.5*dt
  if t2 < tc[1] || t1 > tc[nc]
    # No points lie within the transit:
    ftmp=zeros(T,6+trans.n)
    ftmp[1]=one(T)
  else
    # Loop over points of conjunction:
    tlim = [t1]
    for j=1:nc
      if t1 < tc[j] && tc[j] < t2
        push!(tlim,tc[j])
      end
    end
    push!(tlim,t2)
    ftmp=zero(T)
    for j=1:length(tlim)-1
      fint,nevali,depthmaxi =  integrate_timestep_gradient!(param,trans,tlim[j],tlim[j+1],tol,maxdepth)
#      println("r: ",trans.r," b: ",trans.b," fint: ",fint," dt: ",tlim[j+1]-tlim[j])
      neval_t[i] += nevali
      if depthmaxi > depthmax[i]
        depthmax[i] = depthmaxi
      end
      ftmp += fint
    end
    ftmp *= dtinv
#    println("ftmp: ",ftmp)
  end
  favg1[i]=ftmp
  # Convert from d_n to u_n derivatives:
#  println("i: ",i," t: ",t[i]," result: ",ftmp)
end
return
end

#integrate_lightcurve!(trans,param,t,dt,favg1,nt,1e-4,3,neval1)
#@time integrate_lightcurve!(trans,param,t,dt,favg1,nt,1e-4,8,neval1)

ntol = 12
tol = logarithmspace(-13.0,-2.0,ntol)
K = 50
neval_mean = zeros(K,ntol)
precision = zeros(K,ntol)
timing = zeros(K,ntol)
for k=1:K
  texp = 0.01*10.0^(rand())
  q1,q2 = rand(2)
  u_n=[2sqrt(q1)*q2,q1*(1-2*q2)]
  r=0.01*10.0^rand()
  b=rand()*(1+r)
  trans = transit_init(r,b,u_n,true)
  tau = 1.0 # Duration of transit
  t .= linearspace(-0.5*(tau+texp),0.5*(tau+texp),nt)
  v = 2/tau*sqrt((1+r)^2-b^2)
#  t  .= x/v
  param = [0.0,v,b]   # [t_0,v,b_0]
  t_nobin1 = time_ns()
  compute_lightcurve!(trans,param,t,favg0,nt)
  t_nobin = time_ns()-t_nobin1
  integrate_lightcurve!(trans,param,t,texp,favg2,nt,1e-15,16,neval2,depthmax2)
  for j=1:ntol
    t_bin1 = time_ns()
    integrate_lightcurve!(trans,param,t,texp,favg1,nt,tol[j],16,neval1,depthmax1)
    t_bin = time_ns()-t_bin1
    precision[k,j] += maximum(abs.(favg1-favg2))
    neval_mean[k,j] += mean(neval1)
    println("k: ",k," j: ",j," tol: ",tol[j]," neval: ",mean(neval1)," relative time: ",t_bin/t_nobin/mean(neval1)," max depth: ",maximum(depthmax1))
  end
#  clf()
#  plot(t,favg2[:,1])
#  read(STDIN,Char)
end
precision_tol = zeros(ntol)
neval_tol = zeros(ntol)
for j=1:ntol
  precision_tol[j]=median(precision[:,j])
  neval_tol[j] = median(neval_mean[:,j])
end
clf()
loglog(precision_tol,neval_tol)


text:
# Dynamics model flags
abstract type DynamicsFlag end
struct CR3BP <: DynamicsFlag end

# Homotopy flags
abstract type HomotopyFlag end
struct MEMF <: HomotopyFlag end     # Min. Energy -> Min. Fuel
struct HypTanMF <: HomotopyFlag end # Hyperbolic Tangent Minimum Fuel

# Desired integration and returned solution
# Flag indicating integration to be performed for initialization cost function
abstract type IntegrationFlag end
# Returns full state vector and time to final desired time (for early termination)
# as first and second return values respectively.
abstract type InitializationFlag <: IntegrationFlag end
struct Initialization <: InitializationFlag end
struct InitializationWithIntegralCost <: InitializationFlag end

# Flag indicating integration to be performed for solving BVP via
# shooting.
#
# Returns final full state vector
abstract type SolvingFlag <: IntegrationFlag end
struct Solving <: SolvingFlag end
struct SolvingWithSTM <: SolvingFlag end

# Flag indicating full state history is disired
#
# Returns DifferentialEquations.jl Solution object
abstract type SolutionHistoryFlag <: IntegrationFlag end
struct FullSolutionHistory <: SolutionHistoryFlag end
struct FullSolutionHistoryWithSTM <: SolutionHistoryFlag end

# Flag indicating full state history without
# control is desired
struct FullSolutionHistoryNoControl <: SolutionHistoryFlag end
# struct FullSolutionHistoryNoControlWithSTM end # Do not have this implemented!


text:
export ec2_terminate, ec2_addprocs, ec2_launch, ec2_start, ec2_stop, ec2_show_status, ec2_hostnames, ec2_instances_by_tag
export ec2_mount_snapshot

ec2_basic(env::AWSEnv, action::String, params_in::Dict) = ec2_execute(env, action, flatten_params(params_in))
export ec2_basic

flatten_params(d::Dict) = flatten_params(Array{Tuple}(0), "", d)
function flatten_params(ft::Array, pfx::String, d::Dict)
    for (k, v) in collect(d)
        if isa(v, Array)
            for (i, x) in enumerate(v)
                npfx = pfx * k * "." * string(i)
                if isa(x, Array) || isa(x, Dict)
                    flatten_params(ft, npfx * ".", x)
                else
                    push!(ft, (npfx, aws_string(x)))
                end
            end

        elseif isa(v, Dict)
            flatten_params(ft, pfx*k*".", v)
        else
            push!(ft, (pfx*k, aws_string(v)))
        end

    end
    ft
end

# Have a set of simple APIS for common use cases. Will grow as per feature requests.

CHK_ERR(resp::EC2Response) = (typeof(resp.obj) == EC2Error) ? error(ec2_error_str(resp.obj)) : resp.obj

function ec2_terminate(instances; env=AWSEnv())
    if length(instances) > 0
        req = TerminateInstancesType(instancesSet=instances)
        resp = CHK_ERR(TerminateInstances(env, req))
    end
    instances
end

function check_running(env, chk_instances)
    new_chk_instances = String[]
    resp = CHK_ERR(DescribeInstanceStatus(env, instancesSet=chk_instances, includeAllInstances=true))
    statuses = resp.instanceStatusSet
    for status in statuses
        info("Status of $(status.instanceId) is $(status.instanceState.code):$(status.instanceState.name)")
        if status.instanceState.code != 16
            push!(new_chk_instances, status.instanceId)
        end
    end
    new_chk_instances
end

function ec2_hostnames(instances; env=AWSEnv())
    names = NTuple[]
    resp = CHK_ERR(DescribeInstances(env, instancesSet=instances))
    reservs = resp.reservationSet
    for reserv in reservs
        instancesSet = reserv.instancesSet
        for i in instancesSet
            push!(names, (i.instanceId, i.dnsName, i.privateDnsName))
        end
    end
    names
end

function ec2_launch(ami::String, seckey::String; env=AWSEnv(), insttype::String="m1.small", n::Integer=1,
                    owner::String="julia", clustername::String="julia", launchset::String="")
# "c1.medium"
# "m1.small"

    resp = CHK_ERR(RunInstances(env, RunInstancesType(imageId=ami, instanceType=insttype, minCount=n, maxCount=n, keyName=seckey)))
    instances = String[]
    for inst in resp.instancesSet
        push!(instances, inst.instanceId)
    end
    info("Launched instances : " )
    info(instances)

    if launchset == ""
       launchset =  aws_string(now(Dates.UTC))
    end

    # Tag the instances..also tests ec2_basic API..
    resp = ec2_basic(env, "CreateTags", Dict(
        "ResourceId"=>instances,
        "Tag"=>[
            Dict("Key"=>"LaunchSet","Value"=>launchset),
            Dict("Key"=>"ClusterName","Value"=>clustername),
            Dict("Key"=>"Name","Value"=>clustername), # Just so that it is visible on the AWS Web console
            Dict("Key"=>"Owner","Value"=>owner)]))
    if (typeof(resp.obj) == EC2Error)
        error(ec2_error_str(resp.obj))
    else
        resp.obj = CreateTagsResponseType(resp.pd)
        if resp.obj._return != true
            error("error adding tags!")
        end
    end
    info("Tagged instances")

    # Wait for the instances to come to a running state....
    wait_till_running(env, instances, 600.0)

    info("Lanched LaunchSet $launchset, ClusterName $clustername, Owner $owner" )

    instances
end

function generate_sshnames(hostnames, num_workers, workers_per_instance, use_public_dnsname, sshflags, hostuser="ubuntu")
    num_hosts = length(hostnames)
    idx = use_public_dnsname ? 2 : 3
    sshnames = Dict{String, Any}()
    if num_workers > 0
        if num_workers < num_hosts
            for i in 1:num_workers
                host = hostnames[i][idx]
                sshnames["$(hostuser)@$(host)"] = 1
            end
        else
            # distribute it equally across available hosts
            nwmap = Array{Tuple}(0)
            for i in 1:length(hostnames)
                host = hostnames[i][idx]
                push!(nwmap, (host, round(Int, floor(num_workers/num_hosts))))
            end

            # distribute any remaining equally
            num_rem = num_workers - foldl((x,y)->((k,v)=y; x+v), 0, nwmap)
            for i in 1:num_rem
                idxnwmap = ((i-1) % num_hosts) + 1
                (host,nw) = nwmap[idxnwmap]
                nwmap[idxnwmap] = (host,nw+1)
            end

            for (host,nw) in nwmap
                sshnames["$(hostuser)@$(host)"] = nw
            end
        end

    elseif workers_per_instance > 0
        for i in 1:length(hostnames)
            host = hostnames[i][idx]
            sshnames["$(hostuser)@$(host)"] = workers_per_instance
        end
    else
        for i in 1:length(hostnames)
            host = hostnames[i][idx]
            sshnames["$(hostuser)@$(host)"] = :auto
        end
    end
    [(k,v) for (k,v) in sshnames]
end

function ec2_addprocs(instances, ec2_keyfile::String; env=AWSEnv(), hostuser::String="ubuntu",
                        dir=JULIA_HOME, tunnel=true, use_public_dnsname=true, exename=joinpath(JULIA_HOME,Base.julia_exename()),
                        workers_per_instance=0, num_workers=0, machines_only=false)

    hostnames = ec2_hostnames(instances, env=env)
    sshflags = `-o UserKnownHostsFile=/dev/null -o StrictHostKeyChecking=no -o LogLevel=ERROR -i $(ec2_keyfile)`

    sshnames = generate_sshnames(hostnames, num_workers, workers_per_instance, use_public_dnsname, sshflags, hostuser)

    info(sshnames)

    if machines_only
        return sshnames
    else
        pids = addprocs(sshnames, sshflags=sshflags, dir=dir, tunnel=tunnel, exename=exename)
        info("Added workers $pids")
        return pids
    end
end

function wait_till_running(env, instances, timeout)
    start = time()
    chk_instances = check_running(env, instances)
    wait_time = 5.0
    while (((time() - start) < timeout) && (length(chk_instances) > 0))
        info("All instances not yet in a running state. Trying again in $wait_time seconds....")
        sleep(wait_time)

        chk_instances = check_running(env, chk_instances)
    end

    if (length(chk_instances) > 0)
        info("All instances not yet in a running state. Please check after some time.")
    end

    # Wait till they are pingable, the network interfaces take some time to come up
    info("Testing TCP connects (on port 22) to all newly started hosts...")
    hosts = ec2_hostnames(instances, env=env)
    testhosts = copy(hosts)
    testidx = 1
    while true
        try
            for (i,h) in enumerate(testhosts)
                info("Testing connection to ", h[2])
                s=connect(h[2], 22)
                close(s)
                testidx = i
            end
            break;
        catch
            info("Some newly started hosts are still unreachable. Trying again in 2.0 seconds.")
            sleep(2.0)
            testhosts=testhosts[testidx:end]
        end
    end

    info(hosts)
    :ok
end

function ec2_start(instances; env=AWSEnv())
    # Tag the instances..also tests ec2_basic API..
    if length(instances) > 0
        resp = ec2_basic(env, "StartInstances", Dict("InstanceId"=>instances))
        if (typeof(resp.obj) == EC2Error)
            error(ec2_error_str(resp.obj))
        end

        # Wait for the instances to come to a running state....
        wait_till_running(env, instances, 600.0)
    end
    instances
end

function ec2_stop(instances; env=AWSEnv())
    # Tag the instances..also tests ec2_basic API..
    if length(instances) > 0
        resp = ec2_basic(env, "StopInstances", Dict("InstanceId"=>instances))
        if (typeof(resp.obj) == EC2Error)
            error(ec2_error_str(resp.obj))
        end
    end
    instances
end

function ec2_show_status(instances; env=AWSEnv())
    # Tag the instances..also tests ec2_basic API..
    resp = ec2_basic(env, "DescribeInstanceStatus", Dict("InstanceId"=>instances, "IncludeAllInstances"=>true))
    if (typeof(resp.obj) == EC2Error)
        error(ec2_error_str(resp.obj))
    end

    resp.obj = DescribeInstanceStatusResponseType(resp.pd)

    statuses = resp.obj.instanceStatusSet
    for status in statuses
        info("Status of $(status.instanceId) is $(status.instanceState.code):$(status.instanceState.name)")
    end

    [(status.instanceId, status.instanceState.code, status.instanceState.name) for status in statuses]
end

function ec2_instances_by_tag(tag, tagvalue; env=AWSEnv(), running_only=true)
    tagfilter = FilterType(name="tag:" * tag, valueSet=[tagvalue])

    if running_only
        filterset = [tagfilter, FilterType(name="instance-state-name", valueSet=["running"])]
    else
        filterset = [tagfilter]
    end
    instances = String[]

    req = DescribeInstancesType(filterSet=filterset)
    resp = CHK_ERR(DescribeInstances(env, req))
    reservs = resp.reservationSet
    for reserv in reservs
        instancesSet = reserv.instancesSet
        for i in instancesSet
            push!(instances, i.instanceId)
        end
    end
    instances
end

function ec2_mount_snapshot(instance::String, snapshot::String, mount::String, ec2_keyfile::String; env=AWSEnv(), dev="/dev/xvdh",
                                hostuser::String="ubuntu")
    # first get the availability zone....
    resp = CHK_ERR(DescribeInstances(env, instancesSet=[instance]))
    zone = resp.reservationSet[1].instancesSet[1].placement.availabilityZone
    hostname = resp.reservationSet[1].instancesSet[1].dnsName

    resp = CHK_ERR(CreateVolume(env, snapshotId=snapshot, availabilityZone=zone))
    volumeId = resp.volumeId
    status = resp.status

    info("Created volume $volumeId")

    while (status == "creating")
        resp = CHK_ERR(DescribeVolumes(env, volumeSet=[volumeId]))
        status = resp.volumeSet[1].status
    end

    if (status != "available")
        error("Volume is unavailable. Status $status")
    end

    resp = CHK_ERR(AttachVolume(env, instanceId=instance, volumeId=volumeId, device=dev))
    status = resp.status
    while (status == "attaching")
        resp = CHK_ERR(DescribeVolumes(env, volumeSet=[volumeId]))
        status = resp.volumeSet[1].attachmentSet[1].status
    end

    if (status != "attached")
        error("Unable to attach volume. Status $status")
    end

    info("Attached $volumeId to $instance")

    run(`ssh -i $(ec2_keyfile) -o UserKnownHostsFile=/dev/null -o StrictHostKeyChecking=no  -o LogLevel=ERROR $(hostuser)@$(hostname) "mkdir -p $(mount); sudo mount $dev $mount"`)

    info("Mounted volume $volumeId at $mount")
    info("To login use 'ssh -i $(ec2_keyfile) $(hostuser)@$(hostname)'")

    volumeId
end


text:
using Qsosed, Test

@testset "read parameters" begin
    config_path = String(@__DIR__) * "/../configs/config_example.yaml"
    model = QsosedModel(config_path)
    @test model.corona.bh == model.bh
    @test model.bh.M == 1e8 * M_SUN
    @test model.bh.mdot == 0.5
    @test model.bh.spin == 0.0
    @test model.corona.reprocessing == true
    @test model.corona.hard_xray_fraction == 0.02
    @test model.corona.electron_energy == 100
    @test model.warm.electron_energy == 0.2
    @test model.warm.photon_index == 2.5
end


text:
export
    RtmStart,
    RtmStartResponse

import Base.==

@slackmethod(RtmStart, "rtm.start",
begin
    simple_latest::Nullable{UTF8String}
    no_unreads::Nullable{UTF8String}
    mpim_aware::Nullable{UTF8String}
end,

begin
    url::UTF8String
    self::Self
    team::Team
    users::Vector{User}
    channels::Vector{SlackChannel}
    groups::Vector{Group}
    mpims::Nullable{Vector{Mpim}}
    ims::Vector{Im}
    bots::Vector{Bot}
end)

==(a::RtmStartResponse, b::RtmStartResponse) = a.url == b.url


text:
using CLBLAS
const clblas = CLBLAS

import OpenCL
using OpenCL: cl

function testDirectCall()
        device, ctx, queue = clblas.get_next_compute_context()
        A = rand(Float32, 1)
        B = rand(Float32, 1)
        C = rand(Float32, 1)
        println(A, B, C)
        a = cl.cl_float(6.6)
        b = cl.cl_float(7.2)

        bufA = cl.Buffer(Float32, ctx, (:r, :copy),hostbuf=A)
        bufB = cl.Buffer(Float32, ctx, (:r, :copy),hostbuf=B)
        bufC = cl.Buffer(Float32, ctx, (:rw, :copy),hostbuf=C)

        event = cl.UserEvent(ctx, retain=true)
        ptrEvent = [event.id]

        println("ccalling!!!")
        clblas.clblasSgemm(1, 0, 0,
                    convert(Csize_t, 1), convert(Csize_t, 1), convert(Csize_t, 1),
                    a, bufA.id, convert(Csize_t, 0), convert(Csize_t, 1),
                    bufB.id, convert(Csize_t, 0), convert(Csize_t, 1), b,
                    bufC.id, convert(Csize_t, 0), convert(Csize_t, 1),
                    cl.cl_uint(1), [queue.id],
                    cl.cl_uint(0),
                    C_NULL,
                    ptrEvent)
        println("ccall was ok")
        cl.api.clWaitForEvents(cl.cl_uint(1), ptrEvent)

        println("wait was ok")
        result = Vector{Float32}(1)
        cl.enqueue_read_buffer(queue, bufC, result, Unsigned(0), nothing, true)
        println("read was ok")
        println(result[1])

end

function testRandMul()

     a = float32(3)
     A = float32([3])
     println(A)
     println(eltype(A))
     B = float32([3])
     b = float32(3)
     C = float32([3])

     println("Normal:", (a*A[1]*B[1] + b*C[1]))

     startT = time_ns()
     futBuf = clblas.clblasSgemm(clblas.clblasColumnMajor.n, clblas.clblasNoTrans.n, clblas.clblasNoTrans.n,
                                 1, a, A,
                                 1, B,
                                 1, b, C)
     result = fetch(futBuf)
     endT = time_ns()

     println("RandMul took:", (endT - startT)/1000000)
     println("Result:", result)
end

 clblas.setup(true)
function testRandBuf(x, y)
     #startT = time_ns()
     futBuf = clblas.randBuf(x, y)
     #futBuf = clblas.clblasSscal(futBuf, cl.cl_float(13.13))
     mat    =  clblas.fetch(futBuf)
     #endT = time_ns()

     #println(futBuf.actualEvent[:status])
     #println(futBuf.actualEvent[:profile_duration])
     #println("RandBuf took:", (endT - startT)/1000000)
     #println(mat)
end

function testRand(x, y)
    startT = time_ns()
    test = rand(x, y) * 13.13
    endT = time_ns()
    println("Rand took:", (endT - startT)/1000000)
end

#result = Vector{Float32}(length(futBuf.resultMem))
#cl.enqueue_read_buffer(futBuf.queue, futBuf.resultMem, result, Unsigned(0), nothing, true)
#println(result)


text:
module VectorSpaces

using ..Traits: rawtype, efficientoperations, getcontent, hascontent

import ..Traits: contentnames
import ...Interfaces: dimension, rank

export VectorSpace, EnumerativeVectorSpace, CartesianVectorSpace, NamedVectorSpace, ZipNamedVectorSpace, ProductNamedVectorSpace, shape, ndimshape

"""
    VectorSpace{B} <: AbstractVector{B}

Abstract vector space.
"""
abstract type VectorSpace{B} <: AbstractVector{B} end
@inline Base.:(==)(vs₁::VectorSpace, vs₂::VectorSpace) = ==(efficientoperations, vs₁, vs₂)
@inline Base.isequal(vs₁::VectorSpace, vs₂::VectorSpace) = isequal(efficientoperations, vs₁, vs₂)
@inline Base.size(vs::VectorSpace) = (dimension(vs),)
@inline Base.iterate(vs::VectorSpace, state::Integer=1) = state>length(vs) ? nothing : (vs[state], state+1)
@inline Base.issorted(vs::VectorSpace) = false
@inline function Base.findfirst(basis::B, vs::VectorSpace{B}) where B
    if issorted(vs)
        index = invoke(searchsortedfirst, Tuple{Any, Any}, vs, basis)
        return (0<index<=length(vs) && vs[index]==basis) ? index : nothing
    end
    for i = 1:length(vs)
        vs[i]==basis && return i
    end
end
@inline function Base.searchsortedfirst(vs::VectorSpace{B}, basis::B) where B
    issorted(vs) && return invoke(searchsortedfirst, Tuple{Any, Any}, vs, basis)
    for i = 1:length(vs)
        vs[i]==basis && return i
    end
    return length(vs)+1
end
@inline Base.in(basis::B, vs::VectorSpace{B}) where B = isa(findfirst(basis, vs), Integer)

"""
    EnumerativeVectorSpace{B} <: VectorSpace{B}

Abstract enumerative vector space.
"""
abstract type EnumerativeVectorSpace{B} <: VectorSpace{B} end
@inline contentnames(::Type{<:EnumerativeVectorSpace}) = (:table,)
@inline dimension(vs::EnumerativeVectorSpace) = length(getcontent(vs, :table))
@inline Base.getindex(vs::EnumerativeVectorSpace, i::Int) = getcontent(vs, :table)[i]

"""
    CartesianVectorSpace{B} <: VectorSpace{B}

Abstract Cartesian vector space.
"""
abstract type CartesianVectorSpace{B} <: VectorSpace{B} end
function shape end
@inline ndimshape(vs::CartesianVectorSpace) = ndimshape(typeof(vs))
@inline dimension(vs::CartesianVectorSpace) = mapreduce(length, *, shape(vs))
@inline Base.getindex(vs::CartesianVectorSpace, i::CartesianIndex) = rawtype(eltype(vs))(i, vs)
@inline Base.getindex(vs::CartesianVectorSpace, i::Int) = getindex(vs, CartesianIndices(shape(vs))[i])
@inline Base.issorted(vs::CartesianVectorSpace) = true
function Base.iterate(vs::CartesianVectorSpace)
    idx = CartesianIndices(shape(vs))
    index = iterate(idx)
    isnothing(index) && return nothing
    return vs[index[1]], (idx, index[2])
end
function Base.iterate(vs::CartesianVectorSpace, state)
    index = iterate(state[1], state[2])
    isnothing(index) && return nothing
    return vs[index[1]], (state[1], index[2])
end
@inline function Base.findfirst(basis::B, vs::CartesianVectorSpace{B}) where B
    index, idx = CartesianIndex(basis, vs), CartesianIndices(shape(vs))
    index∉idx && return nothing
    return LinearIndices(shape(vs))[index-first(idx)+oneunit(typeof(index))]
end
@inline function Base.searchsortedfirst(vs::CartesianVectorSpace{B}, basis::B) where B
    index, idx = CartesianIndex(basis, vs), CartesianIndices(shape(vs))
    index∈idx && return LinearIndices(shape(vs))[index-first(idx)+oneunit(typeof(index))]
    return searchsortedfirst(idx, index)
end
@inline Base.in(basis::B, vs::CartesianVectorSpace{B}) where B = CartesianIndex(basis, vs)∈CartesianIndices(shape(vs))

"""
    NamedVectorSpace{M, NS, BS<:Tuple} <: CartesianVectorSpace{NamedTuple{NS, BS}}

Abstract named vector space.

This is a wrapper of Cartesian indexable vector spaces, each of whose indexable dimensions is associated with a name.

It has three type parameters:
* `M`: mode of the named vector space. It specifies how the indexable dimensions are combined to form the bases of the named vector space, and must take one of the following values:
  - `:⊕`: elements from each indexable dimensions are zipped together to form the bases,
  - `:⊗`: elements from each indexable dimensions are direct producted together to form the bases.
For the `:⊕` mode, all the indexable dimensions should have the same number of elements, and the number of formed bases is equal to this number; for the `:⊗` mode, there are no restriction on the numbers of the indexable dimensions, and the number of the final bases is equal to their product.
* `NS::Tuple{Vararg{Symbol}}`: the names of the indexable dimensions
* `BS<:Tuple`: the eltypes of the indexable dimensions

The subtypes must have the following predefined content:
* `:contents::Tuple`: storage of the contents of the indexable dimensions
"""
abstract type NamedVectorSpace{M, NS, BS<:Tuple} <: CartesianVectorSpace{NamedTuple{NS, BS}} end
@inline contentnames(::Type{<:NamedVectorSpace}) = (:contents,)
@inline Base.:(==)(vs₁::NamedVectorSpace, vs₂::NamedVectorSpace) = ==(efficientoperations, vs₁, vs₂) && keys(vs₁)==keys(vs₂)
@inline Base.isequal(vs₁::NamedVectorSpace, vs₂::NamedVectorSpace) = isequal(efficientoperations, vs₁, vs₂) && isequal(keys(vs₁), keys(vs₂))
@inline rank(vs::NamedVectorSpace) = rank(typeof(vs))
@inline rank(::Type{<:NamedVectorSpace{:⊕}}) = 1
@inline rank(::Type{<:NamedVectorSpace{:⊗, NS}}) where NS = length(NS)
@inline shape(nvs::NamedVectorSpace{:⊕}) = (1:length(getcontent(nvs, :contents)[1]),)
@inline ndimshape(::Type{T}) where {T<:NamedVectorSpace} = rank(T)
@inline @generated shape(nvs::NamedVectorSpace{:⊗}) = Expr(:tuple, [:(1:length(getcontent(nvs, :contents)[$i])) for i = 1:rank(nvs)]...)
@inline function Base.CartesianIndex(basis::NamedTuple{NS}, nvs::NamedVectorSpace{:⊕, NS}) where NS
    index = findfirst(isequal(basis[1]), getcontent(nvs, :contents)[1])::Int
    for i = 2:length(NS)
        @assert getcontent(nvs, :contents)[i][index] == basis[i] "CartesianIndex error: input basis out of range."
    end
    return CartesianIndex(index)
end
@inline @generated function Base.CartesianIndex(basis::NamedTuple{NS, BS}, nvs::NamedVectorSpace{:⊗, NS, BS}) where {NS, BS<:Tuple}
    return Expr(:call, :CartesianIndex, [:(findfirst(isequal(basis[$i]), getcontent(nvs, :contents)[$i])::Int) for i = 1:length(NS)]...)
end
@inline @generated function Base.NamedTuple(index::CartesianIndex{1}, nvs::NamedVectorSpace{:⊕, NS}) where NS
    return Expr(:tuple, [:($name=getcontent(nvs, :contents)[$i][index[1]]) for (i, name) in enumerate(NS)]...)
end
@inline @generated function Base.NamedTuple(index::CartesianIndex{N}, nvs::NamedVectorSpace{:⊗, NS}) where {N, NS}
    @assert N == length(NS) "NamedTuple error: mismatched input index and rank of named vector space."
    return Expr(:tuple, [:($name=getcontent(nvs, :contents)[$i][index[$i]]) for (i, name) in enumerate(NS)]...)
end

"""
    keys(nvs::NamedVectorSpace) -> Tuple{Vararg{Symbol}}
    keys(::Type{<:NamedVectorSpace{M, NS} where M}) where NS -> Tuple{Vararg{Symbol}}

Get the names of a named vector space.
"""
@inline Base.keys(nvs::NamedVectorSpace) = keys(typeof(nvs))
@inline Base.keys(::Type{<:NamedVectorSpace{M, NS} where M}) where {NS} = NS

"""
    values(nvs::NamedVectorSpace) -> Tuple{Vararg{AbstractVector}}

Get the contents of a named vector space.
"""
@inline Base.values(nvs::NamedVectorSpace) = getcontent(nvs, :contents)

"""
    pairs(nvs::NamedVectorSpace) -> Base.Iterators.Pairs

Get the name-value pairs of a named vector space.
"""
@inline Base.pairs(nvs::NamedVectorSpace) = Base.Generator(=>, nvs|>keys, nvs|>values)

"""
    eltype(nvs::NamedVectorSpace, i::Int)
    eltype(::Type{<:NamedVectorSpace{M, NS, BS} where {M, NS}}, i::Int) where BS

Get the eltype of the ith indexable dimensions of a named vector space.
"""
@inline Base.eltype(nvs::NamedVectorSpace, i::Int) = eltype(nvs|>typeof, i)
@inline Base.eltype(::Type{<:NamedVectorSpace{M, NS, BS} where {M, NS}}, i::Int) where BS = fieldtype(BS, i)

"""
    ZipNamedVectorSpace{NS, BS<:Tuple, CS<:Tuple} <: NamedVectorSpace{:⊕, NS, BS}

Zipped named vector space.
"""
struct ZipNamedVectorSpace{NS, BS<:Tuple, CS<:Tuple} <: NamedVectorSpace{:⊕, NS, BS}
    contents::CS
end

"""
    ZipNamedVectorSpace{NS}(contents...) where NS

Construct a zipped named vector space.
"""
@generated function ZipNamedVectorSpace{NS}(contents...) where NS
    @assert length(NS)==length(contents) && isa(NS, Tuple{Vararg{Symbol}}) "ZipNamedVectorSpace error: mismatched names and contents."
    BS = Expr(:curly, :Tuple, [eltype(contents[i]) for i = 1:length(NS)]...)
    return quote
        @assert mapreduce(length, ==, contents) "ZipNamedVectorSpace error: mismatched length of contents."
        ZipNamedVectorSpace{NS, $BS, typeof(contents)}(contents)
    end
end

"""
    ProductNamedVectorSpace{NS, BS<:Tuple, CS<:Tuple} <: NamedVectorSpace{:⊕, NS, BS}

Product named vector space.
"""
struct ProductNamedVectorSpace{NS, BS<:Tuple, CS<:Tuple} <: NamedVectorSpace{:⊗, NS, BS}
    contents::CS
end

"""
    ProductNamedVectorSpace{NS}(contents...) where NS

Construct a product named vector space.
"""
@generated function ProductNamedVectorSpace{NS}(contents...) where NS
    @assert length(NS)==length(contents) && isa(NS, Tuple{Vararg{Symbol}}) "ProductNamedVectorSpace error: mismatched names and contents."
    BS = Expr(:curly, :Tuple, [eltype(contents[i]) for i = 1:length(NS)]...)
    return :(ProductNamedVectorSpace{NS, $BS, typeof(contents)}(contents))
end

end # module


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

mutable struct IoK8sApiCoreV1PersistentVolumeClaimStatus <: SwaggerModel
    accessModes::Any # spec type: Union{ Nothing, Vector{String} } # spec name: accessModes
    capacity::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: capacity
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PersistentVolumeClaimCondition} } # spec name: conditions
    phase::Any # spec type: Union{ Nothing, String } # spec name: phase

    function IoK8sApiCoreV1PersistentVolumeClaimStatus(;accessModes=nothing, capacity=nothing, conditions=nothing, phase=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimStatus, Symbol("accessModes"), accessModes)
        setfield!(o, Symbol("accessModes"), accessModes)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimStatus, Symbol("capacity"), capacity)
        setfield!(o, Symbol("capacity"), capacity)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimStatus, Symbol("phase"), phase)
        setfield!(o, Symbol("phase"), phase)
        o
    end
end # type IoK8sApiCoreV1PersistentVolumeClaimStatus

const _property_map_IoK8sApiCoreV1PersistentVolumeClaimStatus = Dict{Symbol,Symbol}(Symbol("accessModes")=>Symbol("accessModes"), Symbol("capacity")=>Symbol("capacity"), Symbol("conditions")=>Symbol("conditions"), Symbol("phase")=>Symbol("phase"))
const _property_types_IoK8sApiCoreV1PersistentVolumeClaimStatus = Dict{Symbol,String}(Symbol("accessModes")=>"Vector{String}", Symbol("capacity")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}", Symbol("conditions")=>"Vector{IoK8sApiCoreV1PersistentVolumeClaimCondition}", Symbol("phase")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1PersistentVolumeClaimStatus }) = collect(keys(_property_map_IoK8sApiCoreV1PersistentVolumeClaimStatus))
Swagger.property_type(::Type{ IoK8sApiCoreV1PersistentVolumeClaimStatus }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1PersistentVolumeClaimStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PersistentVolumeClaimStatus }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PersistentVolumeClaimStatus[property_name]

function check_required(o::IoK8sApiCoreV1PersistentVolumeClaimStatus)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PersistentVolumeClaimStatus }, name::Symbol, val)
end


text:
using MixedModelsPermutations
using Test

@testset "nonparametricbootstrap" begin
    include("nonparametricbootstrap.jl")
end

@testset "permutation" begin
    include("permutation.jl")
end

@testset "ols" begin
    include("ols.jl")
end


text:
module MPCCSolver

using Ipopt, JuMP, LinearAlgebra, MPCC, NLPModels, NLPModelsIpopt, NLPModelsJuMP, Stopping

include("ParamMPCC.jl")

export ParamMPCC

include("fillMstat.jl")

include("direct-nlp.jl")

export solveIpopt

include("regularization_methods/ThetaFct.jl")
include("regularization_methods/Relaxation.jl")
include("regularization_methods/RlxMPCC.jl")
include("regularization_methods/relaxation_algo.jl")

export mpcc_solve

#export RlxMPCC

end # module


text:
# Example for advanced 2D modeling:
# The receiver locations, sampling intervals for source and receiver,
# number of samples, as well as the wavelets can be different for each
# experiment.
#
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using PyCall
using PyPlot
using opesciSLIM.TimeModeling

## Set up model structure
n = (120,100)	# (x,y,z) or (x,z)
d = (10.,10.)
o = (0.,0.)

# Velocity [km/s]
v = ones(Float32,n) + 0.4f0
v0 = ones(Float32,n) + 0.4f0
v[:,Int(round(end/2)):end] = 4.0f0
rho = ones(Float32,n)

# Slowness squared [s^2/km^2]
m = (1./v).^2
m0 = (1./v0).^2
dm = vec(m - m0)

# Setup model structure
nsrc = 4
model = Model(n,d,o,m,rho;nb=100)	# number of pml points on each side as optional argument
model0 = Model(n,d,o,m0,rho;nb=100)

## Set up receiver geometry
nxrec = [120 160 130 115] # no. of receivers for the four experiments
nyrec = 1
xrec = Array{Any}(nsrc)
yrec = Array{Any}(nsrc)
zrec = Array{Any}(nsrc)

# Receiver x position. Each cell contains the coordinates for the corresponding experiment
xrec[1] = linspace(50.,1150.,nxrec[1])	# 101 receivers in x direction for experiment 1
xrec[2] = linspace(100.,1000.,nxrec[2])	# 50 receivers in x direction for experiment 2
xrec[3] = linspace(150.,900.,nxrec[3])	# 25 receivers in x direction for experiment 3
xrec[4] = linspace(250.,800.,nxrec[4])	# 13 receivers in x direction for experiment 4

# Receiver y position
yrec[1] = 0.	# 0 receivers in y direction for all experiments
yrec[2] = 0.
yrec[3] = 0.
yrec[4] = 0.

# Receiver z position
zrec[1] = linspace(100.,100.,nxrec[1])
zrec[2] = linspace(120.,120.,nxrec[2])
zrec[3] = linspace(80.,100.,nxrec[3])
zrec[4] = linspace(80.,140.,nxrec[4])

# receiver sampling and recording time
timeR = Array{Any}(nsrc)
dtR = Array{Any}(nsrc)
timeR[1] = 950.	# receiver recording length for each experiment
timeR[2] = 800.
timeR[3] = 1000.
timeR[4] = 880.
dtR[1] = 4.	# receiver sampling interval for each experiment
dtR[2] = 2.
dtR[3] = 3.
dtR[4] = 4.

# Setup receiver structure
recGeometry = Geometry(xrec,yrec,zrec;dt=dtR,t=timeR)

## Set up source geometry
xsrc = Array{Any}(nsrc)
ysrc = Array{Any}(nsrc)
zsrc = Array{Any}(nsrc)

# Source x position for each of the four experiments
xsrc[1] = [200, 400]	# inject two sources in experiment 1
xsrc[2] = 500
xsrc[3] = 800
xsrc[4] = 600

# Source y position
ysrc[1] = [0.,0.]
ysrc[2] = 0.
ysrc[3] = 0.
ysrc[4] = 0.

# Source z position
zsrc[1] = [50. 65.]	# inject two sources
zsrc[2] = 55.
zsrc[3] = 60.
zsrc[4] = 65.

# Source sampling and number of time steps
timeS = Array{Any}(nsrc)
dtS = Array{Any}(nsrc)
timeS[1] = 1000.	# source length for each experiment
timeS[2] = 900.
timeS[3] = 1200.
timeS[4] = 950.
dtS[1] = 3.	# source sampling interval for each experiment
dtS[2] = 4.
dtS[3] = 2.
dtS[4] = 2.

# Set up source structure
srcGeometry = Geometry(xsrc,ysrc,zsrc;dt=dtS,t=timeS)

# setup wavelets
f0 = [0.01, 0.015, 0.005, 0.02]	# source peak frequencies
wavelet = Array{Any}(nsrc)
wavelet[1] = source(timeS[1],dtS[1],f0[1])
wavelet[1] = [wavelet[1] wavelet[1]]	# first experiments has two sources
wavelet[2] = source(timeS[2],dtS[2],f0[2])
wavelet[3] = source(timeS[3],dtS[3],f0[3])
wavelet[4] = source(timeS[4],dtS[4],f0[4])

# Info structure for linear operators
ntComp = get_computational_nt(srcGeometry,recGeometry,model)	# no. of computational time steps
info = Info(prod(n),nsrc,ntComp)

###################################################################################################

# Setup operators
Pr = joProjection(info, recGeometry)
F = joModeling(info, model)
F0 = joModeling(info, model0)
Ps = joProjection(info, srcGeometry)
q = joData(srcGeometry, wavelet)

#Ffull = joModeling(info,model,srcGeometry,recGeometry)

# Nonlinear modeling
d = Pr*F*Ps'*q
qad = Ps*F'*Pr'*d

# Linearized modeling
J = joJacobian(Pr*F0*Ps',q)
dD = J*dm
rtm = J'*dD

# Plot results
figure(); imshow(d.data[1], vmin=-1e1, vmax=1e1)
figure(); imshow(reshape(rtm,model.n)', ColorMap("gray"))

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
module DependencyParser
  using TextAnalysis
  using WordTokenizers

  include("core/core.jl")
  include("units.jl")
  include("dependency_parsing/dependency_parsing.jl")
  include("preprocess_units/preprocess_units.jl")
  include("pipeline.jl")
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
########################################################################
#################### Helper Functions for Testing ######################
########################################################################

function resource_path()
    if splitdir(pwd())[2] == "TexTables"
        return joinpath("test", "resources")
    else
        return joinpath("resources")
    end
end

"""
Compares the contents of the file found at `fpath` to `fstring` line by
line, testing for equality.
"""
function compare_file(fpath::String, fstring::String)
    open(fpath) do f
        lines  = readlines(f)
        lines2 = split(fstring, "\n")

        for (l1, l2) in zip(lines, lines2)
            @test l1 == l2
        end
    end
end

"""
Generates the path to test table `i`
"""
function test_table(i)
    return joinpath(resource_path(), "test_table$i.txt")
end

function compare_file(i::Int, fstring::String)
    compare_file(test_table(i), fstring)
end

function export_table(i, fstring::String)
    open(test_table(i), "w") do f
        write(f, fstring)
    end
end

function export_table(t::TexTable; kwargs...)

    export_table(to_ascii(t; kwargs...), table_type="ascii")
    export_table(to_tex(t; kwargs...),   table_type="latex")

end

function next_test_table()
    files = readdir(resource_path())
    nums  = map(files) do file
        m = match(r"(?<=(test_table))\d*(?=(.txt))", file).match
        parse(Int, m)
    end
    return maximum(nums) + 1
end

function export_table(fstring::String; table_type="")
    next_num = next_test_table()
    export_table(next_num, fstring)
    println("Exported $table_type to $(test_table(next_num))")
    return next_num
end


text:
# Module to perform optional kabsch rotation to minimise distance between two structures before interpolating between them
module Kabsch

using DelimitedFiles

using LinearAlgebra
using Statistics

export kabsch_rotate

function translate_to_centroid(coord_matrix)
    # Normalises the molecular coordinates by centering them.
    center = mean(coord_matrix,dims=1)
    translated_geom = coord_matrix .- center
    return translated_geom
end

function optimal_rotation_matrix(CCmatrix)
    # Returns 3x3 matrix that can be applied to P to get Q
    # Previous implementation, not general (for cases of non-invertible matrices)
    # ORmatrix = sqrt(transpose(CCmatrix)*CCmatrix)*inv(CCmatrix)
    # Using SVD instead
    u,sing,vt = svd(CCmatrix)
    if sign(det(transpose(vt)*transpose(u))) == 1.0
        ORmatrix = transpose(vt)*transpose(u)
    elseif sign(det(transpose(vt)*transpose(u))) == -1.0
        mat = [1. 0. 0.;0. 1. 0.; 0. 0. -1.]
        ORmatrix = transpose(vt)*mat*transpose(u)
    else
        println("Error: Issue with SVD routine in finding optimal rotation matrix")
    end
    return ORmatrix
end

function xyz2matrix(input_xyz)
    # Imports .xyz file format as data frame, and converts it to an N x 3 Matrix
    raw_xyz=readdlm(input_xyz)
    natoms=size(raw_xyz)[1]
    just_coords=Array{Float64}(raw_xyz[2:natoms,2:4])
    return just_coords
end

function kabsch_rotate(Pgeom,Qgeom)
    # Returns the two geometries in xyz format, one of which has been translated, the
    # other translated and rotated
    # Pgeom = xyz2matrix(p_xyz)
    # Qgeom = xyz2matrix(q_xyz)

    normalisedP = (translate_to_centroid(Pgeom))
    normalisedQ = (translate_to_centroid(Qgeom))

    # Calculate cross covariance matrix
    xcov = transpose(normalisedP)*normalisedQ
    orot = optimal_rotation_matrix(xcov)

    num_atoms = size(normalisedP)[1]
    rotated = zeros(Float64,num_atoms,3)

    for i=1:num_atoms
        rotated[i,:] = orot*normalisedP[i,:]
    end

    RMSD_value = norm(rotated-normalisedQ)

    println("The RMSD between these two structures is ",RMSD_value)

    return rotated, normalisedQ
end

end

max_stars_repo_name: dasha-shchep/LIIC

max_issues_repo_name: dasha-shchep/LIIC

max_forks_repo_name: dasha-shchep/LIIC


text:
# Test derivatives w.r.t control points
#=
comm = MPI.COMM_WORLD
comm_world = MPI.MPI_COMM_WORLD
comm_self = MPI.COMM_SELF
my_rank = MPI.Comm_rank(comm)
comm_size = MPI.Comm_size(comm)
opts = PdePumiInterface.get_defaults()

# 2D mesh
opts["order"] = 1
opts["dimensions"] = 2
opts["use_DG"] = true
opts["operator_type"] = "SBPOmega"
opts["dmg_name"] = ".null"
opts["smb_name"] = "../src/mesh_files/gvortex1.smb"
opts["numBC"] = 2

opts["coloring_distance"] = 2 # 0 For CG Mesh 2 for DG Mesh
opts["jac_type"] = 2
opts["jac_method"] = 2
opts["run_type"] = 5

opts["numBC"] = 4
opts["BC1"] = [0]
opts["BC1_name"] = "FreeStreamBC"
opts["BC2"] = [1]
opts["BC2_name"] = "FreeStreamBC"
opts["BC3"] = [2]
opts["BC3_name"] = "FreeStreamBC"
opts["BC4"] = [3]
opts["BC4_name"] = "noPenetrationBC"

# Create PumiMesh and SBP objects
sbp, mesh, pmesh, Tsol, Tres, Tmsh, mesh_time = createMeshAndOperator(opts, 1)

# Initialize FFD
ndim = mesh.dim
order = [4,4,2]  # Order of B-splines in the 3 directions
nControlPts = [4,4,2]
offset = [0., 0., 0.5] # No offset in the X & Y direction
geom_faces = opts["BC4"]

# Create a mapping object using nonlinear mapping
map = initializeFFD(mesh, sbp, order, nControlPts, offset, false, geom_faces)

fill!(map.work, 0.0)

# Create seed vector
# - Get original wall coordinates
orig_wallCoords = FFD.getUniqueWallCoordsArray(mesh, geom_faces, false)
Xs_bar = ones(3, size(orig_wallCoords,2))
Xs_bar[3,:] = 0.0 # To accurately simulate a 2D mesh
cp_xyz_bar = zeros(map.cp_xyz)
evaldXdControlPointProduct(map, mesh, vec(Xs_bar))

filename = "./testvalues/evaldXdControlPointProduct.dat"
f = open(filename, "w")
for i = 1:length(map.work)
  println(f,map.work[i])
end
close(f)

=#

@testset "---Checking contractWithdGdB ---" begin
#=
  # opts = PdePumiInterface.get_defaults()
  # # 2D mesh
  # opts["order"] = 1
  # opts["dimensions"] = 2
  # opts["use_DG"] = true
  # opts["operator_type"] = "SBPOmega"
  # opts["dmg_name"] = ".null"
  # opts["smb_name"] = "../src/mesh_files/abc.smb"
  # opts["numBC"] = 2

  # opts["coloring_distance"] = 2 # 0 For CG Mesh 2 for DG Mesh
  # opts["jac_type"] = 2
  # opts["jac_method"] = 2
  # opts["run_type"] = 5

  # opts["numBC"] = 4
  # opts["BC1"] = [0]
  # opts["BC1_name"] = "FreeStreamBC"
  # opts["BC2"] = [1]
  # opts["BC2_name"] = "FreeStreamBC"
  # opts["BC3"] = [2]
  # opts["BC3_name"] = "FreeStreamBC"
  # opts["BC4"] = [3]
  # opts["BC4_name"] = "noPenetrationBC"

  # # Create PumiMesh and SBP objects
  # sbp, mesh, pmesh, Tsol, Tres, Tmsh, mesh_time = createMeshAndOperator(opts, 1)
  opts = Dict{ASCIIString, Any}()
  opts["order"] = 1
  opts["dimensions"] = 2
  opts["use_DG"] = true
  opts["Tsbp"] = Float64
  opts["Tmsh"] = Complex128
  opts["operator_type"] = "SBPOmega"
  opts["dmg_name"] = ".null"
  opts["smb_name"] = "../src/mesh_files/abc.smb"

  opts["coloring_distance"] = 2 # 0 For CG Mesh 2 for DG Mesh
  opts["jac_type"] = 2
  opts["jac_method"] = 2
  opts["run_type"] = 5

  opts["numBC"] = 4
  opts["BC1"] = [0]
  opts["BC1_name"] = "FreeStreamBC"
  opts["BC2"] = [1]
  opts["BC2_name"] = "FreeStreamBC"
  opts["BC3"] = [2]
  opts["BC3_name"] = "FreeStreamBC"
  opts["BC4"] = [3]
  opts["BC4_name"] = "noPenetrationBC"

  # Create PumiMesh and SBP objects
  dofpernode = 1
  ref_verts = [-1. 1 -1; -1 -1 1]
  Tsbp = opts["Tsbp"]
  Tmsh = opts["Tmsh"]
  sbp = getTriSBPOmega(degree=opts["order"], Tsbp=opts["Tsbp"])
  sbpface = TriFace{Tsbp}(opts["order"], sbp.cub, ref_verts.')
  topo = 0
  shape_type = 2
#  mesh = PumiMeshDG2{Tmsh}(

  mesh = PumiMeshDG2(Tmsh, sbp, opts, sbpface, dofpernode=dofpernode,
                     shape_type=shape_type)

#  mesh = PumiMeshDG2{Tmsh}(opts["dmg_name"], opts["smb_name"], opts["order"],
#                           sbp, opts, sbpface; dofpernode=dofpernode,
#                           coloring_distance=opts["coloring_distance"],
#                           shape_type=shape_type)
  # geometry faces to be embedded in FFD Box
  geom_faces = opts["BC4"]

  # Free Form deformation parameters
  ndim = mesh.dim
  order = [2,2,2]  # Order of B-splines in the 3 directions
  nControlPts = [2,2,2]

  # Create Mapping object
  map = PumiMapping{Tmsh}(ndim, order, nControlPts, mesh, full_geom=false, geom_faces=geom_faces)

  # Create knot vector
  calcKnot(map)

  # Create Bounding box
  offset = [0., 0., 0.5] # No offset in the X & Y direction
  box = PumiBoundingBox{Tmsh}(map, mesh, sbp, offset)

  # Control points
  controlPoint(map, box)
  # writeControlPointsVTS(map)

  calcParametricMappingNonlinear(map, box, mesh, geom_faces)

  pert = 1e-6 # Finite difference perturbation for tests
  xi = map.xi[1][:,2,1]
  dJdG = rand(Float64, 3)
  dJdG[3] = 0.0

  FFD.contractWithdGdB(map, xi, dJdG)
  cp_xyz_bar = zeros(map.cp_xyz) # Extract the values out of map.work
  for i = 1:size(map.work,4)
    for j = 1:size(map.work,3)
      for k = 1:size(map.work,2)
        cp_xyz_bar[1:3,k,j,i] = map.work[1:3,k,j,i]
      end
    end
  end

  # Check against finite difference
  # - Get original wall coordinates
  orig_wallCoords = FFD.getUniqueWallCoordsArray(mesh, geom_faces)
  multiplying_vec = zeros(length(orig_wallCoords))
  multiplying_vec[4:5] = dJdG[1:2]
  cp_jacobian = zeros(length(orig_wallCoords), length(map.cp_xyz))
  for i = 1:length(map.cp_xyz)
    map.cp_xyz[i] += pert
    vertices = evalSurface(map, mesh)
    commitToPumi(map, mesh, sbp, vertices, opts)
    new_wallCoords = FFD.getUniqueWallCoordsArray(mesh, geom_faces)
    cp_jacobian[:,i] = (vec(new_wallCoords) - vec(orig_wallCoords))/pert
    map.cp_xyz[i] -= pert
  end # End for i = 1:length(map.cp_xyz)

  prod_val = transpose(cp_jacobian)*multiplying_vec

  # Compute the error between reverse mode and finite difference
  error = vec(cp_xyz_bar) - prod_val
  for i = 1:length(error)
    @test isapprox( error[i], 0.0) atol=1e-10
  end
=#
end # End facts("---Checking contractWithdGdB ---")


text:
# Usage: l1emd.jl MAT_IN MAT_OUT
# Computes the Wasserstein distance between pairs of point sets stored in
# MAT files. The output is also stored in a MAT file. The input file must
# contain three variables:
#
# XPts - m-by-2 matrix of points,
# YPts - n-by-2 matrix of points,
# SegY - k-by-2 matrix of segments,
#
# where m, n, and k are positive integers. The output file contains:
#
# emd  - k-by-1 column vector of Wasserstein distance estimates.
#
# Date: January 8, 2018

using Clp
using SteinDiscrepancy: wassersteindiscrete
using MAT

# Command line arguments
matIn = ARGS[1];
matOut = ARGS[2];

# Read in MATLAB variables
file = matopen(matIn);
XPts = read(file, "XPts");
YPts = read(file, "YPts");
SegY = convert(Array{Int}, read(file, "SegY"));
close(file);

# Create solver ahead of time
solver = Clp.ClpSolver(LogLevel = 4);

# Compute EMD between XPts and each segment of YPts
nSeg = size(SegY, 1);
emd = Array{Float64}(nSeg);
for i in 1:nSeg
    (emd[i], ~, ~, ~) = wassersteindiscrete(
        xpoints = XPts,
        ypoints = YPts[colon(SegY[i, 1], SegY[i, 2]), :];
        solver = solver);
end

# Write output to MAT file
file = matopen(matOut, "w");
write(file, "emd", emd);
close(file);


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule GAP_pkg_cvec_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("GAP_pkg_cvec")
JLLWrappers.@generate_main_file("GAP_pkg_cvec", UUID("03de9d45-54dd-5e36-9810-6ca1d577e6b3"))
end  # module GAP_pkg_cvec_jll


text:
recurs_sym_type(ex::Any) =
  (ex==None || typeof(ex)==Symbol || length(ex.args)==1) ? eval(ex) : Expr(ex.head, ex.args[1], recurs_sym_type(ex.args[2]))
macro c(ret_type, func, arg_types, lib)
  local _arg_types = Expr(:tuple, [recurs_sym_type(a) for a in arg_types.args]...)
  local _ret_type = recurs_sym_type(ret_type)
  local _args_in = Any[ symbol(string('a',x)) for x in 1:length(_arg_types.args) ]
  local _lib = eval(lib)
  quote
    $(esc(func))($(_args_in...)) = ccall( ($(string(func)), $(Expr(:quote, _lib)) ), $_ret_type, $_arg_types, $(_args_in...) )
  end
end

macro ctypedef(fake_t,real_t)
  real_t = recurs_sym_type(real_t)
  quote
    typealias $fake_t $real_t
  end
end

# header: /usr/local/include/idas/idas_band.h
@ctypedef IDADlsDenseJacFnB Ptr{:Void}
@ctypedef IDADlsBandJacFnB Ptr{:Void}
@c Int32 IDADlsSetDenseJacFnB (Ptr{:None},:Int32,:IDADlsDenseJacFnB) shlib
@c Int32 IDADlsSetBandJacFnB (Ptr{:None},:Int32,:IDADlsBandJacFnB) shlib
@c Int32 IDABandB (Ptr{:None},:Int32,:Clong,:Clong,:Clong) shlib

# header: /usr/local/include/idas/idas_bbdpre.h
@ctypedef IDABBDLocalFnB Ptr{:Void}
@ctypedef IDABBDCommFnB Ptr{:Void}
@c Int32 IDABBDPrecInitB (Ptr{:None},:Int32,:Clong,:Clong,:Clong,:Clong,:Clong,:realtype,:IDABBDLocalFnB,:IDABBDCommFnB) shlib
@c Int32 IDABBDPrecReInitB (Ptr{:None},:Int32,:Clong,:Clong,:realtype) shlib

# header: /usr/local/include/idas/idas_dense.h
@c Int32 IDADenseB (Ptr{:None},:Int32,:Clong) shlib

# header: /usr/local/include/idas/idas_direct.h

# header: /usr/local/include/idas/idas.h
@ctypedef IDAQuadRhsFn Ptr{:Void}
@ctypedef IDASensResFn Ptr{:Void}
@ctypedef IDAQuadSensRhsFn Ptr{:Void}
@ctypedef IDAResFnB Ptr{:Void}
@ctypedef IDAResFnBS Ptr{:Void}
@ctypedef IDAQuadRhsFnB Ptr{:Void}
@ctypedef IDAQuadRhsFnBS Ptr{:Void}
@c Int32 IDASetQuadErrCon (Ptr{:None},:Int32) shlib
@c Int32 IDAQuadInit (Ptr{:None},:IDAQuadRhsFn,:N_Vector) shlib
@c Int32 IDAQuadReInit (Ptr{:None},:N_Vector) shlib
@c Int32 IDAQuadSStolerances (Ptr{:None},:realtype,:realtype) shlib
@c Int32 IDAQuadSVtolerances (Ptr{:None},:realtype,:N_Vector) shlib
@c Int32 IDASetSensDQMethod (Ptr{:None},:Int32,:realtype) shlib
@c Int32 IDASetSensParams (Ptr{:None},Ptr{:realtype},Ptr{:realtype},Ptr{:Int32}) shlib
@c Int32 IDASetSensErrCon (Ptr{:None},:Int32) shlib
@c Int32 IDASetSensMaxNonlinIters (Ptr{:None},:Int32) shlib
@c Int32 IDASensInit (Ptr{:None},:Int32,:Int32,:IDASensResFn,Ptr{:N_Vector},Ptr{:N_Vector}) shlib
@c Int32 IDASensReInit (Ptr{:None},:Int32,Ptr{:N_Vector},Ptr{:N_Vector}) shlib
@c Int32 IDASensToggleOff (Ptr{:None},) shlib
@c Int32 IDASensSStolerances (Ptr{:None},:realtype,Ptr{:realtype}) shlib
@c Int32 IDASensSVtolerances (Ptr{:None},:realtype,Ptr{:N_Vector}) shlib
@c Int32 IDASensEEtolerances (Ptr{:None},) shlib
@c Int32 IDAQuadSensInit (Ptr{:None},:IDAQuadSensRhsFn,Ptr{:N_Vector}) shlib
@c Int32 IDAQuadSensReInit (Ptr{:None},Ptr{:N_Vector}) shlib
@c Int32 IDAQuadSensSStolerances (Ptr{:None},:realtype,Ptr{:realtype}) shlib
@c Int32 IDAQuadSensSVtolerances (Ptr{:None},:realtype,Ptr{:N_Vector}) shlib
@c Int32 IDAQuadSensEEtolerances (Ptr{:None},) shlib
@c Int32 IDASetQuadSensErrCon (Ptr{:None},:Int32) shlib
@c Int32 IDAGetQuad (Ptr{:None},Ptr{:realtype},:N_Vector) shlib
@c Int32 IDAGetQuadDky (Ptr{:None},:realtype,:Int32,:N_Vector) shlib
@c Int32 IDAGetQuadNumRhsEvals (Ptr{:None},Ptr{:Clong}) shlib
@c Int32 IDAGetQuadNumErrTestFails (Ptr{:None},Ptr{:Clong}) shlib
@c Int32 IDAGetQuadErrWeights (Ptr{:None},:N_Vector) shlib
@c Int32 IDAGetQuadStats (Ptr{:None},Ptr{:Clong},Ptr{:Clong}) shlib
@c Int32 IDAGetSens (Ptr{:None},Ptr{:realtype},Ptr{:N_Vector}) shlib
@c Int32 IDAGetSens1 (Ptr{:None},Ptr{:realtype},:Int32,:N_Vector) shlib
@c Int32 IDAGetSensDky (Ptr{:None},:realtype,:Int32,Ptr{:N_Vector}) shlib
@c Int32 IDAGetSensDky1 (Ptr{:None},:realtype,:Int32,:Int32,:N_Vector) shlib
@c Int32 IDAGetSensConsistentIC (Ptr{:None},Ptr{:N_Vector},Ptr{:N_Vector}) shlib
@c Int32 IDAGetSensNumResEvals (Ptr{:None},Ptr{:Clong}) shlib
@c Int32 IDAGetNumResEvalsSens (Ptr{:None},Ptr{:Clong}) shlib
@c Int32 IDAGetSensNumErrTestFails (Ptr{:None},Ptr{:Clong}) shlib
@c Int32 IDAGetSensNumLinSolvSetups (Ptr{:None},Ptr{:Clong}) shlib
@c Int32 IDAGetSensErrWeights (Ptr{:None},:N_Vector_S) shlib
@c Int32 IDAGetSensStats (Ptr{:None},Ptr{:Clong},Ptr{:Clong},Ptr{:Clong},Ptr{:Clong}) shlib
@c Int32 IDAGetSensNumNonlinSolvIters (Ptr{:None},Ptr{:Clong}) shlib
@c Int32 IDAGetSensNumNonlinSolvConvFails (Ptr{:None},Ptr{:Clong}) shlib
@c Int32 IDAGetSensNonlinSolvStats (Ptr{:None},Ptr{:Clong},Ptr{:Clong}) shlib
@c Int32 IDAGetQuadSensNumRhsEvals (Ptr{:None},Ptr{:Clong}) shlib
@c Int32 IDAGetQuadSensNumErrTestFails (Ptr{:None},Ptr{:Clong}) shlib
@c Int32 IDAGetQuadSensErrWeights (Ptr{:None},Ptr{:N_Vector}) shlib
@c Int32 IDAGetQuadSensStats (Ptr{:None},Ptr{:Clong},Ptr{:Clong}) shlib
@c Int32 IDAGetQuadSens (Ptr{:None},Ptr{:realtype},Ptr{:N_Vector}) shlib
@c Int32 IDAGetQuadSens1 (Ptr{:None},Ptr{:realtype},:Int32,:N_Vector) shlib
@c Int32 IDAGetQuadSensDky (Ptr{:None},:realtype,:Int32,Ptr{:N_Vector}) shlib
@c Int32 IDAGetQuadSensDky1 (Ptr{:None},:realtype,:Int32,:Int32,:N_Vector) shlib
@c None IDAQuadFree (Ptr{:None},) shlib
@c None IDASensFree (Ptr{:None},) shlib
@c None IDAQuadSensFree (Ptr{:None},) shlib
@c Int32 IDAAdjInit (Ptr{:None},:Clong,:Int32) shlib
@c Int32 IDAAdjReInit (Ptr{:None},) shlib
@c None IDAAdjFree (Ptr{:None},) shlib
@c Int32 IDACreateB (Ptr{:None},Ptr{:Int32}) shlib
@c Int32 IDAInitB (Ptr{:None},:Int32,:IDAResFnB,:realtype,:N_Vector,:N_Vector) shlib
@c Int32 IDAInitBS (Ptr{:None},:Int32,:IDAResFnBS,:realtype,:N_Vector,:N_Vector) shlib
@c Int32 IDAReInitB (Ptr{:None},:Int32,:realtype,:N_Vector,:N_Vector) shlib
@c Int32 IDASStolerancesB (Ptr{:None},:Int32,:realtype,:realtype) shlib
@c Int32 IDASVtolerancesB (Ptr{:None},:Int32,:realtype,:N_Vector) shlib
@c Int32 IDAQuadInitB (Ptr{:None},:Int32,:IDAQuadRhsFnB,:N_Vector) shlib
@c Int32 IDAQuadInitBS (Ptr{:None},:Int32,:IDAQuadRhsFnBS,:N_Vector) shlib
@c Int32 IDAQuadReInitB (Ptr{:None},:Int32,:N_Vector) shlib
@c Int32 IDAQuadSStolerancesB (Ptr{:None},:Int32,:realtype,:realtype) shlib
@c Int32 IDAQuadSVtolerancesB (Ptr{:None},:Int32,:realtype,:N_Vector) shlib
@c Int32 IDACalcICB (Ptr{:None},:Int32,:realtype,:N_Vector,:N_Vector) shlib
@c Int32 IDACalcICBS (Ptr{:None},:Int32,:realtype,:N_Vector,:N_Vector,Ptr{:N_Vector},Ptr{:N_Vector}) shlib
@c Int32 IDASolveF (Ptr{:None},:realtype,Ptr{:realtype},:N_Vector,:N_Vector,:Int32,Ptr{:Int32}) shlib
@c Int32 IDASolveB (Ptr{:None},:realtype,:Int32) shlib
@c Int32 IDASetAdjNoSensi (Ptr{:None},) shlib
@c Int32 IDASetUserDataB (Ptr{:None},:Int32,Ptr{:None}) shlib
@c Int32 IDASetMaxOrdB (Ptr{:None},:Int32,:Int32) shlib
@c Int32 IDASetMaxNumStepsB (Ptr{:None},:Int32,:Clong) shlib
@c Int32 IDASetInitStepB (Ptr{:None},:Int32,:realtype) shlib
@c Int32 IDASetMaxStepB (Ptr{:None},:Int32,:realtype) shlib
@c Int32 IDASetSuppressAlgB (Ptr{:None},:Int32,:Int32) shlib
@c Int32 IDASetIdB (Ptr{:None},:Int32,:N_Vector) shlib
@c Int32 IDASetConstraintsB (Ptr{:None},:Int32,:N_Vector) shlib
@c Int32 IDASetQuadErrConB (Ptr{:None},:Int32,:Int32) shlib
@c Int32 IDAGetB (Ptr{:None},:Int32,Ptr{:realtype},:N_Vector,:N_Vector) shlib
@c Int32 IDAGetQuadB (Ptr{:None},:Int32,Ptr{:realtype},:N_Vector) shlib
@c Ptr{:None} IDAGetAdjIDABmem (Ptr{:None},:Int32) shlib
@c Int32 IDAGetConsistentICB (Ptr{:None},:Int32,:N_Vector,:N_Vector) shlib
@c Int32 IDAGetAdjY (Ptr{:None},:realtype,:N_Vector,:N_Vector) shlib
@ctypedef IDAadjCheckPointRec Void
@c Int32 IDAGetAdjCheckPointsInfo (Ptr{:None},Ptr{:IDAadjCheckPointRec}) shlib
@c Int32 IDAGetAdjDataPointHermite (Ptr{:None},:Int32,Ptr{:realtype},:N_Vector,:N_Vector) shlib
@c Int32 IDAGetAdjDataPointPolynomial (Ptr{:None},:Int32,Ptr{:realtype},Ptr{:Int32},:N_Vector) shlib
@c Int32 IDAGetAdjCurrentCheckPoint (Ptr{:None},Ptr{Ptr{:None}}) shlib

# header: /usr/local/include/idas/idas_impl.h
@ctypedef IDAadjMem Ptr{:Void}
@ctypedef IDABMem Ptr{:Void}
@ctypedef IDAAMMallocFn Ptr{:Void}
@ctypedef IDAAMFreeFn Ptr{:Void}
@ctypedef IDAAGetYFn Ptr{:Void}
@ctypedef IDAAStorePntFn Ptr{:Void}
@c Int32 IDASensResDQ (:Int32,:realtype,:N_Vector,:N_Vector,:N_Vector,Ptr{:N_Vector},Ptr{:N_Vector},Ptr{:N_Vector},Ptr{:None},:N_Vector,:N_Vector,:N_Vector) shlib

# header: /usr/local/include/idas/idas_spbcgs.h
@ctypedef IDASpilsPrecSetupFnB Ptr{:Void}
@ctypedef IDASpilsPrecSolveFnB Ptr{:Void}
@ctypedef IDASpilsJacTimesVecFnB Ptr{:Void}
@c Int32 IDASpilsSetGSTypeB (Ptr{:None},:Int32,:Int32) shlib
@c Int32 IDASpilsSetMaxRestartsB (Ptr{:None},:Int32,:Int32) shlib
@c Int32 IDASpilsSetEpsLinB (Ptr{:None},:Int32,:realtype) shlib
@c Int32 IDASpilsSetMaxlB (Ptr{:None},:Int32,:Int32) shlib
@c Int32 IDASpilsSetIncrementFactorB (Ptr{:None},:Int32,:realtype) shlib
@c Int32 IDASpilsSetPreconditionerB (Ptr{:None},:Int32,:IDASpilsPrecSetupFnB,:IDASpilsPrecSolveFnB) shlib
@c Int32 IDASpilsSetJacTimesVecFnB (Ptr{:None},:Int32,:IDASpilsJacTimesVecFnB) shlib
@c Int32 IDASpbcgB (Ptr{:None},:Int32,:Int32) shlib

# header: /usr/local/include/idas/idas_spgmr.h
@c Int32 IDASpgmrB (Ptr{:None},:Int32,:Int32) shlib

# header: /usr/local/include/idas/idas_spils.h

# header: /usr/local/include/idas/idas_sptfqmr.h
@c Int32 IDASptfqmrB (Ptr{:None},:Int32,:Int32) shlib


text:
"""
    identify_condition(df::AbstractDataFrame,cols::Vector{Symbol},codes::Vector{String})

produces a vector of Boolean `true` or `false` depending on whether columns `cols` contain
any of the string values in `codes`. This function is designed to easily identify comorbid
conditions in inpatient or outpatient records that are coded in ICD-9 or ICD-10 codes.
The values in `cols` are compared to the values (`code`) in `codes` only up to the length of
the `code`. For example, when detecting "diabetes", one may use "250.xx" in ICD-9,
which can be specified as identify_condition(df,[:diag1,:diag2,:diag3,:diag4,:diag5],["250"]).
All values that start with "250" in :diag1 - :diag4 (e.g., "25010") will match.
"""
function identify_condition(ip::AbstractDataFrame,vicd::Vector{Symbol},codes::Vector{String})

    retvec = zeros(Bool,nrow(ip))

    for i in 1:nrow(ip)
        for v in vicd
            if ismissing(ip[i,v])
                continue
            end

            vlen=length(ip[i,v])
            for c in codes
                clen = length(c)
                if vlen >= clen && c == ip[i,v][1:length(c)]
                    retvec[i] = true
                    break
                end
            end
            if retvec[i] == true
                break
            end
        end
    end
    return retvec
end

"""
    identify_condition2(df::AbstractDataFrame,cols::Vector{Symbol},codes::Vector{String})

produces a vector of Boolean `true` or `false` depending on whether columns `cols` contain
any of the string values in `codes`. This function is designed to easily identify comorbid
conditions in inpatient or outpatient records that are coded in ICD-9 or ICD-10 codes.
The values in `cols` are compared to the values in `codes` exactly as they are specified.
"""
function identify_condition2(ip::AbstractDataFrame,vicd::Vector{Symbol},codes::Vector{String})

    retvec = zeros(Bool,nrow(ip))

    for i in 1:nrow(ip)
        for v in vicd
            if ismissing(ip[i,v])
                continue
            end

            if ip[i,v] in codes
                retvec[i] = true
            end
        end
    end
    return retvec
end


text:
module Rendering
using Mustache
using Base.Markdown

"""
Defines the proper variables for the template home
"""
function template(template_directory)
    function (app, ctx)
        ctx.data[:template_dir] = abspath(template_directory)
        app(ctx)
    end
end

"""
# render

should render a file as a request

later on, I should deff cache this
"""
function render(ctx, filename, data=Dict())
    template_location = ctx.data[:template_dir]
    f = open(joinpath(template_location, filename))
    content = join(readlines(f), "")
    close(f)
    Mustache.render(content, data)
end

"""
# markdown
loads the filename, and runs it through a markdown application
"""
function markdown(ctx, file, data=Dict(); isfile=true)
    template_location = ctx.data[:template_dir]
    f = open(joinpath(template_location, filename))
    content = join(readlines(f), "")
    close(f)
    md = Mustache.render(content, data)
    md |> Markdown.parse |> Markdown.html
end

function json(ctx, data)
    mime!(ctx, "application/json")
    JSON.json(data)
end

end


text:
"""
ConstantsFS object for FAM with feature selection.
Augments Constants object with new priors:
W_star_prior, p_prior, and omega_prior.
Note that only one of p_prior or omega_prior will be used.
"""
mutable struct ConstantsFS{T <: ContinuousDistribution}
  W_star_prior::Gamma # W*_ik ~ Gamma(a_W, 1)
  p_prior::Beta # p_c ~ Beta(a_p, b_p)
  omega_prior::T # omage_each ~ Normal(m_omega , s_omega)
  constants::Constants
end

function ConstantsFS(c::Constants)
  # ws_prior = Gamma(1. / c.K, 1.)  # Old version.
  ws_prior = Gamma(10., 1.)
  p_prior = Beta(.1, .1)
  omega_prior = Normal(-5, 1.)

  return ConstantsFS(ws_prior, p_prior, omega_prior, c)
end

function printConstants(c::ConstantsFS, preprintln::Bool=true)
  if preprintln
    println("ConstantsFS:")
  end

  _fieldnames = collect(fieldnames(typeof(c)))

  for fname in filter(fn -> fn != :constants, _fieldnames)
    x = getfield(c, fname)
    T = typeof(x)
    if T <: Number
      println("$fname: $x")
    elseif T <: Vector
      N = length(x)
      for i in 1:N
        println("$(fname)[$i]: $(x[i])")
      end
    elseif T <: Dict
      for (k, v) in x
        println("$(fname)[$k]: $v")
      end
    else
      println("$fname: $x")
    end
  end

  printConstants(c.constants, true)
end


text:
using ..Microphysics_0M
using CLIMAParameters.Planet: Omega, e_int_i0, cv_l, cv_i, T_0

export Source,
    Gravity,
    RayleighSponge,
    Subsidence,
    GeostrophicForcing,
    Coriolis,
    RemovePrecipitation,
    CreateClouds

# kept for compatibility
# can be removed if no functions are using this
function atmos_source!(
    f::Function,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    f(atmos, source, state, diffusive, aux, t, direction)
end
function atmos_source!(
    ::Nothing,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
) end
# sources are applied additively
@generated function atmos_source!(
    stuple::Tuple,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    N = fieldcount(stuple)
    return quote
        Base.Cartesian.@nexprs $N i -> atmos_source!(
            stuple[i],
            atmos,
            source,
            state,
            diffusive,
            aux,
            t,
            direction,
        )
        return nothing
    end
end

abstract type Source end

struct Gravity <: Source end
function atmos_source!(
    ::Gravity,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    if atmos.ref_state isa HydrostaticState
        source.ρu -= (state.ρ - aux.ref_state.ρ) * aux.orientation.∇Φ
    else
        source.ρu -= state.ρ * aux.orientation.∇Φ
    end
end

struct Coriolis <: Source end
function atmos_source!(
    ::Coriolis,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    FT = eltype(state)
    _Omega::FT = Omega(atmos.param_set)
    # note: this assumes a SphericalOrientation
    source.ρu -= SVector(0, 0, 2 * _Omega) × state.ρu
end

struct Subsidence{FT} <: Source
    D::FT
end

function atmos_source!(
    subsidence::Subsidence,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    ρ = state.ρ
    z = altitude(atmos, aux)
    w_sub = subsidence_velocity(subsidence, z)
    k̂ = vertical_unit_vector(atmos, aux)

    source.ρe -= ρ * w_sub * dot(k̂, diffusive.∇h_tot)
    source.moisture.ρq_tot -= ρ * w_sub * dot(k̂, diffusive.moisture.∇q_tot)
end

subsidence_velocity(subsidence::Subsidence{FT}, z::FT) where {FT} =
    -subsidence.D * z

struct GeostrophicForcing{FT} <: Source
    f_coriolis::FT
    u_geostrophic::FT
    v_geostrophic::FT
end
function atmos_source!(
    s::GeostrophicForcing,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    u_geo = SVector(s.u_geostrophic, s.v_geostrophic, 0)
    ẑ = vertical_unit_vector(atmos, aux)
    fkvector = s.f_coriolis * ẑ
    source.ρu -= fkvector × (state.ρu .- state.ρ * u_geo)
end

"""
    RayleighSponge{FT} <: Source

Rayleigh Damping (Linear Relaxation) for top wall momentum components
Assumes laterally periodic boundary conditions for LES flows. Momentum components
are relaxed to reference values (zero velocities) at the top boundary.
"""
struct RayleighSponge{FT} <: Source
    "Maximum domain altitude (m)"
    z_max::FT
    "Altitude at with sponge starts (m)"
    z_sponge::FT
    "Sponge Strength 0 ⩽ α_max ⩽ 1"
    α_max::FT
    "Relaxation velocity components"
    u_relaxation::SVector{3, FT}
    "Sponge exponent"
    γ::FT
end
function atmos_source!(
    s::RayleighSponge,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    z = altitude(atmos, aux)
    if z >= s.z_sponge
        r = (z - s.z_sponge) / (s.z_max - s.z_sponge)
        β_sponge = s.α_max * sinpi(r / 2)^s.γ
        source.ρu -= β_sponge * (state.ρu .- state.ρ * s.u_relaxation)
    end
end

"""
    CreateClouds{FT} <: Source

A source/sink to `q_liq` and `q_ice` implemented as a relaxation towards
equilibrium in the Microphysics module.
The default relaxation timescales are defined in CLIMAParameters.jl.
"""
struct CreateClouds <: Source end
function atmos_source!(
    ::CreateClouds,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    # get current temperature and phase partition
    FT = eltype(state)
    ts = recover_thermo_state(atmos, state, aux)
    q = PhasePartition(ts)
    T = air_temperature(ts)

    # phase partition corresponding to the current T and q.tot
    # (this is not the same as phase partition from saturation adjustment)
    ts_eq = PhaseEquil_ρTq(atmos.param_set, state.ρ, T, q.tot)
    q_eq = PhasePartition(ts_eq)

    # cloud condensate as relaxation source terms
    S_q_liq = conv_q_vap_to_q_liq_ice(atmos.param_set.microphys.liq, q_eq, q)
    S_q_ice = conv_q_vap_to_q_liq_ice(atmos.param_set.microphys.ice, q_eq, q)

    source.moisture.ρq_liq += state.ρ * S_q_liq
    source.moisture.ρq_ice += state.ρ * S_q_ice
end

"""
    RemovePrecipitation{FT} <: Source

A sink to `q_tot` when cloud condensate is exceeding a threshold.
The threshold is defined either in terms of condensate or supersaturation.
The removal rate is implemented as a relaxation term
in the Microphysics_0M module.
The default thresholds and timescale are defined in CLIMAParameters.jl.
"""
struct RemovePrecipitation <: Source
    " Set to true if using qc based threshold"
    use_qc_thr::Bool
end
function atmos_source!(
    s::RemovePrecipitation,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    FT = eltype(state)
    ts = recover_thermo_state(atmos, state, aux)
    if has_condensate(ts)

        q = PhasePartition(ts)
        T::FT = air_temperature(ts)
        λ::FT = liquid_fraction(ts)

        _e_int_i0::FT = e_int_i0(atmos.param_set)
        _cv_l::FT = cv_l(atmos.param_set)
        _cv_i::FT = cv_i(atmos.param_set)
        _T_0::FT = T_0(atmos.param_set)

        S_qt::FT = 0
        if s.use_qc_thr
            S_qt = remove_precipitation(atmos.param_set, q)
        else
            q_vap_sat = q_vap_saturation(ts)
            S_qt = remove_precipitation(atmos.param_set, q, q_vap_sat)
        end

        source.moisture.ρq_tot += state.ρ * S_qt

        source.ρ += state.ρ * S_qt

        source.ρe +=
            (
                λ * _cv_l * (T - _T_0) +
                (1 - λ) * (_cv_i * (T - _T_0) - _e_int_i0) +
                gravitational_potential(atmos.orientation, aux)
            ) *
            state.ρ *
            S_qt
    end
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
struct Lens{T, Fields, Field}
end

@generated function runlens(main :: T, new_field, lens::Type{Lens{T, Fields, Field}}) where {T, Fields, Field}
    quote
        $T($([field !== Field ? :(main.$field) : :new_field for field in Fields]...))
    end
end


text:
using DiffPrecTest
using Statistics, StatsBase, LinearAlgebra
using SparseArrays
using ProximalBase, CoordinateDescent, CovSel
using Random, Distributions
using JLD

@show gethostname()

rep   = parse(Int,ARGS[1])
ip    = parse(Int,ARGS[2])
dir = ARGS[3]

pArr = [50, 100, 150]
n = 300

p = pArr[ip]
Random.seed!(54298)

# generate model
Ωx = Matrix{Float64}(I, p, p)
for k=1:div(p, 10)
    for j=10*(k-1)+2:10*(k-1)+10
        Ωx[10*(k-1)+1, j] = 0.5
        Ωx[j, 10*(k-1)+1] = 0.5
    end
end
de = abs(minimum(eigvals(Ωx)))+0.05
Ωx = (Ωx + de*I)/(1+de)

# generate Delta
d = Vector{Float64}(undef, p)
rand!(Uniform(0.5, 2.5), d)
d .= sqrt.(d)
D = Diagonal(d)
Σ = inv(Symmetric(D * Ωx * D))

dist_X = MvNormal(convert(Matrix, Σ))

# generate data
Random.seed!(12346 + rep)
X = rand(dist_X, n)'
Y = rand(dist_X, n)'

Sx = Symmetric( X'X / n )
Sy = Symmetric( Y'Y / n )

@time boot_res, eS = bootstrap(X, Y)

@save "$(dir)/res_$(ip)_$(rep).jld" boot_res eS


text:
## Scatter API
#   - Scatter:
#     - scatter(op, src, idx)
#     - scatter!(op, dst, src, idx)
#   - Scatter destination backpropagation
#     - ∇scatter!_dst
#   - Scatter source backpropagation
#     - ∇scatter_src
#     - ∇scatter!_src
#

typelength(::Type{<:Number}) = 1
typelength(::Type{<:NTuple{M}}) where M = M
typelength(::Type{CartesianIndex{M}}) where M = M

"""
Performs dimensional consistency checks and return the
dimensionality of the scattered objects.
"""
function scatter_dims(X::AbstractArray{Tx,Nx},
                     Y::AbstractArray{Ty,Ny},
                     idx::AbstractArray{Tidx,Nidx}) where {Tx,Ty,Tidx,Nx,Ny,Nidx}
    M = typelength(Tidx)
    dims = scatter_dims(Nx, Ny, M, Nidx)
    size(X)[1:dims] == size(Y)[1:dims] || throw(ArgumentError("Incompatible input shapes."))
    size(Y)[dims+1:end] == size(idx) || throw(ArgumentError("Incompatible input shapes."))
    return dims
end

function scatter_dims(Nx, Ny, M, Nidx)
    @assert Nx - M == Ny - Nidx "Incompatible input shapes of (dst, src, idx) = ($Nx, $Ny, $Nidx)."
    dims = Nx - M
    dims < 0 && throw(ArgumentError("dims must be non-negative but got dims=$dims."))
    return dims
end

_view(X, colons, k) = view(X, colons..., k...)
_view(X, colons, k::Union{Integer, CartesianIndex}) = view(X, colons..., k)

"""
    scatter!(op, dst, src, idx)

Scatter operation, which scatters data in `src` and assigns to `dst` according to `idx`.
A binary reduction operator `op` is applied during the scatter.
For each index `k` in `idx`, accumulates values in `dst` according to

    dst[:, ..., idx[k]...] = (op).(dst[:, ..., idx[k]...], src[:, ..., k...])

# Arguments

- `op`: Operations to be applied on `dst` and `src`, e.g. `+`, `-`, `*`, `/`, `max`, `min` and `mean`.
- `dst`: The destination for `src` to aggregate to. This argument will be mutated.
- `src`: The source data for aggregating.
- `idx`: The mapping for aggregation from source (index) to destination (value).
         The `idx` array can contain either integers or tuples.
"""
function scatter!(op, dst::AbstractArray, src::AbstractArray, idx::AbstractArray)
    dims = scatter_dims(dst, src, idx)
    colons = Base.ntuple(_->Colon(), dims)
    for k in CartesianIndices(idx)
        dst_v = _view(dst, colons, idx[k])
        src_v = _view(src, colons, k)
        dst_v .= (op).(dst_v, src_v)
    end
    dst
end

function scatter!(op::typeof(mean), dst::AbstractArray, src::AbstractArray, idx::AbstractArray)
    Ns = scatter!(+, zero(dst), one.(src), idx)
    dst_ = scatter!(+, zero(dst), src, idx)
    dst .+= safe_div.(dst_, Ns)
    return dst
end

"""
    scatter(op, src, idx; [init])

Scatter operation allocating a destination array `dst` and
calling `scatter!(op, dst, src, idx)` on it.

If `init` is provided, it is used to initialize the content of `dst`.
Otherwise, the init values is inferred from the reduction operator `op`
for some common operators (e.g. `init = 0` for `op = +`).

See [`scatter!`](@ref) for the details.
"""
function scatter(op,
                src::AbstractArray{Tsrc,Nsrc},
                idx::AbstractArray{Tidx,Nidx};
                init = nothing) where {Tsrc,Tidx,Nsrc,Nidx}

    dims = Nsrc - Nidx
    dstsize = (size(src)[1:dims]..., maximum_dims(idx)...)
    dst = similar(src, Tsrc, dstsize)
    xinit = isnothing(init) ? scatter_empty(op, Tsrc) : init
    fill!(dst, xinit)
    scatter!(op, dst, src, idx)
end

scatter_empty(op, T) = Base.reduce_empty(op, T)
scatter_empty(op::typeof(-), T) = zero(T)
scatter_empty(op::typeof(/), T) = one(T)
scatter_empty(op::typeof(min), T) = typemax(T)
scatter_empty(op::typeof(max), T) = typemin(T)
scatter_empty(op::typeof(mean), T) = zero(T)

## Gradients

∇scatter!_src(op, Δ, dst, src, idx) = ∇scatter_src(op, Δ, dst, src, idx)
∇scatter!_dst(op, Δ, dst, y) = Δ

∇scatter!_dst(op::Union{typeof(max),typeof(min)}, Δ, dst_old, dst) =
    (dst_old .== op.(dst_old, dst)) .* Δ

modify_src(::typeof(+), X) = X
modify_src(::typeof(-), X) = -X
modify_src(::typeof(*), X, Y) = X
modify_src(::typeof(/), X, Y) = .-X ./ Y.^2

∇scatter_src(op::Union{typeof(+),typeof(-)}, Δ, dst, src, idx) = modify_src(op, gather(Δ, idx))

∇scatter!_src(op::Union{typeof(*),typeof(/)}, Δ, dst, src, idx) =
    gather(dst, idx) .* ∇scatter_src(op, Δ, dst, src, idx)

function ∇scatter_src(op::Union{typeof(*),typeof(/)}, Δ, dst,
                      src::AbstractArray{Tsrc,Nsrc},
                      idx::AbstractArray{Tidx,Nidx}) where {Tsrc,Tidx,Nsrc,Nidx}
    dims = Nsrc - Nidx
    Δsrc = modify_src(op, gather(Δ, idx), src)
    rev_idx = reverse_indices(idx)
    ax = CartesianIndices(axes(src)[1:dims])
    for k in CartesianIndices(idx)
        inds = filter(x -> x != k, rev_idx[idx[k]])
        for i in ax
            Δsrc[i, k] = op(Δsrc[i, k], prod(j -> src[i, j], inds))
        end
    end
    Δsrc
end

∇scatter_src(::Union{typeof(max),typeof(min)}, Δ, dst, src, idx) = (src .== gather(dst, idx)) .* gather(Δ, idx)

function ∇scatter_src(::typeof(mean), Δ, dst,
                    src::AbstractArray{Tsrc,Nsrc},
                    idx::AbstractArray{Tidx,Nidx}) where {Tsrc,Tidx,Nsrc,Nidx}
    dims = Nsrc - Nidx
    divide_by_counts!(gather(Δ, idx), idx, dims)
end

function rrule(::typeof(scatter!), op, dst::AbstractArray, src::AbstractArray, idx::AbstractArray)
    dst_old = copy(dst)
    scatter!(op, dst, src, idx)
    scatter!_pullback(Δ) = (NoTangent(), NoTangent(), ∇scatter!_dst(op, Δ, dst_old, dst), ∇scatter!_src(op, Δ, dst, src, idx), NoTangent())
    dst, scatter!_pullback
end

function rrule(::typeof(scatter), op, src::AbstractArray, idx::AbstractArray)
    y = scatter(op, src, idx)
    scatter_pullback(Δ) = (NoTangent(), NoTangent(), ∇scatter_src(op, Δ, y, src, idx), NoTangent())
    y, scatter_pullback
end


text:
using Rock
peers = [(Rock.Hostname("localhost"), Rock.Port(8000)),
         (Rock.Hostname("localhost"), Rock.Port(8001)),
         (Rock.Hostname("localhost"), Rock.Port(8002))]
type C <: Rock.ExternalCriticalCommand
    a
end
function t(px, d, c)
@show c
end

i = parse(Int64, ARGS[1])
if length(ARGS) >1
    Rock.command(Rock.client("localhost", peers[i][2]), C(ARGS[2]))
else
    c = Rock.Context(peers[i][2], peers, t, AbstractString("sqlite$i.db"))
    n  =Rock.Instance(c,i)
    Rock.start(n)
end

max_stars_repo_name: JeremyRubin/Rock


text:
#### Redefinition of push! function for stacks
import Base.push!
## Push functions for CompressedStack : push!, push_compressed!, push_explicit!

# Method to compress a block into a signature
function compress!{T}(stack::CompressedStack,
  sign::Nullable{Signature{T}}, block::Block{T})

  if isempty(block)
    sign = Nullable{Signature}()
  else
    if isnull(sign)
      sign = Nullable(Signature(block))

    else
      get(sign).last = block[end].last
    end
  end
  stack.compressed = sign
end

# Function push_explicit to push from input into the explicit blocks
function push_explicit!{T,D}(stack::CompressedStack{T,D}, elt::D)
  data = Data(elt, stack.index)
  if isempty(stack.first_explicit)
    push!(stack.first_explicit, data)
    sign = Signature(stack.index, get(stack.context), stack.pos)
    stack.first_sign = Nullable(sign)
  else
    head = top(stack)
    start_block = head - mod(head - 1, stack.space)
    if stack.index - start_block < stack.space
      push!(stack.first_explicit,data)
      update_signature!(get(stack.first_sign), stack.index)
    else
      if stack.depth == 1 && !isempty(stack.second_explicit)
        if isnull(stack.compressed)
          stack.compressed = stack.second_sign
        else
          update_signature!(get(stack.compressed), get(stack.second_sign).last)
        end
      end
      stack.second_sign = stack.first_sign
      sign = Signature(stack.index, get(stack.context), stack.pos)
      stack.first_sign = Nullable(sign)
      stack.second_explicit = stack.first_explicit
      stack.first_explicit = [data]
    end
  end
end

# Function to push (possibly fully) compressed index of new data from input
function push_compressed!{T,D}(stack::CompressedStack{T,D}, lvl::Int)
  dist_subblock = stack.space^(stack.depth - lvl)
  dist_block = dist_subblock * stack.space

  if isempty(stack, lvl)
    sign = Signature(stack.index, get(stack.context),stack.pos)
    push!(stack.first_partial[lvl], sign)
  else
    head = top(stack.first_partial[lvl])
    start_block = head - mod(head - 1, dist_block)
    # distance of the new index and current block
    δ = stack.index - start_block + 1
    if δ <= dist_block
      # Distance with the current subblock
      start_subblock = head - mod(head - 1, dist_subblock)
      η = stack.index - start_subblock + 1
      # compress new element in the top of the current Block
      if η <= dist_subblock
        subblock = length(stack.first_partial[lvl])
        update_top!(stack.first_partial[lvl], subblock, stack.index)
      else
        sign = Signature(stack.index, get(stack.context), stack.pos)
        push!(stack.first_partial[lvl], sign)
      end
    else
      if lvl == 1
        compress!(stack, stack.compressed, stack.second_partial[1])
      end
      stack.second_partial[lvl] = stack.first_partial[lvl]
      sign = Signature(stack.index, get(stack.context), stack.pos)
      stack.first_partial[lvl] = [sign]
    end
  end
end

# Function push! that push the data in explicit and index in partial/compressed
function push!{T,D}(stack::CompressedStack{T,D}, elt::D)
  stack.push_action(stack, elt)
  # update the explicit Blocks, with possibly shifting first to second
  push_explicit!(stack, elt)
  # update the compressed Blocks at each levels (including the fully compressed)
  for lvl in 1:(stack.depth - 1)
    push_compressed!(stack, lvl)
  end
end

#####################
## push! for NormalStack
function push!{D}(stack::NormalStack, elt::D)
  stack.push_action(stack, elt)
  data = Data(elt, stack.index)
  push!(stack.data, data)
end


text:
const PeriodicMatrixFunction{T} = Union{Constant{<:AbstractMatrix{T}},
                                        FourierSeries{<:AbstractMatrix{T}}}

# Gauges:
# :n => no particular gauge, left and right fixed points completely generic
# :l => left gauge: identity left fixed point, generic right fixed point
# :L => left canonical form: identity left fixed point, diagonal right fixed point
# :r => right gauge
# :R => right canonical form:
# :s => symmetric: left and right fixed point identical
# :S => symmetric: left and right fixed point identical and diagonal

mutable struct InfiniteCMPS{T<:PeriodicMatrixFunction,N} <: AbstractCMPS{T,N}
    Q::T
    Rs::NTuple{N,T}
    gauge::Symbol
    function InfiniteCMPS(Q::T, Rs::NTuple{N,T}; gauge::Symbol = :n) where {T,N}
        for R in Rs
            domain(R) == domain(Q) || throw(DomainMismatch())
        end
        Q0 = Q[0]
        size(Q0, 1) == size(Q0, 2) || throw(DimensionMismatch())
        for R in Rs
            size(R[0]) == size(Q0) || throw(DimensionMismatch())
        end
        return new{T,N}(Q, Rs, gauge)
    end
end
InfiniteCMPS(Q::T, R::T; kwargs...) where T = InfiniteCMPS(Q, (R,); kwargs...)

domain(::InfiniteCMPS) = (-Inf, +Inf)
period(Ψ::InfiniteCMPS) = period(Ψ.Q)

Base.iterate(Ψ::InfiniteCMPS, args...) = iterate((Ψ.Q, Ψ.Rs), args...)

Base.copy(Ψ::InfiniteCMPS) = InfiniteCMPS(copy(Ψ.Q), map(copy, Ψ.Rs); gauge = Ψ.gauge)

virtualdim(Ψ::InfiniteCMPS) = size(Ψ.Q[0], 1)

const UniformCMPS = InfiniteCMPS{<:Constant{<:AbstractMatrix}}
const FourierCMPS = InfiniteCMPS{<:FourierSeries{<:AbstractMatrix}}


text:
qm7b_init() = register(DataDep(
    "QM7b",
    """
    The QM7b dataset from the `"MoleculeNet: A Benchmark for Molecular
    Machine Learning" <https://arxiv.org/abs/1703.00564>`_ paper, consisting of
    7,211 molecules with 14 regression targets.
    """,
    "http://deepchem.io.s3-website-us-west-1.amazonaws.com/datasets/qm7b.mat",
    "e2a9d670d86eba769fa7b5eadeb592184067d2ec12468b1a220bfc38502dda61";
    post_fetch_method=preprocess_qm7b,
))

function preprocess_qm7b(local_path)
    vars = matread(local_path)
    names = vars["names"]
    X = vars["X"]
    T = Matrix{Float32}(vars["T"])
    raw = Dict(:names=>names, :X=>X, :T=>T)

    rawfile = replace(local_path, "qm7b.mat"=>"qm7b.raw.jld2")
    @save rawfile raw
end

struct QM7b <:Dataset end

function rawdata(::QM7b)
    file = datadep"QM7b/qm7b.raw.jld2"
    @load file raw
    raw
end