text:
using Random, LinearAlgebra, SpecialFunctions, PyPlot, KernelDensity, PyCall
stats = pyimport("scipy.stats")

"Generate (m, n) matrix uniformly distributed over Stiefel manifold."
function rando(m::Integer, n::Integer)
    M = randn((m, n))
    Q = Matrix(qr(M).Q)
end

"Generate matrix with exponentially-distributed singular values and singular vectors ~ O(n)."
function randA(n::Integer, scale::Real)
    s = randexp(n)*scale
    V = rando(n, n)
    U = rando(n, n)
    A = U*Diagonal(s)*transpose(V)
    return A, s, U, V
end

"Generate random symmetric matrix with exponentially-distributed eigenvalues and eigenvectors ~ O(n)."
function randS(n::Integer, scale::Real)
    s = randexp(n)*scale
    U = rando(n, n)
    A = U*Diagonal(s)*transpose(U)
    return A, s, U
end

"Generate solution vector from x0 + V N(0, I) + Y G N(0, I)."
function randx(x0::Array{<:Real, 1}, V::Array{<:Real, 2}, Y::Array{<:Real, 2}, G::Array{<:Real, 2})
    v = randn(size(V)[2])
    y = randn(size(Y)[2])
    x = x0 + V*v + Y*G*y
end

"Generate solution vector from x0 + N(0, I)."
function randx(x0::Array{<:Real, 1})
    x = x0 + randn(size(x0))
end

"Generate solution vector from x0 + L N(0, I)."
function randx(x0::Array{<:Real, 1}, L::Array{<:Real, 2})
    x = x0 + L*randn(size(L)[2])
end

"Naive implementation of projection method."
function projection(x0::Array{<:Real, 1}, A::Array{<:Real, 2}, b::Array{<:Real, 1}, V::Array{<:Real, 2}, W::Array{<:Real, 2})
    x = x0 + V*inv(transpose(W)*A*V)*transpose(W)*(b - A*x0)
end

"Krylov subspace K_{n}(A, v)."
function Krylov(A::Array{<:Real, 2}, v::Array{<:Real, 1}, n::Integer)
    V = Array{Float64, 2}(undef, size(v)[1], n)
    V[:, 1] = v/norm(v)
    for i=2:n
        V[:, i] = A*V[:, i-1]
        for j=1:i-1
            V[:, i] = V[:, i] - dot(V[:, i], V[:, j])*V[:, j]
        end
        V[:, i] = V[:, i]/norm(V[:, i])
    end
    return V
end

"PDF of chi-squared distribution."
function chi2PDF(x::Real, n::Integer)
    return x^(n/2-1)*exp(-x/2)/(2^(n/2)*gamma(n/2))
end

"PDF of F-distribution."
function FPDF(x::Real, d1::Integer, d2::Integer)
    return (d1/d2)^(d1/2)*x^(d1/2-1)*(1+d1*x/d2)^(-(d1+d2)/2)/beta(d1/2, d2/2)
end

"Compute orthonormal basis for Null(A), where A is full-rank MxN matrix with N>M."
function NullSpace(A::Array{<:Real, 2})
    svdA = svd(A, full=true)
    return svdA.V[:, (size(A)[1]+1):size(A)[2]]
end

"Approximation of int dx|F1-F2| with central Riemann sum."
function Integrate(F1::Array{<:Real, 1}, F2::Array{<:Real, 1}, h::Real)
    F = abs.(F1 - F2)
    s = h*(sum(F[2:(size(F)[1]-1)]) + (F[1] + F[size(F)[1]])/2)
end

"Symmetric matices, separate estimation of the scale for each run, mean of inverse gamma."
function SymSamplingIn(N_rounds::Integer, scale::Real, m::Integer,  N::Integer=100, N_samples::Integer=5000; randomization="A")
    A, s, U = randS(N, scale)
    L = U*Diagonal(s.^(-1/2))
    inv_A = U*Diagonal(s.^(-1))*transpose(U)
    Statistics_inv = zeros(N_samples)
    Statistics_proj = zeros(N_samples)
    if randomization == "A"
        draw = (n=N, l=L) -> randx(zeros(n), l)
    else
        draw = (n=N) -> randx(zeros(n))
    end

    for i=1:N_samples
        x = draw()
        b = A*x
        V = Krylov(A, b, m)
        y = projection(zeros(N), A, b, V, V)
        sigma = inv_A - V*inv(transpose(V)*A*V)*transpose(V)
        statistics = transpose(x - y)*pinv(sigma)*(x - y)
        Statistics_inv[i] = statistics

        alpha = 0
        for i=1:N_rounds
            z = draw()
            b = A*z
            V = Krylov(A, b, m)
            P = I - V*inv(transpose(V)*A*V)*transpose(V)*A
            sample = P*z
            alpha = alpha + norm(sample)^2
        end
        alpha = alpha/(N_rounds*(N-m)+2)

        statistics = (norm(x - y)^2)/alpha
        Statistics_proj[i] = statistics
    end
    x = LinRange(0, 500, 5000)
    kernel_inv = stats.gaussian_kde(Statistics_inv.*(Statistics_inv.<500).*(Statistics_inv.>0))
    kernel_proj = stats.gaussian_kde(Statistics_proj.*(Statistics_proj.<500))
    stat_inv = kernel_inv(x)
    stat_proj = kernel_proj(x)
    return stat_inv, stat_proj
end

"Symmetric matices, the single estimation of the scale for all run, mean of inverse gamma."
function SymSamplingOut(N_rounds::Integer, scale::Real, m::Integer,  N::Integer=100, N_samples::Integer=5000; randomization="A")
    A, s, U = randS(N, scale)
    L = U*Diagonal(s.^(-1/2))
    inv_A = U*Diagonal(s.^(-1))*transpose(U)
    Statistics_inv = zeros(N_samples)
    Statistics_proj = zeros(N_samples)
    if randomization == "A"
        draw = (n=N, l=L) -> randx(zeros(n), l)
    else
        draw = (n=N) -> randx(zeros(n))
    end

    alpha = 0
    for i=1:N_rounds
        z = draw()
        b = A*z
        V = Krylov(A, b, m)
        P = I - V*inv(transpose(V)*A*V)*transpose(V)*A
        sample = P*z
        alpha = alpha + norm(sample)^2
    end
    alpha = alpha/(N_rounds*(N-m)+2)

    for i=1:N_samples
        x = draw()
        b = A*x
        V = Krylov(A, b, m)
        y = projection(zeros(N), A, b, V, V)
        sigma = inv_A - V*inv(transpose(V)*A*V)*transpose(V)
        statistics = transpose(x - y)*pinv(sigma)*(x - y)
        Statistics_inv[i] = statistics

        statistics = (norm(x - y)^2)/alpha
        Statistics_proj[i] = statistics
    end
    x = LinRange(0, 500, 5000)
    kernel_inv = stats.gaussian_kde(Statistics_inv.*(Statistics_inv.<500).*(Statistics_inv.>0))
    kernel_proj = stats.gaussian_kde(Statistics_proj.*(Statistics_proj.<500))
    stat_inv = kernel_inv(x)
    stat_proj = kernel_proj(x)
    return stat_inv, stat_proj
end

"Symmetric matices, separate estimation of the scale for each run, multivariate Student."
function HSymSamplingIn(N_rounds::Integer, scale::Real, m::Integer,  N::Integer=100, N_samples::Integer=5000; randomization="A")
    A, s, U = randS(N, scale)
    L = U*Diagonal(s.^(-1/2))
    inv_A = U*Diagonal(s.^(-1))*transpose(U)
    Statistics_proj = zeros(N_samples)
    if randomization == "A"
        draw = (n=N, l=L) -> randx(zeros(n), l)
    else
        draw = (n=N) -> randx(zeros(n))
    end

    for i=1:N_samples
        alpha = 0
        for i=1:N_rounds
            z = draw()
            b = A*z
            V = Krylov(A, b, m)
            P = I - V*inv(transpose(V)*A*V)*transpose(V)*A
            sample = P*z
            alpha = alpha + norm(sample)^2
        end
        alpha = alpha/N_rounds

        x = draw()
        b = A*x
        V = Krylov(A, b, m)
        y = projection(zeros(N), A, b, V, V)
        statistics = (norm(x - y)^2)/alpha
        Statistics_proj[i] = statistics
    end

    x = LinRange(0, 10, 5000)
    kernel_proj = stats.gaussian_kde(Statistics_proj.*(Statistics_proj.<10))
    stat_proj = kernel_proj(x)
    return stat_proj
end

"Symmetric matices, the single estimation of the scale for all run, multivariate Student."
function HSymSamplingOut(N_rounds::Integer, scale::Real, m::Integer,  N::Integer=100, N_samples::Integer=5000; randomization="A")
    A, s, U = randS(N, scale)
    L = U*Diagonal(s.^(-1/2))
    inv_A = U*Diagonal(s.^(-1))*transpose(U)
    Statistics_proj = zeros(N_samples)
    if randomization == "A"
        draw = (n=N, l=L) -> randx(zeros(n), l)
    else
        draw = (n=N) -> randx(zeros(n))
    end

    alpha = 0
    for i=1:N_rounds
        z = draw()
        b = A*z
        V = Krylov(A, b, m)
        P = I - V*inv(transpose(V)*A*V)*transpose(V)*A
        sample = P*z
        alpha = alpha + norm(sample)^2
    end
    alpha = alpha/N_rounds

    for i=1:N_samples
        x = draw()
        b = A*x
        V = Krylov(A, b, m)
        y = projection(zeros(N), A, b, V, V)
        statistics = (norm(x - y)^2)/alpha
        Statistics_proj[i] = statistics
    end

    x = LinRange(0, 10, 5000)
    kernel_proj = stats.gaussian_kde(Statistics_proj.*(Statistics_proj.<10))
    stat_proj = kernel_proj(x)
    return stat_proj
end

max_stars_repo_name: VLSF/BayesKrylov

max_issues_repo_name: VLSF/BayesKrylov

max_forks_repo_name: VLSF/BayesKrylov


text:
function Line(arg0::Line)
    return Line((Line,), arg0)
end

function Line(arg0::Vector3D, arg1::Vector3D, arg2::jdouble)
    return Line((Vector3D, Vector3D, jdouble), arg0, arg1, arg2)
end

function closest_point(obj::Line, arg0::Line)
    return jcall(obj, "closestPoint", Vector3D, (Line,), arg0)
end

function contains(obj::Line, arg0::Vector3D)
    return jcall(obj, "contains", jboolean, (Vector3D,), arg0)
end

function distance(obj::Line, arg0::Line)
    return jcall(obj, "distance", jdouble, (Line,), arg0)
end

function distance(obj::Line, arg0::Vector3D)
    return jcall(obj, "distance", jdouble, (Vector3D,), arg0)
end

function from_direction(::Type{Line}, arg0::Vector3D, arg1::Vector3D, arg2::jdouble)
    return jcall(Line, "fromDirection", Line, (Vector3D, Vector3D, jdouble), arg0, arg1, arg2)
end

function get_abscissa(obj::Line, arg0::Vector3D)
    return jcall(obj, "getAbscissa", jdouble, (Vector3D,), arg0)
end

function get_direction(obj::Line)
    return jcall(obj, "getDirection", Vector3D, ())
end

function get_origin(obj::Line)
    return jcall(obj, "getOrigin", Vector3D, ())
end

function get_tolerance(obj::Line)
    return jcall(obj, "getTolerance", jdouble, ())
end

function intersection(obj::Line, arg0::Line)
    return jcall(obj, "intersection", Vector3D, (Line,), arg0)
end

function is_similar_to(obj::Line, arg0::Line)
    return jcall(obj, "isSimilarTo", jboolean, (Line,), arg0)
end

function point_at(obj::Line, arg0::jdouble)
    return jcall(obj, "pointAt", Vector3D, (jdouble,), arg0)
end

function reset(obj::Line, arg0::Vector3D, arg1::Vector3D)
    return jcall(obj, "reset", void, (Vector3D, Vector3D), arg0, arg1)
end

function revert(obj::Line)
    return jcall(obj, "revert", Line, ())
end

function to_space(obj::Line, arg0::JVector)
    return jcall(obj, "toSpace", Vector3D, (JVector,), arg0)
end

function to_space(obj::Line, arg0::Point)
    return jcall(obj, "toSpace", Vector3D, (Point,), arg0)
end

function to_sub_space(obj::Line, arg0::JVector)
    return jcall(obj, "toSubSpace", Vector1D, (JVector,), arg0)
end

function to_sub_space(obj::Line, arg0::Point)
    return jcall(obj, "toSubSpace", Vector1D, (Point,), arg0)
end

function whole_line(obj::Line)
    return jcall(obj, "wholeLine", SubLine, ())
end

ext: jl


text:
# serialize to a file
function save(item, filename::AbstractString)
    open(filename, "w") do fp
        save(item, fp)
    end
end
save(item, fp::IO) = serialize(fp, item)

# restore from a file
function restore(filename::AbstractString)
    open(filename, "r") do fp
        restore(fp)
    end
end
restore(fp::IO) = deserialize(fp)

# partition an array to n parts
function partition{T}(a::Array{T}, n::Integer)
    b = Array{T}[]
    t = floor(Int, length(a) / n)
    cursor = 1
    for i in 1:n
        push!(b, a[cursor : (i == n ? length(a) : cursor + t - 1)])
        cursor += t
    end
    b
end


