text:
# We force specialisation on all arguments to `scan_emit`, otherwise performance can drop
# off to an unacceptable degree when lots of compiler-intensive things like StaticArrays are
# used.
# See https://docs.julialang.org/en/v1/manual/performance-tips/#Be-aware-of-when-Julia-avoids-specializing
# for more info.

"""
    Like Transducers.ScanEmit, but `f` isn't allowed to have internal state, and slightly
    faster in some cases that I care about for this package.

    This function makes some strong assumptions about the nature of `f` and `xs` so that it
    can achieve good performance. In particular, `f` must not involve any globals, and
    the output of `f` should not change type for different elements of `x`.
"""
function scan_emit(f, xs, state, idx)

    # Heuristic Warning: assume all ys have the same type as the 1st.
    (y, state) = f(state, _getindex(xs, idx[1]))
    ys = Vector{typeof(y)}(undef, length(idx))
    ys[idx[1]] = y

    for t in idx[2:end]
        (y, state) = f(state, _getindex(xs, t))
        ys[t] = y
    end

    return (ys, state)
end

function Zygote._pullback(::AContext, ::typeof(scan_emit), f, xs, init_state, idx)

    state = init_state
    (y, state) = f(state, _getindex(xs, idx[1]))

    # Heuristic Warning: assume all ys and states have the same type as the 1st.
    ys = Vector{typeof(y)}(undef, length(xs))
    states = Vector{typeof(state)}(undef, length(xs))

    ys[idx[1]] = y
    states[idx[1]] = state

    for t in idx[2:end]
        (y, state) = f(state, _getindex(xs, t))
        ys[t] = y
        states[t] = state
    end

    function scan_emit_pullback(Δ)

        Δ === nothing && return nothing
        Δys = Δ[1]
        Δstate = Δ[2]

        # This is a hack to handle the case that Δstate=nothing, and the "look at the
        # type of the first thing" heuristic breaks down.
        Δstate = Δ[2] === nothing ? _get_zero_adjoint(states[idx[end]]) : Δ[2]

        T = length(idx)
        if T > 1
            _, Δstate, Δx = step_pb(
                f, states[idx[T-1]], _getindex(xs, idx[T]), Δys[idx[T]], Δstate,
            )
            Δxs = get_adjoint_storage(xs, idx[T], Δx)

            for t in reverse(2:(T - 1))
                a = _getindex(xs, idx[t])
                b = Δys[idx[t]]
                c = states[idx[t-1]]
                _, Δstate, Δx = step_pb(
                    f, c, a, b, Δstate,
                )
                Δxs = _accum_at(Δxs, idx[t], Δx)
            end

            _, Δstate, Δx = step_pb(
                f, init_state, _getindex(xs, idx[1]), Δys[idx[1]], Δstate,
            )
            Δxs = _accum_at(Δxs, idx[1], Δx)

            return (nothing, nothing, Δxs, Δstate, nothing)
        else
            _, Δstate, Δx = step_pb(
                f, init_state, _getindex(xs, idx[1]), Δys[idx[1]], Δstate,
            )
            Δxs = get_adjoint_storage(xs, idx[1], Δx)

            return (nothing, nothing, Δxs, Δstate, nothing)
        end
    end

    return (ys, state), scan_emit_pullback
end

@inline function step_pb(f::Tf, state, x, Δy, Δstate) where {Tf}
    _, pb = _pullback(NoContext(), f, state, x)
    return pb((Δy, Δstate))
end

# Helper functionality for constructing appropriate differentials.

_getindex(x, idx::Int) = getindex(x, idx)

# This really ought to infer, but it doesn't for some unknown reason.
# _getindex(x::Base.Iterators.Zip, idx::Int) = map(x -> _getindex(x, idx), x.is)

# This is a work around for `map` not inferring properly.
_getindex(x::Base.Iterators.Zip, idx::Int) = __getindex(x.is, idx)
__getindex(x::Tuple{Any}, idx::Int) = (_getindex(x[1], idx), )
__getindex(x::Tuple, idx::Int) = (_getindex(x[1], idx), __getindex(Base.tail(x), idx)...)

_get_zero_adjoint(::Any) = nothing
_get_zero_adjoint(x::AbstractArray) = zero(x)

# Vector. In all probability, only one of these methods is necessary.

function get_adjoint_storage(x::Array, n::Int, Δx::T) where {T}
    x̄ = Array{T}(undef, size(x))
    x̄[n] = Δx
    return x̄
end

# function get_adjoint_storage(x::Vector{T}, n::Int, Δx::T) where {T<:Real}
#     x̄ = Vector{T}(undef, length(x))
#     x̄[n] = Δx
#     return x̄
# end

# function get_adjoint_storage(x::Vector, n::Int, init::T) where {T<:AbstractVecOrMat{<:Real}}
#     Δx = Vector{T}(undef, length(x))
#     Δx[n] = init
#     return Δx
# end

# function get_adjoint_storage(x::Vector, n::Int, init::T) where {T<:NamedTuple{(:diag,)}}
#     Δx = Vector{T}(undef, length(x))
#     Δx[n] = init
#     return Δx
# end

# Diagonal type constraint for the compiler's benefit.
@inline function _accum_at(Δxs::Vector{T}, n::Int, Δx::T) where {T}
    Δxs[n] = Δx
    return Δxs
end

# If there's nothing, there's nothing to do.

_accum_at(Δxs::Nothing, n::Int, Δx::Nothing) = nothing

# Zip

function get_adjoint_storage(x::Base.Iterators.Zip, n::Int, Δx::Tuple)
    return (is=map((x_, Δx_) -> get_adjoint_storage(x_, n, Δx_), x.is, Δx),)
end

# function _accum_at(Δxs::NamedTuple{(:is,)}, n::Int, Δx::Tuple)
#     return (is=map((Δxs_, Δx_) -> _accum_at(Δxs_, n, Δx_), Δxs.is, Δx), )
# end

# function _accum_at(Δxs::NamedTuple{(:is,)}, n::Int, Δx::Tuple{Any, Any})
#     return (is=(_accum_at(Δxs[1], n, Δx[1]), _accum_at(Δxs[2], n, Δx[2])), )
#     # return (is=map((Δxs_, Δx_) -> _accum_at(Δxs_, n, Δx_), Δxs.is, Δx), )
# end

# This is a work-around for `map` not inferring for some unknown reason. Very odd...
function _accum_at(Δxs::NamedTuple{(:is, )}, n::Int, Δx::Tuple)
    return (is=__accum_at(Δxs.is, n, Δx), )
end
__accum_at(Δxs::Tuple{Any}, n::Int, Δx::Tuple{Any}) = (_accum_at(Δxs[1], n, Δx[1]), )
function __accum_at(Δxs::Tuple, n::Int, Δx::Tuple)
    return (_accum_at(Δxs[1], n, Δx[1]), __accum_at(Base.tail(Δxs), n, Base.tail(Δx))...)
end

# Fill

get_adjoint_storage(x::Fill, ::Int, init) = (value=init, axes=nothing)

@inline function _accum_at(
    Δxs::NamedTuple{(:value, :axes), Tuple{T, Nothing}}, ::Int, Δx::T,
) where {T}
    return (value=accum(Δxs.value, Δx), axes=nothing)
end

# StructArray

function get_adjoint_storage(x::StructArray, n::Int, Δx::NamedTuple)
    init_arrays = map(
        (x_, Δx_) -> get_adjoint_storage(x_, n, Δx_), getfield(x, :components), Δx,
    )
    return (components = init_arrays, )
end

function _accum_at(Δxs::NamedTuple{(:components,)}, n::Int, Δx::NamedTuple)
    return (components = map((Δy, y) -> _accum_at(Δy, n, y), Δxs.components, Δx), )
end


text:
module GA

    using Statistics, Random

    export
        SimpleGA,
        optimize!

    abstract type EvolutionaryAlgorithm end
    abstract type GeneticAlgorithm <: EvolutionaryAlgorithm end

    mutable struct SimpleGA <: GeneticAlgorithm
        n_iter::Integer
        fitness
        terminal_c::AbstractFloat
        crossover_p::AbstractFloat
        mutation_p::AbstractFloat
        uniform_p::AbstractFloat

        population::Vector{Vector{Bool}}
        best_chro::Vector{Bool}
        iter_n::Integer

        function SimpleGA(;
                n_pop::Integer=100,
                n_genes::Integer=10,
                n_iter::Integer=500,
                fitness,
                terminal_c::AbstractFloat=0.01,
                crossover_p::AbstractFloat=0.5,
                mutation_p::AbstractFloat=0.1,
                uniform_p::AbstractFloat=0.5)
            new(n_iter, fitness, terminal_c, crossover_p, mutation_p, uniform_p,
                [rand(Bool, n_genes) for i ∈ 1:n_pop], zeros(Bool, n_genes), 1)
        end
    end

    function optimize!(ga::GeneticAlgorithm; stop_bybest::Bool=false, stop_bystd::Bool=false)
        for i ∈ ga.iter_n:ga.n_iter
            ga.iter_n = i
            old_pop = shuffle(copy(ga.population))
            new_pop = []
            while !isempty(old_pop)
                ind1 = pop!(old_pop)
                ind2 = pop!(old_pop)
                n_ind1 = copy(ind1)
                n_ind2 = copy(ind2)
                # crossover
                if rand() < ga.crossover_p
                    if rand() < ga.uniform_p
                        for j ∈ 1:length(ind1)
                            if rand() < 0.5
                                n_ind1[j] = ind2[j]
                                n_ind2[j] = ind1[j]
                            end
                        end
                    else
                        c_point = rand(1:length(ind1))
                        n_ind1[c_point:end] = ind2[c_point:end]
                        n_ind2[c_point:end] = ind1[c_point:end]
                    end
                end
                # mutation
                if rand() < ga.mutation_p
                    m_point = rand(1:length(ind1))
                    n_ind1[m_point] = !n_ind1[m_point]
                    n_ind2[m_point] = !n_ind2[m_point]
                end
                push!(new_pop, n_ind1)
                push!(new_pop, n_ind2)
            end
            ga.population = sort([ga.population..., new_pop...], by=ga.fitness)[1:length(ga.population)]
            n_best_chro = ga.population[argmin(broadcast(ga.fitness, ga.population))]
            if ga.fitness(n_best_chro) < ga.fitness(ga.best_chro)
                ga.best_chro = n_best_chro
            elseif stop_bybest
                break
            end
            if stop_bystd && std(ga.fitness.(ga.population)) < ga.terminal_c
                break
            end
        end
    end
end

ext: jl


text:
using Shapefile, Plots, Colors, ColorSchemes, Pkg
using Plots: px
Plots.theme(:default)

# The directory for this project
basedir = dirname(Pkg.project().path)
mkpath(joinpath(basedir, "output"))

# Inlude models and data scripts
include(joinpath(basedir, "src", "rules.jl"))
include(joinpath(basedir, "src", "data.jl"))

# Define rules and simulation data
growth, allee, localdisp, wind, parasitism, localdisp_p, allee_p, carrycap = define_rules()
init_h, init_p, rate_h, rate_p, mask = load_data(basedir, carrycap)

