text:
function remove_isolates(g)
    vs = findall(degree(g) .> 0)
    induced_subgraph(g, vs)
end

function vmax_neighbors(g, v, w, nbfun)
    nb = nbfun(g, v)
    length(nb) == 0 && return nb
    wnb = w[nb]
    return nb[wnb .== maximum(wnb)]
end

function emax_neighbors(g, v, w, nbfun)
    nb = nbfun(g, v)
    length(nb) == 0 && return nb
    if nbfun == inneighbors
        wnb = w[nb, v]
    else
        wnb = w[v, nb]
    end
    nb[wnb .== maximum(wnb)]
end

"""
    bfs_multi(g, source, neighborfn)

Traverse graph `g` via breadth-first search starting at vertices `source` and
choosing neighbors to visit according to `neighborfn`.

This is a modified breadth-first search traversal where instances of
a vertex being visited multiple times by different parents are recorded.
This function was modified from the Graphs bfs algorithm.
"""
function bfs_multi(g::AbstractGraph{T}, source, neighborfn::Function) where T
    n = nv(g)
    visited = falses(n)
    parents = [Vector{T}() for _ in 1:n]
    cur_level = Vector{T}()
    sizehint!(cur_level, n)
    next_level = Vector{T}()
    sizehint!(next_level, n)
    @inbounds for s in source
        visited[s] = true
        push!(cur_level, s)
        parents[s] = [s]
    end
    while !isempty(cur_level)
        @inbounds for v in cur_level
            @inbounds @simd for i in neighborfn(g, v)
                if visited[i]
                    push!(parents[i], v)
                else
                    push!(next_level, i)
                    push!(parents[i], v)
                    visited[i] = true
                end
            end
        end
        empty!(cur_level)
        cur_level, next_level = next_level, cur_level
        sort!(cur_level)
    end
    return parents
end


text:
module Key
    struct lweParams
        n::Int32
        N::Int32
        l::Int32
        Bgbit::Int32
        t::Int32
        basebit::Int32
        α::Float64
        αᵦₖ::Float64
        αₖₛ::Float64
        offset::Int32
        h::Array{Float64}
        decb::Array{Float64}
    end
    function lweParams(n::Int32,N::Int32,α::Float64,l::Int32,Bgbit::Int32,αᵦₖ::Float64,t::Int32,basebit::Int32,αₖₛ::Float64)
        Bg::Float64 = 2.0 ^ Bgbit
        h::Array{Float64} = [Bg^(-(i+1)) for i in 1:l]
        decb::Array{Float64} = [2.0^(-32) * Bg^(i+1) for i in 1:l]
        offset::UInt32 = Bg/2 * sum(2^32 * h)
        return lweParams(n,N,l,Bgbit,t,basebit,α,αᵦₖ,αₖₛ,offset,h,decb)
    end

    struct lweKey
        tlwe::Array{UInt32,1}
        trlwe::Array{UInt32,1}
    end
    function lweKey(n::Int32,N::Int32)
        return lweKey(rand(0:1,n),rand(0:1,N))
    end

    struct SecretKey
        params::lweParams
        key::lweKey
    end
    function SecretKey(n::Int32,N::Int32,alpha::Float64,l::Int32,Bgbit::Int32,bkalpha::Float64,t::Int32,basebit::Int32,ksalpha::Float64)
        return SecretKey(lweParams(n,N,alpha,l,Bgbit,bkalpha,t,basebit,ksalpha),lweKey(n,N))
    end
    const μ = 2.0^-3
end

ext: jl


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

# Decorated intervals, following the IEEE 1758-2015 standard

"""
    DECORATION

Enumeration constant for the types of interval decorations.
The nomenclature of the follows the IEEE-1788 (2015) standard
(sect 11.2):

- `com -> 4`: common: bounded, non-empty
- `dac -> 3`: defined (nonempty) and continuous
- `def -> 2`: defined (nonempty)
- `trv -> 1`: always true (no information)
- `ill -> 0`: nai ("not an interval")
"""
@enum DECORATION ill=0 trv=1 def=2 dac=3 com=4
# Note that `isless`, and hence ``<` and `min`, are automatically defined for enums

# const decorations = Dict("ill" => ill,
#                          "trv" => trv,
#                          "def" => def,
#                          "dac" => dac,
#                          "com" => com)

"""
    DecoratedInterval

A `DecoratedInterval` is an interval, together with a *decoration*, i.e.
a flag that records the status of the interval when thought of as the result
of a previously executed sequence of functions acting on an initial interval.
"""
struct DecoratedInterval{T<:Real} <: AbstractInterval{T}
    interval::Interval{T}
    decoration::DECORATION

    function DecoratedInterval{T}(I::Interval, d::DECORATION) where T
        dd = decoration(I)
        dd <= trv && return new{T}(I, dd)
        d == ill && return new{T}(nai(I), d)
        return new{T}(I, d)
    end
end

DecoratedInterval(I::Interval{T}, d::DECORATION) where T<:AbstractFloat =
    DecoratedInterval{T}(I, d)

function DecoratedInterval(a::T, b::T, d::DECORATION) where T<:Real
    a > b && return DecoratedInterval(nai(T), ill)
    DecoratedInterval(Interval(a,b), d)
end

DecoratedInterval(a::T, d::DECORATION) where T<:Real =
    DecoratedInterval(Interval(a,a), d)

DecoratedInterval(a::Tuple, d::DECORATION) = DecoratedInterval(a..., d)

DecoratedInterval(a::T, b::S, d::DECORATION) where {T<:Real, S<:Real} =
    DecoratedInterval(promote(a,b)..., d)

# Automatic decorations for an interval
DecoratedInterval(I::Interval) = DecoratedInterval(I, decoration(I))

function DecoratedInterval(a::T, b::T) where T<:Real
    a > b && return DecoratedInterval(nai(T), ill)
    DecoratedInterval(Interval(a,b))
end

DecoratedInterval(a::T, b::T) where T<:Integer =
    DecoratedInterval(float(a),float(b))

DecoratedInterval(a::T) where T<:Real = DecoratedInterval(Interval(a,a))

DecoratedInterval(a::T, b::S) where {T<:Real, S<:Real} =
    DecoratedInterval(promote(a, b)...)

DecoratedInterval(a::Tuple) = DecoratedInterval(a...)

DecoratedInterval(I::DecoratedInterval, dec::DECORATION) = DecoratedInterval(I.interval, dec)

interval_part(x::DecoratedInterval) = x.interval

decoration(x::DecoratedInterval) = x.decoration

# Automatic decorations for an Interval
function decoration(I::Interval)
    isnai(I) && return ill           # nai()
    isempty(I) && return trv         # emptyinterval
    isunbounded(I) && return dac     # unbounded
    com                              # common
end

# Promotion and conversion, and other constructors

# promote_rule{T<:Real, N, R<:Real}(::Type{DecoratedInterval{T}},
#     ::Type{ForwardDiff.Dual{N,R}}) = ForwardDiff.Dual{N, DecoratedInterval{promote_type(T,R)}}

promote_rule(::Type{DecoratedInterval{T}}, ::Type{S}) where {T<:Real, S<:Real} =
    DecoratedInterval{promote_type(T, S)}

promote_rule(::Type{DecoratedInterval{T}}, ::Type{DecoratedInterval{S}}) where
    {T<:Real, S<:Real} = DecoratedInterval{promote_type(T, S)}

convert(::Type{DecoratedInterval{T}}, x::S) where {T<:Real, S<:Real} =
    DecoratedInterval( Interval(T(x, RoundDown), T(x, RoundUp)) )

convert(::Type{DecoratedInterval{T}}, x::S) where {T<:Real, S<:Integer} =
    DecoratedInterval( Interval(T(x), T(x)) )
# function convert{T<:AbstractFloat}(::Type{DecoratedInterval{T}}, x::Float64)
#     convert(DecoratedInterval{T}, rationalize(x))
# end
function convert(::Type{DecoratedInterval{T}}, xx::DecoratedInterval) where T<:Real
    x = interval_part(xx)
    x = atomic(Interval{T},x)
    DecoratedInterval( x, decoration(xx) )
end

convert(::Type{DecoratedInterval{T}}, x::AbstractString) where T<:AbstractFloat =
    parse(DecoratedInterval{T}, x)

big(x::DecoratedInterval) = DecoratedInterval(big(interval_part(x)),
                                                decoration(x))

macro decorated(ex...)
    local x

    if length(ex) == 1
        x = :(@interval($(esc(ex[1]))))
    else
        x = :($(esc(ex[1])), $(esc(ex[2])))
    end

    :(DecoratedInterval($x))
end


text:
#tbeg = time()
using LinearAlgebra
BLAS.set_num_threads(4)
using Flux, Statistics # Flux.Data.MNIST
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated, partition
using BSON, HDF5, JLD, Random
using MLDataUtils
using NPZ
include("genDataScripts.jl")
n_sites = 100
periods = [2,3,4,5,6,7,8,9,10,11]
all_signals = gen_signals(n_sites, periods)

neu = 100
#trainingSizes = [50, 100, 250, 500, 1000, 2000, 4000]
#trainSizes = [10000, 20000, 40000]
trainSizes = [100000]

# models are here
fpath = "/data/finite/lalonso/evaluations/"

sigmasAll = collect(0.05:0.05:4)
sigmasAllselect = [0.5,0.75,1,1.25,1.5,2]
indice = parse(Int64, ARGS[1])
#indice = 9
#σv = round(σvt*0.05, digits=3)
#σv = 2.5 # 1.0

grid = []
for ts in trainSizes, σv in sigmasAll
    push!(grid, [ts, σv])
end
tssize = Int64(grid[indice][1])
trainingSizes = [tssize]
σv = grid[indice][2]
iteraciones = 100

# HHG DATA
x_data9k = npzread("./data/data_simulation2_optimized_9669.npy")
labels9k = npzread("./data/labels_simulation2_optimized_9669.npy")
ycold9k = Int64.(labels9k[:,1,1])
nsamples = size(x_data9k)[1]
npts = 100
data_log9k = zeros(npts, nsamples)
for i in 1:nsamples
    input_log = log10.(x_data9k[i,:,2])
    input_log = input_log .+ minimum(input_log)*(-1)
    #input_log = input_log/maximum(input_log)
    data_log9k[:,i] .= input_log
end

accuracy_tcold(x,y,model) = mean(onecold(model(x)) .== y)

function signs100k(σv)
    Random.seed!(1234)
    signals_mas_ruido = zeros(npts,nsamples)
    for i in 1:nsamples
        signals_mas_ruido[:,i] .= signal_plus_noise(data_log9k[:,i], σv, true)
    end
    for _ in 1:9
        signals_mas_ruidox = zeros(npts, nsamples)
        for i in 1:nsamples
            signals_mas_ruidox[:,i] .= signal_plus_noise(data_log9k[:,i], σv, true)
        end
        signals_mas_ruido = hcat(signals_mas_ruido, signals_mas_ruidox)
    end
    signals_mas_ruido
end

ycold9kx = repeat(ycold9k, 10)
ycold9kxh = onehotbatch(ycold9kx, 1:10)
# Synthetic data
function synthData(σv)
    all_signals = gen_signals(n_sites, periods)
    Random.seed!(1234)
    num_test = 10000
    x_test, y_test = gen_data_set(all_signals, σv, num_test,true)
    xtest_temp, ytest = reshape(x_test, (100, num_test, 10)), reshape(y_test, (num_test, 10))
    x_test = reshape(xtest_temp, (n_sites, 10*num_test))
    labels_test = reshape(ytest, (10*num_test))
    x_test, labels_test
end

Random.seed!(1234)
t_size = 1_000_000 # 4_000_000 old
num_train = Int64(round.(t_size *0.8/10)) # number of samples for each class,
#this just works for 10 classes
num_val = Int64(round.(t_size *0.2/10))
num_test = 1000

x_test, y_test = gen_data_set(all_signals, σv, num_test,true)

x_train, y_train = gen_data_set(all_signals, σv, num_train,true)
x_vali, y_vali = gen_data_set(all_signals, σv, num_val,true)

x_test = reshape(x_test, (100, num_test, 10))
y_test = reshape(y_test, (num_test, 10))
x_train = reshape(x_train, (100, num_train, 10))
y_train = reshape(y_train, (num_train, 10))
x_vali = reshape(x_vali, (100, num_val, 10))
y_vali = reshape(y_vali, (num_val, 10))

xtest_temp = x_test[:,1:num_test,:]
ytest_temp = y_test[1:num_test,:]
x_test = reshape(xtest_temp, (n_sites, 10*num_test))
labels_test = reshape(ytest_temp, (10*num_test))

X_test = x_test # |> gpu
Y_test = onehotbatch(labels_test, 1:10); #|> gpu

accuracyT(x, y, model) = mean(onecold(model(x)) .== onecold(y))
accuracybySample(x, y, model) = onecold(model(x)) .== onecold(y)

#t_size = 100_000

trainSizes = [100_000]
t_size=100_000

num_train = Int64(round.(t_size *0.8/10)) #number of samples for each class,
#this just works for 10 classes
num_val = Int64(round.(t_size *0.2/10))

xtrain_temp = x_train[:,1:num_train,:]
ytrain_temp = y_train[1:num_train,:]

xvali_temp = x_vali[:,1:num_val,:]
yvali_temp = y_vali[1:num_val,:]