text:
points = [
    3im
    1+im
    0im
    -1-im
    0-2im
    1.5+0.5im
    2im
    -1.5+0.5im
    0-2im
    1-im
    0im
    -1+im
]

function draw_curve_and_tangents(points)
    newdraw()
    draw_curve(points,color=:black,line=3)

    S = Spline(points,:closed)
    n = npatches(S)

    d = 1/3
    T = 1:d:n+1-d

    for t in T
        tt = t + d/3
        v = S'(tt)
        v /= 2*abs(v)
        draw_vector(v,S(tt),color=:green)
    end

    finish()
end

function tangent_plot(points)
    S = Spline(points,:closed)
    f(t) = angle(S'(t))/(2pi)
    plot(f,0.99,npatches(S)+1.01)
    plot!(legend=false)
end


text:
using DataFrames: DataFrame

export variables_by_node,
    issmooth,
    isdecomposable,
    isstruct_decomposable,
    isdeterministic,
    iscanonical,
    sat_prob,
    model_count,
    prob_equiv_signature,
    infer_vtree

#####################
# circuit evaluation infrastructure
#####################

import Base: foreach # extend

function foreach(node::LogicCircuit, f_con::Function, f_lit::Function,
                f_a::Function, f_o::Function)
    f_leaf(n) = isliteralgate(n) ? f_lit(n) : f_con(n)
    f_inner(n) = is⋀gate(n) ? f_a(n) : f_o(n)
    foreach(node, f_leaf, f_inner)
    nothing # returning nothing helps save some allocations and time
end

import ..Utils: foldup # extend

"""
    foldup(node::LogicCircuit,
        f_con::Function,
        f_lit::Function,
        f_a::Function,
        f_o::Function,
        ::Type{T}; nload = nload, nsave = nsave, reset=true)::T where {T}

Compute a function bottom-up on the circuit.
`f_con` is called on constant gates, `f_lit` is called on literal gates,
`f_a` is called on conjunctions, and `f_o` is called on disjunctions.
Values of type `T` are passed up the circuit and given to `f_a` and `f_o` through a callback from the children.
"""
function foldup(node::LogicCircuit, f_con::Function, f_lit::Function,
                f_a::Function, f_o::Function, ::Type{T};
                nload = nload, nsave = nsave, reset=true)::T where {T}
    f_leaf(n) = isliteralgate(n) ? f_lit(n)::T : f_con(n)::T
    f_inner(n, call) = is⋀gate(n) ? f_a(n, call)::T : f_o(n, call)::T
    foldup(node, f_leaf, f_inner, T; nload, nsave, reset)
end

import ..Utils: foldup_aggregate # extend

"""
    foldup_aggregate(node::LogicCircuit,
        f_con::Function,
        f_lit::Function,
        f_a::Function,
        f_o::Function,
        ::Type{T};  nload = nload, nsave = nsave, reset=true)::T where T

Compute a function bottom-up on the circuit.
`f_con` is called on constant gates, `f_lit` is called on literal gates,
`f_a` is called on conjunctions, and `f_o` is called on disjunctions.
Values of type `T` are passed up the circuit and given to `f_a` and `f_o` in an aggregate vector from the children.
"""
function foldup_aggregate(node::LogicCircuit, f_con::Function, f_lit::Function,
                          f_a::Function, f_o::Function, ::Type{T};
                          nload = nload, nsave = nsave, reset=true)::T where T
    function f_leaf(n)
        isliteralgate(n) ? f_lit(n)::T : f_con(n)::T
    end
    function f_inner(n, cs)
        is⋀gate(n) ? f_a(n, cs)::T : f_o(n, cs)::T
    end
    foldup_aggregate(node, f_leaf::Function, f_inner::Function, T; nload, nsave, reset)
end

#####################
# variable-scope-based queries
#####################

import .Utils.variables # extend

"""
    variables(root::LogicCircuit)::BitSet

Get a bitset of variables mentioned in the circuit `root`.
"""
function variables(root::LogicCircuit)::BitSet
    f_con(_) = BitSet()
    f_lit(n) = BitSet(variable(n))
    f_inner(n, call) = mapreduce(call, union, children(n))
    foldup(root, f_con, f_lit, f_inner, f_inner, BitSet)
end

"Get the variable scope of each node in the circuit"
function variables_by_node(root::LogicCircuit)::Dict{LogicCircuit,BitSet}
    # variables_by_node(linearize(root))
    scope = Dict{Node,BitSet}()
    f_con(n) = scope[n] = BitSet()
    f_lit(n) = scope[n] = BitSet(variable(n))
    f_inner(n, call) = scope[n] = mapreduce(call, union, children(n))
    foldup(root, f_con, f_lit, f_inner, f_inner, BitSet)
    scope
end

"""
    issmooth(root::LogicCircuit)::Bool

Is the circuit smooth?
"""
function issmooth(root::LogicCircuit)::Bool
    result::Bool = true
    f_con(_) = BitSet()
    f_lit(n) = BitSet(variable(n))
    f_a(_, cs) = reduce(union, cs)
    f_o(_, cs) = begin
        scope = reduce(union, cs)
        result = result && all(c -> c == scope, cs)
        scope
    end
    foldup_aggregate(root, f_con, f_lit, f_a, f_o, BitSet)
    result
end

"""
    isdecomposable(root::LogicCircuit)::Bool

Is the circuit decomposable?
"""
function isdecomposable(root::LogicCircuit)::Bool
    result::Bool = true
    f_con(_) = BitSet()
    f_lit(n) = BitSet(variable(n))
    f_a(_, cs) = begin
        result = result && isdisjoint(cs...)
        reduce(union, cs)
    end
    f_o(_, cs) = reduce(union, cs)
    foldup_aggregate(root, f_con, f_lit, f_a, f_o, BitSet)
    result
end

"""
    isstruct_decomposable(root::LogicCircuit)::Bool

Is the circuit structured-decomposable?
"""
function isstruct_decomposable(root::LogicCircuit)::Bool
    result::Bool = true
    f_con(_) = [BitSet()]
    f_lit(n) = [BitSet(variable(n))]
    f_a(_, cs) = begin
        result = result && isdisjoint(vcat(cs...)...)
        map(c -> reduce(union!, c), cs)
    end
    f_o(_, cs) = begin
        result = result && (length(cs) == 0 || all(==(cs[1]), cs))
        [reduce(union, vcat(cs...))]
    end
    foldup_aggregate(root, f_con, f_lit, f_a, f_o, Vector{BitSet})
    result
end

"""
    infer_vtree(root::LogicCircuit)::Vtree
Infer circuits vtree if the circuit is struct-decomposable it.
Otherwise return `nothing`.
"""
function infer_vtree(root::LogicCircuit)::Union{Vtree, Nothing}
    # TODO remove after the `isstruct_decomposable` bug is fixed
    if !issmooth(root)
        throw("Circuit not smooth. Inferring vtree not supported yet!")
    end

    if !isstruct_decomposable(root)
        return nothing # or should we throw error?
    end

    f_con(_) = nothing # can we have constants when there is a vtree?
    f_lit(n) = begin
        PlainVtree(variable(n))
    end
    f_a(n, call) = begin
        @assert  num_children(n) == 2 "And node had $num_children(n) childern. Should be 2"
        left = call(children(n)[1])::Vtree
        right = call(children(n)[2])::Vtree
        PlainVtreeInnerNode(left, right)
    end
    f_o(n, call) = begin
        # Already checked struct-decomposable so just expand on first child
        @assert num_children(n) > 0 "Or node has no children"
        call(children(n)[1])
    end
    foldup(root, f_con, f_lit, f_a, f_o, Vtree)
end

"""
    isdeterministic(root::LogicCircuit)::Bool

Is the circuit determinstic?
Note: this function is generally intractable for large circuits.
"""
function isdeterministic(root::LogicCircuit)::Bool
    mgr = sdd_mgr_for(root)
    result::Bool = true
    f_con(c) = mgr(constant(c))
    f_lit(n) = mgr(literal(n))
    f_a(_, cs) = reduce(&, cs)
    f_o(_, cs) = begin
        for i = 1:length(cs)
            for j = i+1:length(cs)
               result &= isfalse(cs[i] & cs[j])
            end
        end
        reduce(|, cs)
    end
    foldup_aggregate(root, f_con, f_lit, f_a, f_o, Sdd)
    result
end

#####################
# structural properties
#####################

"""
    iscanonical(circuit::LogicCircuit, k::Int; verbose = false)

Does the given circuit have canonical Or gates, as determined by a probabilistic equivalence check?
"""
function iscanonical(circuit::LogicCircuit, k::Int; verbose = false)
   signatures = prob_equiv_signature(circuit, k)
   decision_nodes_by_signature = groupby(n -> signatures[n], ⋁_nodes(circuit))
   for (signature, nodes) in decision_nodes_by_signature
      if length(nodes) > 1
         if verbose
            println("Equivalent Nodes:")
            for node in nodes
               println("  - Node: $node MC: $(model_count(node))")
            end
         end
         return false
      end
   end
   return true
end

#####################
# algebraic model counting queries
#####################

"""
    sat_prob(root::LogicCircuit; varprob::Function)::Rational{BigInt}

Get the probability that a random world satisties the circuit.
Probability of each variable is given by `varprob` Function which defauls to 1/2 for every variable.
"""
function sat_prob(root::LogicCircuit;
                  varprob::Function = v -> BigInt(1) // BigInt(2))::Rational{BigInt}
    f_con(n) = istrue(n) ? one(Rational{BigInt}) : zero(Rational{BigInt})
    f_lit(n) = ispositive(n) ? varprob(variable(n)) : one(Rational{BigInt}) - varprob(variable(n))
    f_a(n, call) = mapreduce(call, *, children(n))
    f_o(n, call) = mapreduce(call, +, children(n))
    foldup(root, f_con, f_lit, f_a, f_o, Rational{BigInt})
end

"""
    model_count(root::LogicCircuit, num_vars_in_scope::Int = num_variables(root))::BigInt

Get the model count of the circuit.
The `num_vars_in_scope` is set to number of variables in the circuit, but sometimes need to set different values,
for example, if not every variable is mentioned in the circuit.
"""
function model_count(root::LogicCircuit, num_vars_in_scope::Int = num_variables(root))::BigInt
    # note that num_vars_in_scope for computing the model count can be more than num_variables(circuit); in particular when variables in the application are missing from the circuit
    BigInt(sat_prob(root) * BigInt(2)^num_vars_in_scope)
end

const Signature = Vector{Rational{BigInt}}

"""
    prob_equiv_signature(circuit::LogicCircuit, k::Int)::Dict{Union{Var,Node},Signature}

Get a signature for each node using probabilistic equivalence checking.
Note that this implentation may not have any formal guarantees as such.
"""
function prob_equiv_signature(circuit::LogicCircuit, k::Int)::Dict{Union{Var,Node},Signature}
    # uses probability instead of integers to circumvent smoothing, no mod though
    signs::Dict{Union{Var,Node},Signature} = Dict{Union{Var,Node},Signature}()
    prime::Int = 7919 #TODO set as smallest prime larger than num_variables
    randprob() = BigInt(1) .// rand(1:prime,k)
    do_signs(v::Var) = get!(randprob, signs, v)
    f_con(n) = (signs[n] = (istrue(n) ? ones(Rational{BigInt}, k) : zeros(Rational{BigInt}, k)))
    f_lit(n) = (signs[n] = (ispositive(n) ? do_signs(variable(n)) : BigInt(1) .- do_signs(variable(n))))
    f_a(n, call) = (signs[n] = (mapreduce(c -> call(c), (x,y) -> (x .* y), children(n))))
    f_o(n, call) = (signs[n] = (mapreduce(c -> call(c), (x,y) -> (x .+ y), children(n))))
    foldup(circuit, f_con, f_lit, f_a, f_o, Signature)
    signs
end


text:
function filtr(model::DynamicDiscreteModel,data::Array{Int,1})
	T=length(data)
	dx,dy=size(model.mu)

	filter=Array(Float64,dx,T)
	filter[:,1]=model.mu[:,data[1]]/sum(model.mu[:,data[1]])
	for t=2:T
		rho=0.0
		for jx=1:dx
			filter[jx,t]=0
			for ix=1:dx
				filter[jx,t]+=model.m[ix,data[t-1],jx,data[t]]*filter[ix,t-1]
			end
			rho+=filter[jx,t]
		end
		for jx=1:dx
			filter[jx,t]/=rho
		end
	end

	smoother=Array(Float64,dx,T)
	smoother[:,T]=1
	for t=T-1:-1:1
		rho=0.0
		for ix=1:dx
			smoother[ix,t]=0
			for jx=1:dx
				smoother[ix,t]+=model.m[ix,data[t],jx,data[t+1]]*smoother[jx,t+1]
			end
			rho+=smoother[ix,t]
		end
		for ix=1:dx
			smoother[ix,t]/=rho
		end
	end

	conditional=Array(Float64,dx,dx,T)
	for t=1:T-1
		tempsum=0.0
		for jx=1:dx
			for ix=1:dx
				conditional[ix,jx,t]=model.m[ix,data[t],jx,data[t+1]]
				conditional[ix,jx,t]*=filter[ix,t]
				conditional[ix,jx,t]*=smoother[jx,t+1]
				tempsum+=conditional[ix,jx,t]
			end
		end
		for jx=1:dx
			for ix=1:dx
				#note that there is no risk of division by zero
				#reason: this is called only on _observed_ data which must thus have non-zero probability
				conditional[ix,jx,t]=conditional[ix,jx,t]/tempsum
			end
		end
	end
	filter,smoother,conditional
end


text:
# This file was generated, do not modify it. # hide
coerce!(df, :isrenovated => OrderedFactor, :has_basement => OrderedFactor);


text:
VERSION < v"0.7.0-" && __precompile__(false)  # Plots issue?

# If can't load into JuliaBox, use this:
#  include(download("https://raw.githubusercontent.com/mth229/MTH229.jl/master/src/MTH229.jl"))
#
# * Run the command `?visualizations` for a description of some interactive features.
msg = """
Loading the `MTH229` package.

* Run the command `?MTH229` for a short description.

"""

using Compat
Compat.@info msg

"""

`MTH229`: helper functions for using `Julia` with MTH229

This module does two things:

* Install other useful packages with one command (`Plots`, `Roots`,  `ForwardDiff`, `QuadGK`, `SpecialFunctions`, ...)

* Add a number of helper functions.

The helper functions include:

- `secant(f, a, b)`: return a function giving the secant line between ``(a,f(a))`` and ``(b,f(b))``.

- `tangent(f, c)`:  return a function giving the tangent line to ``f(x)`` at the point ``(c,f(c))``.

- `bisection(f, a, b)`: A simple implementation of the bisection
  method. The interval ``[a,b]`` should be a bracketing interval. This
  function makes an illustrative graphic. For real use of the bisection method, the `fzero(f,
  a, b)` function, from the `Roots` package, should be used.

- `adjoint`: This allows the derivative of   a function to be found as with math notation: `f'`.  The notation can be used for higher-order derivatives too: `f''`, `f'''`, ... This uses automatic differentiation from the `ForwardDiff` package.

- `plotif(f, g, a, b)`: Plot the function `f` over the interval `[a,b]` and color differently where ``g(x) > 0`` over ``[a,b]``. By passing in `f` for `g` shows where `f` is positive on `[a,b]`; passing in `f'` shows where `f` is increasing on `[a,b]`; and passing in `f''` shows where `f` is concave up on `[a,b]`.

- `riemann(f, a, b, n; method="right")` An implementation of Riemann sums. The method can be "right" or "left" for Riemann sums, or "trapezoid" or "simpsons" for related approximations.

"""
module MTH229
#There are a collection of "demos". Run `?visualizations` for a description.

using Reexport
@reexport using Plots
#@reexport using Roots

@reexport using SpecialFunctions

import ForwardDiff
import QuadGK: quadgk
export quadgk

#using Interact

###
export tangent, secant
export lim,  bisection, riemann
export plotif, trimplot, signchart

" f'(x) will find the derivative of `f` using Automatic Differentation from the `ForwardDiff` package "
Base.adjoint(f::Function) = x -> ForwardDiff.derivative(f, float(x))
D(f, n=1) = n > 1 ? D(D(f), n-1) : x -> ForwardDiff.derivative(f, float(x))
export D

"""
Returns a function describing the tangent line to the graph of f at x=c.

Example. Where does the tangent line intersect the y axis?
```
f(x) = sin(x)
tl(x) = tangent(f, pi/4)(x)  # or tl = tangent(f, pi/3) to use a non-generic function
tl(0)
```

Uses the automatic derivative of `f` to find the slope of the tangent line at `x=c`.

"""
tangent(f,c) = x -> f(c) + f'(c) * (x-c)

"""
Returns a function describing the secant line to the graph of f at x=a and x=b.

Example. Where does the secant line intersect the y axis?
```
f(x) = sin(x)
a, b = pi/4, pi/3
sl(x) = secant(f, a, b)(x)  # or sl = sl(f, a, b) to use a non-generic function
sl(0)
```

"""
secant(f, a, b) = x -> f(a) + (f(b) - f(a)) / (b-a) * (x - a)

"""

`lim(f, c, n, dir="+")`: means to generate numeric table of values of `f` as `h` gets close to `c`.

Example:
```
f(x) = sin(x) / x
lim(f, 0)
```
"""
function lim(f::Function, c::Real; n::Int=6, dir="+")
	 hs = [(1/10)^i for i in 1:n] # close to 0
	 if dir == "+"
	   xs = c .+ hs
	 else
	   xs = c .- hs
	 end
	 ys = map(f, xs)
	 [xs ys]
end

"""

Simple implementation of the bisection method.

Example:

```julia
bisection(sin, 3, 4)
f(x) = x^5 - x^4 - x^3 - x^2 - x - 1
a = bisection(f, 1, 2)
f(a)
```

The display shows a simple graphic illustrating the method's division for the first few steps.

An easier-to-understand alternative to `Roots.find_zero(f, (a,b), Bisection())`.

"""
function bisection(f::Function, a, b)
    a,b = sort([a,b])

    if f(a) * f(b) > 0
        error("[a,b] is not a bracket. A bracket means f(a) and f(b) have different signs!")
    end

    M = a + (b-a) / 2

    i, j = 0, 64
    ss = fill("#", 65)
    ss[i+1]="a"; ss[j+1]="b"
    println("")
    println(join(ss))
    flag = true

    while a < M < b
        if flag && j-i == 1
            ss = fill(" ", 65)
            ss[j:(j+1)] = "⋮"
            println(join(ss))
            println("")
            flag = false
        end

        if f(M) == 0.0
            println("... exact answer found ...")
	    break
        end
        ## update step
	if f(a) * f(M) < 0
	    a, b = a, M

            if flag
                j = div(i + j, 2)
            end

	else
	    a, b = M, b

            if flag
                i = div(i + j, 2)
            end

	end

        if flag
            ss = fill(".", 65)
            ss[i+1]="a"; ss[j+1]="b"; ss[(i+2):j]="#"
            println(join(ss))
        end

        M = a + (b-a) / 2
    end
    M
end

import Roots
import Roots: newton, find_zero, find_zeros
newton(f, fp, x0; kwargs...) = Roots.find_zero((f,fp), x0, Roots.Newton(); kwargs...)
newton(f, x0; kwargs...) = newton(f, D(f), x0; kwargs...)
fzero(f, x0; kwargs...) = Roots.find_zero(f, x0; kwargs...)
fzero(f, a, b; kwargs...) = Roots.find_zero(f, (a, b); kwargs...)
fzeros(f, a, b; kwargs...) = Roots.find_zeros(f, a, b; kwargs...)

export newton, fzero, fzeros

# some plotting utilities

"""
   trimplot(f, a, b, c=20; kwargs...)

Plot f over [a,b] but break graph if it exceeds c in absolute value.
"""
function trimplot(f, a, b, c=20; kwargs...)
  xs = linspace(a, b, 251) #range(a, stop=b, length=251)
  ys = f.(xs)

  us, vs = Real[], Real[]
  p = plot(us, vs, xlim=(a, b), legend=false, kwargs...)
  for (x,y) in zip(xs, ys)
    if abs(y) <= c
       push!(us, x); push!(vs, y)
    else
      length(us) > 0 && plot!(p, us, vs, color=:blue)
      empty!(us); empty!(vs)
    end
 end
 length(us) > 0 && plot!(p, us, vs, color=:blue)
 p
end

"""
    plotif(f, g, a, b)

Plot f colored depending on g < 0 or not.
"""
function plotif(f, g, a, b)
    xs = linspace(a, b, 251)#range(a, stop=b, length=251)
    ys = f.(xs)
    p = plot(xs, ys, color=:blue, linewidth=5, legend=false)
    zs = [g(x) < 0 ? NaN : f(x) for x in xs]
    plot!(p, xs, zs, color=:red, linewidth=5)
    p
end

"""
   signchart(f, a, b)

Plot f over a,b with different color when negative.
"""
function signchart(f, a, b)
    p = plotif(f, f, a, b)
    plot!(p, zero)
    p
end

# visualize newtons method
function newton_vis(f, x0, a=Inf,b=-Inf; steps=5, kwargs...)
    xs = Float64[x0]
    for i in 1:steps
        push!(xs, xs[end] - f(xs[end]) / f'(xs[end]))
    end

    m,M = extrema(xs)
    m = min(m, a)
    M = max(M, b)

    p = plot(f, m, M; linewidth=3, legend=false, kwargs...)
    plot!(p, zero, m, M)
    for i in 1:steps
        plot!(p, [xs[i],xs[i],xs[i+1]], [0,f(xs[i]), 0])
        scatter!(p, xs[i:i],[0])
    end
    scatter!(p, [xs[steps+1]], [0])
    p
end

"""
riemann: compute Riemann sum approximations to a definite integral. As well, implement trapezoid and Simpson's rule.

Example:
```
f(x) = exp(x^2)
riemann(f, 0, 1, 1000)   # default right-Riemann sums
riemann(f, 0, 1, 1000, method="left")       # left sums
riemann(f, 0, 1, 1000, method="trapezoid")  # use trapezoid rule
riemann(f, 0, 1, 1000, method="simpsons")   # use Simpson's rule
```

"""
function riemann(f::Function, a::Real, b::Real, n::Int; method="right")
  if method == "right"
     meth = (f,l,r) -> f(r) * (r-l)
  elseif method == "left"
     meth= (f,l,r) -> f(l) * (r-l)
  elseif method == "trapezoid"
     meth = (f,l,r) -> (1/2) * (f(l) + f(r)) * (r-l)
  elseif method == "simpsons"
     meth = (f,l,r) -> (1/6) * (f(l) + 4*(f((l+r)/2)) + f(r)) * (r-l)
  end

  xs = a + (0:n) * (b-a)/n
  as = [meth(f, l, r) for (l,r) in zip(xs[1:end-1], xs[2:end])]
  sum(as)
end

#import SymPy: real_roots
#real_roots(f; kwargs...) = PolynomialZeros.poly_roots(f, Over.R, kwargs...)

###
#include("demos.jl")

end


text:
"""
Constrained Optimization and Sensitivity for Thermochronology
R. Stoner
2020

See documentation for detailed explanation of COAST
"""
module COAST

## Dependencies
using Ipopt
using JuMP
using Random
using LinearAlgebra
using Interpolations
using Distributions
## Helium diffusion
include("he_preprocessing.jl")
include("he_forward.jl")
include("he_vars_constraints.jl")
include("he_func_registration.jl")
include("FT_setup.jl")
include("FT_forward.jl")
include("zonation.jl")
include("app.jl")

using .App
using .Zoned

# funcs for unit testing
export loaded_COAST

# preprocessing
export decompose_eu, conc_to_atoms_per_volume, UTh_date_to_He_mols

# treating setup
export define_variables!, constraint_func, constrain_upper
export rdaam_define_constraints

# He forward model
export rdaam_forward_diffusion, zrdaam_forward_diffusion, jac_rdaam_forward_diffusion, initialize_JuMP_model
export register_forward_model!, register_objective_function!, set_objective_function!
export fill_u_term2

# FT forward model
export init_track_len
export FT_forward

# constants
export sec_in_yrs,τ38,τ35,τ32

## testing funcs
function loaded_COAST()
  load_success = true
  return load_success
end

end


text:
function write(io::IO, x::Any, verbose=false)
    e = Encoder(io, verbose)
    encode_top_level(e, x)
end

function to_transit(x::Any, verbose=false)
    sprint(io-> write(io, x, verbose))
end


text:
__precompile__()

module CornerPlot

using DataFrames
using Gadfly
using Compose
using Viridis
using Measures

export corner

function corner(samples::Union{Array{Float64}, AbstractDataFrame};
                bins::Union{Int, Array{Int, 1}, Array{Array{Float64, 1}, 1}} = 20,
                hexbins::Union{Int, Array{Int, 1}} = 50,
                range::Union{Array{Tuple{Float64,Float64}, 1}, Void} = nothing,
                varnames::Union{Array{String, 1}, Array{Symbol, 1}, Void} = nothing,
                truthvals::Union{Array{Float64, 1}, Array{Array{Float64, 1}}, Void} = nothing,
                ylabel_orientation::Symbol = :vertical,
                plotsize::Measures.Length{:mm,Float64} = 20cm)

    # If samples are in dataframe, get array of values
    if typeof(samples) <: AbstractDataFrame

        # only take relevant columns if given
        if typeof(varnames) == Array{Symbol, 1}
            samples = Array(samples[:, varnames])
        end

    # in case a 1-d array is passed
    elseif length(size(samples)) == 1
        samples = reshape(samples, (size(samples)[1], 1))
    end

    nsamps = size(samples)[1]
    ndims = size(samples)[2]

    if ndims > nsamps
        throw(ArgumentError("samples array should not have more dimensions than samples, try reshaping/transposing"))
    end

    # now to deal with bins
    if typeof(bins) == Int
        nbins = [bins for a in 1:ndims]

        if range != nothing
            bins = [linspace(r..., nbins[i]) for (i, r) in enumerate(range)]
        else
#            bins = nothing
            bins = [[minimum(samples[:, i]), maximum(samples[:, i])] for i in 1:ndims]
        end

    elseif typeof(bins) == Array{Int, 1}
        nbins = bins

        if range != nothing
            bins = [linspace(r..., nbins[i]) for (i, r) in enumerate(range)]
        else
#            bins = nothing
            bins = [[minimum(samples[:, i]), maximum(samples[:, i])] for i in 1:ndims]
        end

    else
        nbins = [length(a) - 1 for a in bins]
    end

    # and hexbins
    if typeof(hexbins) == Int
        hexbins = [hexbins for a in 1:ndims]
    end

    # set up the plots
    set_default_plot_size(plotsize, plotsize)
    subplots = Array{Context}(ndims, ndims)

    for i in 1:ndims

        # make an array for the layers of the plots
        diagarr = []

        # make a vertical line layer if requested
        if truthvals != nothing

            # turn a single array into an array of arrays
            if typeof(truthvals) == Array{Float64, 1}
                truthvals = [[val] for val in truthvals]
            end

            vline = layer(xintercept = truthvals[i], Geom.vline,
                          Theme(default_color = colorant"red"))
            append!(diagarr, vline)
        end

        # grab only those samples which will be in the plot range
        if bins != nothing
            histsamps = samples[:, i][bins[i][1] .<= samples[:, i] .< bins[i][end]]
        else
            histsamps = samples[:, i]
        end

        # this is the 1-d histogram layer
        append!(diagarr, layer(x = histsamps, Geom.histogram(bincount = nbins[i]),
                               Theme(default_color = Viridis.viridis(0.5))))

        # for last diagonal, add axis label if given
        xlabel = nothing
        ylabel = nothing
        if varnames != nothing
            if i == 1
                ylabel = String(varnames[i])
            elseif i == ndims
                xlabel = String(varnames[i])
            end
        end

        # and plot the diagonal histogram
        subplots[i, i] = render(plot(diagarr...,
                                     Guide.xlabel(nothing),
                                     Guide.ylabel(nothing),))

        # now loop over the subplots down the column
        for j in (i + 1):ndims

            # make layer array for each offdiagonal plot
            offdiag = []

            # add the horizontal and vertical lines
            if truthvals != nothing
                append!(offdiag, vline)
                hline = layer(yintercept = truthvals[j], Geom.hline,
                              Theme(default_color = colorant"red"))
                append!(offdiag, hline)
            end

            x=histsamps
            y=samples[:, j][bins[j][1] .<= samples[:, j] .< bins[j][end]]

            # get a hexbin layer
            append!(offdiag, layer(x=x, y=y,
                                   Geom.hexbin(xbincount = hexbins[i],
                                               ybincount = hexbins[j])))

            # and add the axis labels if given
            if varnames != nothing

                # ylabels only along the side
                if i == 1
                    ylabel = String(varnames[j])
                else
                    ylabel = nothing
                end

                # xlabels only along the bottom
                if j == ndims
                    xlabel = String(varnames[i])
                else
                    xlabel = nothing
                end

            else
                xlabel = nothing
                ylabel = nothing
            end

            # and plot the layers
            subplots[j, i] = render(plot(offdiag...,
                                    Scale.color_continuous(colormap = Viridis.viridis),
                                    Guide.xlabel(xlabel),
                                    Guide.ylabel(ylabel,
                                                 orientation = ylabel_orientation),
                                    Coord.Cartesian(xmin = bins[i][1], xmax = bins[i][end],
                                                    ymin = bins[j][1], ymax = bins[j][end]),
                                    Theme(key_position = :none)))

            # make subplots above diagonal empty
            subplots[i, j] = Compose.context()
        end
    end
    gridstack(subplots)
end

end


text:
using Distributed
@noinline call_no_inline(f,args) = f(args...)
"""
    s = async_signal(f, args...; init = nothing)

Create a signal initialized to `init` whos action `f(args...)` will run asynchronously in a
different task whenever its arguments update.async signals only work in a push based paradigm.
"""
function async_signal(f, args...; init = nothing)
    res = Signal(init)
    Signal(args...) do args...
        @async begin
            try
                res(call_no_inline(f, args))
            catch e
                st = catch_stacktrace()
                return handle_err(e, st)
            end
        end
    end
    res
end
export async_signal

"""
    s = remote_signal(f, args...; init = nothing, procid = first(workers()))

Create a signal initialized to `init` whos action `f(args...)` will run remotely in
a process with id `procid`, whenever its arguments update.remote signals only work
in a push based paradigm.
"""
function remote_signal(f, args...; init = nothing, procid = first(workers()))
    res = Signal(init)
    Signal(args...) do args...
        @async begin
             try
                 x = @fetchfrom procid call_no_inline(f, args)
                 res(x)
             catch e
                 st = [st_tuple[1] for st_tuple in e.captured.processed_bt]
                 err = e.captured.ex
                 return handle_err(err, st)
             end
        end
    end
    res
end
export remote_signal


text:
module TestExamples2DTestExamples2DPart2

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "2D-Part2" begin

# Run basic tests
@testset "Examples 2D" begin
  # Compressible Euler Multicomponent
  include("test_examples_2d_eulermulti.jl")

  # MHD
  include("test_examples_2d_mhd.jl")

  # Lattice-Boltzmann
  include("test_examples_2d_lbm.jl")

  # Acoustic perturbation
  include("test_examples_2d_ape.jl")
end

# Coverage test for all initial conditions
@testset "Tests for initial conditions" begin
  # GLM-MHD
  @testset "elixir_mhd_alfven_wave.jl one step with initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [7.144325530681224e-17, 2.123397983547417e-16, 5.061138912500049e-16, 3.6588423152083e-17, 8.449816179702522e-15, 3.9171737639099993e-16, 2.445565690318772e-16, 3.6588423152083e-17, 9.971153407737885e-17],
      linf = [2.220446049250313e-16, 8.465450562766819e-16, 1.8318679906315083e-15, 1.1102230246251565e-16, 1.4210854715202004e-14, 8.881784197001252e-16, 4.440892098500626e-16, 1.1102230246251565e-16, 4.779017148551244e-16],
      maxiters = 1,
      initial_condition = initial_condition_constant)
  end

  @testset "elixir_mhd_rotor.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"),
      l2   = [1.2429643014570362, 1.7996363685163963, 1.6899889382208215, 0.0, 2.263014495582255, 0.2127948919559293, 0.23341167012961367, 0.0, 0.003341061230142962],
      linf = [10.401662016997985, 14.061042946011094, 15.556382737749237, 0.0, 16.714999099689578, 1.3250449624793987, 1.4148467737020096, 0.0, 0.0878998114279951],
      tspan = (0.0, 0.05))
  end

  @testset "elixir_mhd_blast_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_blast_wave.jl"),
      l2   = [0.17569823349539196, 3.853292514951796, 2.474036084054808, 0.0, 355.36545703316915, 2.3525087027248084, 1.394705056983077, 0.0, 0.029910308624236454],
      linf = [1.5831869462980066, 44.20237543674303, 12.866364462538552, 0.0, 2237.8614138686, 13.066894956593798, 8.984965484247244, 0.0, 0.5226498664960756],
      tspan = (0.0, 0.003))
  end

  # LBM
  @testset "elixir_lbm_couette.jl with initial_condition_couette_steady" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_lbm_couette.jl"),
      l2   = [9.321369073400123e-16, 1.6498793963435488e-6, 5.211495843124065e-16,
              1.6520893954826173e-6, 1.0406056181388841e-5, 8.801606429417205e-6,
              8.801710065560555e-6, 1.040614383799995e-5, 2.6135657178357052e-15],
      linf = [1.4432899320127035e-15, 2.1821189867266e-6, 8.881784197001252e-16,
              2.2481261510165496e-6, 1.0692966335143494e-5, 9.606391697600247e-6,
              9.62138334279633e-6, 1.0725969916147021e-5, 3.3861802251067274e-15],
      initial_condition=initial_condition_couette_steady,
      tspan = (0.0, 1.0))
  end

  @testset "elixir_lbm_lid_driven_cavity.jl with stationary walls" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_lbm_lid_driven_cavity.jl"),
      l2   = [1.7198203373689985e-16, 1.685644347036533e-16, 2.1604974801394525e-16,
              2.1527076266915764e-16, 4.2170298143732604e-17, 5.160156233016299e-17,
              6.167794865198169e-17, 5.24166554417795e-17, 6.694740573885739e-16],
      linf = [5.967448757360216e-16, 6.522560269672795e-16, 6.522560269672795e-16,
              6.245004513516506e-16, 2.1163626406917047e-16, 2.185751579730777e-16,
              2.185751579730777e-16, 2.393918396847994e-16, 1.887379141862766e-15],
      boundary_conditions=boundary_condition_wall_noslip,
      tspan = (0, 0.1))
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end # 2D-Part2

