text:
function _ensure_not_running_on_functor(f, name)
    # if x itself is a Type, then it is a constructor, thus not a functor.
    # This also catchs UnionAll constructors which have a `:var` and `:body` fields
    f isa Type && return

    if fieldcount(typeof(f)) > 0
        throw(ArgumentError(
            "$name cannot be used on closures/functors (such as $f)"
        ))
    end
end

"""
    _wrap_function(f, xs, ignores)

Return a new version of `f`, `fnew`, that ignores some of the arguments `xs`.

# Arguments
- `f`: The function to be wrapped.
- `xs`: Inputs to `f`, such that `y = f(xs...)`.
- `ignores`: Collection of `Bool`s, the same length as `xs`.
  If `ignores[i] === true`, then `xs[i]` is ignored; `∂xs[i] === nothing`.
"""
function _wrap_function(f, xs, ignores)
    function fnew(sigargs...)
        callargs = Any[]
        j = 1

        for (i, (x, ignore)) in enumerate(zip(xs, ignores))
            if ignore
                push!(callargs, x)
            else
                push!(callargs, sigargs[j])
                j += 1
            end
        end
        @assert j == length(sigargs) + 1
        @assert length(callargs) == length(xs)
        return f(callargs...)
    end
    return fnew
end

"""
    _make_j′vp_call(fdm, f, ȳ, xs, ignores) -> Tuple

Call `FiniteDifferences.j′vp`, with the option to ignore certain `xs`.

# Arguments
- `fdm::FiniteDifferenceMethod`: How to numerically differentiate `f`.
- `f`: The function to differentiate.
- `ȳ`: The adjoint w.r.t. output of `f`.
- `xs`: Inputs to `f`, such that `y = f(xs...)`.
- `ignores`: Collection of `Bool`s, the same length as `xs`.
  If `ignores[i] === true`, then `xs[i]` is ignored; `∂xs[i] === nothing`.

# Returns
- `∂xs::Tuple`: Derivatives estimated by finite differencing.
"""
function _make_j′vp_call(fdm, f, ȳ, xs, ignores)
    f2 = _wrap_function(f, xs, ignores)

    ignores = collect(ignores)
    args = Any[nothing for _ in 1:length(xs)]
    all(ignores) && return (args...,)
    sigargs = xs[.!ignores]
    arginds = (1:length(xs))[.!ignores]
    fd = j′vp(fdm, f2, ȳ, sigargs...)
    @assert length(fd) == length(arginds)

    for (dx, ind) in zip(fd, arginds)
        args[ind] = dx
    end
    return (args...,)
end

"""
    _make_jvp_call(fdm, f, xs, ẋs, ignores)

Call `FiniteDifferences.jvp`, with the option to ignore certain `xs`.

# Arguments
- `fdm::FiniteDifferenceMethod`: How to numerically differentiate `f`.
- `f`: The function to differentiate.
- `xs`: Inputs to `f`, such that `y = f(xs...)`.
- `ẋs`: The directional derivatives of `xs` w.r.t. some real number `t`.
- `ignores`: Collection of `Bool`s, the same length as `xs` and `ẋs`.
   If `ignores[i] === true`, then `ẋs[i]` is ignored for derivative estimation.

# Returns
- `Ω̇`: Derivative of output w.r.t. `t` estimated by finite differencing.
"""
function _make_jvp_call(fdm, f, xs, ẋs, ignores)
    f2 = _wrap_function(f, xs, ignores)

    ignores = collect(ignores)
    all(ignores) && return ntuple(_->nothing, length(xs))
    sigargs = zip(xs[.!ignores], ẋs[.!ignores])
    return jvp(fdm, f2, sigargs...)
end

"""
    test_scalar(f, z; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), fkwargs=NamedTuple(), kwargs...)

Given a function `f` with scalar input and scalar output, perform finite differencing checks,
at input point `z` to confirm that there are correct `frule` and `rrule`s provided.

# Arguments
- `f`: Function for which the `frule` and `rrule` should be tested.
- `z`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).

`fkwargs` are passed to `f` as keyword arguments.
All keyword arguments except for `fdm` and `fkwargs` are passed to `isapprox`.
"""
function test_scalar(f, z; rtol=1e-9, atol=1e-9, fdm=_fdm, fkwargs=NamedTuple(), kwargs...)
    _ensure_not_running_on_functor(f, "test_scalar")
    # z = x + im * y
    # Ω = u(x, y) + im * v(x, y)
    Ω = f(z; fkwargs...)

    # test jacobian using forward mode
    Δx = one(z)
    @testset "$f at $z, with tangent $Δx" begin
        # check ∂u_∂x and (if Ω is complex) ∂v_∂x via forward mode
        frule_test(f, (z, Δx); rtol=rtol, atol=atol, fdm=fdm, fkwargs=fkwargs, kwargs...)
        if z isa Complex
            # check that same tangent is produced for tangent 1.0 and 1.0 + 0.0im
            @test isapprox(
                frule((Zero(), real(Δx)), f, z; fkwargs...)[2],
                frule((Zero(), Δx), f, z; fkwargs...)[2],
                rtol=rtol,
                atol=atol,
                kwargs...,
            )
        end
    end
    if z isa Complex
        Δy = one(z) * im
        @testset "$f at $z, with tangent $Δy" begin
            # check ∂u_∂y and (if Ω is complex) ∂v_∂y via forward mode
            frule_test(f, (z, Δy); rtol=rtol, atol=atol, fdm=fdm, fkwargs=fkwargs, kwargs...)
        end
    end

    # test jacobian transpose using reverse mode
    Δu = one(Ω)
    @testset "$f at $z, with cotangent $Δu" begin
        # check ∂u_∂x and (if z is complex) ∂u_∂y via reverse mode
        rrule_test(f, Δu, (z, Δx); rtol=rtol, atol=atol, fdm=fdm, fkwargs=fkwargs, kwargs...)
        if Ω isa Complex
            # check that same cotangent is produced for cotangent 1.0 and 1.0 + 0.0im
            back = rrule(f, z)[2]
            @test isapprox(
                extern(back(real(Δu))[2]),
                extern(back(Δu)[2]),
                rtol=rtol,
                atol=atol,
                kwargs...,
            )
        end
    end
    if Ω isa Complex
        Δv = one(Ω) * im
        @testset "$f at $z, with cotangent $Δv" begin
            # check ∂v_∂x and (if z is complex) ∂v_∂y via reverse mode
            rrule_test(f, Δv, (z, Δx); rtol=rtol, atol=atol, fdm=fdm, fkwargs=fkwargs, kwargs...)
        end
    end
end

"""
    frule_test(f, (x, ẋ)...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), fkwargs=NamedTuple(), kwargs...)

# Arguments
- `f`: Function for which the `frule` should be tested.
- `x`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).
- `ẋ`: differential w.r.t. `x` (should generally be set randomly).

`fkwargs` are passed to `f` as keyword arguments.
All keyword arguments except for `fdm` and `fkwargs` are passed to `isapprox`.
"""
function frule_test(f, xẋs::Tuple{Any, Any}...; rtol=1e-9, atol=1e-9, fdm=_fdm, fkwargs=NamedTuple(), kwargs...)
    _ensure_not_running_on_functor(f, "frule_test")
    xs, ẋs = first.(xẋs), last.(xẋs)
    Ω_ad, dΩ_ad = frule((NO_FIELDS, ẋs...), f, xs...; fkwargs...)
    Ω = f(xs...; fkwargs...)
    # if equality check fails, check approximate equality
    # use collect so can do vector equality
    # TODO: add isapprox replacement that works for more types
    @test Ω_ad == Ω || isapprox(collect(Ω_ad), collect(Ω); rtol=rtol, atol=atol)

    ẋs_is_ignored = ẋs .== nothing
    # Correctness testing via finite differencing.
    dΩ_fd = _make_jvp_call(fdm, (xs...) -> f(xs...; fkwargs...), xs, ẋs, ẋs_is_ignored)
    @test isapprox(
        collect(extern.(dΩ_ad)),  # Use collect so can use vector equality
        collect(dΩ_fd);
        rtol=rtol,
        atol=atol,
        kwargs...
    )

    # No tangent is passed in to test accumlation, so generate one
    # See: https://github.com/JuliaDiff/ChainRulesTestUtils.jl/issues/66
    acc = rand_tangent(Ω)
    _check_add!!_behavour(acc, dΩ_ad; rtol=rtol, atol=atol, kwargs...)
end

"""
    rrule_test(f, ȳ, (x, x̄)...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), fkwargs=NamedTuple(), kwargs...)

# Arguments
- `f`: Function to which rule should be applied.
- `ȳ`: adjoint w.r.t. output of `f` (should generally be set randomly).
  Should be same structure as `f(x)` (so if multiple returns should be a tuple)
- `x`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).
- `x̄`: currently accumulated adjoint (should generally be set randomly).

`fkwargs` are passed to `f` as keyword arguments.
All keyword arguments except for `fdm` and `fkwargs` are passed to `isapprox`.
"""
function rrule_test(f, ȳ, xx̄s::Tuple{Any, Any}...; rtol=1e-9, atol=1e-9, fdm=_fdm, fkwargs=NamedTuple(), kwargs...)
    _ensure_not_running_on_functor(f, "rrule_test")

    # Check correctness of evaluation.
    xs = first.(xx̄s)
    x̄s_acc = last.(xx̄s)
    y_ad, pullback = rrule(f, xs...; fkwargs...)
    y = f(xs...; fkwargs...)
    # if equality check fails, check approximate equality
    # use collect so can do vector equality
    # TODO: add isapprox replacement that works for more types
    @test y_ad == y || isapprox(collect(y_ad), collect(y); rtol=rtol, atol=atol)
    @assert !(isa(ȳ, Thunk))

    ∂s = pullback(ȳ)
    ∂self = ∂s[1]
    x̄s_ad = ∂s[2:end]
    @test ∂self === NO_FIELDS  # No internal fields

    x̄s_is_dne = x̄s_acc .== nothing
    # Correctness testing via finite differencing.
    x̄s_fd = _make_j′vp_call(fdm, (xs...) -> f(xs...; fkwargs...), ȳ, xs, x̄s_is_dne)
    for (x̄_acc, x̄_ad, x̄_fd) in zip(x̄s_acc, x̄s_ad, x̄s_fd)
        if  x̄_acc === nothing
            @assert x̄_fd === nothing  # this is how `_make_j′vp_call` works
            @test x̄_ad isa DoesNotExist  # we said it wasn't differentiable.
        else
            # The main test of the actual deriviative being correct:
            @test isapprox(x̄_ad, x̄_fd; rtol=rtol, atol=atol, kwargs...)

            _check_add!!_behavour(x̄_acc, x̄_ad; rtol=rtol, atol=atol, kwargs...)
        end
    end

    if count(!, x̄s_is_dne) == 1
        # for functions with pullbacks that only produce a single non-DNE adjoint, that
        # single adjoint should not be `Thunk`ed. InplaceableThunk is fine.
        i = findfirst(!, x̄s_is_dne)
        @test !(isa(x̄s_ad[i], Thunk))
    end
end


text:
function swap_pbs!(pbs :: Vector{RelaxPb}, nb_stages)
    if length(pbs) != nb_stages+1
        to_add = (nb_stages + 1) - length(pbs)
        push!(pbs, [copy(pbs[1]) for _ in 1:to_add]... )
    end
    pbs[1], pbs[2:end] = pbs[end], pbs[1:end-1]
end

function imex_bdf1!(pbs :: Vector{RelaxPb}, dt, dx, ε)
    swap_pbs!(pbs,1)

    local Θ = 1.0 / (ε^2 + dt)
    local dx⁻¹ = 1.0 / dx

    ## Compute vⁿ⁺¹
    compute_dxP!(pbs[2], dx⁻¹)
    @. pbs[1].V.data = (1.0 - dt * Θ) * pbs[2].V.data +
                    dt * Θ * (pbs[2].f(pbs[2].U.data) - pbs[2].dxP )

    ## Compute uⁿ⁺¹
    compute_dxV!(pbs[1], dx⁻¹)
    @. pbs[1].U.data = pbs[2].U.data - dt * pbs[1].dxV

    nothing
end

function imex_bdf2!(pbs :: Vector{RelaxPb}, dt, dx, ε)
    if length(pbs) == 1
        push!(pbs, copy(pbs[1]), copy(pbs[1]))
        imex_bdf1!(pbs, dt, dx, ε)
    else
        swap_pbs!(pbs,2)
        local a0, a1 = 1.333333333333333333, -0.333333333333333333
        local b0, b1 = 1.333333333333333333, -0.666666666666666667
        local c₋₁ = 0.666666666666666667

        local Θ = 1.0 / (ε^2 + c₋₁ * dt)
        local dx⁻¹ = 1.0 / dx

        compute_dxP!(pbs[2], dx⁻¹) # assume pbs[3].dxP is correct
        local U1, U2, U3 = pbs[1].U.data, pbs[2].U.data, pbs[3].U.data
        local V1, V2, V3 = pbs[1].V.data, pbs[2].V.data, pbs[3].V.data
        @. V1 = (1.0 - c₋₁ * dt * Θ) * (a0*V2 + a1*V3) +
                dt * Θ * (b0 * (pbs[2].f(U2) - pbs[2].dxP ) +
                            b1 * (pbs[3].f(U3) - pbs[3].dxP ) )

        compute_dxV!(pbs[1], dx⁻¹)
        @. U1 = a0*U2 + a1*U3 - c₋₁ * dt * pbs[1].dxV
    end

    nothing
end

function imex_bdf4!(pbs :: Vector{RelaxPb}, dt, dx, ε)
    if length(pbs) == 1
        imex_bdf1!(pbs, dt, dx, ε)
    elseif length(pbs) == 2
        imex_bdf2!(pbs, dt, dx, ε)
    else
        swap_pbs!(pbs,3)
        local a0, b0 =  1.636363636363636364,  1.636363636363636364
        local a1, b1 = -0.818181818181818182, -1.636363636363636364
        local a2, b2 =  0.181818181818181818,  0.545454545454545455
        local c₋₁ = 0.545454545454545455

        local Θ = 1.0 / (ε^2 + c₋₁ * dt)
        local dx⁻¹ = 1.0 / dx

        compute_dxP!(pbs[2], dx⁻¹) # assume pbs[3+].dxP is correct
        local U1, V1 = pbs[1].U.data, pbs[1].V.data
        local U2, V2 = pbs[2].U.data, pbs[2].V.data
        local U3, V3 = pbs[3].U.data, pbs[3].V.data
        local U4, V4 = pbs[4].U.data, pbs[4].V.data
        @. V1 = (1.0 - c₋₁ * dt * Θ) * (a0*V2 + a1*V3 + a2*V4) +
                dt * Θ * (b0 * (pbs[2].f(U2) - pbs[2].dxP ) +
                            b1 * (pbs[3].f(U3) - pbs[3].dxP ) +
                            b2 * (pbs[4].f(U4) - pbs[4].dxP))

        compute_dxV!(pbs[1], dx⁻¹)
        @. U1 = a0*U2 + a1*U3 + a2*U4 - c₋₁ * dt * pbs[1].dxV
    end

    nothing