x_val = reshape(xvali_temp, (n_sites, 10*num_val))
labels_val = reshape(yvali_temp, (10*num_val))

x_data = reshape(xtrain_temp, (n_sites, 10*num_train))
labels = reshape(ytrain_temp, (10*num_train))

X_val = x_val # |> gpu
Y_val = onehotbatch(labels_val, 1:10); #|> gpu

# Data preparation
X = x_data # |> gpu
Y = onehotbatch(labels, 1:10); #|> gpu

trainErrors = zeros(10,2,100)
valErrors = zeros(10,2,100)
testErrors = zeros(10,2,100)
testErrors_hhg = zeros(10,2,100,length(sigmasAllselect))
testErrors_s = zeros(10,2,100,length(sigmasAllselect))

for i in 1:100
    model = BSON.load(fpath*"modelsN100ESP_ts$(t_size)_σ$(σv)_ite$(i).bson")[:model]

    test_acc = accuracybySample(X_test, Y_test, model)
    train_acc = accuracybySample(X, Y, model)
    val_acc = accuracybySample(X_val, Y_val, model)

    train_fails = findall(iszero, train_acc*1)
    val_fails = findall(iszero, val_acc*1)
    test_fails = findall(iszero, test_acc*1)

    if length(train_fails) > 0
        samples_train = X[:,train_fails]
        pred_train = onecold(model(samples_train))
        true_train =  onecold(Y[:,train_fails])

        u1 = unique(pred_train)
        u2 = unique(true_train)
        d1 = Dict([(j, count(x->x==j,pred_train)) for j in u1])
        d2 = Dict([(j, count(x->x==j,true_train)) for j in u2])
        vals1 = hcat([[key, val] for (key, val) in d1]...)'
        vals2 = hcat([[key, val] for (key, val) in d2]...)'

        trainErrors[vals1[:,1],1,i] = vals1[:,2]
        trainErrors[vals2[:,1],2,i] = vals2[:,2]

    end
    if length(val_fails)>0
        samples_val = X_val[:,val_fails]
        pred_val = onecold(model(samples_val))
        true_val = onecold(Y[:,val_fails])
        u1 = unique(pred_val)
        u2 = unique(true_val)
        d1 = Dict([(j, count(x->x==j,pred_val)) for j in u1])
        d2 = Dict([(j, count(x->x==j,true_val)) for j in u2])
        vals1 = hcat([[key, val] for (key, val) in d1]...)'
        vals2 = hcat([[key, val] for (key, val) in d2]...)'

        valErrors[vals1[:,1],1,i] = vals1[:,2]
        valErrors[vals2[:,1],2,i] = vals2[:,2]

    end
    if length(test_fails)>0
        samples_test = X_test[:,test_fails]
        pred_test = onecold(model(samples_test))
        true_test = onecold(Y[:,test_fails])
        u1 = unique(pred_test)
        u2 = unique(true_test)
        d1 = Dict([(j, count(x->x==j,pred_test)) for j in u1])
        d2 = Dict([(j, count(x->x==j,true_test)) for j in u2])
        vals1 = hcat([[key, val] for (key, val) in d1]...)'
        vals2 = hcat([[key, val] for (key, val) in d2]...)'

        testErrors[vals1[:,1],1,i] = vals1[:,2]
        testErrors[vals2[:,1],2,i] = vals2[:,2]

    end

    for (indj, σvd) in enumerate(sigmasAllselect)
        signals_mas_ruido = signs100k(σvd)
        test_acc_hhg = accuracybySample(signals_mas_ruido, ycold9kxh, model)
        test_fails_hhg = findall(iszero, test_acc_hhg*1)

        if length(test_fails_hhg)>0
            samples_test = signals_mas_ruido[:,test_fails_hhg]
            pred_test = onecold(model(samples_test))
            true_test = onecold(ycold9kxh[:,test_fails_hhg])
            u1 = unique(pred_test)
            u2 = unique(true_test)
            d1 = Dict([(j, count(x->x==j, pred_test)) for j in u1])
            d2 = Dict([(j, count(x->x==j, true_test)) for j in u2])
            vals1 = hcat([[key, val] for (key, val) in d1]...)'
            vals2 = hcat([[key, val] for (key, val) in d2]...)'
            testErrors_hhg[vals1[:,1], 1, i, indj] = vals1[:,2]
            testErrors_hhg[vals2[:,1], 2, i, indj] = vals2[:,2]

        end

    end
    for (indj, σvd) in enumerate(sigmasAllselect)
        x_test_s, labels_test_s = synthData(σvd)
        labels_test_sh = onehotbatch(labels_test_s, 1:10)
        test_acc_s = accuracybySample(x_test_s, labels_test_sh, model)
        test_fails_s = findall(iszero, test_acc_s*1)
        if length(test_fails_s)>0
            samples_test = x_test_s[:,test_fails_s]
            pred_test = onecold(model(samples_test))
            true_test = onecold(labels_test_sh[:,test_fails_s])
            u1 = unique(pred_test)
            u2 = unique(true_test)
            d1 = Dict([(j, count(x->x==j, pred_test)) for j in u1])
            d2 = Dict([(j, count(x->x==j, true_test)) for j in u2])
            vals1 = hcat([[key, val] for (key, val) in d1]...)'
            vals2 = hcat([[key, val] for (key, val) in d2]...)'
            testErrors_s[vals1[:,1], 1, i, indj] = vals1[:,2]
            testErrors_s[vals2[:,1], 2, i, indj] = vals2[:,2]

        end
    end
    #println("iteration: $(i)")
end

npzwrite("./bump/bump_hhg_$(t_size)_sig_$(σv).npy", testErrors_hhg)
npzwrite("./bump/bump_test100k_$(t_size)_sig_$(σv).npy", testErrors_s)
npzwrite("./bump/bump_train_$(t_size)_sig_$(σv).npy", trainErrors)
npzwrite("./bump/bump_val_$(t_size)_sig_$(σv).npy", valErrors)
npzwrite("./bump/bump_test10k_$(t_size)_sig_$(σv).npy", testErrors)

max_stars_repo_name: lazarusA/noisySignals

max_issues_repo_name: lazarusA/noisySignals

max_forks_repo_name: lazarusA/noisySignals


text:
export PlaneWaves, PolarizationStyle, TE, TM, Coupled

struct PlaneWaves{N}
    ω::Float64
    dims::NTuple{N, Tuple{Int64, Float64}}
    x⃗::NTuple{N, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}
    k⃗::NTuple{N, Frequencies{Float64}}
end

"""
    PlaneWaves(ω::Float64, dims::Tuple{Vararg{Tuple{Int64, Float64}}})

This constructor takes a frequency, and an NTuple of Tuple{Int64, Float64}.
Each tuple represents a periodic dimension, and the Int64 represents the number
of grid points per unit cell in that periodic dimension, and the Float64 gives
the lattice constant in that periodic dimension (length of unit cell).
"""
function PlaneWaves(ω::Float64, dims::Tuple{Vararg{Tuple{Int64, Float64}}})
    x⃗ = Tuple(range(0, step=e[2] / e[1], length=e[1]) for e in dims)
    k⃗ = Tuple(2π * fftfreq(e[1], e[1] / e[2]) for e in dims)
    PlaneWaves(ω, dims, x⃗, k⃗)
end

Base.size(pw::PlaneWaves) = Tuple(e[1] for e in pw.dims)
Base.length(pw::PlaneWaves) = prod(size(pw))

abstract type PolarizationStyle end
struct TE <: PolarizationStyle end
struct TM <: PolarizationStyle end
struct Coupled <: PolarizationStyle end

PolarizationStyle(pw::PlaneWaves) = PolarizationStyle(typeof(pw))
PolarizationStyle(::Type{PlaneWaves{1}}) = TM()
PolarizationStyle(::Type{PlaneWaves{2}}) = Coupled()


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct PeriodicVariableSource <: DynamicInjection
        name::String
        R_th::Float64
        X_th::Float64
        internal_voltage_bias::Float64
        internal_voltage_frequencies::Vector{Float64}
        internal_voltage_coefficients::Vector{Tuple{Float64,Float64}}
        internal_angle_bias::Float64
        internal_angle_frequencies::Vector{Float64}
        internal_angle_coefficients::Vector{Tuple{Float64,Float64}}
        base_power::Float64
        states::Vector{Symbol}
        n_states::Int
        ext::Dict{String, Any}
        internal::InfrastructureSystemsInternal
    end

This struct acts as an infinity bus with time varying phasor values magnitude and angle V(t) 	heta(t). Time varying functions are represented using fourier series

# Arguments
- `name::String`
- `R_th::Float64`: Source Thevenin resistance, validation range: `(0, nothing)`
- `X_th::Float64`: Source Thevenin reactance, validation range: `(0, nothing)`
- `internal_voltage_bias::Float64`: a0 term of the Fourier Series for the voltage
- `internal_voltage_frequencies::Vector{Float64}`: Frequencies in radians/s
- `internal_voltage_coefficients::Vector{Tuple{Float64,Float64}}`: Coefficients for terms n > 1. First component corresponds to sin and second component to cos
- `internal_angle_bias::Float64`: a0 term of the Fourier Series for the angle
- `internal_angle_frequencies::Vector{Float64}`: Frequencies in radians/s
- `internal_angle_coefficients::Vector{Tuple{Float64,Float64}}`: Coefficients for terms n > 1. First component corresponds to sin and second component to cos
- `base_power::Float64`: Base power
- `states::Vector{Symbol}`: State for time, voltage and angle
- `n_states::Int`
- `ext::Dict{String, Any}`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct PeriodicVariableSource <: DynamicInjection
    name::String
    "Source Thevenin resistance"
    R_th::Float64
    "Source Thevenin reactance"
    X_th::Float64
    "a0 term of the Fourier Series for the voltage"
    internal_voltage_bias::Float64
    "Frequencies in radians/s"
    internal_voltage_frequencies::Vector{Float64}
    "Coefficients for terms n > 1. First component corresponds to sin and second component to cos"
    internal_voltage_coefficients::Vector{Tuple{Float64,Float64}}
    "a0 term of the Fourier Series for the angle"
    internal_angle_bias::Float64
    "Frequencies in radians/s"
    internal_angle_frequencies::Vector{Float64}
    "Coefficients for terms n > 1. First component corresponds to sin and second component to cos"
    internal_angle_coefficients::Vector{Tuple{Float64,Float64}}
    "Base power"
    base_power::Float64
    "State for time, voltage and angle"
    states::Vector{Symbol}
    n_states::Int
    ext::Dict{String, Any}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function PeriodicVariableSource(name, R_th, X_th, internal_voltage_bias=0.0, internal_voltage_frequencies=[0.0], internal_voltage_coefficients=[(0.0, 0.0)], internal_angle_bias=0.0, internal_angle_frequencies=[0.0], internal_angle_coefficients=[(0.0, 0.0)], base_power=100.0, ext=Dict{String, Any}(), )
    PeriodicVariableSource(name, R_th, X_th, internal_voltage_bias, internal_voltage_frequencies, internal_voltage_coefficients, internal_angle_bias, internal_angle_frequencies, internal_angle_coefficients, base_power, ext, [:Vt, :θt], 2, InfrastructureSystemsInternal(), )
end

function PeriodicVariableSource(; name, R_th, X_th, internal_voltage_bias=0.0, internal_voltage_frequencies=[0.0], internal_voltage_coefficients=[(0.0, 0.0)], internal_angle_bias=0.0, internal_angle_frequencies=[0.0], internal_angle_coefficients=[(0.0, 0.0)], base_power=100.0, states=[:Vt, :θt], n_states=2, ext=Dict{String, Any}(), internal=InfrastructureSystemsInternal(), )
    PeriodicVariableSource(name, R_th, X_th, internal_voltage_bias, internal_voltage_frequencies, internal_voltage_coefficients, internal_angle_bias, internal_angle_frequencies, internal_angle_coefficients, base_power, states, n_states, ext, internal, )
end