end #module


text:
using Test
using Dorothy

@DorothyAll()

@testset "Selection" begin

	@testset "1BTL" begin # Small orthorhombic system
		model = readf("$(Dorothy.datapath)/1BTL.pdb")
		sel1 = view(model, Water)
		@test sel1.ids == 2039:2237
		@test length(view(model, Protein)) == 2032
		@test length(view(model, Water)) == 199
		@test length(view(model, Water | Protein)) == 2032 + 199
		@test length(view(model, Water & Protein)) == 0
		@test length(view(model, Index(1, by=Residue))) == 10
		@test length(view(model, Expand(Index(1), by=Residue))) == 10
		@test length(view(model, Index(1, by=Chain))) == 2236
		@test length(view(model, Expand(Index(1), by=Chain))) == 2236
		sel2 = view(model, Water & Within(5.0, of=Nter))
		@test length(sel2) == 7
		sel3 = view(model, Water & Within(45.0, of=Nter))
		@test length(sel3) == 199
	end

	@testset "1BTL-np" begin # Non-periodic system
		model = readf("$(Dorothy.datapath)/1BTL.pdb")
		delete!(model.header, :cell)
		sel1 = view(model, Water & Within(45.0, of=Nter))
		@test length(sel1) == 188
	end

	@testset "MHC" begin # Rhombododecahedral system
		model = readf("$(Dorothy.datapath)/MHC.pdb")
		sel1 = view(model, Water & Within(5.0, of=Protein))
		@test length(sel1) == 6475
		sel2 = view(model, Cter)
		@test length(sel2) == 56
		@test length(residues(sel2)) == 3
		sel3 = view(model, Nter)
		@test length(sel3) == 57
		@test length(residues(sel3)) == 3
		bitvec = falses(length(model))
		for i = 1:6539
			bitvec[i] = true
		end
		sel4 = view(model, Map(bitvec))
		@test length(sel4) == 6539
		sel5 = view(model, Protein)
		@test sel5 == sel4
		@test length(sel5) == 6539
		sel6 = view(model, Index(1:11))
		@test length(sel6) == 11
		@test all(sel6.resnames .== "GLY")
		sel7 = view(model, Index(1))
		@test sel7[].name == "MN1"
		sel8 = view(model, Index(277, by=Residue))
		@test length(sel8) == 27
		@test all(sel8.resnames .== "ILE")
		sel9 = view(model, Index(3, by=Chain))
		@test length(sel9) == 149
		@test all(sel9.chainids .== "C")
		sel10 = view(model, Index(first, ineach=Chain))
		@test length(sel10) == 6
		Dorothy.Topology.infertopology!(model)
		@test length(fragments(view(model, 1:6539))) == 3
		sel11 = view(model, Index(3, by=Fragment))
		@test sel11 == sel9
		@test length(fragments(view(model, Protein))) == 3
		sel12 = view(model, Expand(Index(6539), by=Chain))
		@test sel12 == sel11
		sel13 = view(model, Expand(Index(4657), by=Residue))
		@test sel13 == sel8
		sel14 = view(model, Restrict(Index(6000:7000), by=Chain))
		@test sel14 == sel11
		sel15 = view(model, Restrict(Index(4650:4690), by=Residue))
		@test sel15 == sel8
		sel16 = view(model, Restrict(Index(6541:6565), by=Fragment))
		@test length(sel16) == 3*7
		sel17 = view(model, Heavy & Water)
		@test length(sel17) == 20579
		sel18 = view(model, AcidResidue & Cα)
		sel19 = view(model, BasicResidue & Cα)
		sel20 = view(model, ChargedResidue & Cα)
		sel21 = view(model, PolarResidue & Cα)
		sel22 = view(model, HydrophobicResidue & Cα)
		@test sum(length.([sel18, sel19])) == length(sel20)
		@test sum(length.([sel20, sel21, sel22])) == 276 + 99 + 9
		sel23 = view(model, Backbone)
		@test length(sel23) == 3 * (276 + 99 + 9)
		sel24 = view(model, MainChain)
		@test length(sel24) == 5 * (276 + 99 + 9) + 2*3 + 1*3 - (16 + 5)
		sel25 = view(model, SideChain & MainChain)
		@test length(sel25) == 0
		sel26 = view(model, SideChain | MainChain)
		@test length(sel26) == length(sel5)
	end

	@testset "PLC" begin # Very large orthorhombic system
		model = readf("$(Dorothy.datapath)/PLC.pdb.bz2")
		sel1 = view(model, Water & Expand(Within(3.0, of=Lipid), by=Residue))
		@test length(sel1) == 30474
		sel2 = view(model, Water & Restrict(Within(3.0, of=Lipid), by=Residue))
		@test length(sel2) == 14271
		sel3 = view(model, Name("CA"))
		@test length(sel3) == 4854
		sel4 = view(model, Id(1:10))
		@test length(sel4) == 10
		sel5 = view(model, ResName("LEU"))
		@test length(sel5) == 11708
		sel6 = view(model, ResName("LEU") & Name("CA"))
		@test length(sel6) == 509
		sel7 = view(model, ResName("LEU") & ResId(732))
		@test length(sel7) == 24
		cache = SelectionCache()
		sel8 = view(model, ChainId("X") & Element("Cl"), cache)
		@test length(sel8) == 679
		sel9 = view(model, ChainId("X") & Element("Na"), cache)
		@test length(sel9) == 813
		sel10 = view(model, ChainId("X") & !Element("Na"), cache)
		@test sel10 == sel8
		sel11 = view(model, ChainId("X") & Element("Na") & Element("Cl"))
		@test length(sel11) == 0
		sel12 = view(model, ChainId("X") & (Element("Na") | Element("Cl")))
		@test length(sel12) == 679 + 813
		sel13 = view(model, ChainId("X") & Element(["Na", "Cl"]))
		@test sel13 == sel12
		sel14 = view(model, Mass(mi -> (mi < 2.0)))
		sel15 = view(model, Hydrogen | VSite)
		@test sel14 == sel15
		sel16 = view(model, Occupancy(1.0))
		@test length(sel16) == length(model)
		sel17 = view(model, R([120.340, 167.910, 179.540]))
		@test length(sel17) == 1
		sel18 = view(model, R(Ri -> (Ri.z > 100.0)))
		@test length(sel18) == 1104250
		# Stride test requires binary dependency
		# sel19 = view(model, SS("H"))
		sel20 = view(model, Ion)
		@test sel20 == sel12
		sel21 = view(model, Lipid & Element("P"))
		@test length(sel21) == 1943
	end