# Define NamedTuples of initialisation and auxilary data for the simulations.
initdata = (; H=init_h, P=init_p)
auxdata = (; rH=rate_h, rP=rate_p)

# Save gifs, as in the paper but with colors
output = GifOutput(initdata;
    filename=joinpath(basedir, "output", "host_parasitoid.gif"),
    tspan=DateTime(2020, 1):Day(7):DateTime(2030, 1),
    aux=(rH=rate_h, rP=rate_p,), mask=mask, # See S3
    minval=(0.0, 0.0), maxval=(carrycap, carrycap),
    scheme=(ColorSchemes.inferno, ColorSchemes.diverging_linear_bjr_30_55_c53_n256),
    text=nothing,
)
sim!(output, (wind, localdisp, allee, growth, localdisp_p, allee_p, parasitism))

# Plot setup
shp = Shapefile.Handle(joinpath(basedir, "data", "ausborder_polyline.shp"))
shape_opts = (lw=0.5, color=RGB(0.3))
plot_opts = (;
    showaxis=false,
    colorbar_formatter=:scientific,
    bottom_margin=0px,
    foreground_color=RGB(0.3),
    tickfontsize=7,
    title="",
    xlims=(112, 155), ylims=(-44, -8),
)

#### Host simulation plots ##########################################
tspan = DateTime(2020, 1):Day(7):DateTime(2030, 1)
kw = (;
    aux=auxdata, mask=mask, tspan=tspan
)
# Define outputs
wind_output = ArrayOutput(initdata; kw...)
localdisp_output = ArrayOutput(initdata; kw...)
combined_output = ArrayOutput(initdata; kw...)
host_para_output = ArrayOutput(initdata; kw...)
# Run simulations
sim!(wind_output, wind, allee, growth)
sim!(localdisp_output, localdisp, allee, growth)
sim!(combined_output, localdisp, wind, allee, growth)
sim!(host_para_output, localdisp, wind, allee, growth, localdisp_p, allee_p, parasitism)

# Default map options
simplot_opts = (;
    color=:inferno,
    legend=:none,
    clims=(0, carrycap),
    showaxis=false,
    ticks=false,
    yguide="",
    plot_opts...
)

# Make simulation plots
t = Ti(Near(DateTime(2022, 2)))
p3 = Plots.plot(localdisp_output[t][:H]; xguide="Local", simplot_opts...)
p4 = Plots.plot(wind_output[t][:H]; xguide="Wind", simplot_opts...)
p5 = Plots.plot(combined_output[t][:H]; xguide="Combined", simplot_opts...)
map(p -> Plots.plot!(p, shp; shape_opts...), (p3, p4, p5))
simplots = Plots.plot(p3, p4, p5; layout=(1, 3), size=(900, 300), dpi=600)

savefig(joinpath(basedir, "output", "simplots.png"))

#### Host/Parasitoid simulation plots #############################

t1 = DateTime(2021, 1)
t2 = DateTime(2025, 1)
t3 = DateTime(2030, 1)
hp_plot_opts = (;
    ytickfontcolor=RGB(1.0),
    simplot_opts...
)
col_p = cgrad(ColorSchemes.diverging_linear_bjr_30_55_c53_n256)
no_x = (xguide="", xtickfontcolor=RGB(1.0))
hp1, hp2, hp3 = map((t1, t2, t3), (no_x, no_x, ())) do t, kw
    hp = Plots.plot(host_para_output[Ti(Near(t))][:H]; xguide="Host", hp_plot_opts..., kw..., yguide=year(t)),
         Plots.plot(host_para_output[Ti(Near(t))][:P]; xguide="Parasitoid", hp_plot_opts..., ytickfontcolor=RGB(1.0), kw..., color=col_p)
    map(x -> Plots.plot!(x, shp; shape_opts...), hp)
end
host_para_plot = Plots.plot(hp1..., hp2..., hp3...; layout=(3, 2), size=(900, 900), dpi=600)

savefig(joinpath(basedir, "output", "hostparaplots.png"))


text:
using SumOfSquares
using DynamicPolynomials
using MosekTools

max_stars_repo_name: sums-of-squares/tools

max_issues_repo_name: sums-of-squares/tools

max_forks_repo_name: sums-of-squares/tools


text:
mutable struct MD5_CTX <: SHA_CTX
    state::Array{UInt32,1}
    bytecount::UInt64
    buffer::Array{UInt8,1}
end

digestlen(::Type{MD5_CTX}) = 16
state_type(::Type{MD5_CTX}) = UInt32
# blocklen is the number of bytes of data processed by the transform!() function at once
blocklen(::Type{MD5_CTX}) = UInt64(64)

MD5_CTX() = MD5_CTX(copy(MD5_initial_hash_value), 0, zeros(UInt8, blocklen(MD5_CTX)))
Base.copy(ctx::T) where {T<:MD5_CTX} = T(copy(ctx.state), ctx.bytecount, copy(ctx.buffer))
Base.show(io::IO, ::MD5_CTX) = write(io, "MD5 hash state")


text:
include("custom_json.jl")

include("widget.jl")
include("observe.jl")

include("gui.jl")
include("gui/index.jl")

include("trees.jl")
include("demo.jl")


text:
### create a directory for saving the results of the city

function createresultsdirectory(path_to_results_folder, city)

    path_to_results = string(path_to_results_folder,city);

    folder_list = readdir(path_to_results_folder);
    for i = 1:length(folder_list)
        folder_name = folder_list[i];
        if (folder_name == city) # catches the case that a results folder already exists for this city
            return path_to_results
        end
    end

    mkdir(path_to_results);
    return path_to_results

end

max_stars_repo_name: ArsamAryandoust/CarParkingMaps

max_issues_repo_name: ArsamAryandoust/CarParkingMaps

max_forks_repo_name: ArsamAryandoust/CarParkingMaps


text:
#
#
# Adapted to from matlab code from BEM
using DataFrames
using Distributions
using Gadfly

df1 = DataFrame(readdlm("tobit_data.txt"))
names!(df1,[:weeks,:D,:D2,:AFQT,:Spouse,:Kids,:Ed])

#nobs = length(weeks);
#const = ones(nobs,1);
#y = weeks;
#x = [const tobit_data(:,4:7)];
#k = size(x,2);
#set prior values

numObs = size(df1,1)
k = 5 # number of prediction variables + the constant
mu_beta = zeros(k,1);
var_beta = 10^2*eye(k);
invarbeta = inv(var_beta);

a = 3;
b = 1/(2*20);
#initialize vectors and begin the Gibbs sampler
iter  = 5500;
burn = 500;
beta_final = zeros(iter-burn,k);
sig_final = zeros(iter-burn,1);

temp1 = array(df1)
D1 = temp1[:,2]
D2 = temp1[:,3]
y = weeks = temp1[:,1]
x = ([ones(size(temp1,1)) temp1[:,4:7]])
bhat = (x' * x)\ (x' * temp1[:,1])