# Constructor for demo purposes; non-functional.
function PeriodicVariableSource(::Nothing)
    PeriodicVariableSource(;
        name="init",
        R_th=0,
        X_th=0,
        internal_voltage_bias=0,
        internal_voltage_frequencies=Any[0],
        internal_voltage_coefficients=[(0.0, 0.0)],
        internal_angle_bias=0,
        internal_angle_frequencies=Any[0],
        internal_angle_coefficients=[(0.0, 0.0)],
        base_power=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`PeriodicVariableSource`](@ref) `name`."""
get_name(value::PeriodicVariableSource) = value.name
"""Get [`PeriodicVariableSource`](@ref) `R_th`."""
get_R_th(value::PeriodicVariableSource) = value.R_th
"""Get [`PeriodicVariableSource`](@ref) `X_th`."""
get_X_th(value::PeriodicVariableSource) = value.X_th
"""Get [`PeriodicVariableSource`](@ref) `internal_voltage_bias`."""
get_internal_voltage_bias(value::PeriodicVariableSource) = value.internal_voltage_bias
"""Get [`PeriodicVariableSource`](@ref) `internal_voltage_frequencies`."""
get_internal_voltage_frequencies(value::PeriodicVariableSource) = value.internal_voltage_frequencies
"""Get [`PeriodicVariableSource`](@ref) `internal_voltage_coefficients`."""
get_internal_voltage_coefficients(value::PeriodicVariableSource) = value.internal_voltage_coefficients
"""Get [`PeriodicVariableSource`](@ref) `internal_angle_bias`."""
get_internal_angle_bias(value::PeriodicVariableSource) = value.internal_angle_bias
"""Get [`PeriodicVariableSource`](@ref) `internal_angle_frequencies`."""
get_internal_angle_frequencies(value::PeriodicVariableSource) = value.internal_angle_frequencies
"""Get [`PeriodicVariableSource`](@ref) `internal_angle_coefficients`."""
get_internal_angle_coefficients(value::PeriodicVariableSource) = value.internal_angle_coefficients
"""Get [`PeriodicVariableSource`](@ref) `base_power`."""
get_base_power(value::PeriodicVariableSource) = value.base_power
"""Get [`PeriodicVariableSource`](@ref) `states`."""
get_states(value::PeriodicVariableSource) = value.states
"""Get [`PeriodicVariableSource`](@ref) `n_states`."""
get_n_states(value::PeriodicVariableSource) = value.n_states
"""Get [`PeriodicVariableSource`](@ref) `ext`."""
get_ext(value::PeriodicVariableSource) = value.ext
"""Get [`PeriodicVariableSource`](@ref) `internal`."""
get_internal(value::PeriodicVariableSource) = value.internal

"""Set [`PeriodicVariableSource`](@ref) `name`."""
set_name!(value::PeriodicVariableSource, val) = value.name = val
"""Set [`PeriodicVariableSource`](@ref) `R_th`."""
set_R_th!(value::PeriodicVariableSource, val) = value.R_th = val
"""Set [`PeriodicVariableSource`](@ref) `X_th`."""
set_X_th!(value::PeriodicVariableSource, val) = value.X_th = val
"""Set [`PeriodicVariableSource`](@ref) `internal_voltage_bias`."""
set_internal_voltage_bias!(value::PeriodicVariableSource, val) = value.internal_voltage_bias = val
"""Set [`PeriodicVariableSource`](@ref) `internal_voltage_frequencies`."""
set_internal_voltage_frequencies!(value::PeriodicVariableSource, val) = value.internal_voltage_frequencies = val
"""Set [`PeriodicVariableSource`](@ref) `internal_voltage_coefficients`."""
set_internal_voltage_coefficients!(value::PeriodicVariableSource, val) = value.internal_voltage_coefficients = val
"""Set [`PeriodicVariableSource`](@ref) `internal_angle_bias`."""
set_internal_angle_bias!(value::PeriodicVariableSource, val) = value.internal_angle_bias = val
"""Set [`PeriodicVariableSource`](@ref) `internal_angle_frequencies`."""
set_internal_angle_frequencies!(value::PeriodicVariableSource, val) = value.internal_angle_frequencies = val
"""Set [`PeriodicVariableSource`](@ref) `internal_angle_coefficients`."""
set_internal_angle_coefficients!(value::PeriodicVariableSource, val) = value.internal_angle_coefficients = val
"""Set [`PeriodicVariableSource`](@ref) `base_power`."""
set_base_power!(value::PeriodicVariableSource, val) = value.base_power = val
"""Set [`PeriodicVariableSource`](@ref) `ext`."""
set_ext!(value::PeriodicVariableSource, val) = value.ext = val


text:
using Distributed
addprocs(Sys.CPU_THREADS)

import PyPlot; const plt = PyPlot
import DSP: conv
import Random: seed!, randperm
import Statistics: std

import KirchMig

dz, dx = 15, 15
x = 0:dx:1000; nx = length(x)
z = 0:dz:1000; nz = length(z)

rho = 1000ones(nz, nx)
rho[div(330,dx):end,:] .+= 1000
rho[div(670,dz):end,:] .-= 500

vel = 2000 .+ zero(rho)

imp = 1e-3rho.*vel;

refl = [zeros(nx)'
    (imp[2:end,:] .- imp[1:end-1,:])./(imp[2:end,:] .+ imp[1:end-1,:])]

refl[div(nz, 2)-1, div(nx, 2)] += 0.2
refl[div(nz, 2)+1, div(nx, 2)] += 0.2;

fig, ax = plt.subplots(1, 3, figsize=(10,2.5))

# Density
cax = ax[1].imshow(rho, extent=[x[1], x[end], z[end], z[1]], aspect="equal", cmap="Blues", interpolation=nothing)
ax[1].set(xlabel="Position [m]", ylabel="Depth [m]")
cbar = fig.colorbar(cax, ax=ax[1])
cbar.ax.set(ylabel="Density [kg/m³]")

# Velocity
cax = ax[2].imshow(vel, extent=[x[1], x[end], z[end], z[1]], aspect="equal", cmap="Blues", interpolation=nothing)
ax[2].set(xlabel="Position [m]", yticks=[])
cbar = fig.colorbar(cax, ax=ax[2])
cbar.ax.set(ylabel="Velocity [m/s]");

cax = ax[3].imshow(imp, extent=[x[1], x[end], z[end], z[1]], aspect="equal", cmap="Blues", interpolation=nothing)
ax[3].set(xlabel="Position [m]", yticks=[])
cbar = fig.colorbar(cax, ax=ax[3])
cbar.ax.set(ylabel="Impedance [kPa·s/m]");
fig.tight_layout()

xran = x[end]-x[1]
mod = reshape(refl, nz, nx)
vmin, vmax = maximum(abs.(mod))*[-1,1]
fig, ax = plt.subplots(figsize=(6,4))
cax = ax.imshow(mod, extent=[x[1], x[end], z[end], z[1]],
    vmin=vmin, vmax=vmax, aspect="equal", cmap="gray", interpolation=nothing)
cbar = fig.colorbar(cax, ax=ax)
cbar.ax.set(ylabel="Amplitude");
ax.plot(xran*mod[:,div(end,2)]./(4vmax) .+ xran./2, z, color="#d62728");
ax.set(xlim=(x[1], x[end]), ylim=(z[end-1], z[1]), xlabel="Position [m]", ylabel="Depth [m]",
    title="Reflectivity")
fig.tight_layout()

nr = 46

rec_z = zeros(nr)
rec_x = range(x[1], stop=x[end], length=nr)

ns = 10
seed!(12)
src_z = zeros(ns)
src_x = sort(rec_x[randperm(nr)][1:ns])

trav_r = KirchMig.eikonal_const_vel([rec_z rec_x], z, x, vel[1]);
trav_s = KirchMig.eikonal_const_vel([src_z src_x], z, x, vel[1]);

fig, ax = plt.subplots(figsize=(6,4))
cax = ax.imshow(trav_s[:,:,div(ns,3)], extent=[x[1], x[end], z[end], z[1]], aspect="equal", cmap="Paired",
    interpolation="bilinear", vmin=0, vmax=maximum(trav_s[:,:,div(ns,3)]))
cbar = fig.colorbar(cax, ax=ax)
cbar.ax.set(ylabel="Traveltime [s]")
cbar.ax.invert_yaxis()
ax.scatter(rec_x[1:3:end], rec_z[1:3:end], color="#1f77b4", marker="v", s=100, clip_on=false, zorder=100)
ax.scatter(src_x, src_z, color="#d62728", alpha=0.4, marker="*", s=200, clip_on=false, zorder=100)
ax.scatter(src_x[div(ns,3)], src_z[div(ns,3)], color="#d62728", marker="*", s=200, clip_on=false, zorder=100)
con = ax.contour(trav_s[:,:,div(ns,3)], extent=[x[1], x[end], z[end], z[1]], origin="upper", colors="k")
plt.clabel(con, fontsize=9, inline=1, fmt="%1.2fs");
ax.set(xlabel="Position [m]", ylabel="Depth [m]", xlim=(x[1], x[end]), ylim=(z[end-1], z[1]))
fig.tight_layout()

t = 0:0.008:1; nt = length(t)
G = KirchMig.KirchMap(t, trav_r, trav_s);

@time Gm = G*refl[:]; # refl is 67×67 a array, refl[:] is a 4489-element vector

ricker(t0, f) = @. (1 - 2pi^2 * f^2 * t0^2) * exp(-pi^2 * f^2 * t0^2)
rick_dtt = ricker(t .- t[div(nt,5)], 15);
@views rick_dtt[2:end-1] = (rick_dtt[1:end-2] - 2.0*rick_dtt[2:end-1] + rick_dtt[3:end])/(t[2] - t[1])^2;
rick_dtt /= maximum(abs.(rick_dtt));

W = KirchMig.ConvMap(rick_dtt, nr, ns, nt);

d = W*Gm;

seed!(1)
n = randn(size(d))
d += n*std(d)/std(n);

data = reshape(d, nr, ns, nt)

rows = 2
fig, ax = plt.subplots(rows, div(ns, rows), figsize=(10,4rows))
idx = reshape(1:rows*div(ns, rows), div(ns, rows), rows)'
for (i, axi) in enumerate(ax)
    csg = Array(data[:,idx[i],:]')
    axi.imshow(csg, extent=[rec_x[1], rec_x[end], t[end], t[1]],
        vmin=-0.5maximum(abs.(csg)), vmax=0.5maximum(abs.(csg)), aspect="auto", cmap="gray", interpolation="none")
    axi.scatter(src_x[idx[i]], t[1], color="#d62728", marker="*", s=200, clip_on=false, zorder=100)
    axi.set(xlim=(rec_x[1],rec_x[end]), ylim=(t[end], t[1]), xlabel="Position [m]",
        xticks = range(minimum(rec_x), stop=maximum(rec_x), length=6))
    if idx[i] != 1 && idx[i] != div(ns,rows)+1
        axi.set(yticks=[])
    else
        axi.set(ylabel="Time [s]")
    end
end
fig.tight_layout()

L = W*G
@time m_mig = L'd;

xran = x[end]-x[1]
mod_mig = reshape(m_mig, nz, nx)
vmin, vmax = maximum(abs.(mod_mig))*[-1,1]
fig, ax = plt.subplots(figsize=(6,4))
cax = ax.imshow(mod_mig, extent=[x[1], x[end], z[end], z[1]],
    vmin=vmin, vmax=vmax, aspect="equal", cmap="gray", interpolation=nothing)
cbar = fig.colorbar(cax, ax=ax)
cbar.ax.set(ylabel="Amplitude");
ax.plot(xran*mod_mig[:,div(end,2)]./(4vmax) .+ xran./2, z, color="#d62728");
ax.set(xlim=(x[1], x[end]), ylim=(z[end], z[1]), xlabel="Position [m]", ylabel="Depth [m]",
    title="Migrated image")
fig.tight_layout()

import IterativeSolvers: cg
@time m_lsm_0, hist_0 = cg(L'L, m_mig, maxiter=20, log=true);

using Arpack
@time λ, ϕ = eigs(L'L; nev=1, maxiter=2); λ = λ[1];

import LinearMaps: UniformScalingMap

Id = UniformScalingMap(1, size(L'L, 1))
@time m_lsm_id, hist_id = cg(L'L + 0.01λ*Id, m_mig, maxiter=20, log=true);

Δ = KirchMig.LaplacianMap(nz, nx);

@time m_lsm_lap, hist_lap = cg(L'L + 0.01λ*Δ'Δ, m_mig, maxiter=20, log=true);

mod_lsm = reshape(m_lsm_0, nz, nx)
mod_lsm_id = reshape(m_lsm_id, nz, nx)
mod_lsm_lap = reshape(m_lsm_lap, nz, nx)

fig, ax = plt.subplots(2, 2, figsize=(6,6))
pclip=1
vmin, vmax = pclip*maximum(abs.(mod_lsm))*[-1,1]
ax[1].imshow(mod_lsm, extent=[x[1], x[end], z[end], z[1]],
    vmin=vmin, vmax=vmax, aspect="equal", cmap="gray", interpolation=nothing)
ax[1].set(xlabel="Position [m]", ylabel="Depth [m]", title="LSM image: no reg.")
ax[1].plot(250mod_lsm[:,50]./vmax .+ 500, z, color="#d62728");

vmin, vmax = pclip*maximum(abs.(mod_lsm_id))*[-1,1]
ax[2].imshow(mod_lsm_id, extent=[x[1], x[end], z[end], z[1]],
    vmin=vmin, vmax=vmax, aspect="equal", cmap="gray", interpolation=nothing)
ax[2].set(xlabel="Position [m]",title="LSM image: min. norm", yticks=[])
ax[2].plot(250mod_lsm_id[:,50]./vmax .+ 500, z, color="#d62728");

vmin, vmax = pclip*maximum(abs.(mod_lsm_lap))*[-1,1]
ax[3].imshow(mod_lsm_lap, extent=[x[1], x[end], z[end], z[1]],
    vmin=vmin, vmax=vmax, aspect="equal", cmap="gray", interpolation=nothing)
ax[3].set(xlabel="Position [m]",title="LSM image: min. curv.", yticks=[])
ax[3].plot(250mod_lsm_lap[:,50]./vmax .+ 500, z, color="#d62728");

vmin, vmax = pclip*maximum(abs.(mod_mig))*[-1,1]
ax[4].imshow(mod_mig, extent=[x[1], x[end], z[end], z[1]],
    vmin=vmin, vmax=vmax, aspect="equal", cmap="gray", interpolation=nothing)
ax[4].set(xlabel="Position [m]",title="Migrated image", yticks=[])
ax[4].plot(250mod_mig[:,50]./vmax .+ 500, z, color="#d62728");
fig.tight_layout()

fig, ax = plt.subplots(figsize=(6,4))
iter = length(hist_0.data[:resnorm])
ax.semilogy(1:iter, hist_0.data[:resnorm], "k", label="No reg.", linewidth=2)

iter = length(hist_id.data[:resnorm])
ax.semilogy(1:iter, hist_id.data[:resnorm], "#1f77b4", label="Min. norm", linewidth=2)

iter = length(hist_lap.data[:resnorm])
ax.semilogy(1:iter, hist_lap.data[:resnorm], "#d62728", label="Min. curv.", linewidth=2)

ax.set(xlabel="Iteration", ylabel="Residual", xlim=(1, iter))
plt.legend(loc=3)
fig.tight_layout()


text:
using LinearAlgebraicRepresentation
using CuboidGrids
using Profile
using BenchmarkTools

using NearestNeighbors
using DataStructures
using IntervalTrees
using Triangle
using SparseArrays
using LinearAlgebra
using Distributed

Lar = CuboidGrids
OlLar = LinearAlgebraicRepresentation

println("=======================================================================")
println("@btime for general execution of the code: ")
    @btime Lar.executeShowVoxels([10,20,10])

println("=======================================================================")
println("@btime for data preparation before visualization(with pmap maethods): ")
    PAR1,PAR2 = @btime Lar.cuboidGrid([10,20,10])
    #PAR1,PAR2 = @btime Lar.cuboidGridT([10,20,10])
println("-----------------------------------------------------------------------")
println("@btime for data preparation before visualization(with map maethods): ")
    PAR1,PAR2 = @btime OlLar.cuboidGrid([10,20,10])

println("=======================================================================")
println("...Adding pmap does not improve the performance in cuboidGrid")
println("To understand why, let's see the single sub functions with the same input")

println("=======================================================================")
println("@btime for larImageVerts(with pmap maethods): ")
    PAR1 = @btime Lar.larImageVerts([10,20,10])
   # PAR1 = @btime Lar.larImageVertsT([10,20,10])
println("-----------------------------------------------------------------------")
println("@btime for larImageVerts(with map maethods): ")
    PAR1 = @btime OlLar.larImageVerts([10,20,10])

println("=======================================================================")
println("@btime for larGridSkeleton(with map maethods): ")
    PAR1 = @btime Lar.larGridSkeleton([10,20,10])
    #PAR1 = @btime Lar.larGridSkeletonT([10,20,10])
println("-----------------------------------------------------------------------")
println("@btime for larGridSkeleton(with map maethods): ")
    PAR1 = @btime OlLar.larGridSkeleton([10,20,10])


text:
module Traces

import ..JSON: json

using ..Utilities, ..TraceUtilities

export Scatter #, Bar, Box, HeatMap, Histogram

# Export only utiliity types that are needed in constructors
# export Font, Title, ColorBar, Gradient, Marker, Line, Error

mutable struct Scatter <: Trace
    type::ChartType
    attributes::Attributes
    stackgroup::String
    groupnorm::GroupNorm
    stackgaps::StackGaps
    mode::Mode
    hoveron::HoverOn
    line::Line
    connectgaps::Bool
    cliponaxis::Bool
    fill::Fill
    fillcolor::Color
    textposition::TextPosition
    textfont::Font
    error_x::Error
    error_y::Error

    function Scatter(;
        stackgroup::String = "",
        groupnorm::String = "",
        stackgaps::String = "infer zero",
        mode::String = "lines+markers",
        hoveron::String = "points",
        line::Line = Line(),
        connectgaps::Bool = false,
        cliponaxis::Bool = true,
        fill::String = "none",
        fillcolor::String = "#444",
        textposition::String = "middle center",
        textfont::Font = Font(),
        error_x::Error = Error(),
        error_y::Error = Error(),kwargs...)
        return new(
            ChartType("scatter"),
            Attributes(;kwargs...),
            stackgroup,
            GroupNorm(groupnorm),
            StackGaps(stackgaps),
            Mode(mode),
            HoverOn(hoveron),
            line,
            connectgaps,
            cliponaxis,
            Fill(fill),
            Color(fillcolor),
            TextPosition(textposition),
            textfont,
            error_x,
            error_y)
    end
end
default["scatter"] = Scatter()

# struct Bar <: Trace
#     type::ChartType
#     attributes::Attributes
#     cliponaxis::Bool
#     textposition::TextPosition
#     textfont::Font
#     insidetextfont::Font
#     outsidetextfont::Font
#     constraintext::ConstrainText
#     base::Float64
#     offset::Union{Float64,Vector{Float64}}
#     width::Union{Size{Float64},Vector{Size{Float64}}}
#     offsetgroup::String
#     alignmentgroup::String
#     error_x::Error
#     error_y::Error

#     function Bar(;
#         cliponaxis::Bool = true,
#         textposition::String = "none",
#         textfont::Font = Font(),

#         insidetextfont::Font = Font(),
#         outsidetextfont::Font = Font(),
#         constraintext::String = "both",
#         base::Float64,
#         offset::Union{Float64,Vector{Float64}},
#         width::Union{Size{Float64},Vector{Size{Float64}}},
#         offsetgroup::String,
#         alignmentgroup::String,
#         error_x::Error,
#         error_y::Error,

#         error_x::Error = Error(),
#         error_y::Error = Error(),kwargs...)
#         return new(
#             ChartType("bar"),
#             Attributes(;kwargs...),
#             cliponaxis,
#             TextPosition(textposition),
#             textfont,
#             )
#     end
# end

# function Base.display(::Figures.Display, scatter::Scatter)
#     println("Try to plot the scatter trace.")
#     s = diff(default["Scatter"],scatter)
#     println(json(s))
#     # id = isempty(Figures.current[]) ? p.plot.divid : Figures.current[]
#     # Pages.broadcast("script","""Plotly.newPlot("$(id)",$(p.plot.data),$(p.plot.layout),$(PlotlyJS.JSON.json(p.options)));""")
#     # println("""Plotly.newPlot("$(id)",$(p.plot.data),$(p.plot.layout),$(PlotlyJS.JSON.json(p.options)));""")
# end

# struct Box <: Trace
#     option::TraceOptions
#     whiskerwidth::RangeValue{Float64,(0,1),(true,true)}
#     notched::Bool
#     notchwidth::RangeValue{Float64,(0,.5),(true,true)}
#     boxpoints::BoxPoints
#     boxmean::BoxMean
#     jitter::RangeValue{Float64,(0,1),(true,true)}
#     pointpos::RangeValue{Float64,(-2,2),(true,true)}
#     hoveron::HoverOn
#     line::Line
#     fillcolor::String
#     width::Union{Size{Float64},Vector{Size{Float64}}}
#     offsetgroup::String
#     alignmentgroup::String
# end

# struct HeatMap <: Trace
#     options::TraceOptions
#     z::Vector{Float64}
#     transpose::Bool
#     xtype::CoordinateType
#     ytype::CoordinateType
#     zsmooth::ZSmooth
#     connectgaps::Bool
#     xgap::Size{Int}
#     ygap::Size{Int}
#     zhoverformat::Union{String,Vector{String}}
#     zauto::Bool
#     zmin::Float64
#     zmax::Float64
#     zmid::Float64
#     colorscale::ColorScale
#     autocolorscale::Bool
#     reversescale::Bool
#     showscale::Bool
#     colorbar::ColorBar
# end

# struct Histogram <: Trace
#     options::TraceOptions
#     histfunc::HistFunc
#     histnorm::Union{HistNorm,Nothing}
#     cumulative::Cumulative
#     nbinsx::Size{Int}
#     xbins::Bins
#     nbinsy::Size{Int}
#     ybins::Bins
#     autobinx::Bool
#     autobiny::Bool
#     offsetgroup::String
#     alignmentgroup::String
#     error_x::Error
#     error_y::Error
# end

end


text:
"""
    asymptotic_monopole_wavenumbers(ω, medium::PhysicalMedium, species; num_wavenumbers = 2)