end # @testset


text:
using Test
using PSDMatrices
using LinearAlgebra
using Suppressor

M_square = [1 1; 2 20]
M_tall = [1 1; 2 20; 3 30]
M_wide = [1 1 1; 2 20 200]
M_neg = -ones(1, 1)
eltypes = (Int64, Float64, BigFloat)

@testset "PSDMatrices.jl" begin
    @testset "eltype=$t | shape=$(size(Mbase))" for t in eltypes,
        Mbase in (M_square, M_tall, M_wide, M_neg)

        M = t.(Mbase)
        S = PSDMatrix(M)
        X = rand(size(S, 1), size(S, 2))

        @testset "Base" begin
            @test eltype(S) == t
            @test size(S) == size(Matrix(S))
            @test S == S
            @test copy(S) == S
            @test !(copy(S) === S)
            @test typeof(similar(S)) == typeof(S)
            @test (S2 = similar(S); copy!(S2, S); S2 == S)
            if size(M, 1) == size(M, 2)
                @test Matrix(inv(S)) ≈ inv(Matrix(S))
            end
            @suppress_out @test isnothing(show(S))
            @suppress_out @test isnothing(show(stdout, MIME("text/plain"), S))
        end

        @testset "LinearAlgebra" begin
            @test diag(S) ≈ diag(Matrix(S))
            if size(M, 1) == size(M, 2)
                @test det(S) ≈ det(Matrix(S))
                @test logdet(S) ≈ logdet(Matrix(S))
            else
                @test_throws MethodError det(S)
                @test_throws MethodError logdet(S)
            end
            if (size(M, 1) >= size(M, 2))
                @test S \ X ≈ Matrix(S) \ X
                @test X / S ≈ X / Matrix(S)
            end
        end

        @testset "Exports" begin
            @test norm(Matrix(S) - M' * M) == 0.0
            @test Matrix(X_A_Xt(S, X)) ≈ X * Matrix(S) * X'
            @test Matrix(X_A_Xt(A=S, X=X)) ≈ X * Matrix(S) * X'
            @test begin
                product_eltype = typeof(one(eltype(X)) * one(eltype(S)))
                S2 = PSDMatrix(zeros(product_eltype, size(S.R)...))
                X_A_Xt!(S2, S, X)
                Matrix(S2) ≈ Matrix(X_A_Xt(S, X))
            end
            @test begin
                product_eltype = typeof(one(eltype(X)) * one(eltype(S)))
                S2 = PSDMatrix(zeros(product_eltype, size(S.R)...))
                X_A_Xt!(S2, A=S, X=X)
                Matrix(S2) ≈ Matrix(X_A_Xt(S, X))
            end
            @test Matrix(add_qr(S, S)) ≈ Matrix(S) + Matrix(S)
            if (size(M, 1) >= size(M, 2))
                @test Matrix(add_cholesky(S, S)) ≈ Matrix(S) + Matrix(S)
                tri = triangularize_factor(S)
                @test tri.R isa UpperTriangular
                @test Matrix(tri) ≈ Matrix(S)
            end
        end
    end
end
nothing


text:
function build_bnf_rules(ast::ASTNode, addwhitespace::Bool, rulenameprefix="")
    assign_rulenames(ast, rulenameprefix)
    rules = Vector{RuleSource}()
    build_bnf_rule(ast, addwhitespace, rules)
    rules
end

function build_bnf_rule(node::ASTNode, addwhitespace, rules::Vector{RuleSource})

  # TODO: reduce number of rules, by not having separate ones for refs and literals? - have own bnf_build_rulename?
  # TODO: could also collapse RHS "choice" into LHS "variable def" during transform?
  if node.func in [:bnf]
	  build_bnf_call(node, rules)
  elseif node.func in [:variabledef]
	  build_bnf_variable(node, rules)
  elseif node.func in [:variableref]
	  # do nothing - we handle this directly in calling node
  elseif node.func in [:concatenation]
	  build_bnf_concatenation(node, addwhitespace, rules)
  elseif node.func in [:alternation]
	  build_bnf_alternation(node, rules)
  elseif node.func in [:optional]
	  build_bnf_optional(node, rules)
  elseif node.func in [:quantifier]
	  build_bnf_quantifier(node, addwhitespace, rules)
  elseif node.func in [:literal]
	  # do nothing - we handle this directly in calling node
  elseif node.func in [:regexp]
	  build_bnf_regexp(node, rules)
  elseif node.func in [:negation]
	  build_bnf_negation(node, rules)
  else
	  error("Unexpected bnf node with function $(node.func)")
  end

  for child in node.children
	  build_bnf_rule(child, addwhitespace, rules)
  end

end

function build_bnf_call(node::ASTNode, rules::Vector{RuleSource})
  rule = build_rule_shortform_start(node)
  @assert length(node.refs)==1
  methodname = build_bnf_called_rulename(node.refs[1])
  push!(rule.source, "$(methodname)")
  build_rule_shortform_end(rule, node)
  push!(rules, rule)
end

function build_bnf_variable(node::ASTNode, rules::Vector{RuleSource})
  rule = build_rule_shortform_start(node)
  @assert length(node.children)==1
  methodname = build_bnf_called_rulename(node.children[1])
  push!(rule.source, "$(methodname)")
  build_rule_shortform_end(rule, node)
  push!(rules, rule)
end

function build_bnf_concatenation(node::ASTNode, addwhitespace, rules::Vector{RuleSource})
  rule = build_rule_start(node)
  push!(rule.source, "  content = Any[]")
  for child in node.children
    methodname = build_bnf_called_rulename(child)
	push!(rule.source, "  push!(content, $(methodname))")
  end
  if addwhitespace
	  push!(rule.source, "  join(content, \" \")")
  else
	  push!(rule.source, "  join(content)")
  end
  build_rule_end(rule, node)
  push!(rules, rule)
end

function build_bnf_alternation(node::ASTNode, rules::Vector{RuleSource}) # this is handled using a rule choice
  for child in node.children
		rule = build_rule_shortform_start(node)
		methodname = build_bnf_called_rulename(child)
		push!(rule.source, "$(methodname)")
		build_rule_shortform_end(rule, node)
	  push!(rules, rule)
  end
end

function build_bnf_optional(node::ASTNode, rules::Vector{RuleSource})
    rule = build_rule_shortform_start(node)
    @assert length(node.children)==1
    methodname = build_bnf_called_rulename(node.children[1])
    push!(rule.source, "choose(Bool) ? $(methodname) : \"\"")
    build_rule_shortform_end(rule, node)
    push!(rules, rule)
end

function build_bnf_quantifier(node::ASTNode, addwhitespace, rules::Vector{RuleSource})
    rule = build_rule_shortform_start(node)
    @assert length(node.children)==1
    methodname = build_bnf_called_rulename(node.children[1])
    if addwhitespace
	    push!(rule.source, "join(reps($(methodname), $(node.args[:min])), \" \")")
    else
	    push!(rule.source, "join(reps($(methodname), $(node.args[:min])))")
    end
    build_rule_shortform_end(rule, node)
    push!(rules, rule)
end

function build_bnf_regexp(node::ASTNode, rules::Vector{RuleSource})
    rule = build_rule_shortform_start(node)
    push!(rule.source, "choose(String, \"$(escape_str(node.args[:value]))\"")
    build_rule_shortform_end(rule, node)
    push!(rules, rule)
end

function build_bnf_negation(node::ASTNode, rules::Vector{RuleSource})
    rule = build_rule_shortform_start(node)
	# currently do not handle: a space is output instead
    push!(rule.source, "\" \"")
    build_rule_shortform_end(rule, node)
    push!(rules, rule)
end

function build_bnf_called_rulename(node::ASTNode)
	if node.func == :literal
	    "\"$(escape_str(node.args[:value]))\""
	elseif node.func == :variableref
	   @assert length(node.refs)==1
	   build_called_rulename(node.refs[1])
	else
	   build_called_rulename(node)
	end
end


text:
# implements a REPL mode

# most of this is lifted from RCall. Thanks guys!

module DouglassPrompt

    using REPL
    import REPL: REPL, LineEdit

    # import ..Douglass:
    #     @d_str,
    #     parse

    function simple_showerror(io::IO, er)
        Base.with_output_color(:red, io) do io
            print(io, "ERROR: ")
            showerror(io, er)
            println(io)
        end
    end

    function parse_status(script::String)
        status = :ok
        try
            #render(script)
        catch ex
            if isa(ex, DParseIncomplete)
                status = :incomplete
            else
                status = :error
            end
        end
        status
    end

    function repl_eval(script::String, stdout::IO, stderr::IO)
        try
#            @d_str(script)
        catch ex
            # TODO check here for specific error types (parse error etc)
            simple_showerror(stderr, ex)
        finally
            return nothing
        end
    end

    function bracketed_paste_callback(s, o...)
        input = LineEdit.bracketed_paste(s)
        sbuffer = LineEdit.buffer(s)
        curspos = position(sbuffer)
        seek(sbuffer, 0)
        shouldeval = (bytesavailable(sbuffer) == curspos && !occursin(UInt8('\n'), sbuffer))
        seek(sbuffer, curspos)
        if curspos == 0
            # if pasting at the beginning, strip leading whitespace
            input = lstrip(input)
        end

        if !shouldeval
            LineEdit.edit_insert(s, input)
            return
        end

        LineEdit.edit_insert(sbuffer, input)
        input = String(take!(sbuffer))

        m = sizeof(input)
        oldpos = 1
        nextpos = 0
        # parse the input line by line
        while nextpos < m
            next_result = findnext("\n", input, nextpos + 1)
            if next_result == nothing
                nextpos = m
            else
                nextpos = next_result[1]
            end
            block = input[oldpos:nextpos]
            status = parse_status(block)

            if status == :error  || (status == :incomplete && nextpos == m) ||
                    (nextpos == m && !endswith(input, '\n'))
                # error / continue and the end / at the end but no new line
                LineEdit.replace_line(s, input[oldpos:end])
                LineEdit.refresh_line(s)
                break
            elseif status == :incomplete && nextpos < m
                continue
            end

            if !isempty(strip(block))
                # put the line on the screen and history
                LineEdit.replace_line(s, strip(block))
                LineEdit.commit_line(s)
                # execute the statement
                terminal = LineEdit.terminal(s)
                REPL.raw!(terminal, false) && LineEdit.disable_bracketed_paste(terminal)
                LineEdit.mode(s).on_done(s, LineEdit.buffer(s), true)
                REPL.raw!(terminal, true) && LineEdit.enable_bracketed_paste(terminal)
            end
            oldpos = nextpos + 1
        end
        LineEdit.refresh_line(s)
    end

    repl_filename(repl, hp::REPL.REPLHistoryProvider) = "REPL[$(max(length(hp.history)-hp.start_idx, 1))]"
    repl_filename(repl, hp) = "REPL"

    function create_d_repl(repl, main)
        d_mode = LineEdit.Prompt("Douglass> ";
            prompt_prefix="\e[38;5;166m",
            prompt_suffix=main.prompt_suffix,
            sticky=true)
        # alternative prompt_prefix: Base.text_colors[:cyan]

        hp = main.hist
        hp.mode_mapping[:d] = d_mode
        d_mode.hist = hp
        # this is the key component: do whatever the main julia repl does, but add the macro call.
        d_mode.on_done = REPL.respond(x->Base.parse_input_line("d\"\"\"" * x * "\"\"\"",filename=repl_filename(repl,hp)), repl, main)

        bracketed_paste_mode_keymap = Dict{Any,Any}(
            "\e[200~" => bracketed_paste_callback
        )

        search_prompt, skeymap = LineEdit.setup_search_keymap(hp)
        prefix_prompt, prefix_keymap = LineEdit.setup_prefix_keymap(hp, d_mode)

        mk = REPL.mode_keymap(main)
        # ^C should not exit prompt
        delete!(mk, "^C")

        b = Dict{Any,Any}[
            bracketed_paste_mode_keymap,
            skeymap, mk, prefix_keymap, LineEdit.history_keymap,
            LineEdit.default_keymap, LineEdit.escape_defaults
        ]
        d_mode.keymap_dict = LineEdit.keymap(b)

        d_mode
    end

    function repl_init(repl)
        mirepl = isdefined(repl,:mi) ? repl.mi : repl
        main_mode = mirepl.interface.modes[1]
        d_mode = create_d_repl(mirepl, main_mode)
        push!(mirepl.interface.modes,d_mode)

        d_prompt_keymap = Dict{Any,Any}(
            '`' => function (s,args...)
                if isempty(s) || position(LineEdit.buffer(s)) == 0
                    buf = copy(LineEdit.buffer(s))
                    LineEdit.transition(s, d_mode) do
                        LineEdit.state(s, d_mode).input_buffer = buf
                    end
                else
                    LineEdit.edit_insert(s, '`')
                end
            end
        )

        main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict, d_prompt_keymap);
        nothing
    end

    function repl_inited(repl)
        mirepl = isdefined(repl,:mi) ? repl.mi : repl
        any(:prompt in fieldnames(typeof(m)) && m.prompt == "Douglass> " for m in mirepl.interface.modes)
    end