betaParam = bhat;
sig = 500;
for ii in 1:iter
    #augmented latent data
    temp2 = x*betaParam
    ztemp_aug1 = [rand(Truncated(Normal(jj,sig),-1000,0)) for jj in temp2]
	#for those with y=0 (D=1), draw latent data to be truncated from above at zero
    #for now, do this for everyone.

    ztemp_aug2 = [rand(Truncated(Normal(jj,sig),52,1000)) for jj in temp2]  #for those with y=52 (D2=1), draw latent data to be truncated from below at 52
                                                                            #for now, do this for everyone.

    zaug = D1.*ztemp_aug1 + D2.*ztemp_aug2 + (1-D1- D2).*y; #Note that D1,D2 are exclusive!

    #regression paramters
    D_beta = inv(x'*x/sig + invarbeta);
    d_beta = x'*zaug/sig + invarbeta*mu_beta;
    H = chol(D_beta);

    betaParam = D_beta*d_beta + H'*randn(k,1);

    #variance parameter

    sig = rand(InverseGamma( (numObs/2) + a, inv(b) + .5*sum( (zaug - x*betaParam).^2 ) ));

    if ii > burn
        beta_final[ii-burn,:] = betaParam';
        sig_final[ii-burn,:] = sig;
    end
end

weeksX,weeksDist  = hist(beta_final[:,1],50)
AFQTX,AFQTDist  = hist(beta_final[:,2],50)
SpouseX,SpouseDist  = hist(beta_final[:,3],50)
KidsX,KidsDist  = hist(beta_final[:,4],50)
EdX,EdDist  = hist(beta_final[:,5],50)

Gadfly.plot(x = weeksX,y = weeksDist,x = AFQTX,y = AFQTDist)

beta_final = convert(DataFrame,beta_final)
names!(beta_final,[:weeks,:AFQT,:Spouse,:Kids,:Ed])
Gadfly.plot( x = beta_final[:weeks], Geom.histogram)

print("Means and Standard Deviations of parameters");
print("const afqt spouse_inc kids ed, STD_DEV");
print([mean(beta_final) mean(sqrt(sig_final))])
print([std(beta_final) std(sqrt(sig_final))])

print('Marginal Effects');
print('const afqt spouse_inc kids ed, variance');
mu_x = mean(x);
probs = normcdf ( (52 - mu_x*beta_final')'./(sqrt(sig_final))) - normcdf((-mu_x*beta_final')'./(sqrt(sig_final)));
[mean(beta_final(:,2).*probs) mean(beta_final(:,3).*probs) mean(beta_final(:,4).*probs) mean(beta_final(:,5).*probs) ]
[std(beta_final(:,2).*probs) std(beta_final(:,3).*probs) std(beta_final(:,4).*probs) std(beta_final(:,5).*probs) ]

max_stars_repo_name: Darf-Ferrara/bayes-scripts

max_issues_repo_name: Darf-Ferrara/bayes-scripts

max_forks_repo_name: Darf-Ferrara/bayes-scripts


text:
### Methods ###

function ThalamusModel(config::Dict)
  numTR = get(config, "numTR", 8)
  numtRN = get(config, "numtRN", 1)
  szTR = get(config, "szTR", 64)
  sztRN = get(config, "sztRN", 64)

  agent = Agent()

  agent_TR = Agent()
  for i = 1:numTR
    addnode!(agent_TR, GenericNeurons(szTR; b=0.25, d=0.05); name="$i")
  end
  agent_tRN = Agent()
  for i = 1:numtRN
    addnode!(agent_tRN, GenericNeurons(sztRN; a=0.1, d=2); name="$i")
  end
  agent["TR"] = agent_TR
  agent["tRN"] = agent_tRN

  # Reciprocally connect TR and tRN
  for i = 1:numTR
    for j = 1:numtRN
      # TODO: Delays
      # TODO: Masking
      addnode!(agent, GenericSynapses(szTR, sztRN; outputMask=-1); name="TR_$(i)_tRN_$(j)")
      addedge!(agent, agent["TR_$(i)_tRN_$(j)"], (
        (agent["TR"]["$i"], :input),
        (agent["tRN"]["$j"], :output)
      ))
      addnode!(agent, GenericSynapses(sztRN, szTR); name="tRN_$(j)_TR_$(i)")
      addedge!(agent, agent["tRN_$(j)_TR_$(i)"], (
        (agent["tRN"]["$j"], :input),
        (agent["TR"]["$i"], :output)
      ))
    end
  end

  return agent
end


text:
info("Running packet tests")

function on_msg(topic, payload)
    @test topic == "abc"
    @test String(payload)== "qwerty"
end

function on_disconnect(reason)
    @test reason == nothing
end

function on_disconnect_ping(reason)
    @test reason.msg == "ping timed out"
end

function on_disconnect_pingresp(reason)
    @test reason.msg == "protocol error"
end

function is_out_correct(filename_expected::AbstractString, actual::Channel{UInt8}, mid::UInt16)
    file_data = read_all_to_arr(filename_expected)
    actual_data = Vector{UInt8}()

    for i in file_data
      append!(actual_data, take!(actual))
    end

    mid_index = get_mid_index(file_data)
    if mid_index > 0
      buffer = PipeBuffer()
      write(buffer, mid)
      converted_mid = take!(buffer)
      file_data[mid_index] = converted_mid[2]
      file_data[mid_index+1] = converted_mid[1]
    end

    correct = true
    i = 1
    while i <= length(file_data)
      if file_data[i] != actual_data[i]
        correct = false
        break
      end
      i += 1
    end
    return correct
end

function is_out_correct(filename_expected::AbstractString, actual::Channel{UInt8})
    file = open(filename_expected, "r")
    correct = true
    while !eof(file)
        if read(file, UInt8) != take!(actual)
            correct = false
            break
        end
    end
    return correct
end

next_id(c) = c.last_id + 0x0001

function test()
    client = Client(on_msg, on_disconnect, 60)
    opts = ConnectOpts(() -> TestFileHandler())
    opts.client_id = "TestID"

    last_id::UInt16 = 0x0001

    #info("Testing connect")
    connect(client, opts)
    tfh::TestFileHandler = client.io
    @test is_out_correct("data/output/connect.dat", tfh.out_channel)
    # CONNACK is automatically being sent in connect call

    #info("Testing subscribe")
    subscribe(client, ("abc", AT_LEAST_ONCE), ("cba", AT_MOST_ONCE), async=true)
    @test is_out_correct("data/output/subreq.dat", tfh.out_channel, next_id(client))
    put_from_file(tfh, "data/input/suback.dat", client.last_id)

    #info("Testing unsubscribe")
    unsubscribe(client, "abc", "cba", async=true)
    @test is_out_correct("data/output/unsubreq.dat", tfh.out_channel, next_id(client))
    put_from_file(tfh, "data/input/unsuback.dat", client.last_id)

    info("Testing receive publish QOS 0")
    put_from_file(tfh, "data/input/qos0pub.dat")

    info("Testing receive publish QOS 1")
    put_from_file(tfh, "data/input/qos1pub.dat", last_id)
    @test is_out_correct("data/output/puback.dat", tfh.out_channel, last_id)
    #last_id += 1

    info("Testing receive publish QOS 2")
    put_from_file(tfh, "data/input/qos2pub.dat", last_id)
    @test is_out_correct("data/output/pubrec.dat", tfh.out_channel, last_id)
    put_from_file(tfh, "data/input/pubrel.dat", last_id)
    @test is_out_correct("data/output/pubcomp.dat", tfh.out_channel, last_id)
    #last_id += 1

    info("Testing send publish QOS 0")
    publish(client, "test1", "QOS_0", qos=AT_MOST_ONCE, async=true)
    @test is_out_correct("data/output/qos0pub.dat", tfh.out_channel)

    info("Testing send publish QOS 1")
    publish(client, "test2", "QOS_1", qos=AT_LEAST_ONCE, async=true)
    @test is_out_correct("data/output/qos1pub.dat", tfh.out_channel, next_id(client))
    put_from_file(tfh, "data/input/puback.dat", client.last_id)

    info("Testing send publish QOS 2")
    f = publish(client, "test3", "test", qos=EXACTLY_ONCE, async=true)
    id = next_id(client)
    @test is_out_correct("data/output/qos2pub.dat", tfh.out_channel, id)
    put_from_file(tfh, "data/input/pubrec.dat", id)
    @test is_out_correct("data/output/pubrel.dat", tfh.out_channel, id)
    put_from_file(tfh, "data/input/pubcomp.dat", id)

    get(f)

    info("Testing disconnect")
    disconnect(client)
    @test is_out_correct("data/output/disco.dat", tfh.out_channel)

    #This has to be in it's own connect flow to not interfere with other messages
    info("Testing keep alive with response")
    client = Client(on_msg, on_disconnect_ping, 1)
    opts = ConnectOpts(() -> TestFileHandler())
    opts.client_id = "TestID"
    opts.keep_alive = 0x0001
    connect(client, opts)
    tfh = client.io

    @test is_out_correct("data/output/connect_keep_alive1s.dat", tfh.out_channel) # Consume output
    @test is_out_correct("data/output/pingreq.dat", tfh.out_channel)
    put_from_file(tfh, "data/input/pingresp.dat")

    info("Testing keep alive without response")
    sleep(1.1)
    @test is_out_correct("data/output/pingreq.dat", tfh.out_channel)

    info("Testing unwanted pingresp")
    client = Client(on_msg, on_disconnect_pingresp, 60)
    opts = ConnectOpts(() -> TestFileHandler())
    opts.client_id = "TestID"
    opts.keep_alive = 0x0001
    connect(client, opts)
    tfh = client.io

    put_from_file(tfh, "data/input/pingresp.dat")
    sleep(2)
end

test()

max_issues_repo_name: kivaari/jlMQTT

max_forks_repo_name: kivaari/jlMQTT


text:
matrixcompletion_insts = [
    [(false, false, false, d, k * d, 0.8) for d in vcat(2, 5:5:max_d)] # includes compile run
    for (k, max_d) in ((10, 45), (20, 30))
    ]

insts = OrderedDict()
insts["nat"] = (nothing, matrixcompletion_insts)
insts["extEP"] = (:ExpPSD, matrixcompletion_insts)
insts["extSEP"] = (:SOCExpPSD, matrixcompletion_insts)
return (MatrixCompletionJuMP, insts)


text:
struct ObservablePair{S, T} <: AbstractObservable{T}
    first::AbstractObservable{S}
    second::AbstractObservable{T}
    f
    g
    excluded::Vector{Function}
    function ObservablePair(first::AbstractObservable{S}, second::AbstractObservable{T}; f = identity, g = identity) where {S, T}
        excluded = Function[]
        first2second = on(first) do val
            setindex!(second, f(val), notify = !in(excluded))
        end
        push!(excluded, first2second)
        second2first = on(second) do val
            setindex!(first, g(val), notify = !in(excluded))
        end
        push!(excluded, second2first)
        new{S, T}(first, second, f, g, excluded)
    end
end

ObservablePair(first::AbstractObservable; f = identity, g = identity) =
    ObservablePair(first, Observable{Any}(f(first[])); f = f, g = g)

observe(o::ObservablePair) = o.second

function off(o::ObservablePair)
    for i in 1:2
        off(o[i], o.excluded[i])
    end
end

Base.iterate(p::ObservablePair, i=1) = i > 2 ? nothing : (getfield(p, i), i + 1)
Base.getindex(p::ObservablePair, i::Int) = getfield(p, i)
Base.firstindex(p::ObservablePair) = 1
Base.lastindex(p::ObservablePair) = 2
Base.length(p::ObservablePair) = 2
Base.first(p::ObservablePair) = p.first
Base.last(p::ObservablePair) = p.second


text:
# Age ranges
ages = ["$age0-$(age0+4)" for age0 = 0:5:95]
push!(ages, "100+")
female = [328, 319, 307, 290, 284, 293, 293, 263, 241, 235,
    218, 191, 161, 135, 95, 66, 47, 26, 11.3, 3.1, 0.4]
male = [350, 341, 325, 310, 303, 310, 305, 270, 246, 239,
    218, 187, 153, 124, 83, 53, 34, 16, 5.7, 1.2, 0.1]
# Plot as bars - with negative values for female
barplot(ages, -female, horizontal=true, color=0xff8080)
hold(true)
barplot(ages, male, horizontal=true, color=0x8080ff)
# Change to absolute figures
xticklabels(x -> abs(x))
# Add labels
title("World population by age in 2019 (https://www.populationpyramid.net/)")
xlabel("Million people")
legend("female", "male")


text:
using SimpleDirectMediaLayer
const SDL2 = SimpleDirectMediaLayer

function pollEvent!()
    #SDL2.Event() = [SDL2.Event(NTuple{56, Uint8}(zeros(56,1)))]
    SDL_Event() = Array{UInt8}(zeros(56))
    e = SDL_Event()
    success = (SDL2.PollEvent(e) != 0)
    return e,success
end

function getEventType(e::Array{UInt8})
    # HAHA This is still pretty janky, but I guess that's all you can do w/ unions.
    bitcat(UInt32, e[4:-1:1])
end

function getEventType(e::SDL2.Event)
    e._Event[1]
end

function bitcat(outType::Type{T}, arr)::T where T<:Number
    out = zero(outType)
    for x in arr
        out = out << (sizeof(x)*8)
        out |= convert(T, x)  # the `convert` prevents signed T from promoting to Int64.
    end
    out
end

function drawinit()
    SDL2.init()

    win = SDL2.CreateWindow("Hello World!", Int32(100), Int32(100), Int32(800), Int32(600),
        UInt32(SDL2.WINDOW_SHOWN))
    SDL2.SetWindowResizable(win,true)

    renderer = SDL2.CreateRenderer(win, Int32(-1),
        UInt32(SDL2.RENDERER_ACCELERATED | SDL2.RENDERER_PRESENTVSYNC))

    return renderer
end

function drawsetcolor(renderer,r,g,b,a=255)
    SDL2.SetRenderDrawColor(renderer,r,g,b,a)
end

function drawclear(renderer)
    SDL2.RenderClear(renderer)
end

function drawline(renderer,x1,y1,x2,y2)
    SDL2.RenderDrawLine(renderer,trunc(Int,x1),trunc(Int,y1),trunc(Int,x2),trunc(Int,y2))
end

function drawloop(renderer, updatefunction)
    running = true

    while running
        SDL2.PumpEvents()

        updatefunction()

        SDL2.RenderPresent(renderer)

        e,_ = pollEvent!()
        t = getEventType(e)
        if t == SDL2.QUIT
            running = false
        end
    end

    SDL2.Quit()

end


text:
using OmegaGrav
using CSV, DataFrames
using Dierckx
using Plots, Plots.PlotMeasures, LaTeXStrings
# %% Compute the matter power spectrum using CLASS
include("compute_pk_class.jl")
# %% Compute `Ωgrav = Ωm * W / 2` and Ωtherm at redshifts of the data points
d = CSV.read("data/d16_Omega_th_data.csv", DataFrame)
nred = size(d)[1]
redshift = zeros(nred + 1)
Ωghalo = zeros(nred + 1)
Ωgpk = zeros(nred + 1)
Ωglin = zeros(nred + 1)
Ωth = zeros(nred + 1)
Ωthu = zeros(nred + 1)
Ωthl = zeros(nred + 1)
for ired = 1:nred+1
   if ired == 1
      z = 0 # Add z=0 which is not included in the data file
   else
      z = d.z[ired-1]
   end
   redshift[ired] = z
   # Define functions to return linear and non-linear power spectra
   # Note: The CLASS code takes wavenumbers in units of 1/Mpc (no h) and
   # return power spectra in units of Mpc^3 (no 1/h^3).
   pkcb_class(kovh) = cosmo.pk_cb_lin(kovh * h0, z) * h0^3
   pknl_class(kovh) = cosmo.pk(kovh * h0, z) * h0^3
   pklin_class(kovh) = cosmo.pk_lin(kovh * h0, z) * h0^3
   # Spline interpolate in log(k)
   lnk = log(1e-4):0.05:log(100)
   pkcb = Spline1D(lnk, pkcb_class.(exp.(lnk)))
   pknl = Spline1D(lnk, pknl_class.(exp.(lnk)))
   pklin = Spline1D(lnk, pklin_class.(exp.(lnk)))
   # %% Compute Ωgrav from non-linear total matter P(k)
   Ωgpk[ired] = ograv_pk(x -> pknl(log(x)), z, Ωm)
   # %% Compute Ωgrav from linear total matter P(k)
   Ωglin[ired] = ograv_pk(x -> pklin(log(x)), z, Ωm)
   # %% Compute Ωgrav from Halos, excluding the neutrino contribution
   Ωghalo[ired] = ograv_halo(x -> pkcb(log(x)), z, Ωm, Ωcb)
   # %% Compute Ωtherm from Halos, excluding the neutrino contribution
   Ωth[ired] = otherm_upp(x -> pkcb(log(x)), z, Ωm, h0, Ωcb)
   # %% Compute Ωtherm for upper and lower 68% confidence level in B
   Ωthl[ired] = otherm_upp(x -> pkcb(log(x)), z, Ωm, h0, Ωcb, massbias = 1.315)
   Ωthu[ired] = otherm_upp(x -> pkcb(log(x)), z, Ωm, h0, Ωcb, massbias = 1.221)
end

#%% Plot results and save to figure1.pdf
fb = params["omega_b"] / (params["omega_cdm"] + params["omega_b"])
Ωb = params["omega_b"] / h0^2
ΔΩg = Ωgpk - Ωglin
y = 0.5 * (Ωghalo + ΔΩg) * (-2fb / 3)
Δyu = ΔΩg * (-2fb / 3) - y
Δyl = y - Ωghalo * (-2fb / 3)
p = plot(
   redshift,
   -2 * Ωgpk,
   c = :black,
   lab = L"-\Omega_W",
   lw = 2,
   xlims = [0, 1.5],
   ylims = [1e-9, 2e-6],
   yaxis = :log,
   xlab = "Redshift, z",
   ylab = "Comoving Energy Density Parameters",
   legend = :topright,
   legendfontsize = 12,
   labelfontsize = 14,
   right_margin = 45px,
)
p = plot!(
   redshift,
   -2 * Ωghalo,
   c = :black,
   ls = :dash,
   lab = L"-\Omega_W^{halo}",
   lw = 2,
)
p = plot!(redshift, y, ribbon = (Δyl, Δyu), lab = L"\Omega_W^{ref}", c = :green)
p = scatter!(
   (d.z, d.Omega_th),
   yerror = (d.Omega_th .- d.Omega_th_low, d.Omega_th_up .- d.Omega_th),
   ms = 5,
   c = 1,
   lab = L"\Omega_{th}~Data",
)
p = plot!(
   redshift,
   Ωth,
   c = 1,
   ribbon = (Ωth - Ωthl, Ωthu - Ωth),
   lab = "Halo Model",
)
p = plot!(
   twinx(),
   redshift,
   0.2 * Ωth / 1.7755e-8 * 0.049 / Ωb,
   c = 1,
   lw = 2,
   ylab = "Density Weighted Temperature [keV]",
   yaxis = :log,
   legend = false,
   xaxis = false,
   ylims = [1e-9, 2e-6] * 0.2 / 1.7755e-8 * 0.049 / Ωb,
   xlims = [0, 1.5],
   yticks = ([1e-2, 1e-1, 1, 10], ["0.01", "0.1", "1", "10"]),
)
savefig("figure1.pdf")
display(p)

# %% Clean CLASS (the equivalent of the struct_free() in the `main`
# of CLASS. This step is primordial when running in a loop over different
# cosmologies, as you will saturate your memory very fast if you ommit
# it.
cosmo.struct_cleanup()
# If you want to change completely the cosmology, you should also
# clean the arguments, otherwise, if you are simply running on a loop
# of different values for the same parameters, this step is not needed
cosmo.empty()


text:
cd(@__DIR__)

@testset "reading negative delta" begin
    midi = readMIDIFile("negative_delta.mid")
    @test midi.tracks[1].events[3].dT == -48
end

@testset "writing negative delta" begin
    stream = IOBuffer()
    @test_throws ErrorException MIDI.writeevent(stream, MIDIEvent(-1, 0x80, [0x00, 0x00]))
    @test_throws ErrorException MIDI.writeevent(stream, MetaEvent(-1, 0x51, [0x00, 0x00, 0x01]))
    @test_throws ErrorException MIDI.writeevent(stream, SysexEvent(-1, [0x01]))
    @test_throws ErrorException MIDI.writevariablelength(stream, -1)
    @test_throws ErrorException MIDI.writevariablelength(stream, Int(0x1FFFFFFF))
    close(stream)
end


text:
module Jib

using DataFrames,
      Sockets,
      TimeZones

include("client.jl")
include("enums.jl")
include("types.jl")
include("types_condition.jl")
include("types_mutable.jl")
include("types_private.jl")
include("wrapper.jl")
include("reader.jl")          ; using .Reader: check_all, start_reader
include("utils.jl")

"""
    Connection()

Hold a connection to IB TWS or IBGateway.
"""
struct Connection
  socket::TCPSocket
  id::Int
  connectOptions::String
  version::Client.Version
  time::ZonedDateTime
  servertz::TimeZone
end
Connection(socket, id, connectOptions, version, time) = Connection(socket,
                                                                   id,
                                                                   connectOptions,
                                                                   version,
                                                                   time,
                                                                   timezone(time))

include("requests.jl")        ; using .Requests

"""
    connect([host, ]port, clientId, connectOptions="", optionalCapabilities="")

Connect to host `host` on port `port` and set client ID `clientId`.

Return a [`Connection`](@ref).
"""
function connect(host, port, clientId, connectOptions::String="", optionalCapabilities::String="")

  s = Sockets.connect(host, port)

  # Init string
  m, M = Client.Version .|> (typemin, typemax) .|> Int

  buf = Client.buffer(true)
  print(buf, m == M ? "v$m" : "v$m..$M")
  isempty(connectOptions) || print(buf, ' ', connectOptions)

  # Handshake
  Client.write_one(s, buf)

  res = Reader.read_msg(s)

  @assert length(res) == 2

  @info "connected" V=res[1] T=res[2]

  v = parse(Int, res[1])
  m ≤ v ≤ M || error("unsupported version")

  ib = Connection(s, clientId, connectOptions, Client.Version(v), from_ibtime(res[2]))

  Requests.startApi(ib, clientId, optionalCapabilities)

  ib
end

function connect(port, clientId, connectOptions::String="", optionalCapabilities::String="")

  localip = getalladdrinfo("localhost")

  connect(localip[1], port, clientId, connectOptions, optionalCapabilities)
end

"""
    disconnect(ib)

Close the connection.
"""
disconnect(ib) = close(ib.socket)

end


text:
using SeisIO, JLD2, Distributed, Sockets

export sac2jld2, sac2jld2_par
#

# Converting sac data to JLD2 format. There are options to save to JLD2 format
# as single file or separate files for each sac file.
"""
Key features:
1. Can be used to convert individual SAC file or a group/list of SAC files.
2. Handles the situation when multiple data files were download for the same station
    and channel. Those files will be merged into one single SeisData object with
    multiple SeisChannel entries.
3. All SAC headers are kept and saved in :misc filed of the SeisData object.

"""

"""

    sac2jld2(infile)

Saves individual sac file (in) to JLD2 file.
The output file name will be the same as the input, with new extension *.jld2.

"""

function sac2jld2(infile::String)
    print("Same file, new extension of jld2.\n")
    print("Converting SAC file: [ ",infile," ] ... \n")
    outfile = join([infile[1:end-3],"jld2"]); #trim the file name to remove .sac
    sacin = SeisIO.read_data("sac",infile,full=true);
    sacin.misc[1]["dlerror"] = 0

    stationlist = String[];
    push!(stationlist,sacin.id[1])  #pass the 1st element of sacin.id[] to the station list
    # stationinfo = Dict(["stationlist" => stationlist,"stationmethod" => "sac2jld2","stationsrc" => "converted"])
    stimetemp = sacin.t[1][3]*1e-6; #convert mseconds to seconds;
    saclength = sacin.t[1][2]/sacin.fs; #convert from number of samples to time duration.
    starttime = u2d(stimetemp[1]); #convert union time to calendar time;
    endtime = u2d(stimetemp[1] + saclength[1]);
    DLtimestamplist = String[];
    push!(DLtimestamplist,replace(string(starttime),"-" => "."));

    #save to jld2
    stemp = joinpath(DLtimestamplist[1],stationlist[1]);
    file = jldopen(outfile, "w");
    file["info/DLtimestamplist"] = DLtimestamplist;
    file["info/stationlist"] = stationlist;
    file["info/starttime"]   = string(starttime)
    file["info/endtime"]     = string(endtime)
    file["info/DL_time_unit"]= saclength[1]; #length of the data in time (duration), not number of samples

    #save the waveform data as SeisData format.
    file[stemp] = SeisData(sacin);

    JLD2.close(file);

    print("Saved JLD2 file to: ", outfile, "\n");
end

"""

    sac2jld2(infile,outfile)

Saves individual sac file (in) to JLD2 file (out), which is explicitly specified.

"""

function sac2jld2(infile::String,outfile::String)
    # print("I convert SAC file to JLD2 file.\n")
    print("Converting SAC file: [ ",infile," ] ... \n")
    sacin = SeisIO.read_data("sac",infile,full=true);
    sacin.misc[1]["dlerror"] = 0

    stationlist = String[];
    push!(stationlist,sacin.id[1])  #pass the 1st element of sacin.id[] to the station list
    # stationinfo = Dict(["stationlist" => stationlist,"stationmethod" => "sac2jld2","stationsrc" => "converted"])
    stimetemp = sacin.t[1][3]*1e-6; #convert mseconds to seconds;
    saclength = sacin.t[1][2]/sacin.fs; #convert from number of samples to time duration.
    starttime = u2d(stimetemp[1]); #convert union time to calendar time;
    endtime = u2d(stimetemp[1] + saclength[1]);
    DLtimestamplist = String[];
    push!(DLtimestamplist,replace(string(starttime),"-" => "."));

    #save to jld2
    stemp = joinpath(DLtimestamplist[1],stationlist[1]);
    file = jldopen(outfile, "w");
    file["info/DLtimestamplist"] = DLtimestamplist;
    file["info/stationlist"] = stationlist;
    file["info/starttime"]   = string(starttime)
    file["info/endtime"]     = string(endtime)
    file["info/DL_time_unit"]= saclength[1]; #length of the data in time (duration), not number of samples

    #save the waveform data as SeisData format.
    file[stemp] = SeisData(sacin);

    JLD2.close(file);

    print("Saved JLD2 file to: ", outfile, "\n");
end

"""

    sac2jld2(filelist,timestamp,outfile,verbose=false)

Saves sac files based on a filelist generated using ls() in Julia to
JLD2 file (out), which is explicitly specified.

    Arguments:
    filelist: An String Array including all sacfiles with the same tiemstamp, e.g., same day
    tiemstamp: this must be specified for the sac fiels, assuming they all share
                the same time stamp. Suggested format: 2019.06.27T00.00.00
    outfile: output file name for the JLD2 file. All data fiels will be readin and
                saved to a single JLD2 file.
    verbose: if true, the code will print out more messages. Default is false.
"""

function sac2jld2(filelist::Array{String,1},timestamp::String,outfile::String,verbose::Bool=false)
    # print("I convert and pack all SAC files to a JLD2 file.\n")

    #read the filenames from the infilelist into a strinng array
    # filelist = readfilelist(infilelist)

    stationlist = String[];
    DLtimestamplist = String[];
    push!(DLtimestamplist,timestamp); #use the specified timestamp for the group.

    #get some intial information from the first sac file.
    sacin0 = SeisIO.read_data("sac",filelist[1]);
    stimemin = sacin0.t[1][3]*1e-6; #initial starttime in seconds;
    ftemp = sacin0.t[1][2]/sacin0.fs; #convert from number of samples to time duration.
    saclengthmax = ftemp[1]
    etimemax = stimemin + saclengthmax;

    file = jldopen(outfile, "w");
    file["info/DLtimestamplist"] = DLtimestamplist;
    stemp_pre = "-"
    count = 0
    # stationlisttemp = String[];
    #To deal with multiple channel issue, I am doing the loop twice.
    # 1. The first loop with collect multiple channels into one SeisData object.
    # 2. The second loop will work through the actual number of SeisData objects, counted by `count`
    Sall = SeisData[];
    Str_all = String[];

    #in the loop for each sac file. we will find the minimum starttime, maximum endtime, and maximum saclength.
    for infile = filelist
        println("Converting SAC file: [ ",infile," ]")

        sacin = SeisIO.read_data("sac",infile,full=true);
        sacin.misc[1]["dlerror"] = 0
        push!(stationlist,sacin.id[1])  #pass the 1st element of sacin.id[] to the station list
        # stationinfo = Dict(["stationlist" => stationlist,"stationmethod" => "sac2jld2","stationsrc" => "converted"])
        stimetemp = sacin.t[1][3]*1e-6; #convert mseconds to seconds;
        if stimetemp < stimemin
            stimemin = stimetemp
        end
        saclength = sacin.t[1][2]/sacin.fs; #convert from number of samples to time duration.
        if saclength[1] > saclengthmax
            saclengthmax = saclength[1]
        end
        etimetemp = stimetemp + saclength[1]
        if etimetemp > etimemax
            etimemax = etimetemp
        end

        #save to jld2
        stemp = joinpath(timestamp,sacin.id[1]);
        if stemp == stemp_pre
            println(stemp," exists. Appending to form multiple channels.")
            append!(Sall[count],SeisData(sacin))
        else
            push!(Sall,SeisData(sacin))
            push!(Str_all,stemp)
            count += 1
            # file[stemp] = SeisData(sacin);
        end
        stemp_pre = stemp
    end
    starttime = u2d(stimemin); #convert union time to calendar time;
    endtime = u2d(etimemax);
    #the following info dat is saved after running through the whole group.
    file["info/stationlist"] = unique(stationlist);
    file["info/starttime"]   = string(starttime)
    file["info/endtime"]     = string(endtime)
    file["info/DL_time_unit"]= saclengthmax; #length of the data in time (duration), not number of samples

    #2. start to save the data to jld2
    for i = 1:count
        file[Str_all[i]] = Sall[i]
    end #end of loop for all sac files within one group/directory

    JLD2.close(file);

    print("Saved JLD2 file to: ", outfile, "\n");
end

"""

    sac2jld2(filelist,timestamplist,outfile,verbose=false)

Saves sac files based on a filelist generated using ls() in Julia to
JLD2 file (out), which is explicitly specified. This function assumes each file in
filelist is assigned a timestamp in timestamplist.

    Arguments:
    filelist: An String Array including all sacfiles with the same tiemstamp, e.g., same day
    tiemstamplist: this must be specified for all the sac files in filelist. This list has to be
                the same length as filelist.
    outfile: output file name for the JLD2 file. All data fiels will be readin and
                saved to a single JLD2 file.
    verbose: if true, the code will print out more messages. Default is false.
"""

function sac2jld2(filelist::Array{String,1};timestamplist::Array{String,1},outfile::String,verbose::Bool=false)
    if length(filelist) != length(timestamplist)
        error("filelist and timestamplist must be the same length!")
    end
    stationlist = String[];

    #get some intial information from the first sac file.
    sacin0 = SeisIO.read_data("sac",filelist[1]);
    stimemin = sacin0.t[1][3]*1e-6; #initial starttime in seconds;
    ftemp = sacin0.t[1][2]/sacin0.fs; #convert from number of samples to time duration.
    saclengthmax = ftemp[1]
    etimemax = stimemin + saclengthmax;

    file = jldopen(outfile, "w");
    file["info/DLtimestamplist"] = unique(timestamplist);
    stemp_pre = "-"
    count = 0
    # stationlisttemp = String[];
    #To deal with multiple channel issue, I am doing the loop twice.
    # 1. The first loop with collect multiple channels into one SeisData object.
    # 2. The second loop will work through the actual number of SeisData objects, counted by `count`
    Sall = SeisData[];
    Str_all = String[];

    #in the loop for each sac file. we will find the minimum starttime, maximum endtime, and maximum saclength.
    for (infile,ts) in zip(filelist,timestamplist)
        println("Converting SAC file: [ ",infile," ]")

        sacin = SeisIO.read_data("sac",infile,full=true);
        sacin.misc[1]["dlerror"] = 0
        push!(stationlist,sacin.id[1])  #pass the 1st element of sacin.id[] to the station list
        # stationinfo = Dict(["stationlist" => stationlist,"stationmethod" => "sac2jld2","stationsrc" => "converted"])
        stimetemp = sacin.t[1][3]*1e-6; #convert mseconds to seconds;
        if stimetemp < stimemin
            stimemin = stimetemp
        end
        saclength = sacin.t[1][2]/sacin.fs; #convert from number of samples to time duration.
        if saclength[1] > saclengthmax
            saclengthmax = saclength[1]
        end
        etimetemp = stimetemp + saclength[1]
        if etimetemp > etimemax
            etimemax = etimetemp
        end

        #save to jld2
        stemp = joinpath(ts,sacin.id[1]);
        if stemp == stemp_pre
            println(stemp," exists. Appending to form multiple channels.")
            append!(Sall[count],SeisData(sacin))
        else
            push!(Sall,SeisData(sacin))
            push!(Str_all,stemp)
            count += 1
            # file[stemp] = SeisData(sacin);
        end
        stemp_pre = stemp
    end
    starttime = u2d(stimemin); #convert union time to calendar time;
    endtime = u2d(etimemax);
    #the following info dat is saved after running through the whole group.
    file["info/stationlist"] = unique(stationlist);
    file["info/starttime"]   = string(starttime)
    file["info/endtime"]     = string(endtime)
    file["info/DL_time_unit"]= saclengthmax; #length of the data in time (duration), not number of samples

    #2. start to save the data to jld2
    for i = 1:count
        file[Str_all[i]] = Sall[i]
    end #end of loop for all sac files within one group/directory

    JLD2.close(file);

    print("Saved JLD2 file to: ", outfile, "\n");
end

"""

    sac2jld2(fhandle::JLDfile, sacdir::String,timestamp::String,verbose::Bool=false)

Saves sac files based in sacdir to JLD2 file handle (fhandle), tagged as timestamp.

    Arguments:
    fhandle: A JLD2 file handle, generated using jldopen(), with write permission.
    sacdir: directory containing the sac files need to convert.
    tiemstamp: this must be specified for the sac fiels, assuming they all share
                the same time stamp. Suggested format: 2019.06.27T00.00.00
    verbose: if true, the code will print out more messages. Default is false.
"""

function sac2jld2(fhandle::JLD2.JLDFile, sacdir::String,timestamp::String,verbose::Bool=false)
    # print("I convert and pack all SAC files to a JLD2 file.\n")

    filelist=ls(joinpath(sacdir,"*.sac"));
    println("Converting SAC dir: ",sacdir,"=> ",length(filelist)," files.")

    stemp_pre = "-"
    count = 0
    #To deal with multiple channel issue, I am doing the loop twice.
    # 1. The first loop with collect multiple channels into one SeisData object.
    # 2. The second loop will work through the actual number of SeisData objects, counted by `count`
    Sall = SeisData[];
    Str_all = String[];
    #in the loop for each sac file. we will find the minimum starttime, maximum endtime, and maximum saclength.
    for infile = filelist
        if verbose == true
            print("Converting SAC file: [ ",infile," ] ...\n")
        end

        sacin = SeisIO.read_data("sac",infile,full=true);
        sacin.misc[1]["dlerror"] = 0

        #save to jld2
        stemp = joinpath(timestamp,sacin.id[1]);
        println(stemp)
        if stemp == stemp_pre
            if verbose == true
                println(stemp," exists. Appending to form multiple channels.")
            end
            append!(Sall[count],SeisData(sacin))
        else
            push!(Sall,SeisData(sacin))
            push!(Str_all,stemp)
            count += 1
            # file[stemp] = SeisData(sacin);
        end
        stemp_pre = stemp
    end

    #2. start to save the data to jld2
    for i = 1:count
        fhandle[Str_all[i]] = Sall[i]
    end #end of loop for all sac files within one group/directory
end

"""

    sac2jld2(sacdirlist,timestamplist,outfile,verbose=false)

Searches for all sac files scanning throgh the sacdirlist (Array{String,1}). The sac files
    in each sac directory are assigned to one time stamp based on timestamplist (Array{String,1}).
    The two Arrays must be equal in size. The list can be generated from a file list (using ls in shell),
    using a utility function in this package: readfilelist(infilelist::String).

All sac files will be saved into one single JLD2 file (out), which is explicitly specified as outfile.

    Arguments:
    sacdirlist: An array of strings storing the directories grouped by time stamp, e.g., one day.
    timestamplist: An array of strings storing the time stamps assigning to each SAC file group.
    outfile: output file name for the JLD2 file. All data fiels will be readin and
                saved to a single JLD2 file.
    verbose: if true, the code will print out more messages. Default is false.
"""

function sac2jld2(sacdirlist::Array{String,1},timestamplist::Array{String,1},outfile::String,verbose::Bool=false)
    if length(sacdirlist) != length(timestamplist)
        error("sacdirlist and timestamplist must be the same length!")
    end

    stationlist = String[];

    # saclengthmax = 0; #initiate the length as zero, will find the maximum value after looping through the files.

    file = jldopen(outfile, "w");
    file["info/DLtimestamplist"] = timestamplist;
    file["info/starttime"] = timestamplist[1];
    #here we use the first value in timestamplist as the starttime for a group of multiple timestamp data.
    file["info/endtime"] = timestamplist[end]; #similar to starttime, here we use the last value in timestamplist as the endtime.

    #in the loop for each sac file in each SAC directory. we will find the maximum saclength.
    # countdir=1;
    t1all=0
    t2all=0
    for (sacdir,ts) in zip(sacdirlist,timestamplist)

        t0=time()
        filelist=ls(joinpath(sacdir,"*.sac"));

        print("Converting for directory: [ ",sacdir," ] ... \n")

        count = 0
        # stationlisttemp = String[];
        #To deal with multiple channel issue, I am doing the loop twice.
        # 1. The first loop with collect multiple channels into one SeisData object.
        # 2. The second loop will work through the actual number of SeisData objects, counted by `count`
        Sall = SeisData[];
        Str_all = String[];
        stemp_pre = "-"
        for infile = filelist
            # t1=@elapsed sacin = SeisIO.read_data("sac",infile,full=true);
            sacin = SeisIO.read_data("sac",infile,full=true);
            sacin.misc[1]["dlerror"] = 0

            if sacdir == sacdirlist[1] && infile == filelist[1]
                # println(sacdir)
                file["info/DL_time_unit"]= sacin.t[1][2]/sacin.fs; #length of the data in time (duration), not number of samples;
                #convert from number of samples to time duration.
            end

            stemp = joinpath(ts,sacin.id[1]);
            # if verbose == true
            #     println("00 current: ",stemp)
            #     println("11 previous: ",stemp_pre)
            # end
            # t2=@elapsed if stemp == stemp_pre
            if stemp == stemp_pre
                if verbose == true
                    println(stemp," exists. Append to form multiple channels.")
                end
                append!(Sall[count],SeisData(sacin))
            else
                push!(Sall,SeisData(sacin))
                push!(Str_all,stemp)
                count += 1
                # file[stemp] = SeisData(sacin);
            end

            stemp_pre = stemp;

        end #end of loop for all sac files within one group/directory

        #2. start to save the data to jld2
        for i = 1:count
            file[Str_all[i]] = Sall[i]
        end #end of loop for all sac files within one group/directory

        stationlist = unique(vcat(stationlist,keys(file[ts])))
        if verbose == true
            print("                ----------------> ",length(filelist)," sac files, time used: ",time() - t0,"\n")
        end

    end #end of loop for all SAC directories.

    # println([t1all, t2all])
    #the following info dat is saved after running through the whole group.
    file["info/stationlist"] = stationlist;

    JLD2.close(file);

    print("Saved JLD2 file to: ", outfile, "\n");
end

"""
    sac2jld2_par(sacdirlist,timestamplist,outfile,verbose=false)

Searches for all sac files scanning throgh the sacdirlist (Array{String,1}). The sac files
    in each sac directory are assigned to one time stamp based on timestamplist (Array{String,1}).
    The two Arrays must be equal in size. The list can be generated from a file list (using ls in shell),
    using a utility function in this package: readfilelist(infilelist::String).

All sac files will be saved into one single JLD2 file (out), which is explicitly specified as outfile.

    Arguments:
    sacdirlist: An array of strings storing the directories grouped by time stamp, e.g., one day.
    timestamplist: An array of strings storing the time stamps assigning to each SAC file group.
    outfile: output file name for the JLD2 file. All data fiels will be readin and
                saved to a single JLD2 file.
    verbose: if true, the code will print out more messages. Default is false.
"""

# For parallel processing, the "info/DL_time_unit" is based on ONLY the first sac file.
function sac2jld2_par(sacdirlist::Array{String,1},timestamplist::Array{String,1},outfile::String,verbose::Bool=false)
    if length(sacdirlist) != length(timestamplist)
        error("sacdirlist and timestamplist must be the same length!")
    end

    stationlist = String[];
    #length of the data in time (duration), not number of samples
    file = jldopen(outfile,"w")

    file["info/DLtimestamplist"] = timestamplist;
    file["info/starttime"] = timestamplist[1];
    #here we use the first value in timestamplist as the starttime for a group of multiple timestamp data.
    file["info/endtime"] = timestamplist[end]; #similar to starttime, here we use the last value in timestamplist as the endtime.

    # function testread(x, filelist)
    #     #println(myid())
    #     SeisIO.read_data("sac",filelist[x],full=true)
    # end

    t1all = 0
    t2all = 0
    t3all = 0

    for (sacdir,ts) in zip(sacdirlist,timestamplist)
        t0=time()
        filelist=ls(joinpath(sacdir,"*.sac"));

        print("Converting for directory: [ ",sacdir," ] ... \n")

        S=pmap(x->SeisIO.read_data("sac",x,full=true),filelist)

        if ts == timestamplist[1]
            file["info/DL_time_unit"]= S[1].t[1][2]/S[1].fs;
        end
        count = 0
        # stationlisttemp = String[];
        #To deal with multiple channel issue, I am doing the loop twice.
        # 1. The first loop with collect multiple channels into one SeisData object.
        # 2. The second loop will work through the actual number of SeisData objects, counted by `count`
        Sall = SeisData[];
        Str_all = String[];

        stemp_pre = "-"
        for Sdata = S
            Sdata.misc[1]["dlerror"] = 0
            stemp = joinpath(ts,Sdata.id[1])
            if stemp == stemp_pre
                println(stemp," exists. Append to form multiple channels.")
                append!(Sall[count],SeisData(Sdata))
            else
                push!(Sall,SeisData(Sdata))
                push!(Str_all,stemp)
                count += 1
                # file[stemp] = SeisData(sacin);
            end
            stemp_pre = stemp;
        end

        #2. start to save the data to jld2
        for i = 1:count
            file[Str_all[i]] = Sall[i]
        end #end of loop for all sac files within one group/directory

        stationlist = unique(vcat(stationlist,keys(file[ts])))
        if verbose == true
            print("                ----------------> ",length(filelist)," sac files, time used: ",time() - t0,"\n")
        end

        # t1all += t1
    end #end of loop for all SAC directories.

    # println([t1all, t2all, t3all])

    #the following info dat is saved after running through the whole group.
    file["info/stationlist"] = stationlist;

    JLD2.close(file);

    print("Saved JLD2 file to: ", outfile, "\n");
    rmprocs(workers())
end


text:
# Example for Variable Projection:
# Author: Ziyi Yin, Rafael Orozco, Philipp Witte, Mathias Louboutin, Gabrio Rizzuti, Felix J. Herrmann
# Date: April, 2020
#

using DrWatson
@quickactivate "ExtSrcImg"

include("gen_geometry.jl")

using PyPlot, FFTW, Images
using JUDI, JOLI, JUDI4Flux
using LinearAlgebra, Random, Statistics
using IterativeSolvers
using JLD

function weight(n, d, xs, zs, nsrc,mode)
	if mode == "medical"
		zs_ave = zeros(1,n[2],nsrc); #hacky way
		zs_ave .= mean(zs)

		xs_ave = zeros(1,n[2],nsrc); #hacky way
		xs_ave .= mean(zs)

		x = zeros(n[1],1,nsrc);
	    x .= reshape((collect(1:n[1]).-1)*d[1], :, 1)
		z = zeros(1,n[2],nsrc);
	    z .= reshape((collect(1:n[2]).-1)*d[2], 1, :)

		W = zeros(Float32,n[1],n[2],nsrc)
		W .= sqrt.((x.-zs_ave).^2f0.+(z.-zs_ave).^2f0)

	elseif mode == "seismic"
		x = zeros(n[1],1,nsrc);
	    x .= reshape((collect(1:n[1]).-1)*d[1], :, 1)
		z = zeros(1,n[2],nsrc);
	    z .= reshape((collect(1:n[2]).-1)*d[2], 1, :)
		W = zeros(Float32,n[1],n[2],nsrc)
		W .= sqrt.((x.-xs).^2f0.+(z.-zs).^2f0)
		#W = sqrt.(r.^2f0.+δ^2f0)./δ
	end

    return reshape(W,n[1],n[2],1,nsrc)
end

function CG(d_obs, q0, m0, W; maxiter=5,λ=1f2)
	# use CG to solve damped least square normal equations
	r1 = d_obs-𝒢(q0, m0);
	r = 𝒢_T(r1,m0)-λ*W.^2f0.*q0;
	p = r;
	r_vec = zeros(maxiter+1);
	println("residual: $(norm(r1)) before updating q")
	for i = 1:maxiter
		if norm(r1)<1f-3
			break
		end
		temp = sum(sum(r.*r,dims=1),dims=2);
		temp1 = 𝒢(p, m0);
		temp2 = 𝒢_T(temp1,m0)+λ*W.^2f0.*p;
		alpha = temp./sum(sum(p.*temp2,dims=1),dims=2);
		q0 = q0 + alpha.* p;
		r1 = r1 - alpha.* temp1;
		r = r - alpha .* temp2;
		beta = sum(sum(r.*r,dims=1),dims=2)./temp;
		p = r + beta.*p;
		r_vec[i+1]=norm(r1)
		println("residual: $(norm(r1)) after updating q")
	end
	return q0, r1, r_vec
end

# Objective function: estimate source and compute residual + function value
function objective(d_obs, q0, m0, W; λ=1f2)

	# Data residual and source residual
  	r_data   = 𝒢(q0, m0) - d_obs
	r_source = q0

	# Function value
	fval = .5f0 *sum(r_data.^2) + 0.5f0*λ*sum((W.*r_source).^2)

	print("Function value: ", fval, "\n")
	return fval
end

function GenSimSourceMulti(xsrc_index, zsrc_index, nsrc, n)
	weights = zeros(Float32, n[1], n[2], 1, nsrc)

	num_simultaneous_sources = size(xsrc_index)[1]
	#check if Simultaneous
	if num_simultaneous_sources > 1 && nsrc == 1
		#for Simultaneous we need to put all the weights in the same slice.
		for j=1:num_simultaneous_sources
			#random weights
			weights[xsrc_index[j], zsrc_index[j], 1, 1] = rand()
        	#weights[xsrc_index[j], zsrc_index[j], 1, 1] = 1f0
    	end
	else
		for j=1:nsrc
        	weights[xsrc_index[j], zsrc_index[j], 1, j] = 1f0
    	end
	end

    return weights
end

function GenLocRec(xrec, zrec, nrec, n)
	weights = zeros(Float32, n[1], n[2], 1, 1)

	#for Simultaneous we need to put all the weights in the same slice.
	for j=1:nrec
        weights[Int(round(xrec[j]/d[1])), Int(round(zrec[j]/d[2])), 1, 1] = 1f0
    end

    return weights
end

function GenSimData(q₀; model=model)   # pass q₀ directly instead of z
    # Simultaneous observed data
    d_obs = F*vec(q₀)
    d_obs = reshape(d_obs, recGeometry.nt[1], nrec, 1, nsrc)  # reshape into tensor
    return d_obs
end

############################################ Models ###############################################

M = load("../data/sigsbee2A_model.jld")

m = M["m"][700:1900,100:800]                     # True model

N = load("../data/background_model.jld")
m0 = (N["m0"])[700:1900,100:800]   	# background model

dm = m - m0
n = size(m)
d = M["d"]
o = M["o"]

model0 = Model(n, d, o, m0)
model = Model(n, d, o, m)

m0 = reshape(model0.m, n[1], n[2], 1, 1)
m = reshape(model.m, n[1], n[2], 1, 1)

extentx = n[1]*d[1] / 1000 #in km
extentz = n[2]*d[2] / 1000 #in km

# receiver sampling and recording time
timeR = 5500f0   # receiver recording time [ms]
dtR = 4f0    # receiver sampling interval
# Source sampling and number of time steps
timeS = 5500f0
dtS = 4f0

mode = "seismic"
#mode = "medical"

if mode == "seismic"
	# Set up receiver structure
	xrec,yrec,zrec,nrec,nsrc =GenRecGeometrySeismic();
	recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

	# Set up source structure
	xsrc_index,zsrc_index = GenSrcGeometrySeismic(nsrc);

	#generate tensor with sources
	q       = GenSimSourceMulti(xsrc_index, zsrc_index, nsrc, model.n);

	#get muting region
	mute_mask = GenMuteSeismic()

elseif mode == "medical"

	# Set up receiver structure
	xrec, yrec, zrec, nrec, nsrc = GenRecGeometryMedical();
	recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

	# Set up source structure
	xsrc_index, zsrc_index, nsrc_simult = GenSrcGeometryMedical();

	#generate tensor with sources
	q           = GenSimSourceMulti(xsrc_index, zsrc_index, nsrc, model.n);

	#get muting region
	mute_mask = GenMuteMedical()

else
	error("Please choose a correct mode.")
end

figure(6); imshow(reshape(m, n[1], n[2])',extent=(0,extentx,extentz,0));
cb = colorbar(fraction=0.038, pad=0.04);
cb[:set_label]("Squared slowness [km^2/s^2]")
#PyPlot.scatter(xrec*d[1]/1000,zrec*d[1]/1000,marker=".",color="black",label="receivers");
PyPlot.scatter(xrec/1000,zrec/1000,marker=".",color="red",label="receivers");
#PyPlot.scatter(xsrc_index*d[1]/1000,zsrc_index*d[2]/1000,label="sources",marker=".",color="white");
legend(loc=1);
xlabel("Lateral position [km]");ylabel("Depth [km]");
fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

PyPlot.rc("font", family="serif"); PyPlot.rc("xtick", labelsize=8); PyPlot.rc("ytick", labelsize=8)
figure(); imshow(reshape(sqrt.(1f0./m), n[1], n[2])',extent=(0,extentx,extentz,0));
#cb = colorbar(fraction=0.027, pad=0.04,format="%.4f");
cb = colorbar(fraction=0.027, pad=0.04);
#cb[:set_label]("Velocity [km/s]", labelsize=20)
#cb[:set_label]("Squared slowness [km^2/s^2]")
PyPlot.scatter(xrec/1000,zrec/1000,marker=".",color="red",label="receivers",s=1);
PyPlot.scatter(xsrc_index*d[1]/1000,zsrc_index*d[2]/1000,label="sources",marker="x",color="white",s=1);
legend(loc=3,fontsize=8);
xlabel("Lateral position [km]", fontsize=8);ylabel("Horizontal position [km]", fontsize=8);
fig = PyPlot.gcf(); fig.set_size_inches(7.2, 4)
#savefig("seismic_setup.png", bbox_inches="tight")
savefig("seismic_setup.png")

# Setup wavelet
f0 = 0.015  # dominant frequency in [kHz]
wavelet = ricker_wavelet(timeS,dtS,f0)

######################################### Set up JUDI operators ###################################

# Set up info structure for linear operators
ntComp = get_computational_nt(recGeometry, model)
info   = Info(prod(model0.n), nsrc, ntComp)

# Write shots as array
opt = Options(return_array = true, subsampling_factor=20, isic=true)

# Modeling operators
Pr  = judiProjection(info, recGeometry)
F0  = judiModeling(info, model0; options=opt)	# modeling operator true model
F  = judiModeling(info, model; options=opt)	# modeling operator initial model
Pw = judiLRWF(info, wavelet)

# Combined operators
F = Pr*F*adjoint(Pw)
F0 = Pr*F0*adjoint(Pw)

# Extended modeling CNN layers
𝒢 = ExtendedQForward(F)
𝒢_T = ExtendedQAdjoint(F)

######################################### Run Forward simulations ###################################

# Generate shots
d_obs   = GenSimData(q); # Data in true model
######################################### Generate 3 gradient images ###################################

#### 1.) RTM
J1=judiJacobian(F0,q)
res_rtm = d_obs-𝒢(q,m0);
grad_rtm=reshape(adjoint(J1)* vec(res_rtm),n[1],n[2])

#### 2.) Extended source gradient - inversion scheme
#Generate weight matrix W. Used in CG for source estimation
if mode == "medical"
	xs = reshape(xsrc_index*d[1],:,1,nsrc_simult);
	zs = reshape(zsrc_index*d[2],:,1,nsrc_simult); #z coords of sources
elseif mode == "seismic"
	xs = reshape(xsrc_index*d[1],:,1,nsrc);
	zs = reshape(zsrc_index*d[2],:,1,nsrc); #z coords of sources
end

W =  weight(n, d, xs, zs, nsrc, mode);

#run cg to invert for estimated source

q_iter = 5
λ      = 100
q_inv = zeros(Float32, n[1], n[2], 1, nsrc);
q_inv, r_end, rvec = CG(d_obs, q_inv, m0, W; maxiter=q_iter,λ=λ)

#get jacobian using estimated source and use it to get gradient
J = judiJacobian(F0, q_inv)
#grad_extend=reshape(adjoint(J)* vec(mask.*r_end),n[1],n[2])
grad_extend=reshape(adjoint(J)* vec(r_end),n[1],n[2])
#### 3.) Extended source gradient - adjoint scheme
q_adj = 𝒢_T(d_obs,m0)
J2=judiJacobian(F0,q_adj) #check with francis on this
res_adj = d_obs-𝒢(q_adj,m0)
grad_adj =reshape(adjoint(J2)* vec(res_adj),n[1],n[2])

######################################### Visualize results   ###################################

#show gradient images

dm = reshape(m-m0,n[1],n[2])
mute_end = zeros(Integer,n[1])
for i = 1:n[1]
	for j = 1:n[2]
		if dm[i,j] >= 5f-3
			mute_end[i] = j
			break
		end
	end
end
D = judiDepthScaling(model)

new_inv=reshape(model_topmute(n,mute_end,5,grad_extend),n[1],n[2]);
final_inv = reshape(D*D*D*D*D*D*vec(deepcopy(new_inv)),n[1],n[2])
PyPlot.rc("font", family="serif"); PyPlot.rc("xtick", labelsize=8); PyPlot.rc("ytick", labelsize=8)
figure();imshow(final_inv'/norm(final_inv',Inf),vmin=-0.1,vmax=0.1,cmap="gray",extent=(0,extentx,extentz,0))
xlabel("Lateral position [km]", fontsize=8);ylabel("Horizontal position [km]", fontsize=8);
fig = PyPlot.gcf(); fig.set_size_inches(7.2, 4)
savefig("seismic_inv_imaging.png")

new_rtm=reshape(model_topmute(n,mute_end,5,grad_rtm),n[1],n[2]);
final_rtm = reshape(D*D*D*D*D*D*vec(deepcopy(new_rtm)),n[1],n[2])
PyPlot.rc("font", family="serif"); PyPlot.rc("xtick", labelsize=8); PyPlot.rc("ytick", labelsize=8)
figure();imshow(final_rtm'/norm(final_rtm',Inf),vmin=-0.1,vmax=0.1,cmap="gray",extent=(0,extentx,extentz,0))
xlabel("Lateral position [km]", fontsize=8);ylabel("Horizontal position [km]", fontsize=8);
fig = PyPlot.gcf(); fig.set_size_inches(7.2, 4)
savefig("seismic_rtm_imaging.png")

new_adj=reshape(model_topmute(n,mute_end,5,grad_adj),n[1],n[2]);
final_adj = reshape(D*D*D*D*D*D*vec(deepcopy(new_adj)),n[1],n[2])
PyPlot.rc("font", family="serif"); PyPlot.rc("xtick", labelsize=8); PyPlot.rc("ytick", labelsize=8)
figure();imshow(final_adj'/norm(final_adj',Inf),vmin=-0.1,vmax=0.1,cmap="gray",extent=(0,extentx,extentz,0))
cb = colorbar(fraction=0.027, pad=0.04);
"Horizontal position [km]", fontsize=8);
fig = PyPlot.gcf(); fig.set_size_inches(5.6, 3.5)
savefig("seismic_adj_imaging.png", bbox_inches="tight")

if mode == "seismic"	# one single extended source in the middle

	q_inv_stack=sum(q_inv,dims=4)[:,:,1,1];
	figure();imshow(q_inv_stack',vmin=-2f-3,vmax=2f-3,extent=(0,extentx,extentz,0));title("Extended Source from Inversion");
	xlabel("Lateral position [samples]");ylabel("Depth [samples]");
	colorbar(fraction=0.028, pad=0.04);
	fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

	q_adj_stack=sum(q_adj,dims=4)[:,:,1,1];
	figure();imshow(q_adj_stack',vmin=-40000,vmax=40000,extent=(0,extentx,extentz,0));title("Extended Source from Adjoint");
	xlabel("Lateral position [samples]");ylabel("Depth [samples]");
	colorbar(fraction=0.028, pad=0.04);
	fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

	if nsrc != 1
		mid_loc = Int(round(nsrc/2)+1)
		q_inv_single=q_inv[:,:,1,mid_loc];
		figure(18);imshow(q_inv_single',vmin=-2f-3,vmax=2f-3,extent=(0,extentx,extentz,0));title("Extended Source from Inversion");
		xlabel("Lateral position [samples]");ylabel("Depth [samples]");
		colorbar(fraction=0.028, pad=0.04);
		fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

		q_adj_single=q_adj[:,:,1,mid_loc];
		figure(19);imshow(q_adj_single',vmin=-40000,vmax=40000,extent=(0,extentx,extentz,0));title("Extended Source from Adjoint");
		xlabel("Lateral position [samples]");ylabel("Depth [samples]");
		colorbar(fraction=0.028, pad=0.04);
		fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)
	end

elseif mode == "medical"
	#show estimated sources
	q_adj_stack=sum(q_adj,dims=4)[:,:,1,1];
	figure(70);imshow(q_adj_stack',extent=(0,extentx,extentz,0),vmin=-40000,vmax=40000);title("Extended Source from Adjoint");
	xlabel("Lateral position [km]");ylabel("Depth [km]");
	colorbar(fraction=0.028, pad=0.04);
	fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

	q_inv_stack=sum(q_inv,dims=4)[:,:,1,1];
	figure(8);imshow(q_inv_stack',extent=(0,extentx,extentz,0),vmin=-2f-3,vmax=2f-3);title("Extended Source from Inversion");
	xlabel("Lateral position [km]");ylabel("Depth [km]");
	colorbar(fraction=0.028, pad=0.04);
	fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

	#show estimated sources
	q_adj_stack=sum(q_adj,dims=4)[:,:,1,1];
	figure(70);imshow(q_adj_stack',extent=(0,extentx,extentz,0));title("Extended Source from Adjoint");
	xlabel("Lateral position [km]");ylabel("Depth [km]");
	colorbar(fraction=0.028, pad=0.04);
	fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

	q_inv_stack=sum(q_inv,dims=4)[:,:,1,1];
	figure(8);imshow(q_inv_stack',extent=(0,extentx,extentz,0));title("Extended Source from Inversion");
	xlabel("Lateral position [km]");ylabel("Depth [km]");
	colorbar(fraction=0.028, pad=0.04);
	fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

	q_stack=sum(q,dims=4)[:,:,1,1];
	figure(19);imshow(q_stack',extent=(0,extentx,extentz,0));title("Real Source");
	xlabel("Lateral position [km]");ylabel("Depth [km]");
	colorbar(fraction=0.028, pad=0.04);
	fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)
end

######################################### Save results   ########################################

#@save "1src_precondition.jld" grad_adj grad_rtm grad_extend q_inv q_adj m m0 mute_mask xrec yrec zrec xsrc_index zsrc_index rvec


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module DynamicModel

export
    AbstractDynamicModel,
    addObserver,

    initialize,
    update,

    initializeDynamicModel,
    simulateDynamicModel,

    SimpleADM,
    SimpleADMInitialState,
    SimpleADMOutputState,
    SimpleADMCommand,

    LLADM,
    LLADMState,
    LLADMOutputState,
    LLADMCommand

import CommonInterfaces: initialize, update

using AbstractDynamicModelImpl

using SimpleADMImpl

using LLADMImpl

end


text:
module RegisterDeformation

using ImageCore, ImageAxes, Interpolations, StaticArrays, HDF5, JLD2, ProgressMeter
using RegisterUtilities, LinearAlgebra, Rotations, Base.Cartesian
using Distributed, Statistics, SharedArrays
using Base: tail
using Interpolations: AbstractInterpolation, AbstractExtrapolation
import ImageTransformations: warp, warp!
# to avoid `scale` conflict with Interpolations, selectively import CoordinateTransformations:
using CoordinateTransformations: AffineMap
import CoordinateTransformations: compose, transform
using OffsetArrays: IdentityUnitRange   # for Julia-version compatibility

export
    # types
    AbstractDeformation,
    GridDeformation,
    WarpedArray,
    # functions
    arraysize,       # TODO: don't export?
    centeraxes,
    compose,
    eachnode,
    extrapolate,
    extrapolate!,
    getindex!,
    griddeformations,
    interpolate,
    interpolate!,
    medfilt,
    nodegrid,
    regrid,
    similarϕ,
    tform2deformation,
    tinterpolate,
    translate,
    vecindex,        # TODO: don't export?
    vecgradient!,    # TODO: don't export?
    warp,
    warp!,
    warpgrid,
    # old AffineTransfroms.jl code (TODO?: remove)
    tformtranslate,
    tformrotate,
    tformeye,
    rotation2,
    rotation3,
    rotationparameters,
    transform!,
    transform

const DimsLike = Union{Vector{Int}, Dims}
const InterpExtrap = Union{AbstractInterpolation,AbstractExtrapolation}

"""
# RegisterDeformation

A deformation (or warp) of space is represented by a function `ϕ(x)`.
For an image, the warped version of the image is specified by "looking
up" the pixel value at a location `ϕ(x) = x + u(x)`.  `u(x)` thus
expresses the displacement, in pixels, at position `x`.  Note that a
constant deformation, `u(x) = x0`, corresponds to a shift of the
*coordinates* by `x0`, and therefore a shift of the *image* in the
opposite direction.

In reality, deformations will be represented on a grid, and
interpolation is implied at locations between grid points. For a
deformation defined directly from an array, make it interpolating
using `ϕi = interpolate(ϕ)`.

The major functions/types exported by RegisterDeformation are:

    - `GridDeformation`: create a deformation
    - `tform2deformation`: convert an `AffineMap` to a deformation
    - `ϕ_old(ϕ_new)` and `compose`: composition of two deformations
    - `warp` and `warp!`: deform an image
    - `WarpedArray`: create a deformed array lazily
    - `warpgrid`: visualize a deformation

"""
RegisterDeformation

abstract type AbstractDeformation{T,N} end
Base.eltype(::Type{AbstractDeformation{T,N}}) where {T,N} = T
Base.ndims(::Type{AbstractDeformation{T,N}}) where {T,N} = N
Base.eltype(::Type{D}) where {D<:AbstractDeformation} = eltype(supertype(D))
Base.ndims(::Type{D}) where {D<:AbstractDeformation} = ndims(supertype(D))
Base.eltype(d::AbstractDeformation) = eltype(typeof(d))
Base.ndims(d::AbstractDeformation) = ndims(typeof(d))

include("griddeformation.jl")
include("utils.jl")
include("timeseries.jl")
include("tformedarrays.jl")
const Extrapolatable{T,N} = Union{TransformedArray{T,N},AbstractExtrapolation{T,N}}
include("warpedarray.jl")
include("warp.jl")
include("visualize.jl")

include("deprecated.jl")

end  # module


text:
using Plots
plotlyjs();

include("../julia/sskernel.jl");
include("../julia/utils.jl");

x_true = readcsv("../data/spikes.csv")[:];
rs_true = readcsv("../data/firing_rate_true.csv")[:];

ws = linspace(0.01,2.5,100);
rs_kernel = sskernel(x_true, ws=ws);

plot(ts, rs_true, label="True")
plot!(ts, rs_kernel, label="Estimated")

xlabel!("Time (s)")
ylabel!("Firing rate (Hz)")

# C = get_costfunc(x_inhomogeneous)

# plot(ws, C.(ws));
# xlabel!("w");
# ylabel!("C(w)");

max_stars_repo_name: ganow/julia-point-process

max_issues_repo_name: ganow/julia-point-process

max_forks_repo_name: ganow/julia-point-process


text:
module DTWs

## A small module to perform Dynamic Time Warping alignment. ##

import StatsBase: fit!

export DTW, fit!, update!, set_template!, backward

# type DTW represents a dynamic time warping.
# TODO(ryuichi) should be more generic
type DTW
    fstep::Int  # forward transition constraint
    bstep::Int  # backward transition constraint
    template::Matrix{Float64}
    costtable::Matrix{Float64}
    backpointer::Matrix{Int}
end

function DTW(;fstep=0, bstep=1)
    DTW(fstep, bstep, zeros(1, 1), zeros(1, 1), zeros(Int, 1, 1))
end

function transition(d::DTW, i::Int, j::Int)
    if j == i+1
        return 0.0
    elseif i == j
        return 1.0
    end

    2.0 # abs(i-j)
end

function observation(d::DTW, v::AbstractVector, i::Int)
    sumabs2(v  - d.template[:,i])
end

# lazy_init! performs state initialization.
function lazy_init!(d::DTW, S::Int)
    d.costtable = reshape(1:S, S, 1)
    d.backpointer = reshape(1:S, S, 1)
end

# lazy_init! performs pre-allocations
function lazy_init!(d::DTW, S::Int, T::Int)
    # plus initial state
    d.costtable = zeros(Float64, S, T+1)
    d.backpointer = ones(Int, S, T+1)

    d.costtable[:,1] = [1:S;]
    d.backpointer[:,1] = [1:S;]
end

function set_template!(d::DTW, template::Matrix{Float64})
    d.template = template
    lazy_init!(d, size(template, 2))
end

# update! updates cost table for a given vector in an on-line manner.
# In off-line situations, please use fit! because it is more efficient
# and faster than this on-line version.
function update!(d::DTW, v::AbstractVector)
    # S: length of tempalte, T: current time length
    S, T = size(d.costtable)
    lastcost = d.costtable[:, T]
    currentcost = zeros(S)
    current_backpointer = zeros(Int, S)

    for i=1:S
        minindex = i
        const obs = observation(d, v, i)
        const trans = transition(d, minindex, i)
        mincost = lastcost[minindex] + obs + trans

        # search minindex
        for j=i-d.bstep:i+d.fstep
            if j < 1 || j > S
                continue
            end
            cost = lastcost[j] + obs + transition(d, j, i)
            if cost < mincost
                mincost, minindex = cost, j
            end
        end
        currentcost[i] = mincost
        current_backpointer[i] = minindex
    end

    d.costtable = [d.costtable currentcost]
    d.backpointer = [d.backpointer current_backpointer]
end

# fit! aligns two sequences using dynamic time warping algorithm.
function fit!(d::DTW, template::Matrix{Float64}, sequence::Matrix{Float64})
    # S: length of template, T: length of target sequence
    S, T = size(template, 2), size(sequence, 2)

    # pre-allocations
    lazy_init!(d, S, T)
    @assert size(d.costtable) == (S, T+1)

    # set template matrix to align
    d.template = template

    for t=1:T
        v = sequence[:,t]
        for i=1:S
            minindex = i
            ocost = observation(d, v, i)
            tcost = transition(d, minindex, i)
            @inbounds mincost = d.costtable[minindex, t] + ocost + tcost

            # search minindex
            for j=i-d.bstep:i+d.fstep
                if j < 1 || j > S
                    continue
                end
                @inbounds cost = d.costtable[j, t] + ocost + transition(d, j, i)
                if cost < mincost
                    mincost, minindex = cost, j
                end
            end
            @inbounds d.costtable[i, t+1] = mincost
            @inbounds d.backpointer[i, t+1] = minindex
        end
    end

    backward(d)
end

fit!(d::DTW, sequence::Matrix{Float64}) = fit!(d, d.template, sequence)

# backward searches the path that minimizes the total cost.
function backward(d::DTW)
    T  = size(d.costtable, 2) - 1 # exclude the initial state
    minpath = zeros(Int, T)

    minpath[end] = indmin(d.costtable[:,T+1])

    # revursive search
    for i=reverse(2:T)
        @inbounds minpath[i-1] = d.backpointer[minpath[i], i+1]
    end

    minpath
end

end # module