Estimates a box in complex k_eff space where all needed effective wavenumbers, for the given tolerance, are inside.
""" # Should really use assymptotics for monopole scatterers to properly estimate box size. Would be more robust.
function asymptotic_monopole_wavenumbers(ω::T, medium::PhysicalMedium, species::Species{T};
        num_wavenumbers = 2
    ) where T

    # eff_medium = effective_medium(medium, species)
    # ko = real(ω / eff_medium.c)
    # if isnan(ko) ko = real(ω / medium.c) end

    ko = real(ω / medium.c)

    # After two wavelengths distance the solution would be smaller than tol
    max_imag = - log(tol) * ko / (T(4) * pi)
    max_real = T(4) * ko

   return [[-max_real,max_real],[min_imag, max_imag]]
end

"""
    box_keff(ω::T, medium::PhysicalMedium, species::Species{T};
            tol::T = 1e-5, min_imag::T = zero(T))

Estimates a box in complex k_eff space where all needed effective wavenumbers, for the given tolerance, are inside.
""" # Should really use assymptotics for monopole scatterers to properly estimate box size. Would be more robust.
function box_keff(ω::T, medium::PhysicalMedium, species::Species{T};
        tol::T = 1e-5, min_imag::T = zero(T)
    ) where T

    # eff_medium = effective_medium(medium, species)
    # ko = real(ω / eff_medium.c)
    # if isnan(ko) ko = real(ω / medium.c) end

    ko = real(ω / medium.c)

    # After two wavelengths distance the solution would be smaller than tol
    max_imag = - log(tol) * ko / (T(4) * pi)
    max_real = T(4) * ko

   return [[-max_real,max_real],[min_imag, max_imag]]
end

struct MySimplexer{T<:AbstractFloat} <: Optim.Simplexer
    dx::T
    dy::T
end

function simplexer(A::MySimplexer, initial_x::Array{T}) where T
    initial_simplex = [initial_x, initial_x + [zero(T), A.dy], initial_x + [A.dx, zero(T)]]
end

NelderMeadparameters(;α = 1.0, β = 1.0, γ = 1.0, δ = 0.5)::Optim.FixedParameters = Optim.FixedParameters(
    α = α, # reflection
    β = β, # expansion 1.0 = does not expand
    γ = γ, # contraction
    δ = δ  # shrink step
)

reduce_kvecs(vec::Vector,tol) = vec

function reduce_kvecs(vecs::Vector{Vector{T}},tol::T) where T<:AbstractFloat
    # All wavenumbers should attenuate
    inds = findall([w[2] < -tol for w in vecs])
    vecs[inds] = - vecs[inds]

    all_inds = collect(eachindex(vecs))
    vecs = map(vecs) do w
        ind_ins = findall([norm(v - w) < tol for v in vecs[all_inds]])
        inds = all_inds[ind_ins]
        deleteat!(all_inds,ind_ins)
        isempty(inds) ? [zero(T),-one(T)] :  mean(vecs[inds])
    end
    vecs = deleteat!(vecs, findall(v-> v == [zero(T),-one(T)], vecs) )

    return vecs
end

function reduce_kvecs(ks::Vector{Complex{T}},tol::T) where T<:AbstractFloat
    # All wavenumbers should attenuate
    inds = findall(imag.(ks) .< -tol)
    ks[inds] = - ks[inds]

    all_inds = collect(eachindex(ks))
    ks = map(ks) do w
        ind_ins = findall([norm(v - w) < tol for v in ks[all_inds]])
        inds = all_inds[ind_ins]
        deleteat!(all_inds,ind_ins)
        isempty(inds) ? zero(T)-one(T)*im :  mean(ks[inds])
    end

    ks = deleteat!(ks, findall(ks .== zero(T)-one(T)*im))

    return ks
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
"""
    Node{T}

Computation node with output of type `T`.

Should always be constructed indirectly through `node!`(@ref) or combinator.
"""
struct Node{T}
    # TODO maybe could have a name (for debug)
    nid
end

"""
    eltype(n::Node{T})

Return type `T`.
"""
eltype(n::Node{T}) where {T} = T

"""
    Dag

Directed Acyclic Graph of [`Node{T}`](@ref).
"""
struct Dag
    nid::Ref{Int}  # node id generator (start at 1)
    nodes # nid => Node  # TODO array instead of Dict

    # TODO current typing is weak ... could remove ... or could replace the
    #      link function with Connector{T} type to make sure Node links would
    #      contain only expected connector type
    #      ... maybe have branch without type, and from there see if we can
    #      build new version to reify connector.
    links  # nid => [(nid,fn)]  # TODO array instead of Dict (or store in Node)

    Dag() = new(Ref(0), Dict(), Dict())
end

function _link!(d::Dag, src_nid, dst_nid, fn)
    # TODO check type fn vs src_nid?
    push!(get!(d.links, src_nid, []), (dst_nid, fn))
end

"""
    onfire!(dag, src_nid, fn)

Register a callback `fn` called when new value for node `src_nid` is fired
during evaluation.
"""
function onfire!(d::Dag, src_nid, fn)
    _link!(d, src_nid, 0, fn)
end

"""
    node!(d::Dag, t::Type)

Add new node `Node{t}` to the graph.

The node is given a identifier `nid`.
"""
function node!(d::Dag, t::Type)
    node!(d, t, [])
end

"""
    node!(d::Dag, t::Type, parents)

When given `parents`, create edge from parents to created node.