end


text:
using Oceananigans.Architectures: device_event
using Oceananigans.Grids: AbstractGrid
using Oceananigans.Operators: ∂xᶠᶜᵃ, ∂yᶜᶠᵃ
using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions

using Adapt

struct ExplicitFreeSurface{E, T}
    η :: E
    gravitational_acceleration :: T
end

ExplicitFreeSurface(; gravitational_acceleration=g_Earth) =
    ExplicitFreeSurface(nothing, gravitational_acceleration)

Adapt.adapt_structure(to, free_surface::ExplicitFreeSurface) =
    ExplicitFreeSurface(Adapt.adapt(to, free_surface.η), free_surface.gravitational_acceleration)

#####
##### Interface to HydrostaticFreeSurfaceModel
#####

function FreeSurface(free_surface::ExplicitFreeSurface{Nothing}, velocities, arch, grid)
    η = FreeSurfaceDisplacementField(velocities, free_surface, arch, grid)
    g = convert(eltype(grid), free_surface.gravitational_acceleration)
    return ExplicitFreeSurface(η, g)
end

#####
##### Kernel functions for HydrostaticFreeSurfaceModel
#####

@inline explicit_barotropic_pressure_x_gradient(i, j, k, grid, free_surface::ExplicitFreeSurface) =
    free_surface.gravitational_acceleration * ∂xᶠᶜᵃ(i, j, k, grid, free_surface.η)