end

function imex_bdf3!(pbs :: Vector{RelaxPb}, dt, dx, ε)
    if length(pbs) == 1
        imex_bdf1!(pbs, dt, dx, ε)
    elseif length(pbs) == 2
        imex_bdf2!(pbs, dt, dx, ε)
    elseif length(pbs) == 3
        imex_bdf3!(pbs, dt, dx, ε)
    else
        swap_pbs!(pbs,4)
        local a0, b0 =  1.92,  1.92
        local a1, b1 = -1.44, -2.88
        local a2, b2 =  0.64,  1.92
        local a3, b3 =  0.12, -0.48
        local c₋₁ = 0.48

        local Θ = 1.0 / (ε^2 + c₋₁ * dt)
        local dx⁻¹ = 1.0 / dx

        compute_dxP!(pbs[2], dx⁻¹) # assume pbs[3+].dxP is correct
        local U1, V1 = pbs[1].U.data, pbs[1].V.data
        local U2, V2 = pbs[2].U.data, pbs[2].V.data
        local U3, V3 = pbs[3].U.data, pbs[3].V.data
        local U4, V4 = pbs[4].U.data, pbs[4].V.data
        local U5, V5 = pbs[5].U.data, pbs[5].V.data
        @. V1 = (1.0 - c₋₁ * dt * Θ) * (a0*V2 + a1*V3 + a2*V4 + a3*V5) +
                dt * Θ * (b0 * (pbs[2].f(U2) - pbs[2].dxP ) +
                            b1 * (pbs[3].f(U3) - pbs[3].dxP ) +
                            b2 * (pbs[4].f(U4) - pbs[4].dxP ) +
                            b3 * (pbs[5].f(U5) - pbs[5].dxP))

        compute_dxV!(pbs[1], dx⁻¹)
        @. U1 = a0*U2 + a1*U3 + a2*U4 - c₋₁ * dt * pbs[1].dxV
    end

    nothing
end

ext: jl

max_stars_repo_name: tremelow/StiffKinetic

max_issues_repo_name: tremelow/StiffKinetic

max_forks_repo_name: tremelow/StiffKinetic


text:
immutable PBCSCarrier <: QuasiRandomCarrier
  c::Float64 # center of underlying circle
  len::Int
end

PBCSCarrier(len::Int) = PBCSCarrier(0.0, len)
PBCSCarrier() = PBCSCarrier(0.0, 0)
PBCSCarrier(; c::Float64=0.0, len::Int=5) = PBCSCarrier(c, len)

function invcdf_pbcs(c::PBCSCarrier, x::Float64)
  if ((0.0<=x) && (x<0.5))
    -sqrt(1-2*x)+c.c
  elseif ((0.5<=x) && (x<=1.0))
    sqrt(2*x-1)+c.c
  else
    error("Input out of domain of inverse cdf of PBCS.")
  end
end

function generate!(c::PBCSCarrier, x::Vector{Float64})
  local a::Int
  dunif = Uniform(0, 1)
  dbern = Bernoulli()

  for i = 1:2:c.len
    x[i] = rand(dunif)
  end

  for i = 2:2:c.len
    a = rand(dbern) == 0 ? -1 : 1
    x[i] = c.c+a*sqrt(1-(x[i-1]-c.c)^2)
  end

  return x
end

generate(c::PBCSCarrier) = generate!(c, Array(Float64, c.len))

mean(c::PBCSCarrier) = c.c

var(c::PBCSCarrier) = 0.5


text:
#Example: Using MPI Manager to do parallel modeling with PIPS-NLP
using MPIClusterManagers
using Distributed

manager = MPIManager(np = 3)
addprocs(manager)

@mpi_do manager begin
    #Run mpi on workers
    using Pkg
    Pkg.activate("../")
    using Plasmo
    using MPI
    using PipsNLP
    comm=MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)

    println("I am $(rank) of $(MPI.Comm_size(comm))")
    MPI.Barrier(comm)
end

@mpi_do manager begin
    #All ranks have optigraph and first stage optinode
    graph = OptiGraph()

    @optinode(graph,n1)
    @variable(n1,x0 >= 0)

    subgraph = OptiGraph()
    add_subgraph!(graph,subgraph)
    @optinode(subgraph,n2)
    @optinode(subgraph,n3)
    @optinode(subgraph,n4)

    #Linked variables using linkconstraints
    @variable(n2,x >= 0)
    @variable(n3,x >= 0)
    @variable(n4,x >= 0)

    #Local models
    if rank == 0
        @variable(n2,x2 >= 2)
        @variable(n2,y2 >= 0)
        @constraint(n2,x2 + y2 == 3)
        @NLconstraint(n2,x2^3 <= 10)
        @constraint(n2,n2[:x] + n2[:x2] <= 10)
        @linkconstraint(graph,n1[:x0] == y2)
        @objective(n2,Min,x2)
    elseif rank == 1
        @variable(n3,x3 >= 2)
        @variable(n3,y3 >= 0)
        @constraint(n3,x3 + y3 == 4)
        @constraint(n3,n3[:x] + n3[:x3] <= 10)
        @linkconstraint(graph,n1[:x0] == y3)
        @objective(n3,Min,x3)
    elseif rank == 2
        @variable(n4,x4 >= 2)
        @variable(n4,y4 >= 0)
        @constraint(n4,x4 + y4 == 5)
        @constraint(n4,n4[:x] + n4[:x4] <= 10)
        @linkconstraint(graph,n1[:x0] == y4)
        @objective(n4,Min,x4)
    end

    #Link constraints between subnodes
    @linkconstraint(graph,n2[:x] == n3[:x])
    @linkconstraint(graph,n3[:x] == n4[:x])
end

#Solve with PIPS-NLP
@mpi_do manager begin
    pipsnlp_solve(graph)
end

#Use Julia to fetch solution values from each rank
solution_values = []
for worker in workers()
    solution_worker = fetch(@spawnat(worker,value.(all_variables(graph))))
    push!(solution_values,solution_worker)
end

#Use MPI Manager to print solution values on each rank
@mpi_do manager begin
    if rank == 0
        println("first stage: ")
        @show (value(n1[:x0]))

        println("second stage: n2")
        @show (value(n2[:y2]))
        @show (value(n2[:x]))
    end
end

@mpi_do manager begin
    if rank == 1
        println("second stage: n3")
        @show (value(n3[:y3]))
        @show (value(n3[:x]))
    end
end

@mpi_do manager begin
    if rank == 2
        println("second stage: n4")
        @show (value(n4[:y4]))
        @show (value(n4[:x]))
    end
end


text:
if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

using IJulia
using Interact

profile = normpath(dirname(@__FILE__), "profile.json")
IJulia.init([profile])
redirect_stdout(IJulia.orig_STDOUT)
redirect_stderr(IJulia.orig_STDERR)

include(Interact.ijulia_setup_path)

sliderWidget = slider(1:5)

@test """Interact.Slider{Int64}(Signal{Int64}(3, nactions=0),"",3,1:5,"d",true)""" == stringmime("text/plain", sliderWidget)
@test "3" == stringmime("text/plain", signal(sliderWidget))

@test "" == stringmime("text/html", sliderWidget)

close(IJulia.ctx)


text:
@testset "Permute parameters" begin
    begin
        # all combinations
        universe = Dict{Any, Vector{Any}}("wnd" => [1,2,3], :mode => ["A", "B"], "coef" => [0.1, 0.5, 1.0])
        combinations = permute_params(universe)

        # check correct number of combinations
        @test length(combinations) == prod(map(length, values(universe)))

        # check all combinations have exactly three parameters
        @test all(map(x -> length(x) == 3, combinations))
    end

    begin
        # filtered combinations for which hold: wnd > 1 if mode == 'A'
        universe = Dict{Any, Vector{Any}}(:wnd => [1,2,3], :mode => ["A", "B"], :coef => [0.1, 0.5, 1.0])
        f = x -> x[:mode] != "A" || x[:wnd] > 1
        combinations = permute_params(universe; filter=f)

        # check correct number of combinations
        @test length(combinations) == 15

        # check all combinations have exactly three parameters
        @test all(map(x -> length(x) == 3, combinations))
    end

    begin
        # filter all (empty result)
        universe = Dict{Any, Vector{Any}}(:wnd => [1,2,3], :mode => ["A", "B"], :coef => [0.1, 0.5, 1.0])
        combinations = permute_params(universe; filter=x -> false)

        # check correct number of combinations
        @test length(combinations) == 0
    end

    begin
        # shuffle
        universe = Dict{Any, Vector{Any}}(:a => [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], :b => ["a", "b", "c"])
        combinations1 = permute_params(universe; shuffle=true)
        combinations2 = permute_params(universe; shuffle=false)

        @test !all(combinations1 .== combinations2)
    end
end


text:
@testset "fusion" begin
    θ = fill(1.0, 5)
    setup = MCMCSetup(
        BSI.Imputation(BSI.NoBlocking(), 0.5, BSI.Vern7()),
        BSI.Imputation(BSI.NoBlocking(), 0.5, BSI.Vern7()),
        BSI.ParamUpdate(BSI.MetropolisHastingsUpdt(), [1], θ,
                        BSI.UniformRandomWalk(0.5, true),
                        BSI.ImproperPosPrior(),
                        BSI.UpdtAuxiliary(BSI.Vern7(), true)),
        BSI.ParamUpdate(BSI.MetropolisHastingsUpdt(), [2], θ,
                        BSI.UniformRandomWalk(0.5, true),
                        BSI.ImproperPosPrior(),
                        BSI.UpdtAuxiliary(BSI.Vern7(), false)),
        BSI.ParamUpdate(BSI.MetropolisHastingsUpdt(), [4], θ,
                        BSI.UniformRandomWalk(0.5, true),
                        BSI.ImproperPosPrior(),
                        BSI.UpdtAuxiliary(Vern7(), false)))
    updt_idx_af = [[1,6],[2,6]]
    schedule = BSI.MCMCSchedule(19, [[1,3],[2,4,5]],
                                (save=5, verbose=4, warm_up=3,
                                readjust=(x->false), fuse=(x->x%10==0)))
    ws = BSI.MCMCWorkspace(setup, schedule, θ)
    cov_mat = rand(5,5)
    ws.cov .= copy(cov_mat)
    for step in schedule
        step.fuse && BSI.fuse!(ws, schedule)
        if step.iter >= 10
            @test step.idx == updt_idx_af[mod1(step.iter, 2)]
        end
    end
    @test length(ws.updates) == 6
    foo = ws.updates[6]
    @test foo.updt_coord == Val{(true, true, false, true, false)}()
    @test foo.t_kernel.gsn_A.Σ == 2.38^2/3.0*Matrix(view(cov_mat, [1,2,4], [1,2,4]))
    @test foo.t_kernel.gsn_B.Σ == 0.1^2/3.0*Matrix{Float64}(I, 3, 3)
    @test all(foo.t_kernel.gsn_A.pos)
    @test all(foo.t_kernel.gsn_B.pos)
    @test foo.priors.priors == (BSI.ImproperPosPrior(), BSI.ImproperPosPrior(),
                               BSI.ImproperPosPrior())
    @test foo.priors.coord_idx == (Val{(true, false, false, false, false)}(),
                                  Val{(false, true, false, false, false)}(),
                                  Val{(false, false, false, true, false)}())
end


text:
module Suri
using Hecke
import Hecke: valuation, divexact, parent_type, elem_type, mul!, addeq!, parent
import Base: +, -, *, ^

#= follows Sebastian Posur's idea
  Start with a pseudo_matrix with ideals A_i and rows alpha_i
  plus an additional pseudo-element      A            alpha
  assume, at least for now, the alpha_i are independent.
  assume sum gamma_i alpha_i + gamma alpha = 0
  (gamma are the dependency)
  Now P :=
  (inv(A_i) gamma_i)
  ( inv(A)  gamma  )
  is a (n+1) x 1 pseudo matrix. Using the (pseudo) HNF we find a matrix
  T in Gl(n+1, K) and ideals C_i s.th.
  T*P = H is a pseudo-HNF. By construction, H has 1. row 1, all other rows 0
  and the coeff. ideals C_i
  Then T is an isomorphism between
  sum inv(A_i) + inv(A) and sum C_i
  as zk-modules
  The inverse-transpose of T is an isomorphism between
  sum inv(C_i) and sum A_i + A
  the dual modules to the ones above

  The idea:
  T is an iso between all modules given by the ideals
  The construction forces inv(T) to have the relation as 1. row/col
  The dual of an abstract module is giving by the inverse ideals
  The corresponding map is transpose(inv(T)) = inv(T)'
  Thus it works...
=#

function mod_lll(a::NfAbsOrdElem, I::NfAbsOrdIdl)
  l = lll(I)[2]
  S = l*basis_matrix(I)
  Si = pseudo_inv(S)
  c = matrix(FlintZZ, 1, nrows(l), coordinates(a)) * Si[1]
  d = matrix(FlintZZ, 1, nrows(l), [round(fmpz, x, Si[2]) for x = c])
  return a - Hecke.parent(a)(collect(d*S))
end

function mod_lll(a::nf_elem, I::Hecke.NfAbsOrdFracIdl)
  O = order(I)
  d = lcm(denominator(a, O), denominator(I))
  return divexact(Hecke.parent(a)(mod_lll(O(d*a), simplify(I*d).num)), d)
end