The `parents` is a sequence of pair of node identifier and function; function
should expect one argument of the parent node output type and should return a
value of type `Nullable{t}`. Note that all the functions should return the same
type.
"""
function node!(d::Dag, t::Type, parents)
    nid = d.nid[] += 1
    n = Node{t}(nid)
    d.nodes[nid] = n
    for (pnode, fn) in parents
        # all fn must have same return type Nullable{t}

        # TODO fn type could be checked through undocumented (not recommended)
        #      function `return_types` but this seems to fail for combine!
        # @assert Base.return_types(fn, (eltype(pnode),)) == [Nullable{t}]
        # TODO can we use functor (link below) to ensure fn(eltype(pnode)) isa Nullable{t}
        # # https://discourse.julialang.org/t/enforcing-function-signatures-by-both-argument-return-types/8174
        # struct Functor{T,R} f end
        # function (ff::Functor{T,R})(x::T)::R where {T,R}
        #     R(ff.f(x))  # not sure if R() is needed
        # end
        # this would raise at call time if f is not ok ... probably inefficient.
        # Another approach is to define new type for each callback and not have
        # f as a field by embedding in type method definition
        # TODO or https://github.com/mauro3/SimpleTraits.jl and enforce user
        # to pass a specific type
        # TODO or https://github.com/yuyichao/FunctionWrappers.jl
        # TODO or just removed Node typing altogether?

        _link!(d, pnode.nid, nid, fn)
    end
    n
end

"""
    return_type_fn1(fn, T)

Return type returned by calling fn with one argument of type T.

fn must be type stable, else AssertionError is raised.

Note: guessing type is not recommended:
      https://groups.google.com/forum/#!topic/julia-users/gb_DR5EzmV4
"""
function return_type_fn1(fn, input_type)
    rets = Base.return_types(fn, (input_type,))
    @assert length(rets) == 1
    return first(rets)
end


text:
"""
    Xtal._is_linearly_independent(vecs::AbstractMatrix{<:Number}) -> Bool
    Xtal._is_linearly_independent(vecs::AbstractVector{<:Number}...) -> Bool

Determines whether a set of vectors is linearly independent.

This function always returns `false` if the first dimension of the matrix is less than the second.
"""
function is_linearly_independent(vecs::AbstractMatrix{<:Number})
    size(vecs)[1] < size(vecs)[2] || return false
    return LinearAlgebra.rank(vecs) == minimum(size(vecs))
end

# Same thing but for sets of vectors
function is_linearly_independent(vecs::AbstractVector{<:Number}...)
    return is_linearly_independent(hcat(vecs...))
end

#=
"""
    _to_vectors(a::AbstractArray{T,N}) where {T,N} -> Vector{...Vector{T}}

Converts an array to a nest of vectors within vectors.
"""
function _to_vectors(a::AbstractArray{T,N}) where {T,N}

end
=#


text:
import Base: ==, -, *, +, append!, convert, copy, delete!, deleteat!,
  firstindex, getindex, hash, in, isempty, isequal, lastindex, length,
  push!, merge, merge!, read, setindex!, show, size, sizeof, sort!,
  sort, summary, write
import SeisIO.Formats: formats, FmtVer, FormatDesc, HistVec
import SeisIO: BUF,
  ChanSpec,
  FloatArray,
  KW,
  NodalLoc,
  TimeSpec,
  checkbuf!,
  checkbuf_strict!,
  code2loctyp,
  code2resptyp,
  code2typ,
  datafields,
  default_fs,
  default_gain,
  do_trace,
  dtchars,
  dtconst,
  fillx_i16_be!,
  flat_resp,
  getbandcode,
  loctyp2code,
  merge_ext!,
  mk_t,
  mkxstr,
  parsetimewin,
  proc_note!,
  prune!,
  read_misc,
  read_string_vec,
  resample!,
  resample,
  resptyp2code,
  showloc_full,
  show_os,
  show_str,
  show_t,
  show_x,
  sμ,
  typ2code,
  write_misc,
  write_string_vec,
  μs


text:
@with_kw type LightDark2DTarget <: AbstractLD2
    min_noise::Float64      = 0.0
    min_noise_loc::Float64  = 5.0
    term_radius::Float64    = 0.05
    init_dist::Any          = SymmetricNormal2([2.0, 2.0], 0.5)
    discount::Float64       = 1.0
    count::Int              = 0
end

reward(p::LightDark2DTarget, s::Vec2) = -1.0
reward(p::LightDark2DTarget, s::Vec2, a::Vec2) = reward(p, s)
reward(p::LightDark2DTarget, s::Vec2, a::Vec2, sp::Vec2) = reward(p, s)


text:
#
# Pagination.
#

TakeQuery(qbase::Query, qn::Query, rev::Bool=false) =
    let q = RecordQuery(qbase, qn)
        q >>
        Query(
            TakeSig(rev),
            Input(domain(output(q))),
            output(qbase) |> setoptional(true))
    end

SkipQuery(qbase::Query, qn::Query) =
    TakeQuery(qbase, qn, true)

immutable TakeSig <: AbstractPrimitive
    rev::Bool
end

ev(sig::TakeSig, dv::DataVector) =
    take_impl(sig.rev, column(dv, 1), column(dv, 2))

function take_impl{OPT,PLU}(rev::Bool, icol::Column{OPT,PLU}, ncol::Column)
    len = length(icol)
    ioffs = offsets(icol)
    ivals = values(icol)
    noffs = offsets(ncol)
    nvals = values(ncol)
    size = 0
    for k = 1:len
        W = ioffs[k+1] - ioffs[k]
        if noffs[k+1] == noffs[k]
            size += W
        else
            w = nvals[noffs[k]]
            size +=
                !rev ?
                    (w >= 0 ? min(w, W) : max(0, W + w)) :
                    (w >= 0 ? max(W - w, 0) : min(W, -w))
        end
    end
    if size == length(ivals)
        return icol
    end
    offs = Vector{Int}(len+1)
    idxs = Vector{Int}(size)
    n = 1
    offs[1] = 1
    for k = 1:len
        L = ioffs[k]
        W = ioffs[k+1] - ioffs[k]
        if noffs[k+1] == noffs[k]
            l = 1
            r = W
        else
            w = nvals[noffs[k]]
            if !rev
                l = 1
                r = w >= 0 ? min(w, W) : max(0, W + w)
            else
                l = w >= 0 ? min(w + 1, W + 1) : max(1, W + w + 1)
                r = W
            end
        end
        for i = (L + l - 1):(L + r - 1)
            idxs[n] = i
            n += 1
        end
        offs[k+1] = n
    end
    return Column{true,PLU}(offs, ivals[idxs])
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
# This file includes digital modulation schemes.

export ASK, PSK, QAM, FSK, PAM, constelsize

abstract type AbstractScheme end

for name in [:PAM, :ASK, :PSK, :QAM, :FSK]
    @eval begin
        struct $name <: AbstractScheme
            "Constellation size"
            M::Int
            function $name(M::Int)
                ispow2(M) || throw(ArgumentError("Expected `M` to be power of 2, got $M instead"))
                new(M)
            end
        end
    end
end

"""
    $SIGNATURES

Returns the alphabet of `scheme`.

# Example
```julia
julia> scheme = PSK(4)
4-PSK

julia> alphabet(scheme)
4-element Array{Array{Float64,1},1}:
 [1.0, 0.0]
 [6.123233995736766e-17, 1.0]
 [-1.0, 1.2246467991473532e-16]
 [-1.8369701987210297e-16, -1.0]
```
"""
function alphabet end

"""
    $SIGNATURES

Returns ``[2m - 1 - M], \\; m = 1, \\ldots, M`` where `M` is constellation size of `scheme`.
"""
function alphabet(scheme::Union{PAM, ASK})
    M = scheme.M
    [[2m - 1 - M + 0im] for m in 1 : M]
end

"""
    $SIGNATURES

Returns ``[cos(\\theta_m), sin(\\theta_m)], \\; m = 1, \\ldots, M`` where ``\\theta_m = \\dfrac{2\\pi(m - 1)}{M}`` and
`M` is the constellation size of the `scheme`.
"""
function alphabet(scheme::PSK)
    M = scheme.M
    [[cos(θ) + 1im * sin(θ)] for θ in 2π * (0 : M - 1) / M]
end

"""
    $SIGNATURES

Returns ``[2m - 1 - M, 2n - 1 - M], \\; m, n = 1, \\ldots, M`` where `M` is the constellation size of the `scheme`.
"""
function alphabet(scheme::QAM)
    M = Int(sqrt(scheme.M))
    Mrange = -(M - 1) : 2 : (M - 1)
    vec([[i + 1im * j] for i in Mrange, j in Mrange])
end

"""
    $SIGNATURES

Returns ``[0, 0, \\ldots, 1, 0, \\ldots, 0], \\; m = 1, \\ldots, M`` where the index of the nonzero element is ``m`` and
`M` is the constellation size of the `scheme`.
"""
function alphabet(scheme::FSK)
    M = scheme.M
    map(i -> setindex!(zeros(M), 1, i), 1 : M)
end

show(io::IO, scheme::T) where T <: AbstractScheme = print(io, "$(scheme.M)-$(T.name)")

"""
    $SIGNATURES

Returns the symbols size of `scheme`.

# Example
```julia
julia> sch = PSK(4);

julia> constelsize(sch)
4
```
"""
constelsize(scheme::AbstractScheme) = scheme.M

"""
    $SIGNATURES

Plots the constellation diagram of the `modulator`.
"""
function constellation(scheme::T) where {T<:AbstractScheme}
    s = alphabet(scheme)
    if T <: ASK || T <: PAM
        ymax = 1
        plt = scatter(vcat(s...), zeros(length(s)), ylims=(-ymax, ymax))
        foreach(item -> annotate!(item[2][1],  0.1, item[1]), enumerate(s))
    elseif T <: PSK || T <: QAM
        ymax = maximum(norm.(s)) * 1.25
        plt = scatter(getindex.(s, 1), getindex.(s, 2), marker=:circle, ylims=(-ymax, ymax))
        foreach(item -> annotate!(item[2][1] * 0.9, item[2][2] * 0.9, item[1]), enumerate(s))
    else
        error("Unknown modulation scheme. Expected `PAM, ASK, PSK, QAM`, got $T instead.")
    end
    plt
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -3.357520876571253 <= q <= 3.7351340364874854)

                     add_NL_constraint(m, :(log(1 + exp(0.618782057809276 + -0.4828154761503023*log(1 + exp(0.8842828756523886 + -0.6089179668773923*log(1 + exp(0.019018219753813348 + 0.8244023518746615*log(1 + exp(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2]))) + 0.3296833580573324*log(1 + exp(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))))) + 0.44138625638024864*log(1 + exp(-0.9126485111761089 + 0.9176394232569307*log(1 + exp(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2]))) + 0.371990627187599*log(1 + exp(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))))))) + 0.3197591788003944*log(1 + exp(0.7251254505916935 + -0.923704860259535*log(1 + exp(0.019018219753813348 + 0.8244023518746615*log(1 + exp(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2]))) + 0.3296833580573324*log(1 + exp(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))))) + 0.28767280168552967*log(1 + exp(-0.9126485111761089 + 0.9176394232569307*log(1 + exp(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2]))) + 0.371990627187599*log(1 + exp(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))))))))) + log(1 + exp(-0.0992210488643468 + 0.772244110024173*log(1 + exp(0.8842828756523886 + -0.6089179668773923*log(1 + exp(0.019018219753813348 + 0.8244023518746615*log(1 + exp(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2]))) + 0.3296833580573324*log(1 + exp(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))))) + 0.44138625638024864*log(1 + exp(-0.9126485111761089 + 0.9176394232569307*log(1 + exp(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2]))) + 0.371990627187599*log(1 + exp(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))))))) + 0.07492471407403789*log(1 + exp(0.7251254505916935 + -0.923704860259535*log(1 + exp(0.019018219753813348 + 0.8244023518746615*log(1 + exp(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2]))) + 0.3296833580573324*log(1 + exp(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))))) + 0.28767280168552967*log(1 + exp(-0.9126485111761089 + 0.9176394232569307*log(1 + exp(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2]))) + 0.371990627187599*log(1 + exp(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# To make composition possible, we implement [`compose`](#), which
# defaults to returning a [`Sequence`](#).

"""
    Sequence(transforms...)

`Transform` that applies multiple `transformations`
after each other.