@inline explicit_barotropic_pressure_y_gradient(i, j, k, grid, free_surface::ExplicitFreeSurface) =
    free_surface.gravitational_acceleration * ∂yᶜᶠᵃ(i, j, k, grid, free_surface.η)

#####
##### Time stepping
#####

ab2_step_free_surface!(free_surface::ExplicitFreeSurface, model, Δt, χ, velocities_update) =
    explicit_ab2_step_free_surface!(free_surface, model, Δt, χ)

explicit_ab2_step_free_surface!(free_surface, model, Δt, χ) =
    launch!(model.architecture, model.grid, :xy,
            _explicit_ab2_step_free_surface!, free_surface.η, Δt, χ,
            model.timestepper.Gⁿ.η, model.timestepper.G⁻.η,
            dependencies = device_event(model.architecture))

#####
##### Kernel
#####

@kernel function _explicit_ab2_step_free_surface!(η, Δt, χ::FT, Gηⁿ, Gη⁻) where FT
    i, j = @index(Global, NTuple)

    @inbounds begin
        η[i, j, 1] += Δt * ((FT(1.5) + χ) * Gηⁿ[i, j, 1] - (FT(0.5) + χ) * Gη⁻[i, j, 1])
    end
end


text:
export read_file

"""
read_file(aws,bicket,path,warn_user)

Read entire SEGY files from object in the AWS S3 bucket.

"""
function read_file(aws::AWSCore.AWSConfig, bucket::String, path::String, warn_user::Bool;
    buffer_size::Int=1024, start_byte::Int=3600, end_byte::Int=0 )

    verbose = false
    debug = false

    # Read and check sizes
    if end_byte==0
        end_byte = s3_get_meta(aws,bucket,path)["Content-Length"]
        end_byte = parse(Int,end_byte)
    end
    end_byte > start_byte || @error "Fatal error: start_byte ($start_byte) < end_byte ($end_byte)"
    data_sz = end_byte - start_byte
    max_buf_sz = buffer_size * (1024^2)

    # Read File Header
    s=IOBuffer(AWSS3SegyIO.my_s3_get_range(aws,bucket,path,0,3599))
    fh = read_fileheader(s)

    # Move to start of block
    seek(s, start_byte)
    debug && println((start_byte,end_byte,length(s.data)))

    # Check datatype of file
    datatype = Float32
    if fh.bfh.DataSampleFormat == 1
        datatype = SegyIO.IBMFloat32
    elseif fh.bfh.DataSampleFormat != 5
        @error "Data type not supported ($(fh.bfh.DataSampleFormat))"
    end

    # Check fixed length trace flag
    (fh.bfh.FixedLengthTraceFlag!=1 & warn_user) && @warn "Fixed length trace flag set in stream: $s"

    ## Check for extended text header

    # Read traces
    trace_byte_len = 240 + fh.bfh.ns*4
    ntraces = Int(data_sz/trace_byte_len)
    max_tr_per_buf = Int(floor(max_buf_sz/trace_byte_len))
    nbuffers = Int(ceil(ntraces/max_tr_per_buf))
    debug && println((ntraces,fh.bfh.ns,trace_byte_len,data_sz,max_buf_sz,max_tr_per_buf,nbuffers))
    buf_sizes, trace_idx = trace_buffer_parts(nbuffers,ntraces)
    debug && println(buf_sizes)
    debug && println(trace_idx)

    # Preallocate memory
    headers = Array{BinaryTraceHeader, 1}(undef, ntraces)
    data = Array{datatype, 2}(undef, fh.bfh.ns, ntraces)
    th_b2s = th_byte2sample()

    # Read buffer
    for buffer=1:nbuffers
        buffer_start = start_byte + (trace_idx[buffer]-1)*trace_byte_len
        buffer_end = start_byte + (trace_idx[buffer+1]-1)*trace_byte_len-1
        s=IOBuffer(AWSS3SegyIO.my_s3_get_range(aws,bucket,path,buffer_start,buffer_end))
        # Read each trace
        for trace in trace_idx[buffer]:trace_idx[buffer+1]-1

            verbose && println("buffer $buffer $buffer_start $buffer_end trace $trace $(position(s)) $(position(s)+trace_byte_len-1)")
            SegyIO.read_trace!(s, fh.bfh, datatype, headers, data, trace, th_b2s)

        end
    end

    return SeisBlock(fh, headers, data)
end

"""
read_file(aws,bicket,path,keys,warn_user)

Read entire SEGY files from object in the AWS S3 bucket,
only reading the header values in 'keys'.

"""
function read_file(aws::AWSCore.AWSConfig, bucket::String, path::String, keys::Array{String,1}, warn_user::Bool;
    buffer_size::Int=1024, start_byte::Int=3600, end_byte::Int=0 )

    verbose = false
    debug = false

    # Read and check sizes
    if end_byte==0
        end_byte = s3_get_meta(aws,bucket,path)["Content-Length"]
        end_byte = parse(Int,end_byte)
    end
    end_byte > start_byte || @error "Fatal error: start_byte ($start_byte) < end_byte ($end_byte)"
    data_sz = end_byte - start_byte
    max_buf_sz = buffer_size * (1024^2)

    # Read File Header
    s=IOBuffer(AWSS3SegyIO.my_s3_get_range(aws,bucket,path,0,3599))
    fh = read_fileheader(s)

    # Move to start of block
    seek(s, start_byte)
    debug && println((start_byte,end_byte,length(s.data)))

    # Check datatype of file
    datatype = Float32
    if fh.bfh.DataSampleFormat == 1
        datatype = SegyIO.IBMFloat32
    elseif fh.bfh.DataSampleFormat != 5
        @error "Data type not supported ($(fh.bfh.DataSampleFormat))"
    end

    # Check fixed length trace flag
    (fh.bfh.FixedLengthTraceFlag!=1 & warn_user) && @warn "Fixed length trace flag set in stream: $s"

    ## Check for extended text header

    # Read traces
    trace_byte_len = 240 + fh.bfh.ns*4
    ntraces = Int(data_sz/trace_byte_len)
    max_tr_per_buf = Int(floor(max_buf_sz/trace_byte_len))
    nbuffers = Int(ceil(ntraces/max_tr_per_buf))
    debug && println((ntraces,fh.bfh.ns,trace_byte_len,data_sz,max_buf_sz,max_tr_per_buf,nbuffers))
    buf_sizes, trace_idx = trace_buffer_parts(nbuffers,ntraces)
    debug && println(buf_sizes)
    debug && println(trace_idx)

    # Preallocate memory
    headers = Array{BinaryTraceHeader, 1}(undef, ntraces)
    data = Array{datatype, 2}(undef, fh.bfh.ns, ntraces)
    th_b2s = th_byte2sample()

    # Read buffer
    for buffer=1:nbuffers
        buffer_start = start_byte + (trace_idx[buffer]-1)*trace_byte_len
        buffer_end = start_byte + (trace_idx[buffer+1]-1)*trace_byte_len-1
        s=IOBuffer(AWSS3SegyIO.my_s3_get_range(aws,bucket,path,buffer_start,buffer_end))
        # Read each trace
        for trace in trace_idx[buffer]:trace_idx[buffer+1]-1

            verbose && println("buffer $buffer $buffer_start $buffer_end trace $trace $(position(s)) $(position(s)+trace_byte_len-1)")
            SegyIO.read_trace!(s, fh.bfh, datatype, headers, data, trace, keys, th_b2s)

        end
    end

    return SeisBlock(fh, headers, data)
end


text:
struct TransmissionRateCache <: AbstractTN
  external::Vector{Union{Nothing, Float64}}
  internal::Array{Union{Nothing, Float64}, 2}

  function TransmissionRateCache(x::Integer)
    return new(Union{Nothing, Float64}[nothing for i = 1:x],
               Union{Nothing, Float64}[nothing for i = 1:x, j = 1:x])
  end
end


text:
module Oulu

using HTTP, DataFrames

greet() = print("Hello World!")

function myfunc(x)
    return x+2
end

function coronavirusupdate()
    r = HTTP.get("https://www.worldometers.info/coronavirus/", cookies=true);
    data = String(r.body)
    println(data[findfirst(r"Coronavirus Update*?(.*) Wuhan China Virus Outbreak"s, data)])
end

coronavirusupdate()

end # module


text:
pkg_list = [
    :BenchmarkTools,
    :StaticArrays,
    :Plots,
    :PlotlyJS,
    :Measures,
    :QuadGK,
    :Unrolled,
    :Setfield,
    :Optim,
    :Traceur,
    :Makie,
    :GLMakie,
    :MakieGallery,
]


text:
# Π(θ) = γ[1] * quality + γ[2] log(y - p)
# C(i;θ) = (β[1] - β[2] * η ) i + β[3] i^2 + β[4] ownQ 1(i>0& ownQ > 0 )

struct Parameter
    γ::Array{Real,1}
    β::Array{Real,1}
    Parameter() = new([1],[0.2, 0.1, 0.2, 0.0])
    Parameter(γ,β) = new(γ,β)
end
#
struct State
    ownQ::Real
    otherQ::Array{Real,1}
    η::Real  #The private shock
    State() = new(sqrt(3) * randn() + 0.5,[],sqrt(2)*randn());
    State(ownQ::Real,otherQ::RealVector,η::Real) = new(ownQ,otherQ,η)
end
#
function profit(s::State,invest::Real,p::Parameter)
    profit = s.ownQ * p.γ[1]
    cost = (invest * (p.β[1] - p.β[2] * s.η) +
    p.β[3] * invest^2 + p.β[4]*(invest>0) * (s.ownQ > 0) * invest * s.ownQ)
    return profit - cost
end
#
mutable struct ProfitFn
    γ::RealVector
    β::RealVector
    p::Parameter
    ProfitFn(p) = new(p.γ, p.β,p) #constructor of a profit function
    (self::ProfitFn)() = 0
    (self::ProfitFn)(k) = println(k)
    (self::ProfitFn)(s,i) = profit(s,i,self.p) #callable struct
    # function profit(s,i) = profit(s,i,p)
end


text:
"""
    Code to compute CCFs
Author: Michael Palumbo
Created: December 2019
Contact: mlp95@psu.edu
Based on code by Alex Wise (aw@psu.edu)
Refactored and optimized by Eric Ford
"""

"""
    `ccf_1D!(ccf_out, λs, fluxes; ccf_plan )`

Compute the cross correlation function of a spectrum with a mask.
    Generalized version that should work with different mask shapes.
# Inputs:
- `ccf_out`: 1-d array of size length(calc_ccf_v_grid(plan)) to store output
- `λs`: 1-d array of wavelengths
- `fluxes`:  1-d array of fluxes
# Optional Arguments:
- `plan`:  parameters for computing ccf (BasicCCFPlan())
"""
function ccf_1D!(ccf_out::A1, λ::A2, flux::A3,
                plan::PlanT = BasicCCFPlan(); projection_workspace::AbstractArray{T4,2} = zeros(length(λ),1),
                assume_sorted::Bool = false ) where {
                T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,1}, T4<:Real,
                PlanT<:AbstractCCFPlan }
    # make sure wavelengths and spectrum are 1D arrays
    @assert ndims(λ) == 1
    @assert ndims(flux) == 1
    @assert length(λ) == length(flux)
    @assert assume_sorted || issorted( plan.line_list.λ )
    v_grid = calc_ccf_v_grid(plan)
    @assert size(ccf_out,1) == length(v_grid)

    if length(plan.line_list) < 1  # If no lines in chunk
        ccf_out .= zero(eltype(A1))
        return ccf_out
    end

    @assert size(projection_workspace,1) == length(λ)

    # loop through each velocity, project the mask and compute ccf at given v
    for i in 1:size(ccf_out,1)
        # project shifted mask on wavelength domain
        doppler_factor = calc_doppler_factor(v_grid[i])
        project_mask!(projection_workspace, λ, plan, shift_factor=doppler_factor)

        # compute the ccf value at the current velocity shift
        ccf_out[i] = plan.allow_nans ? NaNMath.sum(flux .* projection_workspace) : sum(flux .* projection_workspace)
    end
    return ccf_out
end

"""
    `ccf_1D!(ccf_out, ccf_var_out, λs, fluxes, var; ccf_plan )`

Compute the cross correlation function of a spectrum with a mask.
    Generalized version that should work with different mask shapes.
# Inputs:
- `ccf_out`: 1-d array of size length(calc_ccf_v_grid(plan)) to store output
- `ccf_var_out`:  1-d array of size length(calc_ccf_v_grid(plan)) to store output
- `λs`: 1-d array of wavelengths
- `fluxes`:  1-d array of fluxes
- `var`:  1-d array of flux variances
# Optional Arguments:
- `plan`:  parameters for computing ccf (BasicCCFPlan())
"""
function ccf_1D!(ccf_out::A1, ccf_var_out::A1, λ::A2, flux::A3, var::A4,
                plan::PlanT = BasicCCFPlan(); projection_workspace::AbstractArray{T5,2} = zeros(length(λ),1),
                ccf_var_scale::Real = 1.0, assume_sorted::Bool = false ) where {
                T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,1}, T4<:Real, A4<:AbstractArray{T4,1}, T5<:Real,
                PlanT<:AbstractCCFPlan }
    # make sure wavelengths and spectrum are 1D arrays
    @assert ndims(λ) == 1
    @assert ndims(flux) == 1
    @assert length(λ) == length(flux)
    @assert assume_sorted || issorted( plan.line_list.λ )
    v_grid = calc_ccf_v_grid(plan)
    @assert size(ccf_out,1) == size(ccf_var_out,1)
    @assert size(ccf_out,1) == length(v_grid)

    if length(plan.line_list) < 1  # If no lines in chunk
        ccf_out .= zero(eltype(A1))
        ccf_var_out .= zero(eltype(A1))
        return (ccf=ccf_out, ccf_var=ccf_var_out)
    end

    @assert size(projection_workspace,1) == length(λ)
    @assert size(projection_workspace,2) == 1   # Keep until implement CCF w/ multiple masks simultaneously

    # loop through each velocity, project the mask and compute ccf at given v
    for i in 1:size(ccf_out,1)
        # project shifted mask on wavelength domain
        doppler_factor = calc_doppler_factor(v_grid[i])
        project_mask!(projection_workspace, λ, plan, shift_factor=doppler_factor)

        # compute the ccf value at the current velocity shift
        ccf_out[i] = plan.allow_nans ? NaNMath.sum(flux .* projection_workspace) : sum(flux .* projection_workspace)
        ccf_var_out[i] = plan.allow_nans ? NaNMath.sum(var .* projection_workspace.^2) : sum(var .* projection_workspace.^2)
        ccf_var_out[i] *= ccf_var_scale
    end
    return (ccf=ccf_out, ccf_var=ccf_var_out)
end

"""
    `ccf_1D!(ccf_out, ccf_covar_out, λs, fluxes, var; ccf_plan )`

Compute the cross correlation function of a spectrum with a mask.
    Generalized version that should work with different mask shapes.
    WIP:  Generalized version to compute covariance matrix for ccf and optionally to account for LSF
# Inputs:
- `ccf_out`: 1-d array of size length(calc_ccf_v_grid(plan)) to store output
- `ccf_covar_out`:  2-d array of size length(calc_ccf_v_grid(plan))^2 to store output
- `λs`: 1-d array of wavelengths
- `fluxes`:  1-d array of fluxes
- `var`:  1-d array of flux variances
# Optional Arguments:
- `plan`:  parameters for computing ccf (BasicCCFPlan())
"""
function ccf_1D!(ccf_out::A1, ccf_covar_out::A2, λ::A3, flux::A4, var::A5,
                plan::PlanT = BasicCCFPlan(); projection_workspace::AbstractArray{T6,3} = zeros(length(λ),1,size(ccf_out,1)),
                ccf_var_scale::Real = 1.0, assume_sorted::Bool = false, calc_covar_terms::Bool = false ) where {
                T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,2}, T3<:Real, A3<:AbstractArray{T3,1}, T4<:Real, A4<:AbstractArray{T4,1}, T5<:Real, A5<:AbstractArray{T5,1}, T6<:Real,
                PlanT<:AbstractCCFPlan }
    # make sure wavelengths and spectrum are 1D arrays
    @assert ndims(λ) == 1
    @assert ndims(flux) == 1
    @assert length(λ) == length(flux)
    @assert assume_sorted || issorted( plan.line_list.λ )
    v_grid = calc_ccf_v_grid(plan)
    num_vels = length(v_grid)
    @assert size(ccf_out,1) == num_vels
    @assert size(ccf_covar_out,1) == size(ccf_covar_out,2)
    @assert size(ccf_covar_out,1) == num_vels

    if length(plan.line_list) < 1  # If no lines in chunk
        ccf_out .= zero(eltype(A1))
        ccf_covar_out .= zero(eltype(A2))
        return (ccf=ccf_out, ccf_var=ccf_covar_out)
    end

    @assert size(projection_workspace,1) == length(λ)
    @assert size(projection_workspace,2) == 1   # Keep until implement CCF w/ multiple masks simultaneously
    @assert size(projection_workspace,3) == size(ccf_out,1)

    ccf_covar_out .= zero(eltype(A2))
    # loop through each velocity, project the mask and compute ccf at given v
    for i in 1:num_vels
        # project shifted mask on wavelength domain
        doppler_factor = calc_doppler_factor(v_grid[i])
        project_mask!(view(projection_workspace, :, :, i), λ, plan, shift_factor=doppler_factor)

        # compute the ccf value at the current velocity shift
        ccf_out[i] = plan.allow_nans ? NaNMath.sum(flux .* view(projection_workspace,:,1,i) ) : sum(flux .* view(projection_workspace,:,1,i) )
        # computing the variance term of CCF covar's diagonal entries moved to separate loop for memory access locality
        # ccf_covar_out[i,i] = plan.allow_nans ? NaNMath.sum(var .* view(projection_workspace,:,1,i).^2) : sum(var .* view(projection_workspace,:,1,i).^2)
    end
    # Add variance terms for off-diagonal terms of CCF covar
    for i in 1:num_vels
        for j in i:num_vels
            if i == j
                ccf_covar_out[i,i] = plan.allow_nans ? NaNMath.sum(var .* view(projection_workspace,:,1,i).^2) : sum(var .* view(projection_workspace,:,1,i).^2)
            else  # i != j
                covar_term = plan.allow_nans ? NaNMath.sum(var .* view(projection_workspace,:,1,i) .* view(projection_workspace,:,1,j) ) : sum(var .* view(projection_workspace,:,1,i) .* view(projection_workspace,:,1,j) )
                ccf_covar_out[i,j] = covar_term
                ccf_covar_out[j,i] = covar_term
            end  # if i== j ... else ...
        end # for j
    end # for i
    #= =#
    if calc_covar_terms   # If we want to include covar terms due to LSF  # WARNING: SLOW
        pdf_norm(x,σ²) = exp(-0.5*x^2/σ²)/sqrt(2π*σ²)
        lsf_σ = 2400.0
        ccf_off_diag_σ_factor = 2.0

        for k in 1:length(λ)
            for l in (k+1):length(λ)
                #if k == l    continue   end     #  l starts at k+1 anyway; variance terms computed above
                Δv = log(λ[k]/λ[l])*speed_of_light_mps
                if abs(Δv) > (ccf_off_diag_σ_factor * lsf_σ)   continue    end
                flux_covar_norm = pdf_norm(Δv,2*lsf_σ^2)
                ave_var_λ = sqrt(var[k] * var[l])    # TODO: Improve this approximation based on flux var at v_mean?
                prefactor = flux_covar_norm * ave_var_λ
                for i in 1:num_vels
                    for j in i:num_vels
                        covar_term = prefactor * projection_workspace[k,1,i] * projection_workspace[l,1,j]
                        if i == j         # diagnonal terms of covar
                            ccf_covar_out[i,i] += 2*covar_term
                        else # i != j     # off diagonal terms of covar
                            ccf_covar_out[i,j] += 2*covar_term
                            ccf_covar_out[j,i] += 2*covar_term
                        end   # if i==j ... else ...
                    end # for j
                end # for i
            end # for l
        end  # for k
    else  #  !calc_covar_terms
        ccf_covar_out .*= ccf_var_scale
    end # if false
    return (ccf=ccf_out, ccf_covar=ccf_covar_out )
end

"""
    `ccf_1D( λs, fluxes; ccf_plan )`
    Compute the cross correlation function of a spectrum with a mask.
# Inputs:
- `λs`: 1-d array of wavelengths
- `fluxes`:  1-d array of fluxes
# Optional Arguments:
- `ccf_plan`:  parameters for computing ccf (BasicCCFPlan())
# Returns:
- 1-d array of size length(calc_ccf_v_grid(plan))
"""
function ccf_1D(λ::A2, flux::A3,
                plan::PlanT = BasicCCFPlan() ) where {
                T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,1},
                PlanT<:AbstractCCFPlan }
    @assert ndims(λ) == 1
    @assert ndims(flux) == 1
    @assert length(λ) == length(flux)
    len_v_grid = calc_length_ccf_v_grid(plan)
    ccf_out = zeros(len_v_grid)
    if length(plan.line_list) < 1  # If no lines in chunk
        return ccf_out
    end
    ccf_1D!(ccf_out, λ, flux, plan )
    return ccf_out
end

"""
    `ccf_1D( λs, fluxes; ccf_plan )`
    Compute the cross correlation functions of spectra with a mask.
# Inputs:
- `λs`: 1-d array of wavelengths
- `fluxes`:  2-d array of fluxes, individual spectra along first dim
# Optional Arguments:
- `ccf_plan`:  parameters for computing ccf (BasicCCFPlan())
# Returns:
- 2-d array of size (length(calc_ccf_v_grid(plan)), size(flux, 2))
"""
function ccf_1D(λ::A2, flux::A3,
                plan::PlanT = BasicCCFPlan() ) where {
                T2<:Real, A2<:AbstractArray{T2,1},
                T3<:Real, A3<:AbstractArray{T3,2},
                PlanT<:AbstractCCFPlan }
    @assert ndims(λ) == 1
    @assert ndims(flux) == 2
    @assert length(λ) == size(flux, 1)
    len_v_grid = calc_length_ccf_v_grid(plan)
    ccf_out = zeros(len_v_grid, size(flux, 2))
    if length(plan.line_list) < 1  # If no lines in chunk
        return ccf_out
    end

    for i in 1:size(flux, 2)
        ccf_1D!(view(ccf_out, :, i), λ, view(flux, :, i), plan)
    end
    return ccf_out
end

"""
    `ccf_1D( λs, fluxes, vars; ccf_plan )`

Compute the cross correlation function of a spectrum with a mask and its variance^1.
^1 = This version computes the diagonal terms for the ccf variance and neglects covariances due to the line spread function.
Can roughly compensate by scaling up the ccf_var_scale from the default of 1.
WARNING: Still needs more testing.
# Inputs:
- `λs`: 1-d array of wavelengths
- `fluxes`:  1-d array of fluxes
# Optional Arguments:
- `ccf_plan`:  parameters for computing ccf (BasicCCFPlan())
# Returns:
- 1-d array of size length(calc_ccf_v_grid(plan))
"""
function ccf_1D(λ::A2, flux::A3, var::A4,
                plan::PlanT = BasicCCFPlan(); ccf_var_scale::Real = 1.0,
                calc_ccf_covar::Bool = false, calc_covar_terms::Bool = false ) where {
                T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,1}, T4<:Real, A4<:AbstractArray{T4,1},
                PlanT<:AbstractCCFPlan }
    @assert ndims(λ) == 1
    @assert ndims(flux) == 1
    @assert length(λ) == length(flux)
    len_v_grid = calc_length_ccf_v_grid(plan)
    ccf_out = zeros(len_v_grid)
    ccf_covar_out = calc_ccf_covar ? zeros(len_v_grid,len_v_grid) : zeros(len_v_grid)
    if length(plan.line_list) < 1  # If no lines in chunk
        return (ccf=ccf_out, ccf_var=ccf_covar_out)
    end
    if calc_ccf_covar
        ccf_1D!(ccf_out, ccf_covar_out, λ, flux, var, plan, ccf_var_scale=ccf_var_scale, calc_covar_terms=calc_covar_terms )
        return (ccf=ccf_out, ccf_covar=ccf_covar_out)
    else
        ccf_1D!(ccf_out, ccf_covar_out, λ, flux, var, plan, ccf_var_scale=ccf_var_scale )
        return (ccf=ccf_out, ccf_var=ccf_covar_out)
    end
end

""" `project_mask!( output, λs, ccf_plan; shift_factor )`

Compute the projection of the mask onto the 1D array of wavelengths (λs) at a given shift factor (default: 1).
The mask is computed from the ccf_plan, including a line_list and mask_shape (default: tophat).
Assumes plan.line_list is already sorted.

"""
function project_mask!(projection::A2, λs::A1, plan::PlanT ; shift_factor::Real=one(T1)
                            ) where { T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,2}, PlanT<:AbstractCCFPlan }
    num_lines_in_mask = length(plan.line_list)
    @assert num_lines_in_mask >= 1
    @assert size(projection,1) == length(λs)
    @assert size(projection,2) == 1         # Keep until implement CCF w/ multiple masks simultaneously
    projection .= zero(eltype(projection))  # Important to zero out projection before accumulating there!

    # set indexing variables
    p = 1
    p_left_edge_of_current_line = 1
    λsle_cur = λs[p]-0.5*(λs[p+1]-λs[p])   # Current pixel's left edge
    λsre_cur = 0.5*(λs[p]+λs[p+1])   # Current pixel's left edge
    #@assert issorted( plan.line_list.λ )  # Moved outside of inner loop
    m = searchsortedfirst(plan.line_list.λ, λsle_cur/shift_factor)
    while m <= length(plan.line_list.λ) && λ_min(plan.mask_shape,plan.line_list.λ[m]) * shift_factor < λsle_cur   # only includemask entry if it fit entirely in chunk
        m += 1
    end
    if m > length(plan.line_list.λ)   # Return early if no lines fall within chunk's lambda's
        return projection
    else
        #println("# Starting with mask entry at λ=", plan.line_list.λ[m], " for chunk with λ= ",first(λs)," - ",last(λs))
    end
    on_mask = false
    mask_mid = plan.line_list.λ[m] * shift_factor
    mask_hi = λ_max(plan.mask_shape,plan.line_list.λ[m]) * shift_factor   # current line's upper limit
    mask_lo = λ_min(plan.mask_shape,plan.line_list.λ[m]) * shift_factor   # current line's lower limit

    c_mps = speed_of_light_mps
    mask_weight = plan.line_list.weight[m]
    # loop through lines in mask, weight mask by fraction of PSF falling in each wavelength bin and divide by Δlnλ for pixel
    while m <= num_lines_in_mask
        if !on_mask
            if λsre_cur > mask_lo
                if λsre_cur > mask_hi   # Pixel overshot this mask entry, so weight based on mask filling only a portion of the pixel,
                    #For Tophat: projection[p,1] += mask_weight * (mask_hi - mask_lo) / (λsre_cur - λsle_cur)
                    #frac_of_psf_in_v_pixel = 1.0
                    one_over_delta_z_pixel = 0.5*(λsre_cur + λsle_cur) / (λsre_cur - λsle_cur)
                    projection[p,1] += one_over_delta_z_pixel * mask_weight
                    m += 1
                    if m<=length(plan.line_list)      # Move to next line
                        mask_mid = plan.line_list.λ[m] * shift_factor
                        mask_lo = λ_min(plan.mask_shape,plan.line_list.λ[m]) * shift_factor
                        mask_hi = λ_max(plan.mask_shape,plan.line_list.λ[m]) * shift_factor
                        mask_weight = plan.line_list.weight[m]
                    else         # We're out of lines, so exit
                        break
                    end
                else                       # Right edge of current pixel has entered the mask for this line, but left edge hasn't
                    #For Tophat: projection[p,1] += mask_weight * (λsre_cur - mask_lo) / (λsre_cur - λsle_cur)
                    frac_of_psf_in_v_pixel = integrate(plan.mask_shape, (mask_lo-mask_mid)/mask_mid*c_mps, (λsre_cur-mask_mid)/mask_mid*c_mps)
                    one_over_delta_z_pixel = 0.5*(λsre_cur + λsle_cur) / (λsre_cur - λsle_cur)
                    projection[p,1] += frac_of_psf_in_v_pixel * one_over_delta_z_pixel * mask_weight
                    on_mask = true         # Indicate next pixel will hav something to contribute based on the current line
                    p_left_edge_of_current_line = p                 # Mark the starting pixel of the line in case the lines overlap
                    p+=1                   # Move to next pixel
                    λsle_cur = λsre_cur   # Current pixel's left edge
                    λsre_cur = p<size(projection,1) ? 0.5*(λs[p]+λs[p+1]) : λs[p]+0.5*(λs[p]-λs[p-1])    # Current pixel's right edge

                end
            else   # ALl of pixel is still to left of beginning of mask for this line.
                p+=1    # TODO: Opt try to accelerate with something like searchsortedfirst?  Probably not worth it.
                λsle_cur = λsre_cur   # Current pixel's left edge
                λsre_cur = p<size(projection,1) ? 0.5*(λs[p]+λs[p+1]) : λs[p]+0.5*(λs[p]-λs[p-1])    # Current pixel's right edge
            end
        else  # on_mask == true
            if λsre_cur > mask_hi               # Right edge of this pixel moved past the edge of the current line's mask
                #For Tophat: projection[p,1] += mask_weight * (mask_hi - λsle_cur) / (λsre_cur - λsle_cur)
                frac_of_psf_in_v_pixel = integrate(plan.mask_shape, (λsle_cur-mask_mid)*c_mps/mask_mid, (mask_hi-mask_mid)*c_mps/mask_mid)
                one_over_delta_z_pixel = 0.5*(λsre_cur + λsle_cur) / (λsre_cur - λsle_cur)
                projection[p,1] += frac_of_psf_in_v_pixel * one_over_delta_z_pixel * mask_weight
                on_mask = false                 # Indicate that we're done with this line
                m += 1                          # Move to next line
                if m<=length(plan.line_list)
                    mask_mid = plan.line_list.λ[m] * shift_factor
                    mask_lo = λ_min(plan.mask_shape,plan.line_list.λ[m]) * shift_factor
                    mask_hi = λ_max(plan.mask_shape,plan.line_list.λ[m]) * shift_factor
                    mask_weight = plan.line_list.weight[m]
                    if mask_lo < λsle_cur       # If the lines overlap, we may have moved past the left edge of the new line. In that case go back to the left edge of the previous line.
                        p = p_left_edge_of_current_line
                        λsle_cur = p>1 ? 0.5*(λs[p]+λs[p-1]) : λs[p]-0.5*(λs[p+1]-λs[p])   # Current pixel's left edge
                        λsre_cur = 0.5*(λs[p]+λs[p+1])   # Current pixel's right edge
                    end

                else                            # We're done with all lines, can return early
                    break
                end
            else                                # This pixel is entirely within the mask
                # assumes plan.mask_shape is normalized to integrate to unity and flux is constant within pixel
                #For Tophat: projection[p,1] += mask_weight    # Add the full mask weight
                frac_of_psf_in_v_pixel = integrate(plan.mask_shape, (λsle_cur-mask_mid)/mask_mid*c_mps, (λsre_cur-mask_mid)/mask_mid*c_mps)
                one_over_delta_z_pixel = 0.5*(λsre_cur + λsle_cur) / (λsre_cur - λsle_cur)
                projection[p,1] += frac_of_psf_in_v_pixel * one_over_delta_z_pixel * mask_weight
                p += 1                          # move to next pixel
                λsle_cur = λsre_cur   # Current pixel's left edge
                λsre_cur = p<size(projection,1) ? 0.5*(λs[p]+λs[p+1]) : λs[p]+0.5*(λs[p]-λs[p-1])    # Current pixel's right edge
            end
        end
        if p>=size(projection,1) break end       # We're done with all pixels in this chunk.
    end
    projection ./= c_mps
    return projection
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
## Prandtl-Glauert transformation
#======================================================#

prandtl_glauert_scale_coordinates(x, y, z, β) = SVector(x, β * y, β * z)
prandtl_glauert_scale_coordinates(r, β) = @views prandtl_glauert_scale_coordinates(r[1], r[2], r[3], β)

prandtl_glauert_scale_normal(nx, ny, nz, β) = SVector(β * nx, ny, nz)
prandtl_glauert_scale_normal(n, β) = @views prandtl_glauert_scale_normal(n[1], n[2], n[3], β)

prandtl_glauert_scale_coordinates(horseshoe :: Horseshoe, β) =
    setproperties(horseshoe,
                  r1                = prandtl_glauert_scale_coordinates(horseshoe.r1, β),
                  r2                = prandtl_glauert_scale_coordinates(horseshoe.r2, β),
                  collocation_point = prandtl_glauert_scale_coordinates(horseshoe.collocation_point, β),
                  normal            = prandtl_glauert_scale_normal(horseshoe.normal, β),
                )

prandtl_glauert_inverse_scale_velocity(vx, vy, vz, β) = SVector(vx / β^2, vy / β, vz / β)
prandtl_glauert_inverse_scale_velocity(v, β) = @views prandtl_glauert_inverse_scale_velocity(v[1], v[2], v[3], β)


text: println(maximum(filter(x -> string(x) == reverse(string(x)), reduce(vcat, map(x -> map(y -> x * y, 100:999), 100:999), init=[]))))

max_stars_repo_name: seaneshbaugh/rosetta-euler

max_issues_repo_name: seaneshbaugh/rosetta-euler

max_forks_repo_name: seaneshbaugh/rosetta-euler


text:
using Printf
using LinearAlgebra
using Random

include("../common/PWGrid_v02.jl")
include("../common/ortho_gram_schmidt.jl")
include("../common/wrappers_fft.jl")
include("../common/gen_lattice_pwscf.jl")
include("EnergiesT.jl")
include("PotentialsT.jl")
include("op_K.jl")
include("op_V_loc.jl")
include("op_H.jl")
include("calc_rho.jl")
include("calc_grad_v2.jl")
include("calc_Energies.jl")
include("KS_solve_Emin_sd.jl")
include("KS_solve_Emin_cg.jl")
include("Poisson_solve.jl")
include("LDA_VWN.jl")
include("Kprec.jl")

include("../common/gen_dr_center.jl")
include("../common/calc_strfact_v2.jl")
include("../common/calc_ewald_v2.jl")

include("diag_lobpcg.jl")
include("KS_solve_SCF.jl")

function test_main( Ns )

    LatVecs = gen_lattice_sc(16.0)

    pw = PWGrid( Ns, LatVecs )

    Ω  = pw.Ω
    r  = pw.r
    G  = pw.gvec.G
    G2 = pw.gvec.G2
    Npoints = prod(Ns)
    Ngwx = pw.gvecw.Ngwx

    @printf("Ns   = (%d,%d,%d)\n", Ns[1], Ns[2], Ns[3])
    @printf("Ngwx = %d\n", Ngwx)

    actual = Npoints/Ngwx
    theor = 1/(4*pi*0.25^3/3)
    @printf("Compression: actual, theor: %f , %f\n", actual, theor)

    # Helium atom
    Zion = 4.0
    Nstates = 2
    Focc = [2.0, 2.0]

    Xpos = reshape( [1.0, 1.0, 1.0], (3,1) )

    Sf = calc_strfact( Xpos, 1, [1], pw.gvec.G )
    E_nn = calc_ewald( pw, Sf, Xpos, 1, [1], [Zion] )

    # Set up potential, using HGH pseudopotential for H
    # contains only local pseudopotential
    rloc = 0.325000
    c1 = -24.015041
    c2 = 17.204014
    c3 = -3.326390
    c4 = 0.165419
    #
    Vg = zeros(ComplexF64,Npoints)
    pre1 = -4*pi*Zion/Ω
    pre2 = sqrt(8*pi^3)*rloc^3/Ω
    #
    for ig=2:Npoints
        Gr = sqrt(G2[ig])*rloc
        expGr2 = exp(-0.5*Gr^2)
        Vg[ig] = pre1/G2[ig]*expGr2 + pre2*expGr2 * (c1 + c2*(3-Gr^2) +
                 c3*(15 - 10*Gr^2 + Gr^4) + c4*(105 - 105*Gr^2 + 21*Gr^4 - Gr^6) )
    end
    # limiting value
    Vg[1] = 2*pi*Zion*rloc^2 + (2*pi)^1.5 * rloc^3 * (c1 + 3.0*c2 + 15*c3 + 105*c4)
    println("Vg[1] = ", Vg[1])
    V_ionic = real( G_to_R(Ns, Vg .* Sf) ) * Npoints

    # Need to sum up over Nspecies for more than one species type
    # We simply need reshape because we only have one species type here.
    V_ionic = reshape( V_ionic, (Npoints) )

    println("sum(V_ionic) = ", sum(V_ionic))

    psi, Energies, Potentials = KS_solve_Emin_cg( pw, V_ionic, Focc, Nstates, NiterMax=1000, E_NN = E_nn )

    Y = ortho_gram_schmidt(psi)
    mu = Y' * op_H( pw, Potentials, Y )
    evals, evecs = eigen(mu)
    psi = Y*evecs

    #@printf("Solution by self-consistent field method\n")
    #Energies, Potentials, psi, evals = KS_solve_SCF( pw, V_ionic, Focc, Nstates, β=0.3 )

    for ist = 1:Nstates
        @printf("State # %d, Energy = %18.10f\n", ist, real(evals[ist]))
    end

    print_Energies(Energies)

end

@time test_main( [54,54,54] )
#@time test_main( [64,64,64] )
#@time test_main( [100,100,100] )