function _reduce(a::Hecke.PMat, T)
  A = a.coeffs
  M = a.matrix
  for i=2:nrows(M)
    for j=1:i-1
      if iszero(M[j, i])
        continue
      end
      I = A[i]*M[i,i]*inv(A[j])
      c = mod_lll(M[j,i], I)
      @assert (c - M[j,i]) in I
      d = divexact(M[j, i] -c, M[i,i])
      M[j, :] = M[j, :] - d*M[i, :]
      T[j, :] = T[j, :] - d*T[i, :]
      @assert M[j, i] == c
    end
  end
end

function extend(M::Hecke.PMat, b::Generic.MatSpaceElem{nf_elem}, gamma::Generic.MatSpaceElem{nf_elem})

  @assert iszero(hcat(M.matrix', b)*gamma)
  zk = base_ring(M)
  nc = ncols(gamma)
  n = nrows(gamma) - 1
  @assert nc == ncols(b)
  p = pseudo_matrix(hcat(gamma, vcat(identity_matrix(zk, n), zero_matrix(zk, 1, n))), vcat(map(inv, coefficient_ideals(M)), [1*zk for i=1:nc]))
#  @show p

  h, T = Hecke.pseudo_hnf_with_transform(p)
#  @show T*p.matrix == h.matrix
  _reduce(h, T)
#  @show T*p.matrix == h.matrix
  for i=1:nrows(h)
    j, al = Hecke.reduce_ideal(h.coeffs[i])
    T[i, :] *= inv(al)
    h.coeffs[i] = j
    h.matrix[i, :] *= inv(al)
  end
#  @show pseudo_hnf(h) == pseudo_hnf(p)
#  @show T*p.matrix == h.matrix
  #for n x 1 matrices, the transform, especially the inverse can more easily be computed
#  @assert prod(coefficient_ideals(p)) ==  (det(T))*(prod(coefficient_ideals(h)))
#  @assert all(T[j,i] in (coefficient_ideals(p)[i])*inv(coefficient_ideals(h)[j]) for i=1:nrows(M) for j=1:nrows(M))
  Ti = inv(T)
#  @assert all(Ti[i,j] in inv(coefficient_ideals(p)[i])*(coefficient_ideals(h)[j]) for i=1:nrows(M) for j=1:nrows(M))
#there is a transpose problem somewhere...
  e = pseudo_matrix((hcat(M.matrix', b)*Ti)[:, nc+1:nrows(M)+nc]', map(inv, coefficient_ideals(h)[nc+1:nrows(M)+nc]))
  return e
end

function Hecke.denominator(P::Hecke.PMat, M::NfOrd)
  l = fmpz(1)
  p = matrix(P)
  for i=1:nrows(P)
    I = coefficient_ideals(P)[i]
    if isa(I, Hecke.NfAbsOrdFracIdl)
      Hecke.assure_2_normal(I.num)
    else
      Hecke.assure_2_normal(I)
    end
    for j=1:ncols(P)
      l = lcm(l, denominator(simplify(p[i,j]*I)))
    end
  end
  return l
end

function Hecke.simplify(P::Hecke.PMat)
  c = copy(coefficient_ideals(P))
  m = deepcopy(matrix(P))
  for i=1:length(c)
    a, b = Hecke.reduce_ideal(c[i])
    m[i, :] *= inv(b)
    c[i] = fractional_ideal(order(a), a)
  end
  return pseudo_matrix(m, c)
end

function Hecke.dual(P::Hecke.PMat)
  return pseudo_matrix(inv(P.matrix)', map(inv, coefficient_ideals(P)))
end

function Hecke.invmod(A::Generic.MatSpaceElem{nf_elem}, X::fmpz)
  k = base_ring(A)
  zk = maximal_order(k)
  q, mq = quo(zk, X*zk)
  if true
    iA = inv(A)
    id = lcm([denominator(x, zk) for x = iA])
    BX = map_entries(x->mq(zk(id*x)), iA)*invmod(id, X)
  else
    d = lcm([denominator(x, zk) for x= A])
    AX = map_entries(x->mq(zk(d*x)), A)
    BX = inv(AX)*invmod(d, X)
  end
  B = map_entries(x->k(preimage(mq, x)), BX)
  my_mod_sym!(B, X)
  return B
end

function Hecke.invmod(A::fmpz_mat, X::fmpz)
  B0 = map_entries(lift, inv(map_entries(quo(ZZ, X)[1], A)))
  mod_sym!(B0, X)
  return B0
end

function my_mod_sym!(A::fmpz_mat, X::fmpz, ::Any)
  mod_sym!(A, X)
end

function valuation(a::NfAbsOrdElem{AnticNumberField,nf_elem}, X::fmpz)
  v = 0
  first = true
  for x = coordinates(a)
    iszero(x) && continue
    if first
      v = valuation(x, X)
      first = false
    else
      v = min(v, valuation(x, X))
    end
    iszero(v) && return v
  end
  return v
end

function mod_sym(A::NfAbsOrdElem{AnticNumberField,nf_elem}, X::fmpz)
  c = coordinates(A)
  d = map(x->Hecke.mod_sym(x, X), c)
  return parent(A)(d)
end

function my_mod_sym!(A::Generic.MatSpaceElem{nf_elem}, X::fmpz)
  k = base_ring(A)
  zk = maximal_order(k)
  for i=1:nrows(A)
    for j=1:ncols(A)
      A[i,j] = k(mod_sym(zk(A[i,j]), X))
    end
  end
end

#modelled after
#  http://www.ti3.tu-harburg.de/paper/rump/NiRuOi11.pdf
#  On the generation of very ill-conditioned integer matrices
#  Tetsuo Nishi1a, Siegfried M. Rump, and Shin’ichi Oishi
#  Nonlinear Theory and Its Applications, IEICE, vol. 2, no. 2, pp. 226N245
#should generate a matrix where the entries are roughly U^2 (U a range)
#determinant 1, and the inverse should have much larger coefficients.
function bad_mat(R::Hecke.Ring, n::Int, U)
  z = zero_matrix(R, n, n)
  for i=1:n-1
    z[i+1,i] = one(R)
    z[i+1, i+1] = rand(R, U)
  end
  d = one(R)
  for i=n:-1:1
    k = rand(R, U)
    z[1, i] = d+k*z[i,i]
    d = k
  end
  return z
end

function bad_mat_suri(R::Hecke.Ring, n::Int, U)
  m = identity_matrix(R, n)
  for i=1:n
    for j=i+1:n
      m[i,j] = rand(R, U)
    end
  end
  return m
end

mutable struct RRS <: Hecke.Ring
  p::Vector{fmpz}
  P::Vector{fmpz}
  Pi::Vector{fmpz}
  r::fmpz
  N::fmpz
  ce

  function RRS(p::Vector{fmpz})
    s = new()
    s.p = p
    P = prod(p)
    s.P = [divexact(P, x) for x = p]
    s.Pi = [invmod(s.P[i], s.p[i]) for i=1:length(p)]
    s.r = next_prime(2^50)
    s.N = P
    s.ce = Hecke.crt_env(p)
    return s
  end

  function RRS(p::Vector{<:Integer})
    return RRS(fmpz[x for x = p])
  end
end
function Base.show(io::IO, R::RRS)
  print(io::IO, "crt ring with moduli ", R.p)
end

mutable struct RRSelem <: Hecke.RingElem
  x::Vector{fmpz}
  r::fmpz
  R::RRS
  function RRSelem(X::RRS, a::fmpz)
    s = new()
    s.x = [mod(a, p) for p = X.p]
    s.r = mod(a, X.r)
    s.R = X
    return s
  end
  function RRSelem(X::RRS, a::Integer)
    return RRSelem(X, fmpz(a))
  end
  function RRSelem(X::RRS, a::Vector{fmpz}, k::fmpz)
    r = new()
    r.R = X
    r.x = a
    r.r = k
    return r
  end
end

function Base.show(io::IO, a::RRSelem)
  print(io, "crt: ", a.x)
end

elem_type(::RRS) = RRSelem
parent_type(::RRSelem) = RRS
parent_type(::Type{RRSelem}) = RRS

parent(a::RRSelem) = a.R

-(a::RRSelem, b::RRSelem) = RRSelem(a.R, [mod(a.x[i]-b.x[i], a.R.p[i]) for i=1:length(a.x)], mod(a.r-b.r, a.R.r))
+(a::RRSelem, b::RRSelem) = RRSelem(a.R, [mod(a.x[i]+b.x[i], a.R.p[i]) for i=1:length(a.x)], mod(a.r+b.r, a.R.r))
*(a::RRSelem, b::RRSelem) = RRSelem(a.R, [mod(a.x[i]*b.x[i], a.R.p[i]) for i=1:length(a.x)], mod(a.r*b.r, a.R.r))
*(a::Integer, b::RRSelem) = RRSelem(b.R, [mod(a*b.x[i], b.R.p[i]) for i=1:length(b.x)], mod(a*b.r, b.R.r))
divexact(a::RRSelem, b::RRSelem) = RRSelem(a.R, [mod(a.x[i]*invmod(b.x[i], a.R.p[i]), a.R.p[i]) for i=1:length(a.x)], mod(a.r*invmod(b.r, a.R.r), a.R.r))
-(a::RRSelem) = RRSelem(a.R, [mod(-a.x[i], a.R.p[i]) for i=1:length(a.x)], -a.r)
^(a::RRSelem, e::Integer) = RRSelem(a.R, [powermod(a.x[i], e, a.R.p[i]) for i=1:length(a.x)], powermod(a.r, e, a.R.r))
(R::RRS)() = RRSelem(R, fmpz[0 for i=1:length(R.p)], fmpz(0))
(R::RRS)(a::Integer) = RRSelem(R, a)
(R::RRS)(a::RRSelem) = a

function addeq!(a::RRSelem, b::RRSelem)
  for i=1:length(a.x)
    a.x[i] = mod(a.x[i] + b.x[i], a.R.p[i])
    a.r    = mod(a.r    + b.r   , a.R.r)
  end
  return a
end

function mul!(a::RRSelem, b::RRSelem, c::RRSelem)
  for i=1:length(a.x)
    a.x[i] = mod(b.x[i] * c.x[i], a.R.p[i])
    a.r    = mod(b.r    * c.r   , a.R.r)
  end
  return a
end

function check(a::RRSelem)
  z = fmpz(a)
  @assert mod(z, a.R.r) == a.r
end

#given x mod p_i and p_r, find x mod p
function extend(R::RRS, a::RRSelem, p::fmpz)
  k = sum(((a.x[i]*R.Pi[i]) % R.p[i]) * (R.P[i] % R.r) for i=1:length(R.p)) - a.r
  k = (k*invmod(R.N, R.r)) % R.r
  @assert k <= length(R.p)
  return (sum(((a.x[i]*R.Pi[i]) % R.p[i]) * (R.P[i] % p) for i=1:length(R.p)) - k*(R.N % p)%p)%p
end

function mixed_radix(R::RRS, a::RRSelem, li::Int = typemax(Int))
  A = fmpz[]
  for i=1:min(length(a.x), li)
    y = a.x[i]
    for j=1:i-1
      y = mod(((y-A[j])*invmod(R.p[j], R.p[i])), R.p[i])
    end
    push!(A, y)
  end
  return A
  #so a = A[1] + A[2]*p[1] + A[3]*p[1]*p[2] ...s
end

function rss_elem_from_radix(R::RRS, a::Vector{fmpz})
  x = fmpz[]
  for i=1:length(R.p)
    z = a[1]
  end
end

function gen(R::RRS, i::Int)
  p = fmpz[0 for i=1:length(R.p)]
  p[i] = fmpz(1)
  r = mod(R.P[i]*R.Pi[i], R.r)
  return RRSelem(R, p, r)
end

Hecke.fmpz(a::RRSelem) = Hecke.crt(a.x, a.R.ce)

# a random invertable matrix with coeffs in R
#start with identity, do i rounds of random elementary row and column
#operations
function rand_gl(R::Hecke.Ring, n::Int, u, i::Int)
  A = identity_matrix(R, n)
  for j=1:i
    a = rand(R, u)
    mu = rand(1:n)
    nu = rand(1:n)
    if mu == nu
      continue
    end
    A[mu, :] += a*A[nu, :]
    a = rand(R, u)
    mu = rand(1:n)
    nu = rand(1:n)
    if mu == nu
      continue
    end
    A[:, mu] += a*A[:, nu]
  end
  return A
end

function DoublePlusOne(A, X::fmpz, k::Int)
  R = base_ring(A)

  B0 = invmod(A, X)
  R = [divexact(identity_matrix(R, nrows(A))-A*B0, X)]
  M = []

  for i=0:k
    r = R[end]^2
    m = B0*r
    my_mod_sym!(m, X)
    push!(M, m)
    push!(R, divexact(r-A*M[end], X))
    if iszero(R[end])
      break
    end
  end
  return B0, M, R
end

function getB(B0, M, R, X)
  B = [B0]
  I = identity_matrix(base_ring(B0), nrows(B0))
  XX = X
  for i=1:length(M)
    push!(B, B[end]*(I + R[i]*XX) + M[i]*XX^2)
    XX = XX^2*X
    my_mod_sym!(B[end], XX)
  end
  return B
end

end

#= example

Hecke.example("Suri.jl")
Hecke.revise("Suri.jl")

n = 6

k, a = wildanger_field(3, 13)
k, a = quadratic_field(-11)
m = rand(MatrixSpace(k, n, n), 1:10);
m = cat(m,m, dims=(1,2));
b = rand(MatrixSpace(k, 2*n, 1), 1:10);
S = kernel(hcat(m, b));

m1 = Suri.extend(pseudo_matrix(m'), b, S[2]);

norm(det(m))
norm(det(m1))

b = rand(MatrixSpace(k, 2*n, 1), 1:10);
S = kernel(hcat(m1.matrix', b));
m2 = Suri.extend(m1, b, S[2]);

norm(det(m2))

=#


text:
# Setting up initial conditions for a normally incident plane P wave
# Analytical plane wave solutions are based on Peng (1994), PhD Thesis, MIT.
# dependent on: rfft/irfft (FFTW)
#             : Spline1D (Dierckx)
#             : besselj, hankelh1 (SpecialFunctions)

#-------------------------------
# tapering initial field
#-------------------------------
function taper_initial_field!(field,nz,nr,dz,src_iz,wavelength,ntaper)
    #field[nz,nr]
    #flat weight 1 between src_iz-W:src_iz+W
    #ntaper after then

    #final vertical tapering
    #ntaper=10 :given
    tmp_angle=collect(range(0,pi/2,length=ntaper))
    #wavelength_half=Vp_solid/f0/2
    wavelength_half=wavelength/2
    wl_half_iz=Int(round(wavelength_half/dz)+1)+10
    if (src_iz-wl_half_iz-ntaper<=0)
        error("srcdepth is too shallow or src wavelength is too large!")
    elseif (src_iz+wl_half_iz+ntaper>nz)
        error("srcdepth is too deep or src wavelength is too large!")
    end

    weightvec=zeros(1,nz)
    cnt=0
    for iz=src_iz-(ntaper+wl_half_iz)+1:src_iz-wl_half_iz
        cnt=cnt+1
        weightvec[iz]=(sin(tmp_angle[cnt]))^2
    end
    for iz=src_iz-wl_half_iz+1:src_iz+wl_half_iz
        weightvec[iz]=1.0
    end
    cnt=0
    for iz=src_iz+wl_half_iz+1:src_iz+wl_half_iz+ntaper
        cnt=cnt+1
        weightvec[iz]=(sin(tmp_angle[ntaper-cnt+1]))^2
    end

    for iz=1:nz
      tmp_amp=weightvec[iz]
      field[iz,:]=tmp_amp*field[iz,:];
    end

    println("Initial field tapered until iz=",src_iz+wl_half_iz+ntaper)
end #function

#--------------------------------------------------------------
# Peng's analytical solution of borehole pressure and
# elastic field in surrounding media
# (open hole only)
#--------------------------------------------------------------

function Peng_solution(nr,nz,dr,dz,ir_wall,tvec,src_func,srcdepth,f0,ntaper,noffset,
    Vp,Vs,Rho)
   #important: Schoenberg/Peng -> FT Aki-Rechards
   #         : rfft <-> irfft -> FT MATLAB

   dt=tvec[2]-tvec[1]
   T=tvec[end]
   ns=length(tvec)
   fvec=[0:1/T:1/T*(ns-1)]
   fvec=collect(fvec[1])

   fsrc=rfft(src_func[:])
   test_src=irfft(fsrc,ns)
   ns_half=length(fsrc)
   fvec=fvec[1:ns_half]
#   plot(fvec,map(abs,fsrc),xlims=(0,1000))

   #calculating Schoenberg's exact solution for vz, and vr, at fluid and solid phase
   #normal incident P wave only (positive z propagation)
   src_iz=Int(round(srcdepth-dz/2)/dz+1) #src defined at vz cell
   rb=(ir_wall-1)*dr+dr/2.0 # make sure Vp[1:ir_wall]=1500

   max_freq_src=3*f0
   if_max=minimum(findall(x -> x > max_freq_src,fvec)) #to be used to filter high-freq noise

   Vp_fluid=Vp[src_iz,1]
   Rho_fluid=Rho[src_iz,1]
   Vp_solid=Vp[src_iz,ir_wall+1] # make sure Vp[1:ir_wall]=1500
   Vs_solid=Vs[src_iz,ir_wall+1] # make sure Vp[1:ir_wall]=1500
   Rho_solid=Rho[src_iz,ir_wall+1] # make sure Vp[1:ir_wall]=1500

   angle_delta=1E-4 # =0 normal incidence, but give a small number as it diverges due to Hankel function.

   #checking ir_wall value
   flag_error_irwall=0
   if (Vs_solid < 1.0)
      flag_error_irwall=1
   end
   if (Vp_solid == Vp_fluid)
      flag_error_irwall=1
   end

   if (flag_error_irwall==1)
      error("Check ir_wall in Peng_solution! Is ir_wall at the boundary??")
   end

   #output waveforms
   #outputs are vr[ns,nr] and vz[ns,nr]
   vr_Peng=zeros(ns,nr)
   vz_Peng=zeros(ns,nr)

   vphi_Peng=zeros(ns,nr) #zero (axisymmetric)

   trr_Peng=zeros(ns,nr)
   tpp_Peng=zeros(ns,nr)
   tzz_Peng=zeros(ns,nr)
   trp_Peng=zeros(ns,nr) #zero (axisymmetric)
   trz_Peng=zeros(ns,nr)
   tpz_Peng=zeros(ns,nr) #zero (axisymmetric)

   #setting up
   #solid: vz
   fvz_solid=fsrc #independent of r
   maxamp_fvz=maximum(map(abs,fvz_solid))
   max_freq_src=3*f0
   if_max=minimum(findall(x -> x > max_freq_src,fvec)) #to be used to filter high-freq noise

   # Creating -dz/2/Vp_solid time delayed src function for later use.
   # Time rewinds because vr and tii are dz/2 upwards vz, and plane wave propagates downwards.
   # In addition, we consider stress as dt/2 time delayed in the main loop
   # Therefore, this function is used only to create stress BC
   #Option:
   #Default stress is offset by dt/2 from velocity.
   #When enabling pre-update before the main loop, stress is offset by -dt/2 from velocity

   println("===============================")
   println("Output initial stress is +dt/2 from Velocity (default).")
   src_func_forStress=(src_func[1:end]+[src_func[2:end]; 0])/2.0
 #  println("Output initial stress is -dt/2 from Velocity (make sure to enable pre-update!).")
 #  src_func_forStress=(src_func[1:end]+[0; src_func[1:end-1]])/2.0
   println("===============================")

   fsrc_forStress=rfft(src_func_forStress[:])
   tmp_phase_delay=im*2pi*fvec*(dz/2)/Vp_solid #negative phase -> time delays -> MATLAB FT
   fsrc_forStress_delayed=fsrc_forStress.*map(exp,tmp_phase_delay) #It says "delay" but it actually rewinds as exlained above
   fsrc_forStress_delayed[if_max+1:end]=zeros(ns_half-if_max,1)
   src_forStress_delayed=irfft(fsrc_forStress_delayed,ns)

   #repeat every r_now
   for ir=1:nr-noffset

      #Start frequency response
      fvz=complex(zeros(1,ns_half))
      fvr=complex(zeros(1,ns_half))

      ftrr=complex(zeros(1,ns_half))
      ftpp=complex(zeros(1,ns_half))
      ftzz=complex(zeros(1,ns_half))

      ftrz=complex(zeros(1,ns_half)) #other shear component are zero (axisymmetric)

      for iw=2:if_max
         freq=fvec[iw]
         w=2pi*freq
         #Solving B.C. to get Coefficients of potential functions
         A0,B0,C0=SolvePengBC(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)

        if (ir<=ir_wall)
         #fluid:vz
         r_now=(ir-1)*dr
         fvz[iw]=(-im*w)*uz_f(r_now,A0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
         fvz[iw]=conj(fvz[iw]) #FT Aki-Rechards -> FT MATLAB
         fvz[iw]=fvz[iw]*fsrc[iw]/w^2 #/w^2 scaling is due to V(w)=w^2 in Peng (1994).

         #fluid:vr (with delay=-(dz/2)/Vp_solid), time reversed because vr is dz/2 upwards vz, and plane wave downwards
         r_now=(ir-1)*dr+dr/2.0
         fvr[iw]=(-im*w)*ur_f(r_now,A0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
         fvr[iw]=conj(fvr[iw]) #FT Aki-Rechards -> FT MATLAB
         fvr[iw]=fvr[iw]*fsrc[iw]/w^2 #/w^2 scaling is due to V(w)=w^2 in Peng (1994).
            #additinoal delay (positive z propagation)
            tmp_phase_delay=im*2pi*fvec[iw]*(dz/2)/Vp_solid #negative phase -> time advance -> MATLAB FT
         fvr[iw]=fvr[iw]*exp(tmp_phase_delay)

         #fluid: pressure (with delay=-(dz/2)/Vp_solid)
         r_now=(ir-1)*dr
         ftrr[iw]=trr_f(r_now,A0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
         ftrr[iw]=conj(ftrr[iw]) #FT Aki-Rechards -> FT MATLAB
         ftrr[iw]=ftrr[iw]*fsrc_forStress_delayed[iw]/w^2 #/w^2 scaling is due to V(w)=w^2 in Peng (1994).
            #additinoal delay (positive z propagation)
            tmp_phase_delay=im*2pi*fvec[iw]*(dz/2)/Vp_solid #negative phase -> time advance -> MATLAB FT
         ftrr[iw]=ftrr[iw]*exp(tmp_phase_delay)

         ftpp[iw]=ftrr[iw]
         ftzz[iw]=ftrr[iw]
         ftrz[iw]=0.0
      else #(solid)
         #solid:vz
         r_now=(ir-1)*dr
         fvz[iw]=(-im*w)*(
             uz_phi(r_now,B0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)+
             uz_xi(r_now,C0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)+
             uz_P(r_now,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
             )
         fvz[iw]=conj(fvz[iw]) #FT Aki-Rechards -> FT MATLAB
         fvz[iw]=fvz[iw]*fsrc[iw]/w^2 #/w^2 scaling is due to V(w)=w^2 in Peng (1994).

         #solid:vr (with delay=-(dz/2)/Vp_solid)
         r_now=(ir-1)*dr+dr/2.0
         fvr[iw]=(-im*w)*(
             ur_phi(r_now,B0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)+
             ur_xi(r_now,C0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)+
             ur_P(r_now,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
             )
         fvr[iw]=conj(fvr[iw]) #FT Aki-Rechards -> FT MATLAB
         fvr[iw]=fvr[iw]*fsrc[iw]/w^2 #/w^2 scaling is due to V(w)=w^2 in Peng (1994).
            #additinoal delay (positive z propagation)
            tmp_phase_delay=im*2pi*fvec[iw]*(dz/2)/Vp_solid #negative phase -> time advance -> MATLAB FT
         fvr[iw]=fvr[iw]*exp(tmp_phase_delay)

         #solid: trr,tpp,tzz (with delay=-(dz/2)/Vp_solid)
         r_now=(ir-1)*dr
         ftrr[iw]=(
             trr_phi(r_now,B0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)+
             trr_xi(r_now,C0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)+
             trr_P(r_now,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
             )
         ftpp[iw]=(
             tpp_phi(r_now,B0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)+
             tpp_xi(r_now,C0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)+
             tpp_P(r_now,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
             )
         ftzz[iw]=(
             tzz_phi(r_now,B0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)+
             tzz_xi(r_now,C0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)+
             tzz_P(r_now,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
             )
          ftrr[iw]=conj(ftrr[iw]) #FT Aki-Rechards -> FT MATLAB
          ftpp[iw]=conj(ftpp[iw]) #FT Aki-Rechards -> FT MATLAB
          ftzz[iw]=conj(ftzz[iw]) #FT Aki-Rechards -> FT MATLAB
          ftrr[iw]=ftrr[iw]*fsrc_forStress_delayed[iw]/w^2 #/w^2 scaling is due to V(w)=w^2 in Peng (1994).
          ftpp[iw]=ftpp[iw]*fsrc_forStress_delayed[iw]/w^2 #/w^2 scaling is due to V(w)=w^2 in Peng (1994).
          ftzz[iw]=ftzz[iw]*fsrc_forStress_delayed[iw]/w^2 #/w^2 scaling is due to V(w)=w^2 in Peng (1994).
            #additinoal delay (positive z propagation)
            tmp_phase_delay=im*2pi*fvec[iw]*(dz/2)/Vp_solid #negative phase -> time advance -> MATLAB FT
          ftrr[iw]=ftrr[iw]*exp(tmp_phase_delay)
          ftpp[iw]=ftpp[iw]*exp(tmp_phase_delay)
          ftzz[iw]=ftzz[iw]*exp(tmp_phase_delay)

          #solid: trz
          r_now=(ir-1)*dr+dr/2.0
          ftrz[iw]=(
             trz_phi(r_now,B0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)+
             trz_xi(r_now,C0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)+
             trz_P(r_now,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
             )
          ftrz[iw]=conj(ftrz[iw]) #FT Aki-Rechards -> FT MATLAB
          ftrz[iw]=ftrz[iw]*fsrc_forStress_delayed[iw]/w^2 #/w^2 scaling is due to V(w)=w^2 in Peng (1994).
       end #if ir
    end #iw
    vz_tmp=irfft(fvz[:],ns)
    vr_tmp=irfft(fvr[:],ns)
    trr_tmp=irfft(ftrr[:],ns)
    tpp_tmp=irfft(ftpp[:],ns)
    tzz_tmp=irfft(ftzz[:],ns)
    trz_tmp=irfft(ftrz[:],ns)

    vz_Peng[:,ir]=vz_tmp
    vr_Peng[:,ir]=vr_tmp
    trr_Peng[:,ir]=trr_tmp
    tpp_Peng[:,ir]=tpp_tmp
    tzz_Peng[:,ir]=tzz_tmp
    trz_Peng[:,ir]=trz_tmp
  end #ir

    #final tapering
    if (ntaper>1)
    tmp_angle=collect(range(0,pi/2,length=ntaper))
    for ir=1:ntaper
      tmp_amp=(sin(tmp_angle[ir]))^2
      vz_Peng[:,nr-ir+1-noffset]=tmp_amp*vz_Peng[:,nr-ir+1-noffset]
      vr_Peng[:,nr-ir+1-noffset]=tmp_amp*vr_Peng[:,nr-ir+1-noffset]

      trr_Peng[:,nr-ir+1-noffset]=tmp_amp*trr_Peng[:,nr-ir+1-noffset]
      tpp_Peng[:,nr-ir+1-noffset]=tmp_amp*tpp_Peng[:,nr-ir+1-noffset]
      tzz_Peng[:,nr-ir+1-noffset]=tmp_amp*tzz_Peng[:,nr-ir+1-noffset]
      trz_Peng[:,nr-ir+1-noffset]=tmp_amp*trz_Peng[:,nr-ir+1-noffset]

    end
    end

 return vz_Peng,vr_Peng,vphi_Peng,
        trr_Peng,tpp_Peng,tzz_Peng,trp_Peng,trz_Peng,tpz_Peng,
        src_iz
end

#--------------------------------------------------------------
# Converting time-domain solution to
# space domain to construct initial conditions
#--------------------------------------------------------------
#Schoenberg's (or Peng's) solution as an initial condition
#Assuming open hole
function Schoenberg_solution_initialBC(nr,nz,dr,dz,ir_wall,Vp,
         tvec_vel,tvec_stress,
         srcdepth,f0,delay,
         vz_Schoenberg,vr_Schoenberg,vphi_Schoenberg,
         trr_Schoenberg,tpp_Schoenberg,tzz_Schoenberg,trp_Schoenberg,trz_Schoenberg,tpz_Schoenberg,
         ntaper,noffset)

   dt=tvec_vel[2]-tvec_vel[1]
   T=tvec_vel[end]
   ns=length(tvec_vel)

   #calculating Schoenberg's approximate solution for vz, and vr, at fluid and solid phase
   #normal incidence Pwave only (positive z propagation)
   src_iz=Int(round(srcdepth-dz/2)/dz+1) #src defined at vz cell
   rb=(ir_wall-1)*dr+dr/2.0 # make sure Vp[1:ir_wall]=1500
   Vp_solid=Vp[src_iz,ir_wall+1] # make sure Vp[1:ir_wall]=1500

   #Create initial field matrices
   vz_init=zeros(nz,nr)
   vphi_init=zeros(nz,nr)
   vr_init=zeros(nz,nr)

   trr_init=zeros(nz,nr)
   tpp_init=zeros(nz,nr)
   tzz_init=zeros(nz,nr)
   trp_init=zeros(nz,nr)
   trz_init=zeros(nz,nr)
   tpz_init=zeros(nz,nr)

   #At the spatial location of iz=src_iz, XXX_Schoenberg has been assigned in the time domain.
   #To populate the initial field matrices, I will assume that this wave propagates in space.
   #Note, XXX_Schoenberg was assumed as boundary conditions, where stress values are delayed by dt/2.

   # Without interpotaion (fast, but less accurate)
   #==
   is_zero=Int(round(delay/dt)+1)
   for iz=src_iz:nz
      is_tmp=Int(round((iz-src_iz)*dz/Vp_solid/dt))+is_zero
      if (is_tmp <= ns)
      for ir=1:nr
         vz_init[iz,ir]=vz_Schoenberg[is_tmp,ir];
         vr_init[iz,ir]=vr_Schoenberg[is_tmp,ir];
         vphi_init[iz,ir]=vphi_Schoenberg[is_tmp,ir];
         trr_init[iz,ir]=trr_Schoenberg2[is_tmp,ir];
         tpp_init[iz,ir]=tpp_Schoenberg2[is_tmp,ir];
         tzz_init[iz,ir]=tzz_Schoenberg2[is_tmp,ir];
         trp_init[iz,ir]=trp_Schoenberg2[is_tmp,ir];
         trz_init[iz,ir]=trz_Schoenberg2[is_tmp,ir];
         tpz_init[iz,ir]=tpz_Schoenberg2[is_tmp,ir];
      end
      end
   end

   for iz=1:src_iz-1
      is_tmp=Int(round((src_iz-iz)*dz/Vp_solid/dt))+is_zero
      if (is_tmp > 0)
      for ir=1:nr
         vz_init[iz,ir]=vz_Schoenberg[is_tmp,ir];
         vr_init[iz,ir]=vr_Schoenberg[is_tmp,ir];
         vphi_init[iz,ir]=vphi_Schoenberg[is_tmp,ir];
         trr_init[iz,ir]=trr_Schoenberg2[is_tmp,ir];
         tpp_init[iz,ir]=tpp_Schoenberg2[is_tmp,ir];
         tzz_init[iz,ir]=tzz_Schoenberg2[is_tmp,ir];
         trp_init[iz,ir]=trp_Schoenberg2[is_tmp,ir];
         trz_init[iz,ir]=trz_Schoenberg2[is_tmp,ir];
         tpz_init[iz,ir]=tpz_Schoenberg2[is_tmp,ir];
      end
      end
   end
   ==#

   # With interpolation (slow, but accurate)
   itvec_tmp=zeros(1,nz)
   for iz=src_iz:nz
      it_tmp=(iz-src_iz)*dz/Vp_solid+delay
      itvec_tmp[iz]=it_tmp
   end

   for iz=1:src_iz-1
      it_tmp=(src_iz-iz)*dz/Vp_solid+delay
      itvec_tmp[iz]=it_tmp
   end

   for iz=1:nz
#       println("itvec=",itvec_tmp[iz])
       if (itvec_tmp[iz] <= 0.0)
           itvec_tmp[iz]=0.0
       elseif (itvec_tmp[iz] >= T)
           itvec_tmp[iz]=T
       end
   end

   for ir=1:nr
        vz_spl=Spline1D(tvec_vel,vz_Schoenberg[:,ir],k=1)
        vr_spl=Spline1D(tvec_vel,vr_Schoenberg[:,ir],k=1)
        vphi_spl=Spline1D(tvec_vel,vphi_Schoenberg[:,ir],k=1)
        trr_spl=Spline1D(tvec_stress,trr_Schoenberg[:,ir],k=1)
        tpp_spl=Spline1D(tvec_stress,tpp_Schoenberg[:,ir],k=1)
        tzz_spl=Spline1D(tvec_stress,tzz_Schoenberg[:,ir],k=1)
        trp_spl=Spline1D(tvec_stress,trp_Schoenberg[:,ir],k=1)
        trz_spl=Spline1D(tvec_stress,trz_Schoenberg[:,ir],k=1)
        tpz_spl=Spline1D(tvec_stress,tpz_Schoenberg[:,ir],k=1)

        vz_init[:,ir]=vz_spl(itvec_tmp[:]);
        vr_init[:,ir]=vr_spl(itvec_tmp[:]);
        vphi_init[:,ir]=vphi_spl(itvec_tmp[:]);
        trr_init[:,ir]=trr_spl(itvec_tmp[:]);
        tpp_init[:,ir]=tpp_spl(itvec_tmp[:]);
        tzz_init[:,ir]=tzz_spl(itvec_tmp[:]);
        trp_init[:,ir]=trp_spl(itvec_tmp[:]);
        trz_init[:,ir]=trz_spl(itvec_tmp[:]);
        tpz_init[:,ir]=tpz_spl(itvec_tmp[:]);
    end

 return vz_init,vr_init,vphi_init,trr_init,tpp_init,tzz_init,
    trp_init,trz_init,tpz_init
end

#------------------------
# Solving BC problem in Peng (1994)
#------------------------
function SolvePengBC(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)

    kz=w/Vp_solid*cos(angle_delta*pi/180)
    ka=w/Vp_solid
    kb=w/Vs_solid
    kp=sqrt(ka^2-kz^2)
    ks=sqrt(kb^2-kz^2)
    kf=sqrt((w/Vp_fluid)^2-kz^2)

#@show R0phi(rb,Vp_solid,Vs_solid,Rho_solid,w,kp,kz)
    #BC: n=0. D*[A0;B0;C0]=Ep
    D=[U0f(rb,Vp_fluid,kf) -U0phi(rb,Vp_solid,kp) -U0xi(rb,Vs_solid,kb,ks,kz);
       R0f(rb,Vp_fluid,Rho_fluid,w,kf) -R0phi(rb,Vp_solid,Vs_solid,Rho_solid,w,kp,kz) -R0xi(rb,Vp_solid,Vs_solid,Rho_solid,w,ks,kz,kb);
       0.0  -Z0phi(rb,Vp_solid,Vs_solid,Rho_solid,w,kp,kz) -Z0xi(rb,Vs_solid,Rho_solid,w,ks,kz,kb)]
    Ep=[U0P(rb,Vp_solid,kp);
        R0P(rb,Vp_solid,Vs_solid,Rho_solid,w,kp,kz);
        Z0P(rb,Vp_solid,Vs_solid,Rho_solid,w,kp,kz)]

    Coef0=inv(D)*Ep
    A0,B0,C0=Coef0[1],Coef0[2],Coef0[3]

#==Schoenberg's coefficient (A, B, and C)
    cT=sqrt(Vp_fluid^2/(1+Rho_fluid*Vp_fluid^2/(Rho_solid*Vs_solid^2)))
    g=(Vs_solid/Vp_solid)^2
    KF=kf*rb
    KL=kp*rb
    tmp=Vs_solid^2/cT^2*(1-(cT/Vp_solid)^2*(cos(angle_delta*pi/180))^2)
    Coef0_approx=[-Vp_solid/Vp_fluid*(1-2.0*g*(cos(angle_delta*pi/180))^2)/tmp;
                   im*pi/4.0*(KL^2-KF^2*(1-2.0*g*(cos(angle_delta*pi/180))^2)^2/tmp);
                  -im*pi/2.0*KF^2*cos(angle_delta*pi/180)*(1-2.0*g*(cos(angle_delta*pi/180))^2)/tmp]
==#
    return A0,B0,C0
end

#------Peng's solutions----
#---vel---
#---f----
function U0f(r,Vp_fluid,kf)
    Vp_fluid*kf*(-besselj(1,kf*r))
end
function W0f(r,Vp_fluid,kf,kz)
    Vp_fluid*im*kz*besselj(0,kf*r)
end
#---phi---
function U0phi(r,Vp_solid,kp)
    Vp_solid*kp*(-hankelh1(1,kp*r))
end
function W0phi(r,Vp_solid,kp,kz)
    Vp_solid*im*kz*hankelh1(0,kp*r)
end
#---xi---
function U0xi(r,Vs_solid,kb,ks,kz)
    -Vs_solid*kz*ks/kb*(-hankelh1(1,ks*r))
end
function W0xi(r,Vs_solid,kb,ks)
    Vs_solid*im*ks^2/kb*hankelh1(0,ks*r)
end

#---stress---
#---f----
function R0f(r,Vp_fluid,Rho_fluid,w,kf)
    -Rho_fluid*Vp_fluid*w^2*besselj(0,kf*r)
end
#---phi---
#CAUTION
function R0phi(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz)
    -Rho_solid*Vp_solid*(
        (w^2-2.0*Vs_solid^2*kz^2)*hankelh1(0,kp*r)+
        2.0*Vs_solid^2*kp^2/(kp*r)*(-hankelh1(1,kp*r))
        )
end
function Z0phi(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz)
    2.0*im*Rho_solid*Vp_solid*Vs_solid^2*kz*kp*(-hankelh1(1,kp*r))
end
function L0phi(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz,ka)
    -Rho_solid*Vp_solid*(
        w^2-2.0*Vs_solid^2*(ka^2-kz^2)
        )*hankelh1(0,kp*r)
end
function M0phi(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz,ka)
    -Rho_solid*Vp_solid*(
        (Vp_solid^2-2.0*Vs_solid^2)*ka^2*hankelh1(0,kp*r)-
        2.0*Vs_solid^2*kp^2/(kp*r)*(-hankelh1(1,kp*r))
        )
end
#---xi---
function R0xi(r,Vp_solid,Vs_solid,Rho_solid,w,ks,kz,kb)
    2.0*Rho_solid*Vs_solid^3*kz*ks^2/kb*(
        hankelh1(0,ks*r)-1.0/(ks*r)*hankelh1(1,ks*r)
        )
end
function Z0xi(r,Vs_solid,Rho_solid,w,ks,kz,kb)
    -im*Rho_solid*Vs_solid^3*(kz^2-ks^2)*ks/kb*(-hankelh1(1,ks*r))
end
function L0xi(r,Vs_solid,Rho_solid,w,ks,kz,kb)
    -2.0*Rho_solid*Vs_solid^3*kz/kb*(kb^2-kz^2)*hankelh1(0,ks*r)
end
function M0xi(r,Vs_solid,Rho_solid,w,ks,kz,kb)
    -2.0*Rho_solid*Vs_solid^3*kz/kb*(
        -ks^2/(ks*r)*hankelh1(1,ks*r)
        )
end
#--plane P wave--
function U0P(r,Vp_solid,kp)
    Vp_solid*kp*(-besselj(1,kp*r))
end
function W0P(r,Vp_solid,kp,kz)
    Vp_solid*im*kz*besselj(0,kp*r)
end
function R0P(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz)
    -Rho_solid*Vp_solid*(
        (w^2-2.0*Vs_solid^2*kz^2)*besselj(0,kp*r)+
        2.0*Vs_solid^2*kp^2/(kp*r)*(-besselj(1,kp*r))
        )
end
function Z0P(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz)
    2.0*im*Rho_solid*Vp_solid*Vs_solid^2*kz*kp*(-besselj(1,kp*r))
end
function L0P(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz,ka)
    -Rho_solid*Vp_solid*(
        w^2-2.0*Vs_solid^2*(ka^2-kz^2)
        )*besselj(0,kp*r)
end
function M0P(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz,ka)
    -Rho_solid*Vp_solid*(
        (Vp_solid^2-2.0*Vs_solid^2)*ka^2*besselj(0,kp*r)-
        2.0*Vs_solid^2*kp^2/(kp*r)*(-besselj(1,kp*r))
        )
end

#---utilities--
function get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    kz=w/Vp_solid*cos(angle_delta*pi/180)
    ka=w/Vp_solid
    kb=w/Vs_solid
    kp=sqrt(ka^2-kz^2)
    ks=sqrt(kb^2-kz^2)
    kf=sqrt((w/Vp_fluid)^2-kz^2)
    return kz,ka,kb,kp,ks,kf
end
#---Potentials Evaluation--
#---f---
function ur_f(r,A0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return U0f(r,Vp_fluid,kf)*A0
end
function uz_f(r,A0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return W0f(r,Vp_fluid,kf,kz)*A0
end
function trr_f(r,A0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    #note: -p=trr=tpp=tzz
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return R0f(r,Vp_fluid,Rho_fluid,w,kf)*A0
end
#---phi---
function ur_phi(r,B0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return U0phi(r,Vp_solid,kp)*B0
end
function uz_phi(r,B0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return W0phi(r,Vp_solid,kp,kz)*B0
end
function trr_phi(r,B0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return R0phi(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz)*B0
end
function tzz_phi(r,B0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return L0phi(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz,ka)*B0
end
function tpp_phi(r,B0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return M0phi(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz,ka)*B0
end
function trz_phi(r,B0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return Z0phi(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz)*B0
end

#---xi---
function ur_xi(r,C0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return U0xi(r,Vs_solid,kb,ks,kz)*C0
end
function uz_xi(r,C0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return W0xi(r,Vs_solid,kb,ks)*C0
end
function trr_xi(r,C0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return R0xi(r,Vp_solid,Vs_solid,Rho_solid,w,ks,kz,kb)*C0
end
function tzz_xi(r,C0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return L0xi(r,Vs_solid,Rho_solid,w,ks,kz,kb)*C0
end
function tpp_xi(r,C0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return M0xi(r,Vs_solid,Rho_solid,w,ks,kz,kb)*C0
end
function trz_xi(r,C0,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return Z0xi(r,Vs_solid,Rho_solid,w,ks,kz,kb)*C0
end
#---plane P wave---
function ur_P(r,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return U0P(r,Vp_solid,kp)
end
function uz_P(r,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return W0P(r,Vp_solid,kp,kz)
end
function trr_P(r,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return R0P(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz)
end

function tzz_P(r,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return L0P(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz,ka)
end #my eq. --deduction done. check numerically

function tpp_P(r,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return M0P(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz,ka)
end #my eq. no deduction nor numerical

function trz_P(r,w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid,Rho_fluid,Rho_solid)
    kz,ka,kb,kp,ks,kf=get_wavenumbers(w,angle_delta,rb,Vp_fluid,Vp_solid,Vs_solid)
    return Z0P(r,Vp_solid,Vs_solid,Rho_solid,w,kp,kz)
end

#-------------------------------
# Initializing field variables
# (plane wave incidence)
#-------------------------------
function initialize_planewave(nr,nz,dr,dz,ir_wall,Rho,H,G,nt,tvec,dt,src_func,srcdepth,f0,delay)
   #poroelasticity version
   #plane wave version

   #ir_wall: assuming Vp[1:ir_wall]=1500
   println("--------------------------------")
   println("Plane-wave initial conditions...")

#---initializing field matrices
   vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf=init_fields_Por(nz,nr) #
# By using RightBC1D, and complete filling values till Right edge, no taper and offset are required.
   ntaper=0
   noffset=0

# Schoenberg's plane wave solution as boundary conditions
# We create elastic matrices (Vp,Vs,Rho) of the entire modeling domain
# but access only the elastic part (please make sure if this is the case).
   Vp_E=(H./Rho).^0.5
   Vs_E=(G./Rho).^0.5
   Rho_E=zeros(nz,nr)
   mycopy_mat(Rho,Rho_E,nz,nr)

   #
   vz_srcBC,vr_srcBC,vphi_srcBC,
   trr_srcBC,tpp_srcBC,tzz_srcBC,trp_srcBC,trz_srcBC,tpz_srcBC,
   src_iz=Peng_solution(nr,nz,dr,dz,ir_wall,tvec,src_func,srcdepth,f0,ntaper,noffset,Vp_E,Vs_E,Rho_E)

   #make sure tvec of vel and stress are correct
   tvec_vel=copy(tvec)
   tvec_stress=copy(tvec)
   tvec_stress=tvec_stress+dt/2*ones(nt) #default. see also Peng_solution
#   tvec_stress=tvec_stress-dt/2*ones(nt) #Needs pre-update. see also Peng_solution

   vz_init,vr_init,vphi_init,trr_init,tpp_init,tzz_init,
   trp_init,trz_init,tpz_init=Schoenberg_solution_initialBC(nr,nz,dr,dz,ir_wall,Vp_E,
                              tvec_vel,tvec_stress,
                              srcdepth,f0,delay,
                              vz_srcBC,vr_srcBC,vphi_srcBC,
                              trr_srcBC,tpp_srcBC,tzz_srcBC,trp_srcBC,trz_srcBC,tpz_srcBC,
                              ntaper,noffset)

   # Schoenberg's plane wave solution as initial condition
    println("Plane wave source: depth, index")
    @show srcdepth,src_iz
    println("Plane-wave initial conditions...finished")
    println("--------------------------------")

    return vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
         vz_srcBC,vr_srcBC,
         trr_srcBC,tpp_srcBC,tzz_srcBC,trz_srcBC,
         vz_init,vr_init,trr_init,tpp_init,tzz_init,
         trz_init,
         src_iz,noffset,
         f0

end # function init


text:
module BondLengthDistributions

include("bond_length_distribution.jl")

export BondLengthDistribution, AllBondLengthDistributions

end  # module BondLengthDistributions

max_stars_repo_name: IanAWatson/smu

max_issues_repo_name: IanAWatson/smu

max_forks_repo_name: IanAWatson/smu


text:
export GLFont, TextField, StyledTextSegment, TextCursor, FontProperties

immutable StyledTextSegment
	segment::UnitRange{Int}
	style::Dict{ASCIIString, Any}
end
type TextCursor
	x::Float32
	y::Float32
	intend::Float32
end

type TextField
	id::Symbol
	text::UTF8String
	newLineIndexes::Array{UnitRange{Int}, 1}
	#words::Array{SubString, 1}
	styles::Array{StyledTextSegment, 1}
	selection::UnitRange{Int}
	x::Float32
	y::Float32
	area::Shape
	hasFocus::Bool
	function TextField(id::Symbol, text::String, x::Real, y::Real, area::Shape)
		defaultStyle 	= @compat Dict{ASCIIString, Any}("textColor" => Float32[1,1,1,1], "backgroundColor" => Float32[0,0,0,0])
		styles 			= [StyledTextSegment(1:length(text), defaultStyle)]
		new(id, utf8(text), build_line_indexes(text), styles, length(text) : length(text)-1, float32(x), float32(y), area, false)
	end

	function TextField(id::Symbol,text::String, styles::Array{StyledTextSegment, 1}, x::Real, y::Real, area::Shape)
		new(id, utf8(text), build_line_indexes(text), styles, length(text):length(text)-1, float32(x), float32(y), area, false)
	end

	function TextField(id::Symbol,text::String,
			styles::Array{StyledTextSegment, 1},
			selection::Range,
			x::Real, y::Real, area::Shape)
		new(id, utf8(text), build_line_indexes(text), styles, selection, float32(x), float32(y), area, false)
	end
end

type FontProperties
	lineHeight::GLfloat
    advance::GLfloat
end

immutable GLFont
    data::Dict{Symbol, Any}
	props::Vector{Vec2}
	function GLFont(name::String)
	    flStream            = open("$(name).txt")
	    width               = int(readline(flStream))
	    height              = int(readline(flStream))
	    lineHeight::GLfloat = int(readline(flStream))
	    advance::GLfloat 	= 0

	    values                   = split(readline(flStream))
        charCodet                = char(int(values[1]))
        advancet                 = int(values[2])
        xt::GLfloat              = int(values[3]) / width
        x2t::GLfloat             = (int(values[3]) + advance) / width
        yt::GLfloat              =  ((height - int(values[4])  - lineHeight) / height)
        texLineHeightt::GLfloat  = lineHeight / height

        tmp_uv 	  = Array(Vec2, 1, 5)
        tmp_uv[1] = Vec2(advancet, lineHeight)
        tmp_uv[2] = Vec2(xt,  yt)
        tmp_uv[3] = Vec2(xt,  yt + texLineHeightt)
        tmp_uv[4] = Vec2(x2t, yt + texLineHeightt)
        tmp_uv[5] = Vec2(x2t, yt)

	    uv = vcat(tmp_uv)
	    i = 1
	    for line in eachline(flStream)
	        values                  = split(line)
	        charCode                = char(int(values[1]))
	        advance                 = int(values[2])
	        x::GLfloat              = int(values[3]) / width
	        x2::GLfloat             = (int(values[3]) + advance) / width
	        y::GLfloat              = ((height - int(values[4]) - lineHeight) / height)
	        texLineHeight::GLfloat  = lineHeight / height
	        tmp_uv[1] = Vec2(advance, lineHeight)
	        tmp_uv[2] = Vec2(x,  y)
	        tmp_uv[3] = Vec2(x,  y + texLineHeight)
	        tmp_uv[4] = Vec2(x2, y + texLineHeight)
	        tmp_uv[5] = Vec2(x2, y)
	        uv = vcat(uv, tmp_uv)
	        i += 1
	    end
	    close(flStream)
		data = @compat Dict{Symbol, Any}(
			:dontdelete_uv_index 		=> GLBuffer(GLint[1:6;], 1),
			:dontdelete_indexes 		=> indexbuffer(GLuint[0:5;]),
			:dontdelete_uv 				=> Texture(uv),
			:dontdelete_font_texture 	=> Texture("$(name).bmp")
		)
	    new(data, vec(uv[1,:]))
	end
end


text:
"""
    TF(r, Z)

Approximate solution of the [Thomas-Fermi equation](https://en.wikipedia.org/wiki/Thomas-Fermi_model)
for neutral atom with nuclear charge `Z`. Returns density at the distance `r`.

See Moliere, G. (1947). Theorie der streuung schneller geladener teilchen i. einzelstreuung am abgeschirmten coulomb-feld.
Zeitschrift für Naturforschung A, 2(3), 133-145.

## Example
For a density function ``\\rho(r)`` of ``N`` particles the following equation holds.
```math
N = 4\\pi\\int dr \\rho(r)r^2
```
On the grid ``r_i = \\exp(x_i)`` it is
```math
N = 4\\pi\\int dx \\rho(r(x))r^3
```
```jldoctest
julia> x = -30:0.1:20
-30.0:0.1:20.0

julia> r, n, dx = exp.(x), length(x), step(x);

julia> isapprox(4π * dx * sum(TF.(r, 18) .* r .^3), 18, atol=1e-10)
true
```
"""
function TF(r, Z)
  b = 1/2*(3π/4)^(2/3)*Z^(-1/3)
  x = r/b

  B₁, B₂, B₃ = 0.1, 0.55, 0.35
  β₁, β₂, β₃ = 6.0, 1.20, 0.30

  Z/(4π*b^3*x)*(B₁*β₁^2*exp(-β₁*x) + B₂*β₂^2*exp(-β₂*x) + B₃*β₃^2*exp(-β₃*x))
end


text:
const pyjlbasetype = pynew()

_pyjl_getvalue(x) = @autopy x C.PyJuliaValue_GetValue(getptr(x_))

_pyjl_setvalue!(x, v) = @autopy x C.PyJuliaValue_SetValue(getptr(x_), v)

pyjl(t, v) = pynew(errcheck(@autopy t C.PyJuliaValue_New(getptr(t_), v)))
export pyjl

"""
    pyisjl(x)

Test whether `x` is a wrapped Julia value, namely an instance of `juliacall.ValueBase`.
"""
pyisjl(x) = pytypecheck(x, pyjlbasetype)
export pyisjl

pyjlisnull(x) = @autopy x begin
    if pyisjl(x_)
        C.PyJuliaValue_IsNull(getptr(x_))
    else
        error("Expecting a 'juliacall.ValueBase', got a '$(pytype(x_).__name__)'")
    end
end

"""
    pyjlvalue(x)

Extract the value from the wrapped Julia value `x`.
"""
pyjlvalue(x) = @autopy x begin
    if pyjlisnull(x_)
        error("Julia value is NULL")
    else
        _pyjl_getvalue(x_)
    end
end
export pyjlvalue

function init_jlwrap_base()
    setptr!(pyjlbasetype, incref(C.POINTERS.PyJuliaBase_Type))
    pyjuliacallmodule.ValueBase = pyjlbasetype
end

pyconvert_rule_jlvalue(::Type{T}, x::Py) where {T} = pyconvert_tryconvert(T, _pyjl_getvalue(x))

function C._pyjl_callmethod(f, self_::C.PyPtr, args_::C.PyPtr, nargs::C.Py_ssize_t)
    @nospecialize f
    if C.PyJuliaValue_IsNull(self_)
        errset(pybuiltins.TypeError, "Julia object is NULL")
        return C.PyNULL
    end
    in_f = false
    self = C.PyJuliaValue_GetValue(self_)
    try
        if nargs == 1
            in_f = true
            ans = f(self)::Py
            in_f = false
        elseif nargs == 2
            arg1 = pynew(incref(C.PyTuple_GetItem(args_, 1)))
            in_f = true
            ans = f(self, arg1)::Py
            in_f = false
            pydel!(arg1)
        elseif nargs == 3
            arg1 = pynew(incref(C.PyTuple_GetItem(args_, 1)))
            arg2 = pynew(incref(C.PyTuple_GetItem(args_, 2)))
            in_f = true
            ans = f(self, arg1, arg2)::Py
            in_f = false
            pydel!(arg1)
            pydel!(arg2)
        elseif nargs == 4
            arg1 = pynew(incref(C.PyTuple_GetItem(args_, 1)))
            arg2 = pynew(incref(C.PyTuple_GetItem(args_, 2)))
            arg3 = pynew(incref(C.PyTuple_GetItem(args_, 3)))
            in_f = true
            ans = f(self, arg1, arg2, arg3)::Py
            in_f = false
            pydel!(arg1)
            pydel!(arg2)
            pydel!(arg3)
        else
            errset(pybuiltins.NotImplementedError, "__jl_callmethod not implemented for this many arguments")
        end
        ptr = getptr(ans)
        pystolen!(ans)
        return ptr
    catch exc
        if exc isa PyException
            C.PyErr_Restore(incref(getptr(exc._t)), incref(getptr(exc._v)), incref(getptr(exc._b)))
            return C.PyNULL
        else
            try
                if in_f
                    return pyjl_handle_error(f, self, exc)
                else
                    errset(pyJuliaError, pytuple((pyjlraw(exc), pyjlraw(catch_backtrace()))))
                    return C.PyNULL
                end
            catch
                errset(pyJuliaError, "an error occurred while setting an error")
                return C.PyNULL
            end
        end
    end
end

function pyjl_handle_error(f, self, exc)
    @nospecialize f self exc
    t = pyjl_handle_error_type(f, self, exc)::Py
    if ispynull(t)
        # NULL => raise JuliaError
        errset(pyJuliaError, pytuple((pyjlraw(exc), pyjlraw(catch_backtrace()))))
        return C.PyNULL
    elseif pyistype(t)
        # Exception type => raise this type of error
        errset(t, string("Julia: ", Py(sprint(showerror, exc))))
        return C.PyNULL
    else
        # Otherwise, return the given object (e.g. NotImplemented)
        return incref(getptr(t))
    end
end

function pyjl_methodnum(f)
    @nospecialize f
    C.PyJulia_MethodNum(f)
end

function pyjl_handle_error_type(f, self, exc)
    @nospecialize f self exc
    PyNULL
end


text:
__precompile__()

module mbNewton
export newton, centralDifference, forwardDifference, broyden, bisection

# base function
"""
    newton(homotopy, jacobian, v₀, pred[, init, callback, useOpt])
Newton's method. For example: `v1 = newton(H, J, v0, predEps(0.001))`.
"""
function newton(H::Function, J::Function, v₀::Vector{Float64}, pred::Function;
	init=Void, callback=Void, useOpt::Bool=false)

	local v = deepcopy(v₀)
	local tmpH, tmpJ
	tmpH, tmpJ = H(v), J(v)

	init≠Void && init()
	callback≠Void && callback(v, tmpH, tmpJ)

	if useOpt
		local optv, optH, tmpH2
		optv, optH = v, norm(tmpH)
	end

	while pred(tmpH)
		# v -= size(tmpJ,1)==size(tmpJ,2) ? tmpJ\tmpH : tmpJ'*inv(tmpJ*tmpJ')*tmpH
		v -= tmpJ\tmpH

		tmpH, tmpJ = H(v), J(v)

		callback≠Void && callback(v, tmpH, tmpJ)

		if useOpt && (tmpH2 = norm(tmpH)) < optH
			optv, optH = v, tmpH2
		end
	end

  return useOpt ? optv : v
end

"""
    centralDifference(homotopy, v, epsilon)
	centralDifference(homotopy, epsilon)
Same as `forwardDifference`, but with the symmetric difference quotient.
"""
function centralDifference(H::Function, v::Vector{Float64}; ϵ::Float64=1e-4)
	local J = cell(length(v))
	local w = deepcopy(v)
	for i in 1:length(v)
		w[i] = v[i]+ϵ; 	J[i] = H(w)
		w[i] = v[i]-ϵ; 	J[i] -= H(w)
	 	J[i] /= 2ϵ;		w[i] = v[i]
	end
	return reduce(hcat, J)
end
centralDifference(H::Function; ϵ=1e-4) = v -> centralDifference(H, v; ϵ)

"""
    forwardDifference(homotopy, v, epsilon)
	forwardDifference(homotopy, epsilon)
Returns the difference quotient in `v` or an approximate jacobian using this method.
"""
function forwardDifference(H::Function, v::Vector{Float64}; ϵ::Float64=1e-4)
	local J = cell(length(v))
	local w = deepcopy(v)
	local tmpH = H(v)
	for i in 1:length(v)
	  w[i] = v[i]+ϵ;	J[i] = (H(w) - tmpH) ./ ϵ
	  w[i] = v[i]
	end
	return reduce(hcat, J)
end
forwardDifference(H::Function; ϵ=1e-4) = v -> forwardDifference(H, v; ϵ=ϵ)

# WARNING inefficient, H is evaluated which is not necessary (at least in Newton context...
#		asymptotic behavior is unchanged though) however, this allows to construct a pseudo-
#		Jacobian with the usual one-parameter-signature, instead of requiring to adapt any functions.
# can be countered through rudimentary caching in H
# WARNING mind that this function has an internal state depending on the preceding
#		evaluations. the 'Jacobian' at a point is thus not well defined.
# TODO sequence of evaluation not optimal
"""
    broyden(homotopy, jacobian)
Creates an approximate jacobian based on broyden-updates and maintainance of an internal
state. Consider this for performance improvements.
"""
function broyden(H, J)
	local H₀,J₀,v₀							# closure
	local init = true

	return function (v₁)
		if init
			init = false
			H₀,J₀,v₀ = H(v₁),J(v₁),v₁
			return J₀
		end

		local H₁, dv, dH
		H₁ = H(v₁) 								# ← hmm...
		dv, dH = v₁-v₀, H₁-H₀
		v₀, H₀ = v₁, H₁

		norm(dv)>1e-10 && (J₀ += (dH-J₀*dv) ./ norm(dv)^2 * dv')

		return J₀
	end
end

end


text:
#Linear ARD Covariance Function

"""
    LinArd <: Kernel

ARD linear kernel (covariance)
```math
k(x,x') = xᵀL⁻²x'
```
with length scale ``ℓ = (ℓ₁, ℓ₂, …)`` and ``L = diag(ℓ₁, ℓ₂, …)``.
"""
mutable struct LinArd{T<:Real} <: Kernel
    "Length scale"
    ℓ::Vector{T}
    "Priors for kernel parameters"
    priors::Array
end

"""
Linear ARD Covariance Function

    LinArd(ll::Vector{T})

# Arguments
  - `ll::Vector{Real}`: vector of length scales (given on log scale)
"""
LinArd(ll::Vector{T}) where T = LinArd{T}(exp.(ll), [])

cov(lin::LinArd, x::AbstractVector, y::AbstractVector) = dot(x./lin.ℓ, y./lin.ℓ)

struct LinArdData{D<:AbstractArray} <: KernelData
    XtX_d::D
end

function KernelData(k::LinArd, X1::AbstractMatrix, X2::AbstractMatrix)
    dim1, nobs1 = size(X1)
    dim2, nobs2 = size(X2)
	@assert dim1==dim2
	dim = dim1
    XtX_d = Array{eltype(X2)}(undef, nobs1, nobs2, dim)
    @inbounds @simd for d in 1:dim
        for i in 1:nobs1
            for j in 1:nobs2
                XtX_d[i, j, d] = X1[d, i] * X2[d, j]
            end
        end
    end
    LinArdData(XtX_d)
end
kernel_data_key(k::LinArd, X1::AbstractMatrix, X2::AbstractMatrix) = "LinArdData"
function cov(lin::LinArd, X::AbstractMatrix)
    K = (X./lin.ℓ)' * (X./lin.ℓ)
    LinearAlgebra.copytri!(K, 'U')
    return K
end
# function cov!(cK::AbstractMatrix, lin::LinArd, X::AbstractMatrix, data::LinArdData)
    # dim, nobs = size(X)
    # fill!(cK, 0)
    # for d in 1:dim
        # LinearAlgebra.axpy!(1/lin.ℓ[d]^2, view(data.XtX_d,1:nobs, 1:nobs ,d), cK)
    # end
    # return cK
# end
# function cov(lin::LinArd, X::AbstractMatrix, data::LinArdData)
    # nobs = size(X,2)
    # K = Array{eltype(X)}(undef, nobs, nobs)
    # cov!(K, lin, X, data)
# end
@inline @inbounds function cov_ij(lin::LinArd, X1::AbstractMatrix, X2::AbstractMatrix, i::Int, j::Int, dim::Int)
    ck = 0.0
    for d in 1:dim
        ck += dotijp(X1,X2,i,j,d) * 1/lin.ℓ[d]^2
    end
    return ck
end
@inline @inbounds function cov_ij(lin::LinArd, X1::AbstractMatrix, X2::AbstractMatrix, data::LinArdData, i::Int, j::Int, dim::Int)
    ck = 0.0
    for d in 1:dim
        ck += data.XtX_d[i,j,d] * 1/lin.ℓ[d]^2
    end
    return ck
end
get_params(lin::LinArd) = log.(lin.ℓ)
get_param_names(lin::LinArd) = get_param_names(lin.ℓ, :ll)
num_params(lin::LinArd) = length(lin.ℓ)

function set_params!(lin::LinArd, hyp::AbstractVector)
    length(hyp) == num_params(lin) || throw(ArgumentError("Linear ARD kernel has $(num_params(lin)) parameters"))
    @. lin.ℓ = exp(hyp)
end

@inline dk_dll(lin::LinArd, xy::Real, d::Int) = -2 * xy / lin.ℓ[d]^2
@inline function dKij_dθp(lin::LinArd, X1::AbstractMatrix, X2::AbstractMatrix, i::Int, j::Int, p::Int, dim::Int)
    if p<=dim
        return dk_dll(lin, dotijp(X1,X2,i,j,p), p)
    else
        return NaN
    end
end
@inline function dKij_dθp(lin::LinArd, X1::AbstractMatrix, X2::AbstractMatrix, data::LinArdData, i::Int, j::Int, p::Int, dim::Int)
    if p <= dim
        return dk_dll(lin, data.XtX_d[i,j,p],p)
    else
        return NaN
    end
end


text:
if basename(pwd()) == "aoc"
    cd("2021/1")
end

function part1(filename::AbstractString)
    count(x > 0 for x in diff(parse.(Int, eachline(filename))))
end

@assert part1("example.txt") == 7

function part2(filename::AbstractString)
    xs = parse.(Int, eachline(filename))
    count(x > 0 for x in xs[begin+3:end] .- xs[begin:end-3])
end

@assert part2("example.txt") == 5

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using BinaryBuilder
using GitHub
using Base64

# This is a global github authentication token that is set the first time
# we authenticate and then reused
const _github_auth = Ref{GitHub.Authorization}()
function github_auth()
    if !isassigned(_github_auth) || _github_auth[] isa GitHub.AnonymousAuth
        # If the user is feeding us a GITHUB_AUTH token, use it
        auth = get(ENV, "GITHUB_AUTH", nothing)
        _github_auth[] = (auth === nothing ? GitHub.AnonymousAuth() : GitHub.authenticate(auth))
    end
    return _github_auth[]
end

function get_julia_version(ref::String="master")
    file = GitHub.file("JuliaLang/julia", "VERSION";
                       params = Dict("ref" => ref), auth=github_auth())
    @assert file.encoding == "base64" # GitHub says this will always be the case
    str = String(base64decode(chomp(file.content)))
    return VersionNumber(chomp(str))
end

"""
    version_id = build_julia(ref::String="master"; binarybuilder_args::Vector{String}=String["--verbose"])

Download and build julia at git reference `ref` using BinaryBuilder. Return the `version_id` (what other functions use
to identify this build).
"""
function build_julia(ref::String="master"; binarybuilder_args::Vector{String}=String["--verbose"])
    # This errors if `ref` cannot be found, error message is pretty ok
    version = get_julia_version(ref)

    reference = GitHub.reference("JuliaLang/julia", "heads/$(ref)"; handle_error=false, auth=github_auth())
    if reference.object == nothing
        reference = GitHub.reference("JuliaLang/julia", "tags/$(ref)"; handle_error=false, auth=github_auth())
    end
    if reference.object == nothing
        commit_hash = GitHub.commit("JuliaLang/julia", ref; auth=github_auth()).sha
    else
        commit_hash = reference.object["sha"]
    end

    if version.prerelease != ()
        contrib_path = joinpath(Sys.BINDIR, "..", "..", "contrib")
        found_buildname = false
        if isdir(contrib_path)
            try
                cd(contrib_path) do
                    version = VersionNumber(read(`./commit-name.sh $commit_hash`, String))
                end
                found_buildname = true
            catch e
                @error "failed to get build number" exception=e
            end
        end
        if !found_buildname
            version = VersionNumber(string(version) * string("-", commit_hash[1:6]))
        end
    end

    # Collection of sources required to build julia
    sources = [
        "https://github.com/JuliaLang/julia.git" => commit_hash,
    ]

    # Bash recipe for building across all platforms
    script = raw"""
    cd $WORKSPACE/srcdir
    mount -t devpts -o newinstance jrunpts /dev/pts
    mount -o bind /dev/pts/ptmx /dev/ptmx
    cd julia
    cat > Make.user <<EOF
    LLVM_ASSERTIONS=1
    LIBSSH2_ENABLE_TESTS=0
    JULIA_CPU_TARGET=generic;sandybridge,-xsaveopt,clone_all;haswell,-rdrnd,base(1)
    EOF
    make -j${nproc}
    make install
    cp LICENSE.md ${prefix}
    contrib/fixup-libgfortran.sh ${prefix}/lib/julia
    contrib/fixup-libstdc++.sh ${prefix}/lib ${prefix}/lib/julia
    """

    # These are the platforms we will build for by default, unless further
    # platforms are passed in on the command line
    platforms = [
        Linux(:x86_64, libc=:glibc)
    ]

    # The products that we will ensure are always built
    products = [
        ExecutableProduct("julia", :julia)
    ]

    # Dependencies that must be installed before this package can be built
    dependencies = []

    # Build the tarballs, and possibly a `build.jl` as well.
    product_hashes = cd(joinpath(@__DIR__, "..", "deps")) do
        build_tarballs(binarybuilder_args, "julia", version, sources, script, platforms, products, dependencies, preferred_gcc_version=v"7", skip_audit=true)
    end
    tarball, hash = product_hashes[platforms[1]]

    version_stanza = """
    ["$version"]
    file = "$tarball"
    sha = "$hash"
    """
    download_path = joinpath(@__DIR__, "..", "deps", "downloads")
    if isfile(joinpath(download_path, tarball))
        @warn "$tarball already exists in deps/downloads folder. Not copying."
        println("You may manually copy the file from products/ and add the following stanza to Versions.toml:")
        println(version_stanza)
    else
        mkpath(download_path)
        cp(joinpath(@__DIR__, "..", "deps", "products", tarball), joinpath(download_path, tarball))
        open(versions_file(); append=true) do f
            println(f, version_stanza)
        end
    end
    return version
end


text:
module SurvivalSignature

using LinearAlgebra,
    IterTools,
    ProgressMeter,
    Statistics,
    Random,
    Distributions,
    Distributed,
    BigCombinatorics

export survivalsignature,
    exactentry,
    approximateentry,
    gridnetwork,
    percolation,
    s_t_connectivity,
    reliability,
    percolation_preprocessor!

include("percolation.jl")
include("preprocessors.jl")
include("reliability.jl")
include("signature.jl")
include("structurefunctions.jl")
include("util.jl")

end # module


text:
"""
The **pcl_keypoints** library contains implementations of two point cloud
keypoint detection algorithms

http://docs.pointclouds.org/trunk/group__keypoints.html

## Exports

$(EXPORTS)
"""
module PCLKeyPoints

using DocStringExtensions
using LibPCL
using PCLCommon
using Cxx

const libpcl_keypoints = LibPCL.find_library_e("libpcl_keypoints")
try
    Libdl.dlopen(libpcl_keypoints, Libdl.RTLD_GLOBAL)
catch e
    warn("You might need to set DYLD_LIBRARY_PATH to load dependencies proeprty.")
    rethrow(e)
end

end # module


text:
import Base: sum, next, start, done, IndexStyle

abstract type AbstractDeviceArray{T, N} <: AbstractArray{T, N} end

IndexStyle(::AbstractDeviceArray) = IndexLinear()
start(x::AbstractDeviceArray) = UInt32(1)
next(x::AbstractDeviceArray, state::UInt32) = x[state], state + UInt32(1)
done(x::AbstractDeviceArray, state::UInt32) = state > length(x)

function sum(A::AbstractDeviceArray{T}) where T
    acc = zero(T)
    for elem in A
        acc += elem
    end
    acc
end

const shmem_counter = Ref{Int}(0)

"""
Creates a local static memory shared inside one block.
"""
macro LocalMemory(state, T, N)
    id = (shmem_counter[] += 1)
    quote
        lémem = LocalMemory($(esc(state)), $(esc(T)), Val{$(esc(N))}(), Val{$id}())
        AbstractDeviceArray(lémem, $(esc(N)))
    end
end

export @LocalMemory

"""
Creates a block local array pointer with `T` being the element type
and `N` the length. Both T and N need to be static! C is a counter for
approriately get the correct Local mem id in CUDAnative.
This is an internal method which needs to be overloaded by the GPU Array backends
"""
function LocalMemory(state, ::Type{T}, ::Val{N}, ::Val{C}) where {N, T, C}
    error("Not implemented")
end


text:
using LiquidSDR
using PyPlot

fs = 1000; # sampling frequency

# test signal generation
t = 0:1.0/fs:1;
s1 = sin(2*pi*t*10);  # 10Hz sinusoid
s2 = sin(2*pi*t*100); # 100Hz sinusoid
s = s1+s2; # mixed signal

# let's build butterworth filter
fc = 50.0 / fs; # normalized cut-off frequency (50Hz)
f0 = 51.0 / fs; # f0 is not considered for Low-Pass filters( so using a dummy value)
b,a = butterworth(fc, f0);

filt = createiirfilter(b,a);

# let's try to remove 100Hz component!
sig_result = applyfilter(filt, s);
plot(t, sig_result);


text:
"""
    LinearContact{T,N} <: Contact{T,N}

    contact object for impact and friction with a linearized friction cone

    friction_coefficient: value of friction coefficient
    surface_projector: mapping from world frame to surface tangent frame
    surface_normal_projector: inverse/complement of surface_projector
    contact_point: position of contact on Body relative to center of mass
    offset: position of contact relative to contact_point
"""
mutable struct LinearContact{T,N} <: Contact{T,N}
    friction_coefficient::T
    surface_projector::SMatrix{4,3,T,12}
    surface_normal_projector::Adjoint{T,SVector{3,T}} # inverse matrix
    contact_point::SVector{3,T}
    offset::SVector{3,T}

    function LinearContact(body::Body{T}, normal::AbstractVector, friction_coefficient;
        contact_point=szeros(T, 3),
        offset::AbstractVector=szeros(T, 3)) where T
        V1, V2, V3 = orthogonal_columns(normal)
        A = [V1 V2 V3]
        Ainv = inv(A)
        surface_normal_projector = Ainv[3, SA[1; 2; 3]]'
        surface_projector = SA{T}[
             1  0  0
            -1  0  0
             0  1  0
             0 -1  0
        ]
        new{Float64,12}(friction_coefficient, surface_projector, surface_normal_projector, contact_point, offset)
    end
end

function constraint(mechanism, contact::ContactConstraint{T,N,Nc,Cs}) where {T,N,Nc,Cs<:LinearContact{T,N}}
    model = contact.model
    body = get_body(mechanism, contact.parent_id)
    x2, v25, q2, ϕ25 = current_configuration_velocity(body.state)
    x3, q3 = next_configuration(body.state, mechanism.timestep)

    # transforms the velocities of the origin of the link into velocities along all 4 axes of the friction pyramid
    # vp = V(cp, B / W)_w velocity of the contact point cp, attached to body B wrt world frame, expressed in the world frame.
    vp = v25 + skew(vector_rotate(ϕ25, q3)) * (vector_rotate(model.contact_point, q3) - model.offset)
    γ = contact.impulses[2][1]
    sγ = contact.impulses_dual[2][1]
    ψ = contact.impulses[2][2]
    sψ = contact.impulses_dual[2][2]
    β = contact.impulses[2][@SVector [3,4,5,6]]
    sβ = contact.impulses_dual[2][@SVector [3,4,5,6]]
    SVector{6,T}(
        model.surface_normal_projector * (x3 + vector_rotate(model.contact_point,q3) - model.offset) - sγ,
        model.friction_coefficient * γ - sum(β) - sψ,
        (model.surface_projector * vp + ψ * sones(4) - sβ)...)
end

function constraint_jacobian_configuration(model::LinearContact, x3::AbstractVector, q3::Quaternion,
    x2::AbstractVector, v25::AbstractVector, q2::Quaternion, ϕ25::AbstractVector,
    λ, timestep)
    V = [model.surface_normal_projector;
         szeros(1,3);
         szeros(4,3)]
    ∂v∂q3 = skew(vector_rotate(ϕ25, q3)) * ∂vector_rotate∂q(model.contact_point, q3)
    ∂v∂q3 += skew(model.offset - vector_rotate(model.contact_point, q3)) * ∂vector_rotate∂q(ϕ25, q3)
    Ω = [model.surface_normal_projector * ∂vector_rotate∂q(model.contact_point, q3);
        szeros(1,4);
        model.surface_projector * ∂v∂q3]
    return [V Ω]
end

function constraint_jacobian_velocity(model::LinearContact, x3::AbstractVector, q3::Quaternion,
    x2::AbstractVector, v25::AbstractVector, q2::Quaternion, ϕ25::AbstractVector,
    λ, timestep)
    V = [model.surface_normal_projector * timestep;
         szeros(1,3);
         model.surface_projector]
    ∂v∂q3 = skew(vector_rotate(ϕ25, q3)) * ∂vector_rotate∂q(model.contact_point, q3)
    ∂v∂q3 += skew(model.offset - vector_rotate(model.contact_point, q3)) * ∂vector_rotate∂q(ϕ25, q3)
    ∂v∂ϕ25 = skew(model.offset - vector_rotate(model.contact_point, q3)) * ∂vector_rotate∂p(ϕ25, q3)
    Ω = [model.surface_normal_projector * ∂vector_rotate∂q(model.contact_point, q3) * rotational_integrator_jacobian_velocity(q2, ϕ25, timestep);
        szeros(1,3);
        model.surface_projector * (∂v∂ϕ25 + ∂v∂q3 * rotational_integrator_jacobian_velocity(q2, ϕ25, timestep))]
    return [V Ω]
end

function set_matrix_vector_entries!(mechanism::Mechanism, matrix_entry::Entry, vector_entry::Entry,
    contact::ContactConstraint{T,N,Nc,Cs,N½}) where {T,N,Nc,Cs<:LinearContact{T,N},N½}
    # ∇impulses[impulses .* impulses - μ; g - s] = [diag(impulses); -diag(0,1,1)]
    # ∇impulses[impulses .* impulses - μ; g - s] = [diag(impulses); -diag(1,0,0)]
    # (friction_coefficient γ - ψ) dependent of ψ = impulses[2][1:1]
    # B(z) * zdot - sβ dependent of sβ = impulses[2][2:end]
    friction_coefficient = contact.model.friction_coefficient
    γ = contact.impulses[2] + REG * neutral_vector(contact.model)
    s = contact.impulses_dual[2] + REG * neutral_vector(contact.model)

    ∇s1 = Diagonal(γ) # 6x6
    ∇s2 = Diagonal(-sones(T,6))
    ∇s = vcat(∇s1, ∇s2) # 12x6

    ∇γ1 = Diagonal(s) # 6x6
    ∇γ2 = @SMatrix[ 0  0  0  0  0  0;
                   friction_coefficient  0 -1 -1 -1 -1;
                    0  1  0  0  0  0;
                    0  1  0  0  0  0;
                    0  1  0  0  0  0;
                    0  1  0  0  0  0;]
    ∇γ = vcat(∇γ1, ∇γ2) # 12x6
    matrix_entry.value = hcat(∇s, ∇γ)

    # [-impulses .* impulses + μ; -g + s]
    vector_entry.value = vcat(-complementarityμ(mechanism, contact), -constraint(mechanism, contact))
    return
end


text:
println("===========================================================================")
println("Beginning execution of script GMHPhotoReceptor.jl/scripts/PhotoReceptor1.jl")
println("===========================================================================")

using GeneralizedMetropolisHastings
using GMHPhotoReceptor

###Values of the model
numvilli1 = 30000

#specify the values that determine the priors on the parameters
latencylocation = (2.0,3.5) #uniform distribution with (low,high) values
latencyscale = (0.2,0.7) #uniform distribution with (low,high) values
refractorylocation = (4.0,6.0) #uniform distribution with (low,high) values
refractoryscale = (1.5,2.5) #uniform distribution with (low,high) values
bumpamplitude = (3.0,5.0) #uniform distribution with (low,high) values
bumpshape = (log(3.0),0.3) #lognormal distribution with (location,scale) values; variable can vary roughly between 2.0 and 4.0, but becomes increasingly penalized outside
bumpscale = (log(2.5),0.3) #lognormal distribution with (location,scale) values; variable can vary roughly between 1.5 and 3.5, but becomes increasingly penalized outside

photonfilename = joinpath(splitdir(@__FILE__())[1],"../data/naturallight.jld")
photons1 = photonsequence(photonfilename)
current1 = lightinducedcurrent(photonfilename)

modelpolicy1 = policy(:photoreceptor) #4-parameter model with stochastic lognormal latency and refractory parameters and fixed bump parameters
params1 = parameters(:photoreceptor,modelpolicy1,latencylocation,latencyscale,refractorylocation,refractoryscale)

####Variance for the noise model, estimated from previous runs
variance1 = [3600.0]

println("==========================================")
println("Simulation parameters defined successfully")
println("==========================================")

###Create a PhotoReceptor model
model1 = model(:photoreceptor,params1,photons1,current1,variance1,numvilli1,modelpolicy1)

###Show the model
println("==========================")
println("Model defined successfully")
println("==========================")
show(model1)

###Create a Metropolis sampler with a Normal proposal density
propcov = [0.1,0.01,0.1,0.01]
sampler1 = sampler(:mh,:normal,0.1,propcov)
println("============================")
println("Sampler defined successfully")
println("============================")
show(sampler1)

#MCMC iteration specifications
niterations = 20
npropbase1 = 5000
npropbase2 = 1000
ntunerperiod = npropbase1*niterations

###Create a tuner that just reports the time
tuner1 = tuner(:monitor,ntunerperiod)
println("==========================")
println("Tuner defined successfully")
println("==========================")
show(tuner1)

println("Running a simulation with Standard MH")
runnerpolicy0 = policy(:mh,1;initialize=:prior)
runner0 = runner(runnerpolicy0,npropbase2*niterations)
show(runner0)
@time chain0 = run!(runner0,model1,sampler1,tuner1)
smhruntime = chain0.runtime
smhproposed = chain0.proposed

println("Running a simulation with Generalized MH in 1 process")
runnerpolicy1 = policy(:mh,npropbase2;initialize=:prior)
runner1 = runner(runnerpolicy1,niterations,npropbase2)
@time chain1 = run!(runner1,model1,sampler1,tuner1)
gmhruntime = chain1.runtime
gmhproposed = chain1.proposed

NWORKERS = [1,2,4,8,12,16,24,32]
proposalstrategy = ["total","segment"]
nproposals = [npropbase1*ones(Int,length(NWORKERS)) npropbase2*max(NWORKERS,1)]
parallelruntimes = zeros(length(NWORKERS),2)
parallelproposed = zeros(Int,length(NWORKERS),2)

for (j,s) in enumerate(proposalstrategy)
    println("=======================================")
    println("Testing fixed number of $s proposals")
    println("=======================================")
    for (i,n) in enumerate(NWORKERS)
        try
            addprocs(n)

            @everywhere using GeneralizedMetropolisHastings
            @everywhere using GMHPhotoReceptor

            println("Running simulation with $(nworkers()) workers: ")

            ###Create a Generalized Metropolis-Hastings runner
            runnerpolicy2 = policy(:mh,nproposals[i,j];initialize=:prior)
            runner2 = runner(runnerpolicy2,niterations,nproposals[i,j])
            @time chain2 = run!(runner2,model1,sampler1,tuner1)

            ###Show the result of the simulations
            show(chain2)

            parallelruntimes[i,j] = chain2.runtime
            parallelproposed[i,j] = chain2.proposed

        catch e
            println(e)
        end
        rmprocs(workers())
    end
end

println("Strategy \t Runtime \t Proposals \t Proposal Rate")
println("=========================================================================")
println("Standard MH \t $smhruntime \t $smhproposed \t $(smhproposed/smhruntime)")
println("GMH No Workers \t $gmhruntime \t $gmhproposed \t $(gmhproposed/gmhruntime)")
for i=1:length(NWORKERS)
    println("TOT $(NWORKERS[i]) Workers \t $(parallelruntimes[i,1]) \t $(parallelproposed[i,1]) \t $(parallelproposed[i,1]/parallelruntimes[i,1])")
end
for i=1:length(NWORKERS)
    println("SEG $(NWORKERS[i]) Workers \t $(parallelruntimes[i,2]) \t $(parallelproposed[i,2]) \t $(parallelproposed[i,2]/parallelruntimes[i,2])")
end


text:
"""
    ArrayAllocators.zeros(T=Float64, dims::Integer...)
    ArrayAllocators.zeros(T=Float64, dims::Dims)

Return an Array with element type `T` with size `dims` filled by `0`s via
`calloc`. Depending on the libc implementation, the operating system may lazily
wait for a page fault before obtaining memory initialized to `0`.
This is an alternative to `Base.zeros` that always fills the array with `0`
eagerly.

# Examples
```julia
julia> @time ArrayAllocators.zeros(Int, 3200, 3200);
  0.000026 seconds (4 allocations: 78.125 MiB)

julia> @time Base.zeros(Int, 3200, 3200);
  0.133595 seconds (2 allocations: 78.125 MiB, 67.36% gc time)

julia> ArrayAllocators.zeros(Int, 256, 256) == Base.zeros(Int, 256, 256)
true
```
"""
function zeros(::Type{T}, dims::Integer...) where T
    return Array{T}(calloc, dims...)
end
function zeros(::Type{T}, dims::Dims) where T
    return Array{T}(calloc, dims)
end
zeros(dims::Integer...) = zeros(Float64, dims)
zeros(dims::Dims) = zeros(Float64, dims)


text:
using CSV
using DataFrames
using Plots
using StatsBase

print("hello World")

function load_df(file_path::AbstractString)
    return DataFrame(CSV.File(file_path))
end

df = load_df("data/preprocessed/yaw-1.csv")

print(df)