You should not use this explicitly. Instead use [`compose`](#).
"""
struct Sequence{T<:Tuple} <: Transform
    transforms::T
end

Sequence(tfms...) = Sequence{typeof(tfms)}(tfms)
Sequence(tfm::Transform) = tfm

getrandstate(seq::Sequence) = getrandstate.(seq.transforms)

compose(tfm1::Transform, tfm2::Transform) = Sequence(tfm1, tfm2)
compose(seq::Sequence, tfm::Transform) = Sequence(seq.transforms..., tfm)
compose(tfm::Transform, seq::Sequence) = compose(tfm, seq.transforms...)
compose(::Identity, seq::Sequence) = seq
compose(seq::Sequence, ::Identity) = seq
compose(seq1::Sequence, seq2::Sequence) = compose(seq1.transforms..., seq2.transforms...)

function apply(seq::Sequence, items::Tuple; randstate = getrandstate(seq))
    for (tfm, r) in zip(seq.transforms, randstate)
        items = apply(tfm, items; randstate = r)
    end
    return items
end

apply(seq::Sequence, item::Item; randstate = getrandstate(seq)) =
    apply(seq, (item,); randstate = randstate) |> only

function makebuffer(pipeline::Sequence, items)
    buffers = []
    for tfm in pipeline.transforms
        push!(buffers, makebuffer(tfm, items))
        items = apply(tfm, items)
    end
    return buffers
end

function apply!(buffers::Vector, pipeline::Sequence, items::Item; randstate = getrandstate(pipeline))
    @assert length(buffers) == length(pipeline.transforms)
    for (tfm, buffer, r) in zip(pipeline.transforms, buffers, randstate)
        items = apply!(buffer, tfm, items; randstate = r)
    end
    return items
end


text:
"""
# TauPy

Use the Python package ObsPy to calculate seismic travel times and ray paths for
teleseismic phases.

## Phases defined only by event depth and distance

For this case, functions return vectors of `Phase`s, which only contain details of each
phase in terms of distance from the event and radius.  Two methods are available:

## Phases defined by geographic coordinates

Here, functions give vectors of `PhaseGeog`s, which include information about the
geographic location within the Earth of the event, station, and potentially ray
path between them.

## Functions

In the case of either 1D or geographical phases, two functions are available:

- `travel_time` returns `Phase` objects which contain only time and slowness information.
- `path` returns `Phase`s which also contain the ray path.

Each have two methods, one where only depth and distance are specified, and the other
for when one wants also to include geographic information.

There is also the convenience function `turning_depth`, which gives the turning
depth of a `Phase` or `PhaseGeog` when the path has been calculated, and throws
an error when it has not.
"""
module TauPy

export
    AbstractPhase,
    Phase,
    PhaseGeog,
    path,
    travel_time,
    turning_depth

using PyCall

function __init__()
    # This is just to make sure geographiclib is installed for geographic ray paths
    pyimport_conda("geographiclib", "geographiclib", "conda-forge")
    copy!(Taup, pyimport_conda("obspy.taup", "obspy", "conda-forge"))
    for m in AVAILABLE_MODELS
        global MODEL[m] = Taup.TauPyModel(m)
        global RADIUS[m] = MODEL[m].model.radius_of_planet
    end
end

const Taup = PyNULL()

const TauPyFloat = Float64
"All available models for us in TauPy"
const AVAILABLE_MODELS = ("1066a", "1066b", "ak135", # "ak135f", # FIXME: Doesn't work
                          "ak135f_no_mud", "herrin",
                          "iasp91", "jb", "prem", "pwdk", "sp6")
"Default model used by TauPy"
const DEFAULT_MODEL = "ak135"
"Dictionary of pre-instantiated models, where keys are the model names"
const MODEL = Dict{String,PyObject}()
"Dictionary of the radius (km) of models"
const RADIUS = Dict{String,TauPyFloat}()

"""
    AbstractPhase

Abstract type from which `Phase` and `PhaseGeog` subtype.  Methods which may take
either `Phase`s or `PhaseGeog`s should specify their arguments as `AbstractPhase`s.
"""
abstract type AbstractPhase end

function Base.:(==)(p1::AbstractPhase, p2::AbstractPhase)
    p1 === p2 && return true
    typeof(p1) == typeof(p2) || return false
    for f in fieldnames(typeof(p1))
        if getfield(p1, f) isa AbstractVector
            length(getfield(p1, f)) == length(getfield(p2, f)) || return false
            all((a,b) -> a==b, zip(getfield(p1, f), getfield(p2, f))) ||
                return false
        else
            getfield(p1, f) == getfield(p2, f) || return false
        end
    end
    true
end

include("phase.jl")
include("phase_geog.jl")
include("memo.jl")

"""
    available_models()

Return the models which are available for use in TauPy.
"""
available_models() = AVAILABLE_MODELS

"""
    turning_depth(p::AbstractPhase) -> depth

Return the turning `depth` (km) for a phase which has had its ray path calculated.
Throw an error if the path has not been calculated.
"""
turning_depth(p::AbstractPhase) = length(p.radius) > 0 ? RADIUS[p.model] - minimum(p.radius) :
    error("No path calculated for phase")

"""
    _call_taup(taup_function, model, args...; cache=true) -> arrivals

Generic routine to call the Obspy taup module and process the arguments into a vector
of `Phase`s or `PhaseGeog`s.

`taup_function` is a `Symbol` giving the name of the obspy.taup function to call.
"""
function _call_taup(func::Symbol, model, args...; cache=true)
    key = cache_args_key(func, model, args...)
    arr = if cache && haskey(RAY_CACHE, key)
        get_cache(key)
    else
        arr = _phases_from_arrivals(getproperty(MODEL[model], func)(args...), model, args...)
        cache && update_cache!(key, arr)
        arr
    end
end

end # module

max_stars_repo_name: anowacki/TauPy


text:
struct M1 <: MCMC.Model end
MCMC.make_init_state(m::M1) = (mu=0.0, sigma=0)

struct M2 <: MCMC.Model end
MCMC.make_init_state(m::M2) = (mu=0.0, sigma=0, eta=zeros(3))

struct M3 <: MCMC.Model end
MCMC.make_init_state(m::M3) = (mu=0.0, sigma=0.0, eta=zeros(3))

struct M4 <: MCMC.Model end
MCMC.make_init_state(m::M4) = (theta=0.1, mu=.2)

struct M5 <: MCMC.Model K::Int end
MCMC.make_init_state(m::M5) = (theta=zeros(m.K), mu=.2)


text:
function jt09()
#
#  test polynomial suggested by Jenkins and Traub
#
    p1 = [1;zeros(9);-1e-20]
    p2 = [1;zeros(9);1e20]
    p = conv(p1,p2)
    z1 = [-0.01000000000000;
          -0.00809016994375 + 0.00587785252292im;
          -0.00809016994375 - 0.00587785252292im;
          -0.00309016994375 + 0.00951056516295im;
          -0.00309016994375 - 0.00951056516295im;
          0.00309016994375 + 0.00951056516295im;
          0.00309016994375 - 0.00951056516295im;
          0.01000000000000;
          0.00809016994375 + 0.00587785252292im;
          0.00809016994375 - 0.00587785252292im;
         ]

    z2 = [-0.95105651629515 + 0.30901699437495im;
          -0.95105651629515 - 0.30901699437495im;
          -0.58778525229247 + 0.80901699437495im;
          -0.58778525229247 - 0.80901699437495im;
          -0.00000000000000 + 1.00000000000000im;
          -0.00000000000000 - 1.00000000000000im;
          0.95105651629515 + 0.30901699437495im;
          0.95105651629515 - 0.30901699437495im;
          0.58778525229247 + 0.80901699437495im;
          0.58778525229247 - 0.80901699437495im;
          ] * 100
    z = [[z1; z2] ones(20)]
    p, PolyZeros(z)
end


text:
using DemeNet: Notary, Cypher, DemeSpec, Deme, Signer, Certificate, datadir, save, DemeID, ID

using PeaceVote.BraidChains.Types: Proposal, Vote
using PeaceVote.BraidChains.Ledgers: Ledger, load
using PeaceVote.BraidChains.Braiders: BraiderConfig, Braider, Mixer, braid!
using PeaceVote.BraidChains.Analysis: members, proposals, attest, voters
using PeaceVote.BraidChains.Recorders: RecorderConfig, Recorder, record

using PeaceCypher

for dir in [homedir() * "/.demenet/"]
    isdir(dir) && rm(dir,recursive=true)
end

### Configuration ###

demespec = DemeSpec("PeaceDeme",:default,:PeaceCypher,:default,:PeaceCypher,:PeaceVote)
save(demespec) ### Necessary to connect with Mixer
uuid = demespec.uuid
deme = Deme(demespec)

maintainer = Signer(deme,"maintainer")

# Somewhere far far away
mixer = Signer(deme,"mixer")
mixerserver = Mixer(2001,deme,mixer)

server = Signer(deme,"server")

MIXER_ID = mixer.id
SERVER_ID = server.id

braiderconfig = BraiderConfig(2000,2001,UInt8(3),UInt8(64),SERVER_ID,DemeID(uuid,MIXER_ID))
recorderconfig = RecorderConfig([maintainer.id,],server.id,2002,2003,2004)

braider = Braider(braiderconfig,deme,server)

ledger = Ledger(uuid)
recorder = Recorder(recorderconfig,deme,ledger,braider,server)

### Testing the server ###

for i in 1:3
    account = "account$i"
    member = Signer(deme,account * "/member")
    cert = Certificate(member.id,maintainer)
    @show record(recorderconfig,cert)
end

pmember = Signer(deme,"account2" * "/member")
proposal = Proposal("Found peace for a change?",["yes","no","maybe"])
cert = Certificate(proposal,pmember)

record(recorderconfig,cert);

@sync for i in 1:3
    @async begin
        account = "account$i"
        member = Signer(deme,account * "/member")
        voter = Signer(deme,account * "/voters/$(string(member.id))")
        braid!(braiderconfig,deme,voter,member)
    end
end

### Now I can work on reading in parsing the ledger to a BraidChain

sleep(1)

loadedledger = load(ledger)
messages = attest(loadedledger,deme.notary)
index = proposals(messages)[1]
proposal = messages[index]

for i in 1:3
    account = "account$i"

    member = Signer(deme,account * "/member")
    voter = Signer(deme,account * "/voters/$(string(member.id))")

    option = Vote(index,rand(1:length(proposal.document.options)))
    cert = Certificate(option,voter)

    record(recorderconfig,cert)
end


text:
using SnoopCompile

# using runtests:
@snoopiBot "Juno"


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Container data structure (structure-of-arrays style) for DG elements on curved unstructured mesh
struct UnstructuredElementContainer2D{RealT<:Real, uEltype<:Real}
  node_coordinates     ::Array{RealT, 4}   # [ndims, nnodes, nnodes, nelement]
  jacobian_matrix      ::Array{RealT, 5}   # [ndims, ndims, nnodes, nnodes, nelement]
  inverse_jacobian     ::Array{RealT, 3}   # [nnodes, nnodes, nelement]
  contravariant_vectors::Array{RealT, 5}   # [ndims, ndims, nnodes, nnodes, nelement]
  normal_directions    ::Array{RealT, 4}   # [ndims, nnodes, local sides, nelement]
  surface_flux_values  ::Array{uEltype, 4} # [variables, nnodes, local sides, elements]
end

# construct an empty curved element container to be filled later with geometries in the
# unstructured mesh constructor
function UnstructuredElementContainer2D{RealT, uEltype}(capacity::Integer, n_variables, n_nodes) where {RealT<:Real, uEltype<:Real}
  nan_RealT = convert(RealT, NaN)
  nan_uEltype = convert(uEltype, NaN)

  node_coordinates      = fill(nan_RealT, (2, n_nodes, n_nodes, capacity))
  jacobian_matrix       = fill(nan_RealT, (2, 2, n_nodes, n_nodes, capacity))
  inverse_jacobian      = fill(nan_RealT, (n_nodes, n_nodes, capacity))
  contravariant_vectors = fill(nan_RealT, (2, 2, n_nodes, n_nodes, capacity))
  normal_directions     = fill(nan_RealT, (2, n_nodes, 4, capacity))
  surface_flux_values   = fill(nan_uEltype, (n_variables, n_nodes, 4, capacity))

  return UnstructuredElementContainer2D{RealT, uEltype}(node_coordinates,
                                                        jacobian_matrix,
                                                        inverse_jacobian,
                                                        contravariant_vectors,
                                                        normal_directions,
                                                        surface_flux_values)
end

@inline nelements(elements::UnstructuredElementContainer2D) = size(elements.surface_flux_values, 4)
@inline eachelement(elements::UnstructuredElementContainer2D) = Base.OneTo(nelements(elements))

@inline nvariables(elements::UnstructuredElementContainer2D) = size(elements.surface_flux_values, 1)
@inline nnodes(elements::UnstructuredElementContainer2D) = size(elements.surface_flux_values, 2)

Base.real(elements::UnstructuredElementContainer2D) = eltype(elements.node_coordinates)
Base.eltype(elements::UnstructuredElementContainer2D) = eltype(elements.surface_flux_values)

@inline function get_surface_normal(vec, indices...)
  # way to extract the normal vector at the surfaces without allocating
  surface_vector = SVector(ntuple(j -> vec[j, indices...], 2))
  return surface_vector
end

function init_elements(mesh::UnstructuredMesh2D, equations, basis, RealT, uEltype)
  elements = UnstructuredElementContainer2D{RealT, uEltype}(
    mesh.n_elements, nvariables(equations), nnodes(basis))
  init_elements!(elements, mesh, basis)
  return elements
end

function init_elements!(elements::UnstructuredElementContainer2D, mesh, basis)
  four_corners = zeros(eltype(mesh.corners), 4, 2)

  # loop through elements and call the correct constructor based on whether the element is curved
  for element in eachelement(elements)
    if mesh.element_is_curved[element]
      init_element!(elements, element, basis.nodes, view(mesh.surface_curves, :, element))
    else # straight sided element
      for i in 1:4, j in 1:2
        # pull the (x,y) values of these corners out of the global corners array
        four_corners[i, j] = mesh.corners[j, mesh.element_node_ids[i, element]]
      end
      init_element!(elements, element, basis.nodes, four_corners)
    end
  end
end

# initialize all the values in the container of a general element (either straight sided or curved)
function init_element!(elements, element, nodes, corners_or_surface_curves)

  calc_node_coordinates!(elements.node_coordinates, element, nodes, corners_or_surface_curves)

  calc_metric_terms!(elements.jacobian_matrix, element, nodes, corners_or_surface_curves)

  calc_inverse_jacobian!(elements.inverse_jacobian, element, elements.jacobian_matrix)

  calc_contravariant_vectors!(elements.contravariant_vectors, element, elements.jacobian_matrix)

  calc_normal_directions!(elements.normal_directions, element, nodes, corners_or_surface_curves)

  return elements
end

# generic container for the interior interfaces of an unstructured mesh
struct UnstructuredInterfaceContainer2D{uEltype<:Real}
  u                ::Array{uEltype, 4} # [primary/secondary, variables, i, interfaces]
  start_index      ::Vector{Int}       # [interfaces]
  index_increment  ::Vector{Int}       # [interfaces]
  element_ids      ::Array{Int, 2}     # [primary/secondary, interfaces]
  element_side_ids ::Array{Int, 2}     # [primary/secondary, interfaces]
end

# Construct an empty curved interface container to be filled later with neighbour
# information in the unstructured mesh constructor
function UnstructuredInterfaceContainer2D{uEltype}(capacity::Integer, n_variables, n_nodes) where {uEltype<:Real}

  nan_uEltype = convert(uEltype, NaN)

  u                = fill(nan_uEltype, (2, n_variables, n_nodes, capacity))
  start_index      = fill(typemin(Int), capacity)
  index_increment  = fill(typemin(Int), capacity)
  element_ids      = fill(typemin(Int), (2, capacity))
  element_side_ids = fill(typemin(Int), (2, capacity))

  return UnstructuredInterfaceContainer2D{uEltype}(
    u, start_index, index_increment, element_ids, element_side_ids)
end

@inline ninterfaces(interfaces::UnstructuredInterfaceContainer2D) = length(interfaces.start_index)
@inline nnodes(interfaces::UnstructuredInterfaceContainer2D) = size(interfaces.u, 3)

function init_interfaces(mesh::UnstructuredMesh2D, elements::UnstructuredElementContainer2D)

  interfaces = UnstructuredInterfaceContainer2D{eltype(elements)}(
    mesh.n_interfaces, nvariables(elements), nnodes(elements))

  # extract and save the appropriate neighbour information from the mesh skeleton
  init_interfaces!(interfaces, mesh.neighbour_information, mesh.boundary_names,
                   mesh.n_elements, Val(isperiodic(mesh)))

  return interfaces
end

function init_interfaces!(interfaces, edge_information, boundary_names, n_elements,
                          periodic::Val{false})

  n_nodes = nnodes(interfaces)
  n_surfaces = size(edge_information, 2)
  intr_count = 1
  for j in 1:n_surfaces
    if edge_information[4,j] > 0
      # get the primary/secondary element information and coupling for an interior interface
      interfaces.element_ids[1,intr_count]      = edge_information[3,j]      # primary element id
      interfaces.element_ids[2,intr_count]      = edge_information[4,j]      # secondary element id
      interfaces.element_side_ids[1,intr_count] = edge_information[5,j]      # primary side id
      interfaces.element_side_ids[2,intr_count] = abs(edge_information[6,j]) # secondary side id
      # default the start and increment indexing
      interfaces.start_index[intr_count] = 1
      interfaces.index_increment[intr_count] = 1
      if edge_information[6,j] < 0
        # coordinate system in the secondary element is "flipped" compared to the primary element.
        # Adjust the start and increment indexes such that the secondary element coordinate system
        # can match the primary neighbour when surface coupling is computed
        interfaces.start_index[intr_count] = n_nodes
        interfaces.index_increment[intr_count] = -1
      end
      intr_count += 1
    end
  end

  return nothing
end

function init_interfaces!(interfaces, edge_information, boundary_names, n_elements,
                          periodic::Val{true})

  n_nodes = nnodes(interfaces)
  n_surfaces = size(edge_information, 2)
  # for now this set a fully periodic domain
  #   TODO: possibly adjust to be able to set periodic in only the x or y direction
  for j in 1:n_surfaces
    if edge_information[4,j] > 0
      # get the primary/secondary element information and coupling for an interior interface
      interfaces.element_ids[1,j]      = edge_information[3,j]      # primary element id
      interfaces.element_ids[2,j]      = edge_information[4,j]      # secondary element id
      interfaces.element_side_ids[1,j] = edge_information[5,j]      # primary side id
      interfaces.element_side_ids[2,j] = abs(edge_information[6,j]) # secondary side id
      # default the start and increment indexing
      interfaces.start_index[j] = 1
      interfaces.index_increment[j] = 1
      if edge_information[6,j] < 0
        # coordinate system in the secondary element is "flipped" compared to the primary element.
        # Adjust the start and increment indexes such that the secondary element coordinate system
        # can match the primary neighbour when surface coupling is computed
        interfaces.start_index[j] = n_nodes
        interfaces.index_increment[j] = -1
      end
    else
      # way to set periodic BCs where we are assuming to have a structured mesh with internal curves
      primary_side = edge_information[5,j]
      primary_element = edge_information[3,j]
      # Note: This is a way to get the neighbour element number and local side from a square
      #       structured mesh where the element local surface numbering is right-handed
      if boundary_names[primary_side, primary_element] === :Bottom
        secondary_element = primary_element + (n_elements - convert(Int, sqrt(n_elements)))
        secondary_side    = 3
      elseif boundary_names[primary_side, primary_element] === :Top
        secondary_element = primary_element - (n_elements - convert(Int, sqrt(n_elements)))
        secondary_side    = 1
      elseif boundary_names[primary_side, primary_element] === :Left
        secondary_element = primary_element + (convert(Int, sqrt(n_elements)) - 1)
        secondary_side    = 2
      elseif boundary_names[primary_side, primary_element] === :Right
        secondary_element = primary_element - (convert(Int, sqrt(n_elements)) - 1)
        secondary_side    = 4
      end
      interfaces.element_ids[1,j]      = primary_element
      interfaces.element_ids[2,j]      = secondary_element
      interfaces.element_side_ids[1,j] = primary_side
      interfaces.element_side_ids[2,j] = secondary_side
      # set the start and increment indexing
      #  Note! We assume that the periodic mesh has no flipped element coordinate systems
      interfaces.start_index[j] = 1
      interfaces.index_increment[j] = 1
    end
  end

  return nothing
end

# TODO: Clean-up meshes. Find a better name since it's also used for other meshes
# generic container for the boundary interfaces of an unstructured mesh
struct UnstructuredBoundaryContainer2D{RealT<:Real, uEltype<:Real}
  u               ::Array{uEltype, 3} # [variables, i, boundaries]
  element_id      ::Vector{Int}       # [boundaries]
  element_side_id ::Vector{Int}       # [boundaries]
  node_coordinates::Array{RealT, 3}   # [ndims, nnodes, boundaries]
  name            ::Vector{Symbol}    # [boundaries]
end

# construct an empty curved boundary container to be filled later with neighbour
# information in the unstructured mesh constructor
function UnstructuredBoundaryContainer2D{RealT, uEltype}(capacity::Integer, n_variables, n_nodes) where {RealT<:Real, uEltype<:Real}

  nan_RealT = convert(RealT, NaN)
  nan_uEltype = convert(uEltype, NaN)

  u                = fill(nan_uEltype, (n_variables, n_nodes, capacity))
  element_id       = fill(typemin(Int), capacity)
  element_side_id  = fill(typemin(Int), capacity)
  node_coordinates = fill(nan_RealT, (2, n_nodes, capacity))
  name             = fill(:empty, capacity)

  return UnstructuredBoundaryContainer2D{RealT, uEltype}(
    u, element_id, element_side_id, node_coordinates, name)
end

@inline nboundaries(boundaries::UnstructuredBoundaryContainer2D) = length(boundaries.name)

function init_boundaries(mesh::UnstructuredMesh2D, elements::UnstructuredElementContainer2D)

  boundaries = UnstructuredBoundaryContainer2D{real(elements), eltype(elements)}(
    mesh.n_boundaries, nvariables(elements), nnodes(elements))

  # extract and save the appropriate boundary information provided any physical boundaries exist
  if mesh.n_boundaries > 0
    init_boundaries!(boundaries, mesh.neighbour_information, mesh.boundary_names, elements)
  end
  return boundaries
end

function init_boundaries!(boundaries::UnstructuredBoundaryContainer2D, edge_information,
                          boundary_names, elements)

  n_surfaces = size(edge_information,2)
  bndy_count = 1
  for j in 1:n_surfaces
    if edge_information[4,j] == 0
      # get the primary element information at a boundary interface
      primary_element = edge_information[3,j]
      primary_side    = edge_information[5,j]
      boundaries.element_id[bndy_count]      = primary_element
      boundaries.element_side_id[bndy_count] = primary_side

      # extract the physical boundary's name from the global list
      boundaries.name[bndy_count] = boundary_names[primary_side, primary_element]

      # Store copy of the (x,y) node coordinates on the physical boundary
      enc = elements.node_coordinates
      if primary_side == 1
        boundaries.node_coordinates[:, :, bndy_count] .= enc[:, :,    1, primary_element]
      elseif primary_side == 2
        boundaries.node_coordinates[:, :, bndy_count] .= enc[:, end,  :, primary_element]
      elseif primary_side == 3
        boundaries.node_coordinates[:, :, bndy_count] .= enc[:, :,  end, primary_element]
      else # primary_side == 4
        boundaries.node_coordinates[:, :, bndy_count] .= enc[:, 1,    :, primary_element]
      end
      bndy_count += 1
    end
  end

  return nothing
end

end # @muladd


text:
include("Change.jl")

mutable struct Variable

	name::String

	lowerBound::Vector{Int}
	upperBound::Vector{Int}

	cardinalInf::Int
	cardinalSup::Int

	linkedConstraint::Set{Constraint}

	isFixed::Bool# <=> upperBound == lowerBound

end

function isFixed(var::Variable)
	return (var.cardinalInf == length(var.upperBound) || var.cardinalSup == length(var.lowerBound))
end

function fix!(var::Variable, change::Change)

	if (var.cardinalInf == length(var.upperBound))
		add = setdiff(var.upperBound, var.lowerBound)

		for elt in add
			push!(change.added, elt)
		end
		union!(var.lowerBound, add)
		change.cardRemoved += var.cardinalInf - var.cardinalSup
		var.cardinalSup = var.cardinalInf
	else
		rem = setdiff(var.upperBound, var.lowerBound)
		for elt in rem
			push!(change.removed, elt)
		end
		setdiff!(var.upperBound, rem)
		change.cardAdded += var.cardinalSup - var.cardinalInf
		var.cardinalInf = var.cardinalSup
	end

	#println("On fixe $var")

	var.isFixed = true
	change.fixed = true

	var
end

function forced!(var::Variable, value::Union{Int, Nothing}, model)

	@assert var.isFixed == false "Aled"

	change = Change(var)

	if value == nothing
		cardSup = var.cardinalInf - var.cardinalSup
		var.cardinalSup += cardSup
		change.cardRemoved += cardSup

		rem = setdiff(var.upperBound, var.lowerBound)
		for elt in rem
			push!(change.removed, elt)
		end
		setdiff!(var.upperBound, rem)

		change.fixed = true
		var.isFixed = true
		filter!(var-> !(var.isFixed), model.varsNotFixed)

	else

		union!(var.lowerBound, [value])
		push!(change.added, value)
		cardInf = max(length(var.lowerBound), var.cardinalInf)

		change.cardAdded += cardInf - var.cardinalInf
		var.cardinalInf += change.cardAdded

		if isFixed(var)
			fix!(var, change)
			filter!(var-> !(var.isFixed), model.varsNotFixed)
		end
	end
	return change
end

function unforced!(var::Variable, change::Change, model)

	if change.fixed
		var.isFixed = false
		union!(model.varsNotFixed, [var])
	end

	union!(var.upperBound, change.removed)
	setdiff!(var.lowerBound, change.added)

	var.cardinalInf -= change.cardAdded
	var.cardinalSup -= change.cardRemoved

	var

end

import Base.show
function Base.show(io::IO, var::Variable)
	print(io, var.name)
end

function Variable(model, domain, name::String = "X")
	var = Variable(domain, name)
	push!(model.varsInter, var)
	push!(model.varsNotFixed, var)

	var
end

function Variable(domain::UnitRange{Int}, name::String = "X")
	return Variable(collect(domain), name)
end

function Variable(domain::Vector{Int}, name::String = "X")

	lb = Vector{Int}()
	ub = domain[:]

	cardinalInf = 0
	cardinalSup = length(domain)

	linkedConstraint = Set{Constraint}()

	return Variable(name,
					lb,
					ub,
					cardinalInf,
					cardinalSup,
					linkedConstraint,
					false)

end


text:
"""
Provides the [`walk`](@ref) function.
"""
module Walkers

import ..Documenter:
    Anchors,
    Builder,
    Documents,
    Expanders,
    Formats,
    Utilities

using Compat, DocStringExtensions

"""
$(SIGNATURES)

Calls `f` on `element` and any of its child elements. `meta` is a `Dict` containing metadata
such as current module.
"""
walk(f, meta, element) = (f(element); nothing)

# Change to the docstring's defining module if it has one. Change back afterwards.
function walk(f, meta, block::Markdown.MD)
    tmp = get(meta, :CurrentModule, nothing)
    mod = get(block.meta, :module, nothing)
    mod ≡ nothing || (meta[:CurrentModule] = mod)
    f(block) && walk(f, meta, block.content)
    tmp ≡ nothing ? delete!(meta, :CurrentModule) : (meta[:CurrentModule] = tmp)
    nothing
end

function walk(f, meta, block::Vector)
    for each in block
        walk(f, meta, each)
    end
end

const MDContentElements = Union{
    Markdown.BlockQuote,
    Markdown.Paragraph,
    Markdown.MD,
}
walk(f, meta, block::MDContentElements) = f(block) ? walk(f, meta, block.content) : nothing

walk(f, meta, block::Anchors.Anchor) = walk(f, meta, block.object)

walk(f, meta, block::Expanders.DocsNodes) = walk(f, meta, block.nodes)
walk(f, meta, block::Expanders.DocsNode)  = walk(f, meta, block.docstr)

walk(f, meta, block::Expanders.EvalNode) = walk(f, meta, block.result)

walk(f, meta, block::Documents.RawHTML) = nothing

walk(f, meta, block::Expanders.MetaNode) = (merge!(meta, block.dict); nothing)

const MDTextElements = Union{
    Markdown.Bold,
    Markdown.Header,
    Markdown.Italic,
}
walk(f, meta, block::MDTextElements) = f(block) ? walk(f, meta, block.text)  : nothing

if isdefined(Base.Markdown, :Footnote)
    walk(f, meta, block::Markdown.Footnote) = f(block) ? walk(f, meta, block.text) : nothing
end

if isdefined(Base.Markdown, :Admonition)
    walk(f, meta, block::Markdown.Admonition) = f(block) ? walk(f, meta, block.content) : nothing
end

walk(f, meta, block::Markdown.Image) = f(block) ? walk(f, meta, block.alt)   : nothing
walk(f, meta, block::Markdown.Table) = f(block) ? walk(f, meta, block.rows)  : nothing
walk(f, meta, block::Markdown.List)  = f(block) ? walk(f, meta, block.items) : nothing
walk(f, meta, block::Markdown.Link)  = f(block) ? walk(f, meta, block.text)  : nothing

end


text:
using LinearMaps, SparseArrays, PositiveFactorizations, Roots
using .AbstractOperator, .GP

function makeregR0(F::Operator1D)
    n = length(F.ρ) - F.nfixed
    LinearMap(I, n)
end

function makeregR1(F::Operator1D)
    n = length(F.ρ) - F.nfixed
    LinearMap(R1Dop, Rt1Dop, n)
end

function R1Dop(x::Vector)
    vcat(0, diff(x))
end

function Rt1Dop(y::Vector)
    x = vcat(-diff(y),y[end])
    x[1] = -y[2]
    x
end

function makereg(r::Symbol, F::Operator)
    r == :R0 && return sparse(makeregR0(F))
    r == :R1 && return sparse(makeregR1(F))
    r == :R2 && return sparse(makeregR1(F)*makeregR1(F))
    error("unknown regularization")
end

function pushback(m, lo, hi, reflect=true)
    for i in eachindex(m)
        while (m[i]<lo) || (m[i]>hi)
                (m[i]<lo) && (m[i] = reflect ? 2*lo - m[i] : lo)
                (m[i]>hi) && (m[i] = reflect ? 2*hi - m[i] : hi)
        end
    end
end

function newtonstep(m::AbstractVector, m0::AbstractVector, F::Operator, λ²::Float64, β²::Float64, R::SparseMatrixCSC;
                    regularizeupdate=true)
    JtW, Wr = F.J'*F.W, F.W*F.res
    H = (JtW*(JtW)' + λ²*R'R + λ²*β²*I)
    U = cholesky(Positive, H, Val{false}).U
    if !regularizeupdate # regularize model
        -U\(U'\(JtW*Wr + λ²*R'R*m + λ²*β²*(m - m0)))
    else
        -U\(U'\(JtW*Wr))
    end
end

function occamstep(m::AbstractVector, m0::AbstractVector, Δm::AbstractVector, mnew::Vector{Vector{Float64}},
                χ²::Vector{Float64}, λ²::Vector{Vector{Float64}}, F::Operator, R::SparseMatrixCSC, target,
                lo, hi, λ²min, λ²max, β², ntries; knownvalue=NaN, regularizeupdate = false)

    χ²₀ = getχ²(F, m, computeJ=true)
    knownvalue *= χ²₀
    α = 1
    count = 0
    countmax = 6
    while true
        for (i, l²) in enumerate(10 .^reverse(LinRange(λ²min, λ²max, ntries)))
            count == 0 && (Δm[i] = newtonstep(m, m0, F, l², β², R, regularizeupdate=regularizeupdate))
            mnew[i] = m + α*Δm[i]
            pushback(mnew[i], lo, hi)
            chi2 = getχ²(F, mnew[i])
            nu = [l²; α]
            count == 0 && (push!(χ², chi2); push!(λ², nu))
            χ²[i] = chi2
            λ²[i] = nu
            χ²[i] <= knownvalue && (count = countmax; break)
        end
        all(χ² .>= χ²₀) && (α = α/2)
        count += 1
        count > countmax && break
    end
    idx, foundroot = -1, false
    if all(χ² .> target)
        idx = argmin(χ²)
    else
        idx = findlast(χ² .<= target)
        a = log10(λ²[idx][1])
        b = λ²max # isa(bidx, Nothing) ? λ²max : log10(λ²[bidx][1]) # maybe just set default λ²max
        dophase2(m, m0, F, R, regularizeupdate, lo, hi, target, a, b, mnew,  χ², λ², β², idx)
        foundroot = true
    end
    idx, foundroot
end

function dophase2(m, m0, F, R, regularizeupdate, lo, hi, target, a, b, mnew,  χ², λ², β², idx, reflect=true)
    α = λ²[idx][2]
    if (fλ²(α, m, m0, F, 10. ^b, β², R, regularizeupdate, lo, hi) <= target)
        λ²[idx][1] = 10. ^b
    else
        λ²[idx][1] = 10^find_zero(l² -> fλ²(α, m, m0, F, 10^l², β², R, regularizeupdate, lo, hi) - target, (a, b))
    end
    mnew[idx] = m + α*newtonstep(m, m0, F, λ²[idx][1], β², R, regularizeupdate=regularizeupdate)
    pushback(mnew[idx], lo, hi, reflect)
    χ²[idx] = getχ²(F, mnew[idx])
    nothing
end

function fλ²(α, m, m0, F, l², β², R, regularizeupdate, lo, hi, reflect=true)
    mnew = m + α*newtonstep(m, m0, F, l², β², R, regularizeupdate=regularizeupdate)
    pushback(mnew, lo, hi, reflect)
    getχ²(F, mnew)
end

function getχ²(F, m; computeJ=false)
    r, W = F.res, F.W
    r₀ = copy(r)
    getresidual(F, m, computeJ=computeJ)
    χ² = norm(W*r)^2
    computeJ || (r .= r₀)
    return χ²
end

function initbo(λ²min, λ²max, λ²frac, ntestdivsλ², αmin, αmax, αfrac, ntestdivsα)
    l2 = LinRange(λ²min, λ²max, ntestdivsλ²) # test range for surrogate
    α = LinRange(αmin, αmax, ntestdivsα) # test range for surrogate
    λ²GP = [((λ²max-λ²min)/λ²frac)^2; (abs(α[end]-α[1])/αfrac)^2] # length scale square for surrogate
    X1, X2 = [x1 for x1 in l2, x2 in α], [x2 for x1 in l2, x2 in α]
    t = [X1[:]';X2[:]']
    t, λ²GP
end

function bostep(m::AbstractVector, m0::AbstractVector, mnew::Vector{Vector{Float64}}, χ²::Vector{Float64}, λ²sampled::Vector{Vector{Float64}},
                    t::Array{Float64, 2}, β²::Float64, λ²GP::Array{Float64, 1}, F::Operator, R::SparseMatrixCSC, target, lo, hi;
                   regularizeupdate = false,

                   ## GP stuff
                   demean = true,
                   κ = GP.Mat52(),
                   δtry = 1e-2,
                   acqfun = GP.EI(),
                   ntries = 6,
                   firstvalue = :last,
                   knownvalue = NaN,
                   breakonknown = false)

    χ²₀ = getχ²(F, m, computeJ=true)
    knownvalue *= χ²₀

    ttrain = zeros(size(t,1), 0)
    for i = 1:ntries
        nextpos, = getBOsample(κ, χ²', ttrain, t, λ²GP, δtry, demean, knownvalue, firstvalue, acqfun)
        push!(λ²sampled, [10^t[1, nextpos]; 2^t[2,nextpos]])
        mnew[i] = m + 2^t[2,nextpos]*newtonstep(m, m0, F, 10^t[1,nextpos], β², R, regularizeupdate=regularizeupdate)
        pushback(mnew[i], lo, hi)
        push!(χ², getχ²(F, mnew[i])) # next training value
        ttrain = hcat(ttrain, t[:,nextpos]) # next training location
        (χ²[i] <= knownvalue && breakonknown) && break
    end
    idx, foundroot = -1, false
    if all(χ² .> target)
        idx = argmin(χ²)
    else
        sortedλ²idx = sortperm(vec(ttrain[1,:]))
        sortedχ²idx = findlast(χ²[sortedλ²idx] .<= target)
        idx = sortedλ²idx[sortedχ²idx]
        a = log10(λ²sampled[idx][1])
        b = maximum(t[1,:]) # λ²max in log10 units
        dophase2(m, m0, F, R, regularizeupdate, lo, hi, target, a, b, mnew,  χ², λ²sampled, β², idx)
        foundroot = true
    end
    idx, foundroot
end

function getBOsample(κ, χ², ttrain, t, λ²GP, δtry, demean, knownvalue, firstvalue, acqfun)
    # χ², ttrain, t are row major
    ntrain = length(ttrain)
    if ntrain > 0
        ytest, σ2, = GP.GPfit(κ, χ², ttrain, t, λ²GP, δtry, demean=demean)
        diagσ2 = diag(σ2)
        nextpos = argmax(GP.getAF(acqfun, vec(χ²), vec(ytest), diagσ2, findmin=true, knownvalue=knownvalue))
    else
        nextpos = AFoneiter(firstvalue, size(t, 2))
    end
    nextpos
end

function AFoneiter(firstvalue::Symbol, n::Int)
    if firstvalue == :first
        return 1
    elseif firstvalue == :last
        return n
    elseif firstvalue == :middle
        return round(Int, middle(1:n))
    elseif firstvalue == :random
        return rand(1:n)
    end
end

function gradientinv(   m::AbstractVector,
                        m0::AbstractVector,
                        F::Operator;
                        regtype=:R1,
                        nstepsmax = 10,
                        ntries = 6,
                        target = nothing,
                        lo = -3.,
                        hi = 1.,
                        λ²min = 0,
                        λ²max = 8,
                        λ²frac=5,
                        β² = 0.,
                        ntestdivsλ²=50,
                        αmin=-4,
                        αmax=0,
                        αfrac=4,
                        ntestdivsα=32,
                        regularizeupdate = false,
                        knownvalue=0.7,
                        firstvalue=:last,
                        κ = GP.Mat52(),
                        breakonknown=false,
                        dobo = false,
                        fname="")
    R = makereg(regtype, F)
    ndata = length(F.res)
    isnothing(target) && (target = ndata)
    target₀ = target
    mnew = [[similar(m) for i in 1:ntries] for j in 1:nstepsmax]
    χ²   = [Vector{Float64}(undef, 0) for j in 1:nstepsmax]
    λ² = [Vector{Vector{Float64}}(undef, 0) for j in 1:nstepsmax]
    oidx = zeros(Int, nstepsmax)
    ndata = length(F.res)
    t, λ²GP = initbo(λ²min, λ²max, λ²frac, ntestdivsλ², αmin, αmax, αfrac, ntestdivsα)
    istep = 1
    io = open_history(fname)
    if !dobo
        Δm = [similar(m) for i in 1:ntries]
    end
    while true
        if dobo
            idx, foundroot = bostep(m, m0, mnew[istep], χ²[istep], λ²[istep], t, β², λ²GP, F, R, target, lo, hi,
            regularizeupdate=regularizeupdate, ntries=ntries, κ = κ,
            knownvalue=knownvalue, firstvalue=firstvalue, breakonknown=breakonknown)
        else
            idx, foundroot = occamstep(m, m0, Δm, mnew[istep], χ²[istep], λ²[istep], F, R, target,
                lo, hi, λ²min, λ²max, β², ntries, knownvalue=knownvalue, regularizeupdate = regularizeupdate)
        end
        prefix = isempty(fname) ? fname : fname*" : "
        @info prefix*"iteration: $istep χ²: $(χ²[istep][idx]) target: $target"
        m = mnew[istep][idx]
        oidx[istep] = idx
        isa(io, Nothing) || write_history(io, [istep; χ²[istep][idx]/target₀; vec(m)])
        foundroot && break
        istep == nstepsmax - 1 && (target = χ²[istep][idx]) # exit with smoothest
        istep += 1
        istep > nstepsmax && break
    end
    isa(io, Nothing) || begin @info "Finished "*fname; close(io) end
    return mnew, χ², λ², oidx
end

function open_history(fname)
    @assert !isfile(fname) "$fname exists"
    if !isempty(fname)
        io = open(fname, "w")
    end
end

function write_history(io, v::Vector)
    for el in v
        msg = @sprintf("%e\t", el)
        write(io, msg)
    end
    write(io, "\n")
end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA