text:
using POMDPs
import POMDPs: create_state, discount, reward
using POMDPToolbox

type ImageMDP <: MDP{Matrix{Int},Int}
    size::Tuple{Int, Int}
end

create_state(mdp::ImageMDP) = Array(Int, mdp.size[1], mdp.size[2])
discount(::ImageMDP) = 0.9
function generate_sr(mdp::ImageMDP, s::Matrix{Int}, a::Int, rng::AbstractRNG, sp::Matrix{Int})
    copy!(sp, s)
    r = sp[a,a]
    i = rand(rng, 1:mdp.size[1])
    j = rand(rng, 1:mdp.size[2])
    sp[i,j] += a
    return sp, r
end
function generate_sr(mdp::ImageMDP, s::Matrix{Int}, a::Int, rng::AbstractRNG)
    sp = copy(s)
    r = sp[a,a]
    i = rand(rng, 1:mdp.size[1])
    j = rand(rng, 1:mdp.size[2])
    sp[i,j] += a
    return sp, r
end
function step!(s::Matrix{Int}, mdp::ImageMDP, a::Int, rng::AbstractRNG)
    i = rand(rng, 1:mdp.size[1])
    j = rand(rng, 1:mdp.size[2])
    s[i,j] += a
    return s
end
function reward(mdp::ImageMDP,s::Matrix{Int},a::Int)
    return s[a,a]
end

type IPolicy <: Policy{Matrix{Int}}
    addend::Int
end
action(p::IPolicy, s::Matrix{Int}, a::Int=1) = p.addend

function simulate_no_alloc{S,A}(sim::RolloutSimulator, mdp::MDP{S,A}, policy::Policy, initial_state::S)

    eps = get(sim.eps, 0.0)
    max_steps = get(sim.max_steps, typemax(Int))

    s = initial_state

    disc = 1.0
    r_total = 0.0
    step = 1

    while disc > eps && !isterminal(mdp, s) && step <= max_steps
        a = action(policy, s)

        sp, r = generate_sr(mdp, s, a, sim.rng)

        r_total += disc*r

        s = sp

        disc *= discount(mdp)
        step += 1
    end

    return r_total
end

function simulate_alloc{S,A}(sim::RolloutSimulator, mdp::MDP{S,A}, policy::Policy, initial_state::S)

    eps = get(sim.eps, 0.0)
    max_steps = get(sim.max_steps, typemax(Int))

    s = deepcopy(initial_state)
    sp = create_state(mdp)

    a = create_action(mdp)

    disc = 1.0
    r_total = 0.0
    step = 1

    while disc > eps && !isterminal(mdp, s) && step <= max_steps
        a = action(policy, s, a)

        sp, r = generate_sr(mdp, s, a, sim.rng, sp)

        r_total += disc*r

        # alternates using the memory allocated for s and sp so nothing new has to be allocated
        tmp = s
        s = sp
        sp = tmp

        disc *= discount(mdp)
        step += 1
    end

    return r_total
end

function simulate_step{S,A}(sim::RolloutSimulator, mdp::MDP{S,A}, policy::Policy, initial_state::S)

    eps = get(sim.eps, 0.0)
    max_steps = get(sim.max_steps, typemax(Int))

    s = deepcopy(initial_state)
    sp = create_state(mdp)

    a = create_action(mdp)

    disc = 1.0
    r_total = 0.0
    step = 1

    while disc > eps && !isterminal(mdp, s) && step <= max_steps
        a = action(policy, s, a)

        r = reward(mdp, s, a)
        step!(s, mdp, a, rng)

        r_total += disc*r

        disc *= discount(mdp)
        step += 1
    end

    return r_total
end

mdp = ImageMDP((500,500))
policy = IPolicy(147)

rng = MersenneTwister(123)
is = rand(rng, Int, mdp.size...)
@show simulate_alloc(RolloutSimulator(max_steps=500, rng=rng), mdp, policy, is)
rng = MersenneTwister(123)
is = rand(rng, Int, mdp.size...)
@show simulate_no_alloc(RolloutSimulator(max_steps=500, rng=rng), mdp, policy, is)
rng = MersenneTwister(123)
is = rand(rng, Int, mdp.size...)
@show simulate_step(RolloutSimulator(max_steps=500, rng=rng), mdp, policy, is)

println("Step")
rtot=0.0
@time for i = 1:100
    rng = MersenneTwister(i)
    is = rand(rng, Int, mdp.size...)
    rtot += simulate_step(RolloutSimulator(max_steps=500, rng=rng), mdp, policy, is)
end
@show rtot

println("Pre Allocation:")
rtot=0.0
@time for i = 1:100
    rng = MersenneTwister(i)
    is = rand(rng, Int, mdp.size...)
    rtot += simulate_alloc(RolloutSimulator(max_steps=500, rng=rng), mdp, policy, is)
end
@show rtot

println("No Allocation")
rtot=0.0
@time for i = 1:100
    rng = MersenneTwister(i)
    is = rand(rng, Int, mdp.size...)
    rtot += simulate_no_alloc(RolloutSimulator(max_steps=500, rng=rng), mdp, policy, is)
end
@show rtot


text:
# TermStructures module
using Dates
# module TermStructures
# core TermStructure methods
function check_range(ts::TermStructure, date::Date)
  date < ts.referenceDate && "Date $date before reference_date $(ts.referenceDate)"
end

function check_range(ts::TermStructure, time_frac::Float64)
  time_frac < 0.0 && "Negative time given: $time"
end

max_date(ts::TermStructure) = ts.referenceDate + Date.Year(100)

time_from_reference(ts::TermStructure, date::Date) = year_fraction(ts.dc, reference_date(ts), date, Date(0), Date(0))

function reference_date(ts::TermStructure)
  if ts.referenceDate == Date(0)
    ts.referenceDate = advance(Dates.Day(ts.settlementDays), ts.calendar, settings.evaluation_date)
  end

  return ts.referenceDate
end

## Null TS ##
struct NullTermStructure <: TermStructure end
struct NullDefaultProbabilityTermStructure <: AbstractDefaultProbabilityTermStructure end
# end

ext: jl


text:
module Util

using Revise, Nullables

"""
    add_quotes(str::String) :: String

Adds quotes around `str` and escapes any previously existing quotes.
"""
function add_quotes(str::String) :: String
  if ! startswith(str, "\"")
    str = "\"$str"
  end
  if ! endswith(str, "\"")
    str = "$str\""
  end

  str
end

"""
    strip_quotes(str::String) :: String

Unquotes `str`.
"""
function strip_quotes(str::String) :: String
  if is_quoted(str)
    str[2:end-1]
  else
    str
  end
end

"""
    is_quoted(str::String) :: Bool

Checks weather or not `str` is quoted.
"""
function is_quoted(str::String) :: Bool
  startswith(str, "\"") && endswith(str, "\"")
end

"""
    function walk_dir(dir, paths = String[]; only_extensions = ["jl"], only_files = true, only_dirs = false) :: Vector{String}

Recursively walks dir and `produce`s non directories. If `only_files`, directories will be skipped. If `only_dirs`, files will be skipped.
"""
function walk_dir(dir, paths = String[]; only_extensions = ["jl"], only_files = true, only_dirs = false) :: Vector{String}
  f = readdir(abspath(dir))
  for i in f
    full_path = joinpath(dir, i)
    if isdir(full_path)
      ! only_files || only_dirs && push!(paths, full_path)
      walk_dir(full_path, paths)
    else
      only_dirs && continue

      (last(split(i, ['.'])) in only_extensions) || isempty(only_extensions) && push!(paths, full_path)
    end
  end

  paths
end

"""
    expand_nullable{T}(value::Nullable{T}, default::T) :: T

Returns `value` if it is not `null` - otherwise `default`.
"""
function expand_nullable(value::T)::T where T
  value
end
function expand_nullable(value::Nullable{T}, default::T)::T where T
  if isnull(value)
    default
  else
    Base.get(value)
  end
end

end


text:
using Sockets
import InteractiveUtils

function global_err_handler(e, bt, vscode_pipe_name, cloudRole)
    @warn "Some Julia code in the VS Code extension crashed with" e

    st = stacktrace(bt)
    pipe_to_vscode = connect(vscode_pipe_name)
    try
        # Send cloudRole as one line
        println(pipe_to_vscode, cloudRole)
        # Send error type as one line
        println(pipe_to_vscode, typeof(e))

        # Send error message
        temp_io = IOBuffer()
        InteractiveUtils.versioninfo(temp_io, verbose = false)
        println(temp_io)
        println(temp_io)
        showerror(temp_io, e)
        error_message_str = chomp(String(take!(temp_io)))
        n = count(i->i == '\n', error_message_str) + 1
        println(pipe_to_vscode, n)
        println(pipe_to_vscode, error_message_str)

        # Send stack trace, one frame per line
        # Note that stack frames need to be formatted in Node.js style
        for s in st
            print(pipe_to_vscode, " at ")
            Base.StackTraces.show_spec_linfo(pipe_to_vscode, s)

            filename = string(s.file)

            # Now we need to sanitize the filename so that we don't transmit
            # things like a username in the path name
            filename = normpath(filename)
            if isabspath(filename)
                root_path_of_extension = normpath(joinpath(@__DIR__, "..", ".."))
                if startswith(filename, root_path_of_extension)
                    filename = joinpath(".", filename[lastindex(root_path_of_extension) + 1:end])
                else
                    filename = basename(filename)
                end
            else
                filename = basename(filename)
            end

            # Use a line number of "0" as a proxy for unknown line number
            print(pipe_to_vscode, " (", filename, ":", s.line >= 0 ? s.line : "0", ":1)")

            # TODO Unclear how we can fit this into the Node.js format
            # if s.inlined
            #     print(pipe_to_vscode, " [inlined]")
            # end

            println(pipe_to_vscode)
        end
    finally
        close(pipe_to_vscode)
    end
    Base.display_error(e, bt)
    exit(1)
end


text:
__precompile__()

module DelayDiffEq

using Reexport
@reexport using OrdinaryDiffEq

using DataStructures, RecursiveArrayTools, MuladdMacro,
      DiffEqDiffTools, ForwardDiff, NLSolversBase, Roots

import OrdinaryDiffEq: initialize!, perform_step!, loopfooter!, loopheader!, alg_maximum_order,
                       handle_tstop!, ODEIntegrator, savevalues!, postamble!,
                       handle_callback_modifiers!,handle_discontinuities!,initialize_callbacks!,
                       tstop_saveat_disc_handling

import DiffEqBase: solve!, resize!, u_cache, user_cache, du_cache, full_cache,
                   deleteat!, terminate!, u_modified!, get_proposed_dt, set_proposed_dt!,
                   has_reinit, reinit!, auto_dt_reset!, AbstractContinuousCallback,
                   AbstractDDEAlgorithm, AbstractDDEIntegrator, DEIntegrator, AbstractDDEProblem,
                   find_callback_time

using OrdinaryDiffEq: Rosenbrock23Cache, Rosenbrock32Cache, ImplicitEulerCache,
    TrapezoidCache

using DiffEqBase: check_error!, apply_callback!, addsteps!

include("discontinuity_type.jl")
include("integrator_type.jl")
include("integrator_utils.jl")
include("integrator_interface.jl")
include("interpolants.jl")
include("history_function.jl")
include("algorithms.jl")
include("callbacks.jl")
include("alg_utils.jl")
include("solve.jl")
include("utils.jl")

export Discontinuity, MethodOfSteps

end # module


text:
"""
Provides two functions, [`search_wdi`](@ref) and [`wdi`](@ref), for searching and fetching World Development Indicators data from the World Bank.
"""
module WorldBankData

using HTTP
using JSON
using DataFrames

export wdi, search_wdi

function download_parse_json(url::String; verbose::Bool=false)
    if verbose
        println("download: ", url)
    end
    request = HTTP.get(url)
    if request.status != 200
        error("download failed")
    end
    JSON.parse(String(request.body))
end

# convert json from worldbank for an indicator to dataframe
function parse_indicator(json::Array{Any,1})::DataFrame
    indicator_val = String[]
    name_val = String[]
    description_val = String[]
    source_database_val = String[]
    source_organization_val = String[]

    for d in json[2]
        append!(indicator_val, [d["id"]])
        append!(name_val, [d["name"]])
        append!(description_val, [d["sourceNote"]])
        append!(source_database_val, [d["source"]["value"]])
        append!(source_organization_val, [d["sourceOrganization"]])
    end

    DataFrame(indicator=indicator_val, name=name_val,
              description=description_val, source_database=source_database_val,
              source_organization=source_organization_val)
end

function tofloat(f::AbstractString)::Union{Missing,Float64}
    x = tryparse(Float64, f)
    x isa Nothing ? missing : x
end

# convert country json to DataFrame
function parse_country(json::Array{Any,1})::DataFrame
    iso3c_val = String[]
    iso2c_val = String[]
    name_val = String[]
    region_val = String[]
    capital_val = String[]
    longitude_val = String[]
    latitude_val = String[]
    income_val = String[]
    lending_val = String[]

    for d in json[2]
        append!(iso3c_val, [d["id"]])
        append!(iso2c_val, [d["iso2Code"]])
        append!(name_val, [d["name"]])
        append!(region_val, [d["region"]["value"]])
        append!(capital_val, [d["capitalCity"]])
        append!(longitude_val, [d["longitude"]])
        append!(latitude_val, [d["latitude"]])
        append!(income_val, [d["incomeLevel"]["value"]])
        append!(lending_val, [d["lendingType"]["value"]])
    end

    longitude_val = tofloat.(longitude_val)
    latitude_val = tofloat.(latitude_val)

    DataFrame(iso3c=iso3c_val, iso2c=iso2c_val, name=name_val,
              region=region_val, capital=capital_val, longitude=longitude_val,
              latitude=latitude_val, income=income_val, lending=lending_val)
end

function download_indicators(;verbose::Bool=false)::DataFrame
    dat = download_parse_json("https://api.worldbank.org/v2/indicators?per_page=25000&format=json", verbose=verbose)

    parse_indicator(dat)
end

function download_countries(;verbose::Bool=false)::DataFrame
    dat = download_parse_json("https://api.worldbank.org/v2/countries/all?per_page=25000&format=json", verbose=verbose)

    parse_country(dat)
end

country_cache = false
indicator_cache = false

function reset_country_cache()
    global country_cache = false
end

function reset_indicator_cache()
    global indicator_cache = false
end

function set_country_cache(df::AbstractDataFrame)
    global country_cache = df
end

function set_indicator_cache(df::AbstractDataFrame)
    global indicator_cache = df
end

function get_countries(;verbose::Bool=false)
    if country_cache == false
        set_country_cache(download_countries(verbose=verbose))
    end
    country_cache
end

function get_indicators(;verbose::Bool=false)
    if indicator_cache == false
        set_indicator_cache(download_indicators(verbose=verbose))
    end
    indicator_cache
end

# The "." character is illegal in symbol, but used a lot in WDI. replace by "_".
# example: NY.GNP.PCAP.CD becomes NY_GNP_PCAP_CD
function make_symbol(x::String)::Symbol
    Symbol(replace(x, "." => "_"))
end

make_symbol(x::Symbol) = x

df_match(df::AbstractDataFrame, entry::String, regex::Regex)::DataFrame = df[occursin.(Ref(regex), df[!, make_symbol(entry)]),:]

function country_match(entry::String, regex::Regex)::DataFrame
    df = get_countries()
    df_match(df, entry, regex)
end

function indicator_match(entry::String, regex::Regex)::DataFrame
    df = get_indicators()
    df_match(df, entry, regex)
end

function search_countries(entry::String, regx::Regex)::DataFrame
    entries = ["name","region","capital","iso2c","iso3c","income","lending"]
    if !(entry in entries)
        error("unsupported country entry: \"", entry, "\". supported are:\n", entries)
    end
    country_match(entry, regx)
end

function search_indicators(entry::String, regx::Regex)::DataFrame
    entries = ["name","description","topics","source_database","source_organization"]
    if !(entry in entries)
        error("unsupported indicator entry: \"", entry, "\". supported are\n", entries)
    end
    indicator_match(entry, regx)
end

"""
search_wdi(data::String, entry::String, regx::Regex)::DataFrame

Search World Development Indicators for countries or indicators.

https://datacatalog.worldbank.org/dataset/world-development-indicators

**Arguments**

* `data` : data to search for: "indicators" or "countries"
* `entry` : entry to lookup
  for countries: `name`,`region`,`capital`,`iso2c`,`iso3c`,`income`,`lending`
  for indicators: `name`, `description`, `topics`, `source_database`, `source_organization`
* `regex` : regular expression to find

# Examples
```julia
search_wdi("countries", "name", r"united"i)
search_wdi("indicators", "description", r"gross national"i)
```
"""
function search_wdi(data::String, entry::String, regx::Regex)::DataFrame
    if data == "countries"
        return search_countries(entry, regx)
    elseif data == "indicators"
        return search_indicators(entry, regx)
    else
        error("unsupported data source:", data, ". supported are: \"countries\" or \"indicators\"")
    end
end

function clean_entry(x::Union{AbstractString,Nothing})
    if typeof(x) == Nothing
        return "NA"
    else
        return x
    end
end

function clean_append!(vals::Union{Array{String,1},Array{String,1}}, val::Union{String,String,Nothing})
    append!(vals, [clean_entry(val)])
end

function parse_wdi(indicators::Array{String,1}, jsondata::Array{Any,1}, startyear::Integer, endyear::Integer, returnlong::Bool)::DataFrame
    country_id = String[]
    country_name = String[]
    value = Union{Float64,Missing}[]
    date = String[]
    indicator_id = String[]

    for d in jsondata
        clean_append!(country_id, d["country"]["id"])
        clean_append!(country_name, d["country"]["value"])
        clean_append!(indicator_id, d["indicator"]["id"])
        push!(value, d["value"] isa Nothing ? missing : d["value"])
        clean_append!(date, d["date"])
    end

    date = tofloat.(date)

    dflong = DataFrame(iso2c=country_id, country=country_name, year=date, indicator=indicator_id, value=value)

    if returnlong
        return dflong
    end

    dflong[!,:indicator] = map(x -> replace(x, "." => "_"), dflong[!,:indicator])
    noindcols = filter(x -> (x != "indicator") && (x != "value"), names(dflong))
    dfwide = unstack(dflong, noindcols, :indicator, :value)

    dfwide
end

function array_to_querystr(values::Array{String,1})::String
    res = ""
    for (i, v) in enumerate(values)
        if i == 1
            res = v
        else
            res = string(res, ";", v)
        end
    end
    res
end

function get_url(indicators::Union{String,Array{String,1}}, countries::Union{String,Array{String,1}}, startyear::Integer, endyear::Integer, sourceid::Integer; verbose::Bool=false)::String
    datestr = ""
    if startyear == -1 && endyear == -1
        datestr = ""
    elseif startyear == -1 && endyear != -1
        datestr = string("&date=1800:", endyear)
    elseif startyear != -1 && endyear == -1
        error("need to also set endyear if startyear is given")
    elseif startyear != -1 && endyear != -1
        datestr = string("&date=", startyear, ":", endyear)
    else
        error("internal error. should never get here")
    end

    countriesstr = ""
    if typeof(countries) == String
        countriesstr = countries
    elseif typeof(countries) == Array{String,1}
        countriesstr = array_to_querystr(countries)
    end

    indicatorstr = ""
    if typeof(indicators) == String
        indicatorstr = indicators
    elseif typeof(indicators) == Array{String,1}
        indicatorstr = array_to_querystr(indicators)
    end

    url = string("https://api.worldbank.org/v2/countries/", countriesstr, "/indicators/", indicatorstr,
                 "?format=json&per_page=25000", datestr, "&source=", sourceid)

    url
end

function wdi_download(indicators::Union{String,Array{String,1}}, countries::Union{String,Array{String,1}}, startyear::Integer, endyear::Integer, sourceid::Integer, dflong::Bool; verbose::Bool=false)::DataFrame
    url = get_url(indicators, countries, startyear, endyear, sourceid, verbose=verbose)
    jsondata = download_parse_json(url, verbose=verbose)

    if length(jsondata) == 1
        d = jsondata[1]
        if haskey(d, "message")
            msg = d["message"][1]
            error("request error. response key=\"", msg["key"], "\". value=\"", msg["value"], "\"")
        end
        error("response json data problem:", jsondata)
    end

    if length(jsondata) != 2
        error("wrong length json reply:", jsondata)
    end

    parse_wdi(indicators, jsondata[2], startyear, endyear, dflong)
end

"""
function wdi(indicators::Union{String,Array{String,1}}, countries::Union{String,Array{String,1}}, startyear::Integer=-1, endyear::Integer=-1; extra::Bool=false, sourceid::Integer=2, dflong::Bool=false, verbose::Bool=false)::DataFrame

Download data from World Development Indicators (WDI) Data Catalog of the World Bank.

https://datacatalog.worldbank.org/dataset/world-development-indicators

# Arguments

- `indicators` : indicator name or array of indicators
- `countries` : string or string array of ISO 2 or ISO 3 letter country codes or `all` for all countries.
- `startyear` : first year to include
- `endyear` : last year to include (required if startyear is set)
- `extra` : if `true` additional country data should be included (region, capital, longitude, latitude, income, lending)
- `sourceid` : source number, see https://api.worldbank.org/v2/sources
- `dflong` : long dataframe format. default is wide. with many indicators `long` might be easier to work with.
- `verbose` : if `true` print URLs downloaded

# Examples
```julia
df = wdi("SP.POP.TOTL", "US")
df = wdi("SP.POP.TOTL", "US", 1980, 2012)
df = wdi("SP.POP.TOTL", "USA", 1980, 2012, extra=true)
df = wdi("SP.POP.TOTL", "all", 2000, 2000, verbose=true, extra=true)
df = wdi("SP.POP.TOTL", ["US","BR"], 1980, 2012, extra=true)
df = wdi(["SP.POP.TOTL", "NY.GDP.MKTP.CD"], ["US","BR"], 1980, 2012, extra=true)
df = wdi(["SP.POP.TOTL", "NY.GDP.MKTP.CD"], ["US","BR"], 1980, 2012, dflong=true)
```
"""
function wdi(indicators::Union{String,Array{String,1}}, countries::Union{String,Array{String,1}}, startyear::Integer=-1, endyear::Integer=-1; extra::Bool=false, sourceid::Integer=2, dflong::Bool=false, verbose::Bool=false)::DataFrame
    if typeof(countries) == String
        countries = [countries]
    end

    if ! (startyear <= endyear)
        error("startyear has to be < endyear. startyear=", startyear, ". endyear=", endyear)
    end

    if typeof(indicators) == String
        indicators = [indicators]
    end

    df = wdi_download(indicators, countries, startyear, endyear, sourceid, dflong, verbose=verbose)

    if extra
        cntdat = get_countries(verbose=verbose)
        df = innerjoin(df, cntdat, on=:iso2c)
    end

    sort!(df, [order(:iso2c), order(:year)])

    df
end

end


text:
function ah_time(t::Array{Int32,1}, ts::Float32)
  broadcast!(bswap, t, t)
  y = t[1]
  ts = Float64(bswap(ts))
  return y2μs(y) + Int64(md2j(y, t[2], t[3])-one(Int32))*86400000000 +
    Int64(t[4])*3600000000 + Int64(t[5])*60000000 + round(Int64, ts*sμ)
end

function read_ah_str(io::IO, j::Int64)
  n = bswap(read(io, Int32))
  r = rem(n, 4)
  # these are all very short; a "for" loop is fastest
  for i = 1:n
    j += 1
    BUF.buf[j] = read(io, UInt8)
  end
  r > 0 && skip(io, 4-r)
  return j
end

function read_comm!(io::IO, buf::Array{UInt8,1}, full::Bool)
  n = bswap(read(io, Int32))
  r = rem(n, 4)
  j = 0
  k = 0
  if full
    checkbuf!(buf, n)
    readbytes!(io, buf, n)
    while k < n
      k += 1
      c = getindex(buf, k)
      if c != 0x00
        j += 1
        BUF.buf[j] = c
      end
    end
    (r > 0) && skip(io, 4-r)
  else
    (r > 0) && (n += 4-r)
    skip(io, n)
  end
  return j
end

# obnoxiously, our sample AH files use both 0x00 and 0x20 as string spacers
function mk_ah_id(js::Int64, jc::Int64, jn::Int64)
  fill!(BUF.id, 0x00)
  j = 0
  k = jc+1
  while j < 2 && k < jn
    c = getindex(BUF.buf, k)
    if c != 0x00 && c != 0x20
      j += 1
      BUF.id[j] = c
    end
    k += 1
  end
  j += 1
  J = j+5
  k = 1
  BUF.id[j] = 0x2e
  while j < J && k < js
    c = getindex(BUF.buf, k)
    if c != 0x00 && c != 0x20
      j += 1
      BUF.id[j] = c
    end
    k += 1
  end
  BUF.id[j] = 0x2e
  BUF.id[j+1] = 0x2e
  j += 1
  J = j+3
  k = js+1
  while j < J && k < jc
    c = getindex(BUF.buf, k)
    if c != 0x00 && c != 0x20
      j += 1
      BUF.id[j] = c
    end
    k += 1
  end
  return j
end

# AH-2
function read_ah2!(S::GphysData, ahfile::String;
  v::Int64=SeisIO.KW.v,
  full::Bool=false
  )

  io = open(ahfile, "r")
  str = getfield(BUF, :sac_cv)
  ti = BUF.date_buf
  resize!(ti, 5)

  while !eof(io)
    loc = GeoLoc()

    ver = bswap(read(io, Int32))
    ver == 1100 || error("Not a valid AH-2 file!")
    len = bswap(read(io, UInt32))

    # Station header =========================================================
    skip(io, 4)
    js = read_ah_str(io, 0)
    skip(io, 4)
    jc = read_ah_str(io, js)
    skip(io, 4)
    jn = read_ah_str(io, jc)
    jrec = read_ah_str(io, jn)
    jsen = read_ah_str(io, jrec)
    j = mk_ah_id(js, jc, jn)
    id = unsafe_string(pointer(BUF.id), j)

    # Fill SeisChannel header
    setfield!(loc, :az, Float64(bswap(read(io, Float32))))
    setfield!(loc, :inc, 90.0-bswap(read(io, Float32)))
    setfield!(loc, :lat, bswap(read(io, Float64)))
    setfield!(loc, :lon, bswap(read(io, Float64)))
    setfield!(loc, :el, Float64(bswap(read(io, Float32))))
    gain = bswap(read(io, Float32))

    C = SeisChannel(id = id,
                    loc = loc,
                    gain = Float64(gain)
                    )

    (v > 1) && println("id = ", C.id)
    if full
      C.misc["recorder"] = unsafe_string(pointer(BUF.buf, jn+1), jrec-jn)
      C.misc["sensor"] = unsafe_string(pointer(BUF.buf, jrec+1), jsen-jrec)
    end
    A0 = bswap(read(io, Float32))
    NP = bswap(read(io, Int32))
    P = zeros(Complex{Float32}, NP)
    read!(io, P)
    NZ = bswap(read(io, Int32))
    Z = zeros(Complex{Float32}, NZ)
    read!(io, Z)
    setfield!(C, :resp, PZResp(a0 = A0, p = P, z = Z))
    j = read_comm!(io, str, full)
    if full
      C.misc["sta_comment"] = unsafe_string(pointer(BUF.buf), j)
    end

    # Event header ==========================================================
    if full
      C.misc["ev_lat"] = bswap(read(io, Float64))
      C.misc["ev_lon"] = bswap(read(io, Float64))
      C.misc["ev_dep"] = bswap(read(io, Float32))
      read!(io, ti)
      t_s = read(io, Float32)
      skip(io, 4)
      j = read_comm!(io, str, full)
      C.misc["event_comment"] = unsafe_string(pointer(BUF.buf), j)
      C.misc["ot"] =  ah_time(ti, t_s)
      (v > 1) && println("ev_lat = ", C.misc["ev_lat"], ", ev_lon = ", C.misc["ev_lon"])
    else
      skip(io, 48)
      read_comm!(io, str, full)
    end

    # Data header ===========================================================
    fmt   = bswap(read(io, Int32))
    nx    = bswap(read(io, UInt32))
    dt    = bswap(read(io, Float32))
    Amax  = bswap(read(io, Float32))
    read!(io, ti)
    t_s   = read(io, Float32)
    n     = bswap(read(io, Int32))
    C.units = n > 0 ? unsafe_string(pointer(read(io, n))) : ""
    n = bswap(read(io, Int32))
    inunits = read(io, n)
    n = bswap(read(io, Int32))
    outunits = read(io, n)
    skip(io, 4)
    j = read_comm!(io, str, full)
    if full
      C.misc["data_comment"] = unsafe_string(pointer(BUF.buf), j)
      C.misc["units_in"] = unsafe_string(pointer(inunits))
      C.misc["units_out"] = unsafe_string(pointer(outunits))
      C.misc["Amax"] = Amax
      C.misc["ti"] = bswap.(copy(ti))
    end
    skip(io, 4)
    k = read_comm!(io, str, full)
    (v > 1) && println("nx = ", nx, ", dt = ", dt, ", Amax = ", Amax, ", k = ", k)
    if full
      # Log all processing to C.notes
      i = 0
      j = 0
      while j < k
        j += 1
        c = BUF.buf[j]
        if c == 0x00
          i = j
        elseif c == 0x3b && j-i > 1
          note!(C, unsafe_string(pointer(BUF.buf[i+1:j-1])) * ", recorded in .ah file log")
          i = j
        end
      end
    end

    nattr = bswap(read(io, Int32))
    if nattr > 0
      UA = Dict{String, String}()
      for i = 1:nattr
        j1 = read_ah_str(io, 0)
        j2 = read_ah_str(io, j)
        k = unsafe_string(pointer(BUF.buf), j1)
        V = unsafe_string(pointer(BUF.buf, j1+1), j2-j1)
        UA[k] = V
      end
      v > 1 && println("UA = ", UA)
      merge!(C.misc, UA)
    end

    # Create C.t
    t0 = ah_time(ti, t_s) #) + round(Int64, t_s*1.0e6)
    t = Array{Int64, 2}(undef, 2, 2)
    setindex!(t, one(Int64), 1)
    setindex!(t, nx, 2)
    setindex!(t, t0, 3)
    setindex!(t, zero(Int64), 4)
    setfield!(C, :t, t)

    # Set C.fs
    setfield!(C, :fs, 1.0/dt)

    # Data
    T = fmt == 1 ? Float32 : Float64
    x = Array{T, 1}(undef, nx)
    read!(io, x)
    x .= bswap.(x)
    v > 2 && println("x = [", x[1], ", ", x[2], ", ", x[3], ", ", x[4], ", ... ", x[nx-3], ", ", x[nx-2], ", ", x[nx-1], ", ", x[nx], "]")

    # Cleanup
    setfield!(C, :x, x)
    push!(S, C)
  end
  close(io)
  resize!(str, 192)
  resize!(ti, 7)
  if length(BUF.buf) != 65535
    resize!(BUF.buf, 65535)
  end
  return S
end

function read_ah1!(S::GphysData, ahfile::String;
  v     ::Int64 = SeisIO.KW.v,
  full  ::Bool  = false
  )

  fullfile = realpath(ahfile)
  io = open(ahfile, "r")
  str = getfield(BUF, :sac_cv)
  ti  = getfield(BUF, :date_buf)
  pz_buf = getfield(BUF, :x)
  if full
    stamp = timestamp() * ": "
  end
  resize!(ti, 5)
  resize!(pz_buf, 125)

  while !eof(io)
    # Create SeisChannel, location container
    misc = Dict{String,Any}()
    notes = Array{String,1}(undef,0)

    # Station header =========================================================
    js = read_ah_str(io, 0)
    jc = read_ah_str(io, js)
    jn = read_ah_str(io, jc)
    j  = mk_ah_id(js, jc, jn)
    id = unsafe_string(pointer(BUF.id), j)
    read!(io, pz_buf)

    # Event header ==========================================================
    if full
      misc["ev_lat"] = bswap(read(io, Float32))
      misc["ev_lon"] = bswap(read(io, Float32))
      misc["ev_dep"] = bswap(read(io, Float32))
      read!(io, ti)
      t_s = read(io, Float32)
      j   = read_comm!(io, str, full)
      misc["event_comment"] = unsafe_string(pointer(BUF.buf), j)
      misc["ot"] =  ah_time(ti, t_s)
      (v > 1) && printstyled("ev_lat = ", misc["ev_lat"], ", ev_lon = ", misc["ev_lon"], ", ev_dep = ", misc["ev_dep"], ", ")
    else
      skip(io, 36)
      read_comm!(io, str, full)
    end

    # Data header ===========================================================
    fmt   = bswap(read(io, Int32))
    nx    = bswap(read(io, UInt32))
    dt    = bswap(read(io, Float32))
    (v > 1) && println("nx = ", nx, ", dt = ", dt)
    if full
      misc["Amax"]  = bswap(read(io, Float32))
      read!(io, ti)
      t_s           = read(io, Float32)
      misc["xmin"]  = bswap(read(io, Float32))
      j             = read_comm!(io, str, full)

      # Log all processing to :notes
      misc["data_comment"] = unsafe_string(pointer(BUF.buf), j)
      k = read_comm!(io, str, full)

      # Log all processing to C.notes
      i = 0
      j = 1
      while j ≤ length(str)
        c = BUF.buf[j]
        if c == 0x00
          i = j
        elseif c == 0x3b && j-i > 1
          push!(notes, stamp * unsafe_string(pointer(BUF.buf, i+1), j-i-1) * ", recorded in .ah file log")
          i = j
        end
        j += 1
      end

      # Set all "extras" in :misc
      ne = bswap(read(io, UInt32))
      if ne > 0
        misc["extras"] = zeros(Float32, ne)
        read!(io, misc["extras"])
        broadcast!(bswap, misc["extras"], misc["extras"])
      end
    else
      skip(io, 4)
      read!(io, ti)
      t_s = read(io, Float32)
      skip(io, 4)
      read_comm!(io, str, full)
      read_comm!(io, str, full)
      ne = bswap(read(io, UInt32))
      if ne > 0
        skip(io, 0x00000004*ne)
      end
    end

    # Set fields in C =======================================================
    broadcast!(bswap, pz_buf, pz_buf)

    # :resp
    NP = round(Int32, pz_buf[6])
    NZ = round(Int32, pz_buf[8])
    P = zeros(Complex{Float32}, NP)
    Z = zeros(Complex{Float32}, NZ)
    NC = max(NP, NZ)
    k = 10
    for i = 1:NC
      (i > NP) || (P[i] = complex(pz_buf[k], pz_buf[k+1]))
      (i > NZ) || (Z[i] = complex(pz_buf[k+2], pz_buf[k+3]))
      k += 4
    end

    j = mk_ah_id(js, jc, jn)
    C = SeisChannel(id,
                    "",
                    GeoLoc("",
                           Float64(pz_buf[1]),
                           Float64(pz_buf[2]),
                           Float64(pz_buf[3]),
                           0.0,
                           0.0,
                           0.0),
                    1.0/dt,
                    Float64(pz_buf[4]),
                    PZResp(pz_buf[5], 0.0f0, P, Z),
                    "",
                    fullfile,
                    misc,
                    notes,
                    mk_t(nx, ah_time(ti, t_s)),
                    Array{fmt == one(Int32) ? Float32 : Float64, 1}(undef, nx))

    x = getfield(C, :x)
    read!(io, x)
    broadcast!(bswap, x, x)
    (v > 2) && println("x[1:5] = ", x[1:5])

    # Cleanup
    push!(S, C)
  end
  close(io)
  resize!(BUF.sac_cv, 192)
  resize!(BUF.date_buf, 7)
  resize!(BUF.x, 65535)
  if length(BUF.buf) != 65535
    resize!(BUF.buf, 65535)
  end
  return S
end


text:
@testset "Haar" begin

    haar_udwt = UDWT_Filter_Haar{Float64}()

    @test ϕ_filter(haar_udwt) ≈ [1/2, 1/2]
    @test tildeϕ_filter(haar_udwt) ≈ [1/2, 1/2]

    @test ψ_filter(haar_udwt) ≈ [-1/2, 1/2]
    @test tildeψ_filter(haar_udwt) ≈ [-1/2, 1/2]

end

@testset "Starck2" begin

    starck2_udwt = UDWT_Filter_Starck2{Float64}()

    @test fcoef(tildeϕ_filter(starck2_udwt)) ≈ [1]
end

@testset "UDWT Transform" begin

    signal = rand(20)

    for filter in [UDWT_Filter_Haar{Float64}()
                   UDWT_Filter_Starck2{Float64}()]

        s = 4
        m = udwt(signal,filter,scale=s)

        @test size(m.W) == (length(signal),s)
        @test size(m.V) == (length(signal),)
        @test scale(m) == s

        signal_from_inv = inverse_udwt(m)

        @test signal ≈ signal_from_inv
    end
end

max_stars_repo_name: vincent-picaud/DirectConvolution

max_issues_repo_name: vincent-picaud/DirectConvolution

max_forks_repo_name: vincent-picaud/DirectConvolution


text:
#simplified concept
#	two components: continuation & system
#	no interface requirements: components call stuff and it either exists or it doesn't
#		exceptions: continuationExec, systemExec must exist
#		types only for datatypes (if at all), generalization is user responsibility
#		global Proj variable is accessed by both components
#			never simultaneously, no sync, no communication (in general, otherwise: user responsibility)
#		no contraints on what components are allowed to do: user responsibility

#	this file specifies the whole interface...

#TODO overarching  Session  type, encapsulating Project, SystemCore, Continuation, Visualization
#	make activeSolution a property of session? removes dependency between SystemCore and Project

#TODO ifft for trigoninterp

# addprocs(1)

const MASTER = splitdir(@__FILE__())[1]
push!(LOAD_PATH, "$MASTER/lib")
using Gtk.ShortNames
using mbNewton, mbPred, mbUtil, mbInterpolate, mbObserve
map(include, [
	"lib/ncmprojMKCONTROLGRID.jl"; "lib/ncmprojRelay.jl";
	"type/Solution.jl"; "type/Branch.jl"; "type/Project.jl";
	"type/SystemCore.jl"; "type/ContinuationMethod.jl"; "type/Visualization.jl";
	"type/Session.jl";
	"type/PC.jl"; "type/Galerkin.jl"; "type/GalerkinViz.jl"
	"lib/ncmprojSAVEHELPER.jl"
	])

#(do menu stuff, other global GUI stuff, saving, ...)

#select continuation method, select system
# include(open_dialog("Select continuation method.", Gtk.GtkNullContainer(), ASCIIString[]))
# include(open_dialog("Select system.", Gtk.GtkNullContainer(), ASCIIString[]))

# @everywhere include("system/lorenz.jl")
# sesL = create(lorenz, lorenz′, lorenzProjection)
# sesL = load("lorenz350nodoublings", lorenz, lorenz′, lorenzProjection)
# save("lorenz", sesL, overwrite=true)

# @everywhere include("system/roessler.jl")
# ses = create(roessler, Hroessler, Jroessler, roesslerProjection)
# ses = create(roessler, roessler′, roesslerProjection)
# ses = load("...", roessler, Hroessler, Jroessler, roesslerProjection)
# ses = load("roessler4612126", roessler, Hroessler, Jroessler, roesslerProjection)
# save("roessler4612126", ses, overwrite=true)

#=
using mbRK
global data
rk4((t,x)->ses.core.f(t,[x;100.5]), .0, rand(3), .01, predCount(10000), init=()->(global data=Array{Float64}(0,3)), callback=(t,y,f)->(global data = [data;y']))

ion()
f = figure()
axes(projection="3d")
plot(data[end÷2:end,1], data[end÷2:end,2], data[end÷2:end,3])
=#

#=
using ProfileView
Profile.clear()
@profile for i in 1:10
	step(ses.cont)
end
ProfileView.view()
=#

# homokline orbits
# thompsen steward
# index fixpoint -> perturbation
# proj T


text:
using Breadth: x, A, B, C, result, result_wp, ⊂, TRUE, FALSE, ⟹, ∧, ∨, P, Q, rewrites, rewrites_wp
using Test

@testset "Baseline" begin
    # check for no errors basically for all functions
    println(result((~P ∨ ~Q) ∨ Q, ~P ∨ (Q ∨ ~Q), 1))
    println(result_wp((~P ∨ ~Q) ∨ Q, ~P ∨ (Q ∨ ~Q), 1))
end
Breadth.Node[(~P ∨ ~Q) ∨ Q,
~P ∨ (~Q ∨ Q),
Q ∨ (~P ∨ ~Q),
~(~(~P ∨ ~Q) ∧ ~Q),
~(~((~P ∨ ~Q) ∨ Q)),
(~Q ∨ ~P) ∨ Q,
~(~(~P) ∧ ~(~Q)) ∨ Q,
~(~(~P ∨ ~Q)) ∨ Q,
(P ⟹ ~Q) ∨ Q, (~(~(~P)) ∨ ~Q) ∨ Q,
(~P ∨ ~(~(~Q))) ∨ Q, (~P ∨ ~Q) ∨ ~(~Q)]

Breadth.Node[(~P ∨ ~Q) ∨ Q,
~P ∨ (~Q ∨ Q),
Q ∨ (~P ∨ ~Q),
~(~((~P ∨ ~Q) ∨ Q)), (~Q ∨ ~P) ∨ Q, ~(~(~P) ∧ ~(~Q)) ∨ Q, ~(~(~P ∨ ~Q)) ∨ Q, (P ⟹ ~Q) ∨ Q, (~(~(~P)) ∨ ~Q) ∨ Q, (~P ∨ ~(~(~Q))) ∨ Q, (~P ∨ ~Q) ∨ ~(~Q)]

Breadth.Node[(~P ∨ ~Q) ∨ Q, ~P ∨ (~Q ∨ Q), Q ∨ (~P ∨ ~Q), ~(~(~P ∨ ~Q) ∧ ~Q), ~(~((~P ∨ ~Q) ∨ Q)), (~P ∨ ~Q) ∨ Q, (~Q ∨ ~P) ∨ Q, ~(~(~P) ∧ ~(~Q)) ∨ Q, ~(~(~P ∨ ~Q)) ∨ Q, (P ⟹ ~Q) ∨ Q, (~P ∨ ~Q) ∨ Q, (~(~(~P)) ∨ ~Q) ∨ Q, (~P ∨ ~Q) ∨ Q, (~(~(~P)) ∨ ~Q) ∨ Q, (~P ∨ ~Q) ∨ Q, (~P ∨ ~(~(~Q))) ∨ Q, (~P ∨ ~Q) ∨ Q, (~P ∨ ~(~(~Q))) ∨ Q, (~P ∨ ~Q) ∨ Q, (~P ∨ ~Q) ∨ ~(~Q)]

Breadth.Node[(~P ∨ ~Q) ∨ Q, ~P ∨ (~Q ∨ Q), Q ∨ (~P ∨ ~Q), ~(~(~P ∨ ~Q) ∧ ~Q), ~(~((~P ∨ ~Q) ∨ Q)), (~P ∨ ~Q) ∨ Q, (~Q ∨ ~P) ∨ Q, ~(~(~P) ∧ ~(~Q)) ∨ Q, ~(~(~P ∨ ~Q)) ∨ Q, (P ⟹ ~Q) ∨ Q, (~P ∨ ~Q) ∨ Q, (~(~(~P)) ∨ ~Q) ∨ Q, (~P ∨ ~Q) ∨ Q, (~(~(~P)) ∨ ~Q) ∨ Q, (~P ∨ ~Q) ∨ Q, (~P ∨ ~(~(~Q))) ∨ Q, (~P ∨ ~Q) ∨ Q, (~P ∨ ~(~(~Q))) ∨ Q, (~P ∨ ~Q) ∨ Q, (~P ∨ ~Q) ∨ ~(~Q)]
Breadth.Node[(~P ∨ ~Q) ∨ Q, ~P ∨ (~Q ∨ Q), Q ∨ (~P ∨ ~Q), ~(~(~P ∨ ~Q) ∧ ~Q), ~(~((~P ∨ ~Q) ∨ Q)), (~P ∨ ~Q) ∨ Q, (~Q ∨ ~P) ∨ Q, ~(~(~P) ∧ ~(~Q)) ∨ Q, ~(~(~P ∨ ~Q)) ∨ Q, (P ⟹ ~Q) ∨ Q, (~P ∨ ~Q) ∨ Q, (~(~(~P)) ∨ ~Q) ∨ Q, (~P ∨ ~Q) ∨ Q, (~(~(~P)) ∨ ~Q) ∨ Q, (~P ∨ ~(~(~Q))) ∨ Q, (~P ∨ ~Q) ∨ Q, (~P ∨ ~(~(~Q))) ∨ Q, (~P ∨ ~Q) ∨ Q, (~P ∨ ~Q) ∨ ~(~Q)]
Breadth.Node[(~P ∨ ~Q) ∨ Q, ~P ∨ (~Q ∨ Q), Q ∨ (~P ∨ ~Q), ~(~(~P ∨ ~Q) ∧ ~Q), ~(~((~P ∨ ~Q) ∨ Q)), (~P ∨ ~Q) ∨ Q, (~Q ∨ ~P) ∨ Q, ~(~(~P) ∧ ~(~Q)) ∨ Q, ~(~(~P ∨ ~Q)) ∨ Q, (P ⟹ ~Q) ∨ Q, (~P ∨ ~Q) ∨ Q, (~(~(~P)) ∨ ~Q) ∨ Q, (~P ∨ ~Q) ∨ Q, (~(~(~P)) ∨ ~Q) ∨ Q, (~P ∨ ~Q) ∨ Q, (~P ∨ ~(~(~Q))) ∨ Q, (~P ∨ ~Q) ∨ Q, (~P ∨ ~(~(~Q))) ∨ Q, (~P ∨ ~Q) ∨ Q, (~P ∨ ~Q) ∨ ~(~Q)]
@testset "BUG" begin
    new_values = Any[Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ Q, 1, [])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), (~P ∨ (~Q ∨ Q), 20, [])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), (Q ∨ (~P ∨ ~Q), 22, [])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), (~(~(~P ∨ ~Q) ∧ ~Q), 24, [])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), (~(~((~P ∨ ~Q) ∨ Q)), 35, [])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ Q, 1, [1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~Q ∨ ~P) ∨ Q, 22, [1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), (~(~(~P) ∧ ~(~Q)) ∨ Q, 24, [1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), (~(~(~P ∨ ~Q)) ∨ Q, 35, [1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((P ⟹ ~Q) ∨ Q, 37, [1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ Q, 1, [1, 1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~(~(~P)) ∨ ~Q) ∨ Q, 35, [1, 1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ Q, 1, [1, 1, 1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~(~(~P)) ∨ ~Q) ∨ Q, 35, [1, 1, 1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ Q, 1, [1, 2])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~(~(~Q))) ∨ Q, 35, [1, 2])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ Q, 1, [1, 2, 1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~(~(~Q))) ∨ Q, 35, [1, 2, 1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ Q, 1, [2])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ ~(~Q), 35, [2])])]
    new_values = Any[Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ Q, 1, [])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), (~P ∨ (~Q ∨ Q), 20, [])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), (Q ∨ (~P ∨ ~Q), 22, [])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), (~(~(~P ∨ ~Q) ∧ ~Q), 24, [])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), (~(~((~P ∨ ~Q) ∨ Q)), 35, [])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ Q, 1, [1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~Q ∨ ~P) ∨ Q, 22, [1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), (~(~(~P) ∧ ~(~Q)) ∨ Q, 24, [1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), (~(~(~P ∨ ~Q)) ∨ Q, 35, [1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((P ⟹ ~Q) ∨ Q, 37, [1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ Q, 1, [1, 1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~(~(~P)) ∨ ~Q) ∨ Q, 35, [1, 1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ Q, 1, [1, 1, 1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~(~(~P)) ∨ ~Q) ∨ Q, 35, [1, 1, 1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ Q, 1, [1, 2])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~(~(~Q))) ∨ Q, 35, [1, 2])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ Q, 1, [1, 2, 1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~(~(~Q))) ∨ Q, 35, [1, 2, 1])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ Q, 1, [2])]), Breadth.Proof(Tuple{Breadth.Node,Int64,Array{Int64,N} where N}[((~P ∨ ~Q) ∨ Q, 1, []), ((~P ∨ ~Q) ∨ ~(~Q), 35, [2])])]
    last_steps = [steps(x)[end] for x in new_values]
    println(length(unique(new_values)))
    println(length(unique(last_steps)))
    print(setdiff(new_values, last_steps))
end


text:
using AES

#AES CBC PoC IV reuse

function xorUIntVector(a::Vector{UInt8}, b::Vector{UInt8})
        min_len = min(length(a), length(b))
        result = Vector{UInt8}(undef, min_len)
        for i = 1:min_len
                result[i] = xor(a[i],b[i])
        end
        return result
end

#usually this has nothing to do with passphrases and salts
key = rand(UInt8, div(256, 8))
iv = rand(UInt8, 16)

println("Key: ", key)
println("IV: ", iv)

plain1 = Vector{UInt8}("This plaintext is known and at least as long as the flag.")
plain2 = Vector{UInt8}("I am the flag")

enc1 = AESOFB(plain1, key, iv)
println("Encrypted 1: ", enc1)
println("IV: ", iv)

enc2 = AESOFB(plain2, key, iv)
println("Encrypted 2: ",  enc2)
println("IV: ", iv)

#only using the plain1 to get plain2

keystream = xorUIntVector(Vector{UInt8}(plain1), enc1)
hacked_text = xorUIntVector(keystream, enc2)

println(Vector{UInt8}(hacked_text))
println(Vector{UInt8}(plain2))


text:
# y: data
# yhat_generator: function that samples data given the parameters
# algo parameters: Namedtuple of algorithm parameters
# N_samples: number of samples desired
function ABC_MCMC(y,yhat_generator,algo_parameters,N_samples)
    CPUtic()

    proposal = algo_parameters[:proposal]
    proposalRatio = algo_parameters[:proposalRatio]
    sd = algo_parameters[:sd]
    eta_y = algo_parameters[:eta](y)
    size_y = size(y)
    q = length(algo_parameters[:prior])
    thinning_interval = algo_parameters[:thinning]
    N_samples = N_samples + algo_parameters[:burn_in]
    result = zeros(N_samples,q)
    dist = Inf
    p_prev=Tuple(zeros(q))

    # initial point using ABC
    while (dist > algo_parameters[:epsilon])
        p_prev = rand.(algo_parameters[:prior])
        yhat = yhat_generator(p_prev)
        if size_y != size(yhat)
            continue
        end
        dist = algo_parameters[:d](eta_y, algo_parameters[:eta](yhat))
    end

    naccept=0
    k = 1
    i = 0
    p_cand = Tuple(zeros(q))
    while i <= (N_samples*thinning_interval)
        # propose candidate parameters in log space
        valid = true
        while valid
            p_cand = proposal(p_prev, sd)
            valid = any(p_cand .< 0)
        end
        # p_cand = proposal(p_prev, sd)

        # p_cand = kernel(p_prev, sd)[1]
        #generate data
        yhat = yhat_generator(p_cand)
        # yhat = yhat_generator( p_cand)
        if size_y != size(yhat)
            continue
        end
        i += 1

        dist = algo_parameters[:d](eta_y, algo_parameters[:eta](yhat))

        u = rand(Uniform(0,1))

        proposal_ratio = proposalRatio_Normal(p_prev,p_cand,sd)

        log_prior_prev = sum(logpdf.(algo_parameters[:prior], p_prev))
        log_prior_cand = sum(logpdf.(algo_parameters[:prior], p_cand))

        log_alpha = log_prior_cand + sum(log.(proposal_ratio)) - log_prior_prev
        if (i % thinning_interval==0)
            if (log(u) < log_alpha && dist < algo_parameters[:epsilon])
                result[k,:] .= p_cand
                p_prev = p_cand
                naccept += 1
            else
                result[k,:] .= p_prev
            end
            k += 1
        end
    end
    if algo_parameters[:verbose]
        @printf("acceptance rate=%f\n", naccept/N_samples)
    end
    if algo_parameters[:burn_in] > 0
        result = result[(1+algo_parameters[:burn_in]):end,:]
    end
    t_end = CPUtoc_modified(false)

    return result,t_end
end


text:
## ----------  SymbolTable

struct SymbolTable
    d :: Dict{Symbol, SymbolInfo}
end

Base.show(io :: IO, s :: SymbolTable) =
    Base.print(io, "SymbolTable of length ", length(s));

Base.isempty(s :: SymbolTable) = Base.isempty(s.d);
Base.length(s :: SymbolTable) = Base.length(s.d);

"""
	$(SIGNATURES)

Erase the SymbolTable in place.
test this +++
"""
function erase!(s :: SymbolTable)
    for k in keys(s.d)
        delete!(s.d, k);
    end
end

"""
	$(SIGNATURES)

Constructor for empty SymbolTable.
"""
SymbolTable() = SymbolTable(Dict{Symbol, SymbolInfo}());

"""
	$(SIGNATURES)

Constructor from String array with names, latex symbols, descriptions, groups.

# Example
```
SymbolTable(["capShare"  "\\alpha"  "Capital share"  "Technology";
    "tfp"  "A"  "Total factor productivity"  "Technology"]);
```
"""
function SymbolTable(m :: Matrix{String})
    @assert size(m, 2) == 4  "Invalid no of columns: $(size(m))"
    s = SymbolTable();
    for j = 1 : size(m, 1)
        add_symbol!(s, SymbolInfo(m[j,:]...));
    end
    return s
end

"""
	$(SIGNATURES)

Read a file that defines notation. Comma delimited.

Format:
- header row: Group | Name | Latex | Description
- group rows: only the group entry is filled
- other rows: only the other entries are filled

# Example
```
Group,Name,Latex,Description
G1,,,
,nameG1,lG1,DescrG1
G2,,,
,nameG2,lG2,DescrG2
```
"""
function SymbolTable(fPath :: AbstractString; delimChar = ',')
    s = SymbolTable();
    load_from_csv!(s, fPath; delimChar = delimChar);
    return s
end

"""
	$(SIGNATURES)

Load a comma delimited file into an existing SymbolTable.
Does not erase existing content.
"""
function load_from_csv!(s :: SymbolTable, fPath :: AbstractString; delimChar = ',')
    m = readdlm(fPath, delimChar, String, skipstart = 1);
    grp = "";
    for j = 1 : size(m, 1)
        if !isempty(m[j,1])
            grp = m[j,1];
        else
            add_symbol!(s, SymbolInfo(m[j,2], m[j,3], m[j,4], grp));
        end
    end
    return nothing
end

## ----------  Adding symbols

Base.setindex!(s :: SymbolTable, si :: SymbolInfo) =
    s.d[si.name] = si;

"""
	$(SIGNATURES)

Add a new symbol. Overwrite existing if `replaceExisting == true`.
"""
function add_symbol!(s :: SymbolTable, si :: SymbolInfo;
    replaceExisting :: Bool = false)

    if !replaceExisting
        @assert !has_symbol(s, si.name)  "$(si.name) already exists"
    end
    s.d[si.name] = si;
    return nothing
end

"""
	$(SIGNATURES)

Does symbol exist?
"""
has_symbol(s :: SymbolTable, name) = haskey(s.d, name);

# ---------------  Retrieving symbols

"""
	$(SIGNATURES)

Retrieve a symbol by name. Returns a `SymbolInfo`.

# Example
```
s[:alpha]
```
"""
function Base.getindex(s :: SymbolTable, name :: Symbol)
    @assert has_symbol(s, name)  "Not found: $name"
    return s.d[name]
end

"""
	$(SIGNATURES)

Return description of a symbol.
"""
description(s :: SymbolTable, name :: Symbol) =
    description(s[name]);

"""
	$(SIGNATURES)

Latex symbol for a symbol.
"""
latex(s :: SymbolTable, name :: Symbol) =
    latex(s[name]);

"""
	$(SIGNATURES)

Group for a symbol.
"""
group(s :: SymbolTable, name :: Symbol) =
    group(s[name]);

function group_list(s :: SymbolTable)
    n = length(s);
    gl = Vector{String}();
    if n > 0
        for (name, si) ∈ s.d
            if !(group(si) ∈ gl)
                push!(gl, group(si));
            end
        end
    end
    return gl
end

## -------------  Preamble file

"""
	$(SIGNATURES)

Write a set of symbols to a preamble. The idea is that this defines notation in a paper to be consistent with notation in program generated tables and figures.
"""
function write_preamble(io :: IO, s :: SymbolTable)
    if !isempty(s)
        groupV = sort(group_list(s));
        for group ∈ groupV
            write_preamble_one_group(io, s, group);
        end
    end
    return nothing
end

function write_preamble_one_group(io :: IO, s :: SymbolTable, grp)
    println(io, "% =====  ", grp, ":");
    for (name, si) ∈ s.d
        if group(si) == grp
            println(io, "% ", description(si));
            println(io, newcommand(si));
        end
    end
    println(io, "% ----------");
    return nothing
end

## -----------  Latex doc with notation, by group

"""
	$(SIGNATURES)

Write a tex file with notation, by group. The output looks like:

    Endowments
    * \$\\alpha\$: some endowment
    * \$\\beta\$: another endowment
    Preferences
    * [...]

This writes a document section. It is meant to be included in the Notation section of an existing document.

In `Lyx`: Insert | File | Child document. Use the `input` format (not `include`).
"""
function write_notation_tex(io :: IO, s :: SymbolTable)
    # println(io, "\\Section{Notation}")
    groupV = sort(group_list(s));
    for group ∈ groupV
        write_notation_tex_one_group(io, s, group);
    end
end

function write_notation_tex_one_group(io :: IO, s :: SymbolTable, grp :: String)
    println(io,  grp, "\n");
    # println(io,  "\\begin{itemize}");
    println(io, "\\begin{tabular}{lll}")
    println(io, "\\hline");
    # println(io, "Symbol & Abbrev. & Description", "\\tabularnewline");
    # println(io, "\\hline");

    listV = Vector{String}();
    for (name, si) ∈ s.d
        if group(si) == grp
            push!(listV, notation_line(si));
        end
    end
    listV = sort(listV);
    for line ∈ listV
        println(io, line);
    end
    # println(io, "\\end{itemize}");
    println(io, "\\hline");
    println(io, "\\end{tabular}")
    println(io, " ")
end

# Line in latex file of the form
# `* \\alpha:  Capital share (capShare)`
# where `capShare` is the `newcommand` defined for the symbol.
function notation_line(si :: SymbolInfo)
    return "\$" * latex(si) * "\$ & $(si.name) & " * description(si) *
        "\\tabularnewline";
        # println(io, "\\item \$", latex(si), "\$:  ", description(si),
        # "  ($(si.name))");
end

## --------  For testing

function test_symbol_table(nGroups, nSyms)
    s = SymbolTable();
    for ig = nGroups : -1 : 1
        gName = "$group$ig";
        # Reverse, so that sorting can be checked in output files.
        for j = nSyms : -1 : 1
            sStr = "s_$(ig)_$j";
            sName = Symbol(sStr);
            si = SymbolInfo(sName, "\\alpha_{$ig,$j}", "Description $sStr in $gName", gName);
            add_symbol!(s, si);
        end
    end
    return s
end

# -------------


text:
export get_L,
       get_Q,
       get_B,
       get_Sel,
       get_Exp

"""
## Create FDTD Laplatian matrix
* usage: L = get_L(g::CuboidRoom)
"""
function get_L(g::CuboidRoom)

	Dx = spdiagm((-2*ones(g.Nx), ones(g.Nx-1), ones(g.Nx-1)),(0,1,-1),g.Nx,g.Nx )
	Dy = spdiagm((-2*ones(g.Ny), ones(g.Ny-1), ones(g.Ny-1)),(0,1,-1),g.Ny,g.Ny )
	Dz = spdiagm((-2*ones(g.Nz), ones(g.Nz-1), ones(g.Nz-1)),(0,1,-1),g.Nz,g.Nz )

	#perfect reflection in phase (forward diff)
	Dx[1,2],Dy[1,2],Dz[1,2] =  1, 1, 1
	Dx[1,1],Dy[1,1],Dz[1,1] =  -1, -1, -1
	Dx[end, end-1], Dy[end, end-1], Dz[end, end-1] =   1, 1, 1
	Dx[end, end], Dy[end, end], Dz[end, end] =   -1, -1, -1

	Dxx = kron(speye(g.Nz),kron(speye(g.Ny), Dx))
	Dyy = kron(speye(g.Nz),kron(Dy, speye(g.Nx)))
	Dzz = kron(Dz, kron(speye(g.Ny),speye(g.Nx)))

	L = Dxx+Dyy+Dzz+g.env.scheme.a*(Dxx*Dyy+Dxx*Dzz+Dyy*Dzz)+g.env.scheme.b*(Dxx*Dyy*Dzz)

	return L
end

"""
## Create FDTD loss matrix Q
* usage: Q = get_Q(g::CuboidRoom)
"""
function get_Q(g::CuboidRoom)

	qx = spdiagm([1/(g.ξ[1]*2);zeros(g.Nx-2);1/(g.ξ[2]*2)]) #factor 2 due to fwd diff
	qy = spdiagm([1/(g.ξ[3]*2);zeros(g.Ny-2);1/(g.ξ[4]*2)])
	qz = spdiagm([1/(g.ξ[5]*2);zeros(g.Nz-2);1/(g.ξ[6]*2)])

	Qx = kron(speye(g.Nz),kron(speye(g.Ny),qx))
	Qy = kron(speye(g.Nz),kron(qy,speye(g.Nx)))
	Qz = kron(qz,kron(speye(g.Ny),speye(g.Nx)))

	return Qx+Qy+Qz
end

"""
## Create FDTD matrix B
* usage: B = get_B(g::CuboidRoom, Nt::Int64)
"""
function get_B(g::CuboidRoom, Nt::Int64)

	L = get_L(g)
	Q = get_Q(g)

	λ2 = (g.env.scheme.λ)^2
	Nxyz = size(L,1)

	Qm,A,Qp = get_QmAQp(g)

	B =(kron(speye(Nt+2),spdiagm(Qp))
	+kron(spdiagm(ones(Nt+2-1),-1,Nt+2,Nt+2),-A)
	+kron(spdiagm(ones(Nt+2-2),-2,Nt+2,Nt+2),-spdiagm(Qm)))

	B[1:2*Nxyz,1:2*Nxyz] = speye(2*Nxyz) #initial conditions

	return B

end

"""
## Create selection matrix
* usage: F = get_Sel(g::CuboidRoom,Nt::Int,xr::Vector{Tuple{Int,Int,Int}})
"""
function get_Sel(g::CuboidRoom,Nt::Int, xr::Vector{Tuple{Int,Int,Int}})

	Ks = length(xr)   #Ks number of mics

	indxr = [sub2ind((g.Nx,g.Ny,g.Nz),x...) for x in xr]

	Nxyz = g.Nx*g.Ny*g.Nz

	indposvec = zeros(Int64,Ks*Nt)
	for k =1:Ks
		for i = 1:Nt
			indposvec[i+(k-1)*Nt] = indxr[k]
			indxr[k] += Nxyz
		end
	end
	indposvec = indposvec+2*Nxyz
	F =  sparse((1:Ks*Nt),indposvec,ones(Int64,Ks*Nt),Ks*Nt,Nxyz*(Nt+2))

	return F

end

"""
## Create selection matrix
* usage: E = get_Exp(g::CuboidRoom,Nt::Int64,pos::Vector{Tuple{Int,Int,Int}})
"""
function get_Exp(g::CuboidRoom,Nt::Int,pos::Vector{Tuple{Int,Int,Int}})

	Ks = length(pos)   #Ks number of sources

	indposm = [sub2ind((g.Nx,g.Ny,g.Nz) ,x...) for x in pos]

	Nxyz = g.Nx*g.Ny*g.Nz

	indposvec = zeros(Int64,Ks*Nt)
	for k =1:Ks
		for i = 1:Nt
			indposvec[i+(k-1)*(Nt)] = indposm[k]
			indposm[k] = indposm[k]+Nxyz
		end
	end
	indposvec = indposvec +2*Nxyz #remove ICs
	E =  sparse(indposvec,(1:Ks*Nt),ones(Int64,Ks*Nt),Nxyz*(Nt+2),Ks*Nt)

	return E

end

function get_QmAQp(geo::CuboidRoom)

	λ2 = (geo.env.scheme.λ)^2
	Q = get_Q(geo)
	L = get_L(geo)

	Qm = geo.env.scheme.λ*diag(Q)-1
	A  = λ2*L+2*speye(geo.Nx*geo.Ny*geo.Nz)
	Qp = geo.env.scheme.λ*diag(Q)+1

	return Qm,A,Qp
end

max_stars_repo_name: nantonel/JuFdtd

max_issues_repo_name: nantonel/JuFdtd

max_forks_repo_name: nantonel/JuFdtd


text:
"""
    function optimize(f, sol, range, method::PopulationBase;
        rng=nothing,
        iters=2_000,
        pop_size=30,
        kwargs...
    ) -> OptimizationResults

The `optimize` constitutes an interface to quickly dispatch over all the possible implementations.

In particular, this function dispatches over all `PopulationBase` implementations.

# Arguments

- `f`: The function to optimize.
- `sol`: A candidate solution, should be an `AbstractArray`. The dimension of this array will determine the dimension of the problem.
- `range`: A collection that contains just two elements, regarding the global bounds of the search space.
- `method`: A `PopulationBase`, e.g. `PSO` for Particle Swarm Optimization.

# Keyword arguments

- `rng`: An `AbstractRNG` object. This will be used for creating the population of solutions as well as for the implementation itself. If you wish to enforce reproducibility of results, create an RNG object with a set seed. If you do not pass an RNG object, a `Xorshifts.Xoroshiro128Plus` RNG will be created, using a truly random seed taken from the system.
- `iters`: The total number of iterations for the main loop in the implementation.
- `pop_size`: The total size of the population. The larger the population, the longer it will take to converge. A good population size is somewhere between 20-35.
- `kwargs`: Keyword arguments passed on to the particular implementation. Each implementation has its own keyword arguments so check the documentation for a specific implementation in order to modify these arguments.
"""
function optimize(f, sol, range, method::PopulationBase;
    rng=nothing,
    iters=2_000,
    pop_size=35,
    kwargs...
)
    @assert length(range) == 2 "No more than two global bounds."
    # If no RNG is passsed, create a randomly seeded one
    if isnothing(rng)
        rng = Xorshifts.Xoroshiro128Plus()
    end

    # Create a population
    dim = length(sol)
    lower, upper = range
    pops = Population(pop_size, dim, lower, upper, rng)
    for p in pops
        p.x_best = copy(sol)
    end

    # Dispatch to the PopulationBase method
    val = optimize(f, pops, iters, rng, method; kwargs...)

    return val
end

"""
    function optimize(f, sol, range, method::TrajectoryBase;
        rng=nothing,
        iters=2_000,
        kwargs...
    ) -> OptimizationResults

The `optimize` constitutes an interface to quickly dispatch over all the possible implementations.

In particular, this function dispatches over all `TrajectoryBase` implementations.

# Arguments

- `f`: The function to optimize.
- `sol`: A candidate solution, should be an `AbstractArray`. The dimension of this array will determine the dimension of the problem.
- `range`: A collection that contains just two elements, regarding the global bounds of the search space.
- `method`: A `TrajectoryBase` object, e.g. `SimulatedAnnealing` for the simulated annealing implementation.

# Keyword arguments

- `rng`: An `AbstractRNG` object. This will be used for creating the population of solutions as well as for the implementation itself. If you wish to enforce reproducibility of results, create an RNG object with a set seed. If you do not pass an RNG object, a `Xorshifts.Xoroshiro128Plus` RNG will be created, using a truly random seed taken from the system.
- `iters`: The total number of iterations for the main loop in the implementation.
- `kwargs`: Keyword arguments passed on to the particular implementation. Each implementation has its own keyword arguments so check the documentation for a specific implementation in order to modify these arguments.
"""
function optimize(f, sol, range, method::TrajectoryBase;
    rng=nothing,
    iters=2_000,
    pop_size=35,
    kwargs...
)
    @assert length(range) == 2 "No more than two global bounds."
    # If no RNG is passsed, create a randomly seeded one
    if isnothing(rng)
        rng = Xorshifts.Xoroshiro128Plus()
    end

    # Dispatch to the TrajectoryBase method
    val = optimize(f, range, length(sol), iters, rng, method; kwargs...)

    return val
end


text:
module ErrorTypes

using SumTypes

@sum_type Result{O, E} begin
    Ok{O, E}(::O)
    Err{O, E}(::E)
end

"""
    Result{O, E}

A sum type of either `Ok{O}` or `Err{E}`. Used as return value of functions that
can error with an informative error object of type `E`.
"""
Result

@sum_type Option{T} begin
    Thing{T}(::T)
    None{T}()
end

"""
    Option{T}

A sum type of either `Thing{T}` or `None`. Used instead of `Result` as return
value of functions whose error state need not carry any information.
"""
Option

struct ResultConstructor{R,T}
    x::T
end

"""
    Ok

The success state of a `Result{O, E}`, carrying an object of type `O`.
For convenience, `Ok(x)` creates a dummy value that can be converted to the
appropriate `Result type`.
"""
Ok(x) = ResultConstructor{Ok,typeof(x)}(x)

"""
    Err

The error state of a `Result{O, E}`, carrying an object of type `E`.
For convenience, `Err(x)` creates a dummy value that can be converted to the
appropriate `Result type`.
"""
Err(x) = ResultConstructor{Err,typeof(x)}(x)

function Base.convert(::Type{<:Result{O, E}}, x::ResultConstructor{Ok, O2}
) where {O, O2 <: O, E}
    Ok{O, E}(x.x)
end

function Base.convert(::Type{<:Result{O, E}}, x::ResultConstructor{Err, E2}
) where {O, E, E2 <: E}
    Err{O, E}(x.x)
end

"""
    None{T}

The variant of `Option` signifying absence of a value of type `T`. The singleton
`none` is the instance of `None{Nothing}` (NOT an `Option`), and can be freely
converted to any option value containing a `None`.
"""
None

const none = None{Nothing}().data
Base.:(==)(::None{A}, ::None{B}) where {A, B} = true
Base.convert(::Type{<:Option{T}}, ::None{Nothing}) where T = None{T}()

"""
    Thing(x::T)

Construct a `Thing{T}(x)`, variant of `Option` signifying the presence of a
value `x`.
"""
Thing(x::T) where T = Thing{T}(x)

is_error(x::Result{O, E}) where {O, E} = x.data isa Err{O, E}
is_none(x::Option{T}) where T = x.data isa None{T}

"""
    expect(x::Union{Result, Option}, s::AbstractString

If `x` is of the associated error type, error with message `s`. Else, return
the contained result type.
"""
function expect end

function expect(x::Option, s::AbstractString)
    data = x.data
    data isa Thing ? data._1 : error(s)
end

function expect(x::Result, s::AbstractString)
    data = x.data
    data isa Ok ? data._1 : error(s)
end

"""
    expect_none(x::Option, s::AbstractString)

If `x` contains a None, return nothing. Else, throw an error with message `s`.
"""
expect_none(x::Option, s::AbstractString) = is_none(x) ? nothing : error(s)

@noinline unwrap_err() = error("unwrap on unexpected type")

"""
    unwrap(x::Union{Result, Option})

If `x` is of the associated error type, throw an error. Else, return the contained
result type.
"""
function unwrap end

function unwrap(x::Result)
    data = x.data
    data isa Ok ? data._1 : unwrap_err()
end

function unwrap(x::Option)
    data = x.data
    data isa Thing ? data._1 : unwrap_err()
end

"""
    unwrap_none(x::Option, s::AbstractString)

If `x` contains a None, return nothing. Else, throw an error.
"""
unwrap_none(x::Option) = is_none(x) ? nothing : unwrap_err()

"""
    @?(expr)

Propagate a `Result` or `Option` error value to the outer function.
Evaluate `expr`, which should return a `Result` or `Option`. If it contains
an error value, return the error value from the outer function. Else, the macro
evaluates to the result value of `expr`.

# Example
julia> (f(x::Option{T})::Option{T}) where T = Thing(@?(x) + one(T))

julia> f(Thing(1.0)), f(None{Int}())
(Option{Float64}: Thing(2.0), Option{Int64}: None())
"""
macro var"?"(expr)
    sym = gensym()
    quote
        $(sym)::Union{Result, Option} = $(esc(expr))
        if $sym isa Option
            if ($sym).data isa Thing
                ($sym).data._1
            else
                return none
            end
        else
            if ($sym).data isa Ok
                ($sym).data._1
            else
                return $sym
            end
        end
    end
end

export Result, Option,
    Thing, none, Ok, Err,
    is_error, is_none,
    expect, expect_none,
    unwrap, unwrap_none,
    @?

end # module


text:
# TargetOptions
TargetOptions() = TargetOptions(create_target_options())

"""
    create_target_options() -> CXTargetOptions
Return a pointer to a `clang::TargetOptions` object.
"""
function create_target_options()
    status = Ref{CXInit_Error}(CXInit_NoError)
    opts = clang_TargetOptions_create(status)
    @assert status[] == CXInit_NoError
    return opts
end

function setTriple(x::TargetOptions, triple::String)
    @check_ptrs x
    clang_TargetOptions_setTriple(x, triple, length(triple))
    return nothing
end

function PrintStats(x::TargetOptions)
    @check_ptrs x
    return clang_TargetOptions_PrintStats(x)
end


text:
export run_acdcopf_bf

""
function run_acdcopf_bf(file::String, model_type::Type{T}, solver; kwargs...) where T <: _PM.AbstractBFModel
    data = _PM.parse_file(file)
    PowerModelsACDC.process_additional_data!(data)
    return run_acdcopf_bf(data, model_type, solver; ref_extensions = [add_ref_dcgrid!], kwargs...)
end

""
function run_acdcopf_bf(data::Dict{String,Any}, model_type::Type{T}, solver; kwargs...) where T <: _PM.AbstractBFModel
    return _PM.run_model(data, model_type, solver, post_acdcopf_bf; ref_extensions = [add_ref_dcgrid!], kwargs...)
end

function post_acdcopf_bf(pm::_PM.AbstractPowerModel)
    _PM.variable_bus_voltage(pm)
    _PM.variable_gen_power(pm)
    _PM.variable_branch_power(pm)
    _PM.variable_branch_current(pm)

    variable_active_dcbranch_flow(pm)
    variable_dcbranch_current(pm)
    variable_dc_converter(pm)
    variable_dcgrid_voltage_magnitude(pm)

    _PM.objective_min_fuel_cost(pm)

    _PM.constraint_model_current(pm)
    constraint_voltage_dc(pm)

    for i in _PM.ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i)
    end

    for i in _PM.ids(pm, :bus)
        constraint_power_balance_ac(pm, i)
    end

    for i in _PM.ids(pm, :branch)

        _PM.constraint_power_losses(pm, i)
        _PM.constraint_voltage_magnitude_difference(pm, i)

        _PM.constraint_voltage_angle_difference(pm, i)
        _PM.constraint_thermal_limit_from(pm, i)
        _PM.constraint_thermal_limit_to(pm, i)
    end
    for i in _PM.ids(pm, :busdc)
        constraint_power_balance_dc(pm, i)
    end
    for i in _PM.ids(pm, :branchdc)
        constraint_ohms_dc_branch(pm, i)
        constraint_dc_branch_current(pm, i)
    end
    for i in _PM.ids(pm, :convdc)
        constraint_converter_losses(pm, i)
        constraint_converter_current(pm, i)
        constraint_conv_transformer(pm, i)
        constraint_conv_reactor(pm, i)
        constraint_conv_filter(pm, i)
        if pm.ref[:nw][pm.cnw][:convdc][i]["islcc"] == 1
            constraint_conv_firing_angle(pm, i)
        end
    end
end


text:
# define the interface that unums should follow:

abstract AbstractUnum{B,ESS,FSS} <: Real

isnumber(u::AbstractUnum) = true
isreal(u::AbstractUnum) = true

# NOTE: the definitions below aren't strictly necessary, but I'll keep them here for now so that there's a clear reference
#       of the functions I should implement vs other functions that should just fail.
#       A better solution is likely to use traits... I may switch to that eventually...

for func in [:typemin, :typemax]
  @eval $func{T<:AbstractUnum}(::Type{T}) = error("Not implemented for type $T")
end

for func in [:isapprox, :isnan, :isnull, :isinf, :isfinite, :isposinf, :isneginf,
              :isnegative, :iszero, :ispositive, :-, :isinteger, :issubnormal]
  @eval $func(u::AbstractUnum) = error("Not implemented for $u\!")
end

for op in [:+, :-, :*, :/, :(==), :<, :>, :<=, :>=, :isless, :isequal]
  @eval $op(u1::AbstractUnum, u2::AbstractUnum) = error("Operation not implemented for $u1 and $u2\!")
  for t in [:Bool, :Integer, :AbstractFloat]
    @eval $op(u::AbstractUnum, r::$t) = error("Operation not implemented for Unum $u and ", $t, " $r\!")
    @eval $op(r::$t, u::AbstractUnum) = error("Operation not implemented for ", $t, " $r and Unum $u\!")
  end
end

# ---------------------------------------------------------------------------------------

# B is the base of the exponent
# ESS is the "exponent size size"
# FSS is the "fraction size size"
bitstype 64 FixedUnum64{B,ESS,FSS} <: AbstractUnum{B,ESS,FSS}
typealias BinaryUnum64{ESS,FSS}   FixedUnum64{2,ESS,FSS}
typealias Unum64                  BinaryUnum64{4,5}

bitstype 16 FixedUnum16{B,ESS,FSS} <: AbstractUnum{B,ESS,FSS}
typealias BinaryUnum16{ESS,FSS}   FixedUnum16{2,ESS,FSS}
typealias Unum16                  BinaryUnum16{1,3}  # note: could also be 2,2

# ---------------------------------------------------------------------------------------

# build a unum from the component values
@generated function call{U<:AbstractUnum}(::Type{U}, isneg::Bool, ubit::Bool, es::Int, fs::Int, e::Int, f::Int)
  c = unumConstants(U)
  quote
    # shift values into the correct section, right to left
    u = $(c.UINT)(fs - 1)
    u = u | ($(c.UINT)(es - 1) << $(c.fsizepos))
    if ubit
      u = u | (one($(c.UINT)) << $(c.esizepos))
    end
    u = u | ($(c.UINT)(f) << $(c.ubitpos))
    u = u | ($(c.UINT)(e) << $(c.fpos))
    if isneg
      u = u | (one($(c.UINT)) << $(c.epos))
    end
    reinterpret(U,u)
  end
end

# some helpful aliases

# typealias Unum{E, U} FixedUnum{2, E, U}
# typealias DecimalUnum{E, U} FixedUnum{10, E, U}

# typealias Unum16 Unum{3, UData16}
# typealias Unum32 Unum{7, UData32}
# typealias Unum64 Unum{16, UData64}
# typealias Unum128 Unum{20, UData128}

# typealias DecimalUnum16 DecimalUnum{3, UData16}
# typealias DecimalUnum32 DecimalUnum{7, UData32}
# typealias DecimalUnum64 DecimalUnum{16, UData64}
# typealias DecimalUnum128 DecimalUnum{20, UData128}


text:
# MonoidAlternative
# =================

a = [1,2,3]
b = [4,5,6]

@test a ⊕ b == [a; b]

# FunctorApplicativeMonad
# =======================

@test eltype(Vector{Int}) == Int

@test map(x -> x*x, a) == [1, 4, 9]

product_ab = [5, 6, 7, 6, 7, 8, 7, 8, 9]
@test mapn(a, b) do x, y
  x + y
end == product_ab

h = @syntax_flatmap begin
  x = a
  y = b
  @pure x + y
end
@test h == product_ab

@test pure(Vector{Int}, 3) == [3]

# flatten combinations with other Monads
# --------------------------------------

@test flatten([Identity(3), Identity(4)]) == [3, 4]
@test flatten([Option(3), Option(nothing), Option(4)]) == [3, 4]
@test flatten([Try(4), (@Try error("hi")), Try(5)]) == [4, 5]
@test flatten([Either{String}(4), either("left", false, 3), Either{Int, String}("right")]) == [4, "right"]

h1 = @syntax_flatmap begin
 a, b = [(1, 2), (2, 3), (3, 4), (4, 5), (5, 6)]
 iftrue(a % 2 == 0) do
   a + b
 end
end
@test h1 == [5, 9]

# FlipTypes
# ========

v = [Option(:a),
     Option(:b),
     Option(:c)]

@test flip_types(v) == Option([:a, :b, :c])


text:
type SimplePopulationDefinition <: PopulationDefinition
  size::Int
  offspring::Int
  species::SpeciesDefinition
  breeder::BreederDefinition

  SimplePopulationDefinition() =
    new(80, 0)
  SimplePopulationDefinition(size::Int,
    offspring::Int,
    species::species.SpeciesDefinition,
    breeder::breeder.BreederDefinition
  ) =
    new(size, offspring, species, breeder)
end

"""
Composes a simple population from a provided definition.
"""
function compose!(def::SimplePopulationDefinition)
  dd = deme(def.size, def.offspring, def.species, def.breeder)
  dm = _deme_.compose!(dd)
  Population(Deme[dm])
end

"""
Simple populations consist of a single fixed-size deme containing individuals of
the same species.

**Properties:**

* `size::Int`, the number of individuals within the population.
* `offspring::Int`, the number of offspring that should be produced at each
  generation.
* `species::Species`, the species to which all individuals within this
  population belong.
* `breeder::Breeder`, the breeder used to generate the offspring for this
  population at each generation.
"""
function simple(def::Function)
  d = SimplePopulationDefinition()
  def(d)
  d
end


text:
using MatFac, Test, Flux, Zygote, CSV, DataFrames

MF = MatFac

logistic(x) = 1 / (1 + exp(-x))

function noise_model_tests()

    @testset "Noise models" begin
        M = 20
        N = 40

        # Normal noise model
        normal_data = ones(M,N)
        normal_Z = zeros(M,N)
        nn = MF.NormalNoise()
        @test MF.invlink(nn, normal_Z) == normal_Z
        @test MF.loss(nn, MF.invlink(nn, normal_Z), normal_data) == 0.5.*ones(M,N)
        @test Flux.gradient(x->sum(MF.loss(nn, MF.invlink(nn, x),normal_data)), normal_Z)[1] == -ones(M,N)
        @test Flux.gradient(x->MF.invlinkloss(nn, x,normal_data), normal_Z)[1] == -ones(M,N)

        # Logistic noise model
        logistic_data = ones(M,N)
        logistic_Z = zeros(M,N)
        ln = MF.BernoulliNoise()
        @test MF.invlink(ln, logistic_Z) == 0.5 .* ones(M,N)
        @test Flux.gradient(x->sum(MF.invlink(ln, x)), logistic_Z)[1] == MF.invlink(ln,logistic_Z).*(1 .- MF.invlink(ln, logistic_Z))
        @test isapprox(MF.loss(ln, MF.invlink(ln, logistic_Z), logistic_data), log(2.0).*ones(M,N), atol=1e-6)
        logistic_A = MF.invlink(ln, logistic_Z)
        @test Flux.gradient(x->sum(MF.loss(ln, x, logistic_data)), logistic_A)[1] == -2.0.*ones(M,N)
        @test Flux.gradient(x->MF.invlinkloss(ln, x, logistic_data), logistic_Z)[1] == (logistic_A .- logistic_data)

        # Poisson noise model
        poisson_data = ones(M,N)
        poisson_Z = zeros(M,N)
        pn = MF.PoissonNoise()
        @test MF.invlink(pn, poisson_Z) == ones(M,N)
        @test isapprox(MF.loss(pn, MF.invlink(pn, poisson_Z), poisson_data), ones(M,N), atol=1e-6)
        @test Flux.gradient(x->sum(MF.loss(pn, MF.invlink(pn, x), logistic_data)), poisson_Z)[1] == zeros(M,N)
        @test Flux.gradient(x->MF.invlinkloss(pn, x, logistic_data), poisson_Z)[1] == zeros(M,N)

        # Ordinal noise model
        ordinal_data = [1 2 3;
                        1 2 3]
        ordinal_Z = [0 0 0;
                     0 0 0]
        on = MF.OrdinalNoise([-Inf, -1, 1, Inf])
        @test MF.invlink(on, ordinal_Z) == ordinal_Z
        @test isapprox(MF.loss(on, ordinal_Z, ordinal_data), [-log.(logistic(-1)-logistic(-Inf)) -log.(logistic(1)-logistic(-1)) -log.(logistic(Inf)-logistic(1));
                                                               -log.(logistic(-1)-logistic(-Inf)) -log.(logistic(1)-logistic(-1)) -log.(logistic(Inf)-logistic(1))],
                       atol=1e-6)
        thresh_grad, z_grad = Flux.gradient((noise, x)->sum(MF.loss(noise, MF.invlink(noise, x), ordinal_data)), on, ordinal_Z)

        lgrad(lt,rt,z) = logistic(lt-z)*(1 - logistic(lt-z))/(logistic(rt-z)-logistic(lt-z))
        rgrad(lt,rt,z) = -logistic(rt-z)*(1 - logistic(rt-z))/(logistic(rt-z)-logistic(lt-z))
        test_thresh_grad = [2*lgrad(-Inf,-1,0),
                              2*(rgrad(-Inf,-1,0)+lgrad(-1,1,0)),
                              2*(rgrad(-1,1,0)+lgrad(1,Inf,0)),
                              2*rgrad(1,Inf,0)]
        test_z_grad = [1-logistic(-1)-logistic(-Inf) 1-logistic(1)-logistic(-1) 1-logistic(Inf)-logistic(1);
                         1-logistic(-1)-logistic(-Inf) 1-logistic(1)-logistic(-1) 1-logistic(Inf)-logistic(1)]

        @test thresh_grad.ext_thresholds == test_thresh_grad
        @test z_grad == test_z_grad

        thresh_grad, z_grad = Flux.gradient((noise,x)->MF.invlinkloss(noise, x, ordinal_data), on, ordinal_Z)
        @test thresh_grad.ext_thresholds == test_thresh_grad
        @test z_grad == test_z_grad

        # Composite noise models
        composite_Z = zeros(M,N)
        composite_data = zeros(M,N)
        composite_data[:,21:30] .= 1
        composite_data[:,31:40] .= 2

        noise_model = repeat(["normal","bernoulli","poisson", "ordinal3"], inner=10)
        cn = MF.CompositeNoise(noise_model)

        composite_A = MF.invlink(cn, composite_Z)
        @test composite_A[:,1:10] == MF.invlink(nn, composite_Z[:,1:10])
        @test composite_A[:,11:20]== MF.invlink(ln, composite_Z[:,11:20])
        @test composite_A[:,21:30]== MF.invlink(pn, composite_Z[:,21:30])
        @test composite_A[:,31:40]== MF.invlink(cn.noises[4], composite_Z[:,31:40])

        composite_l = MF.loss(cn, composite_A, composite_data)
        @test composite_l[:,1:10] == MF.loss(nn, composite_A[:,1:10] , composite_data[:,1:10])
        @test composite_l[:,11:20]== MF.loss(ln, composite_A[:,11:20], composite_data[:,11:20])
        @test composite_l[:,21:30]== MF.loss(pn, composite_A[:,21:30], composite_data[:,21:30])
        @test composite_l[:,31:40]== MF.loss(cn.noises[4], composite_A[:,31:40], composite_data[:,31:40])

        grad_cn, grad_Z = Flux.gradient((noise,x)->MF.invlinkloss(noise, x, composite_data),
                                        cn, composite_Z)
        test_grad_Z = zeros(M,N)
        test_grad_Z[:,11:20] .= 0.5
        test_grad_ordinal = [0.0, 200.0*logistic(-0.5)*(1 - logistic(-0.5))/(logistic(0.5) - logistic(-0.5)),
                             -200.0*logistic(0.5)*(1 - logistic(0.5))/(logistic(0.5) - logistic(-0.5)) ,0.0]
        @test isapprox(grad_Z, test_grad_Z)
        @test isapprox(grad_cn.noises[4].ext_thresholds, test_grad_ordinal)

        # GPU tests
        cn_d = gpu(cn)
        composite_Z_d = gpu(composite_Z)
        composite_data_d = gpu(composite_data)
        test_grad_Z_d = gpu(test_grad_Z)
        test_grad_ordinal_d = gpu(test_grad_ordinal)
        grad_cn_d, grad_Z_d = Flux.gradient((noise,x)->MF.invlinkloss(noise, x, composite_data_d),
                                            cn_d, composite_Z_d)
        @test isapprox(grad_Z_d, test_grad_Z_d)
        @test isapprox(grad_cn_d.noises[4].ext_thresholds, test_grad_ordinal_d)
    end

end

function model_tests()

    @testset "Model" begin

        M = 20
        N = 40
        K = 5
        n_loss_types = 4
        #n_row_batches = 4
        n_logistic = div(N,4)
        n_ordinal = div(N,4)
        n_poisson = div(N,4)
        n_normal = N - (n_logistic + n_ordinal + n_poisson)

        col_losses = [repeat(["bernoulli"], n_logistic);
                      repeat(["ordinal5"], n_ordinal);
                      repeat(["poisson"], n_poisson);
                      repeat(["normal"], n_normal)]

        model = MF.MatFacModel(M,N,K,col_losses)

        @test size(model.X) == (K,M)
        @test map(typeof, model.noise_model.noises) == (MF.BernoulliNoise, MF.OrdinalNoise, MF.PoissonNoise, MF.NormalNoise)
        @test size(model()) == (M,N)
        @test size(model) == (M,N)

        # GPU
        model_d = gpu(model)
        @test isapprox(transpose(model_d.X)*model_d.Y, gpu(transpose(model.X)*model.Y))
        @test isapprox(model_d(), gpu(model()))
    end
end

function update_tests()

    @testset "Update" begin

        # Array update
        A = ones(10,10)
        start_A = deepcopy(A)

        f = x -> sum(x.^2)
        grads = Zygote.gradient(f, A)
        opt = Flux.Optimise.ADAGrad()

        Flux.Optimise.update!(opt, A, grads[1])
        @test !isapprox(A, start_A)

        # Tuple update
        B = (ones(10,10), ones(5,5))
        start_B = deepcopy(B)
        B_grad = (-ones(10,10), -ones(5,5))

        Flux.Optimise.update!(opt, B, B_grad)
        @test !isapprox(B[1], start_B[1])
        @test !isapprox(B[2], start_B[2])

        # NamedTuple update
        C = (cat=ones(10,10), dog=ones(5,5))
        start_C = deepcopy(C)
        C_grad = (cat=-ones(10,10), dog=-ones(5,5), fish=nothing)

        Flux.Optimise.update!(opt, C, C_grad)
        @test !isapprox(C[1], start_C[1])
        @test !isapprox(C[2], start_C[2])

    end
end

function fit_tests()

    @testset "Fit" begin

        M = 20
        N = 40
        K = 3
        n_loss_types = 4
        n_logistic = div(N,n_loss_types)
        n_ordinal = div(N,n_loss_types)
        n_poisson = div(N,n_loss_types)
        n_normal = div(N,n_loss_types)

        col_losses = [repeat(["bernoulli"], n_logistic);
                      repeat(["normal"], n_normal);
                      repeat(["poisson"], n_poisson);
                      repeat(["ordinal5"], n_ordinal)];

        composite_data = zeros(M,N)
        composite_data[:,21:30] .= 1
        composite_data[:,31:40] .= 3

        #################################
        # CPU TESTS
        model = MF.MatFacModel(M,N,K, col_losses)
        X_start = deepcopy(model.X)
        Y_start = deepcopy(model.Y)
        thresholds_start = deepcopy(model.noise_model.noises[4].ext_thresholds)

        # test whether the fit! function can run to
        # completion (under max_epochs condition)
        fit!(model, composite_data; verbosity=1, lr=0.05, max_epochs=10)
        @test true

        # test whether the parameters were modified
        @test !isapprox(model.X, X_start)
        @test !isapprox(model.Y, Y_start)
        @test !isapprox(model.noise_model.noises[4].ext_thresholds, thresholds_start)

        #################################
        # GPU TESTS
        model = MF.MatFacModel(M,N,K, col_losses)
        model = gpu(model)
        composite_data = gpu(composite_data)
        X_start = deepcopy(model.X)
        Y_start = deepcopy(model.Y)
        thresholds_start = deepcopy(model.noise_model.noises[4].ext_thresholds)

        # test whether the fit! function can run to
        # completion (under max_epochs condition)
        fit!(model, composite_data; verbosity=1, lr=0.05, max_epochs=10)
        @test true

        # test whether the parameters were modified
        @test !isapprox(model.X, X_start)
        @test !isapprox(model.Y, Y_start)
        @test !isapprox(cpu(model.noise_model.noises[4].ext_thresholds), cpu(thresholds_start))

    end
end

function io_tests()

    @testset "IO tests" begin

        test_bson_path = "model_test.bson"

        M = 20
        N = 40
        K = 3
        n_loss_types = 4
        n_logistic = div(N,n_loss_types)
        n_ordinal = div(N,n_loss_types)
        n_poisson = div(N,n_loss_types)
        n_normal = div(N,n_loss_types)

        col_losses = [repeat(["bernoulli"], n_logistic);
                      repeat(["normal"], n_normal);
                      repeat(["poisson"], n_poisson);
                      repeat(["ordinal5"], n_ordinal)];
        model = MF.MatFacModel(M,N,K, col_losses)

        MF.save_model(test_bson_path, model)

        new_model = MF.load_model(test_bson_path)

        @test new_model.X == model.X
        @test new_model.Y == model.Y
        @test new_model.noise_model.noises[4].ext_thresholds == model.noise_model.noises[4].ext_thresholds

        rm(test_bson_path)
    end
end

function main()

    noise_model_tests()
    model_tests()
    update_tests()
    fit_tests()
    io_tests()

end

main()

max_stars_repo_name: dpmerrell/BatchMatFac

max_issues_repo_name: dpmerrell/BatchMatFac

max_forks_repo_name: dpmerrell/BatchMatFac


text:
module MyFFT

export myfft, myfftClearCache, myifft, myirealfft, myrealfft

_myBluesteinCache = Dict{(Type, Integer), AbstractArray}()

function _myBluestein{F<:Real}(x::AbstractArray{Complex{F}, 1}, d::Integer)
    n = size(x)[1]
    m = nextpow2((n << 1) - 1)
    a = zeros(Complex{F}, m)
    b = zeros(Complex{F}, m)
    c = Array(Complex{F}, m)
    c[1] = 1 + 0.0im

    for i = 1:m-1
        c[i+1] = cis(convert(F, 1) * i * i * pi / d)
    end

    a[1:n] = x[1:n] .* c[1:n]

    if haskey(_myBluesteinCache, (Complex{F}, d))
        b_fft = _myBluesteinCache[(Complex{F}, d)]
    else
        b[1:n] = conj(c[1:n])
        i = 2
        j = m

        while i <= n
            b[j] = b[i]
            i = i + 1
            j = j - 1
        end

        b_fft = myfft(b)
        _myBluesteinCache[(Complex{F}, d)] = b_fft
    end

    copy((c .* myifft(myfft(a) .* b_fft))[1:n])
end

function _myfftTask!{F<:Real, R<:Range, Rr<:Range, I<:Integer}(
    x::AbstractArray{Complex{F}, 1}, r::R,
    result::AbstractArray{Complex{F}, 1}, rResult::Rr, n::I,
    twiddleBasis::Complex{F})

    if n > 2
        nHalf = n >> 1

        # call sub-tasks
        s = step(r)
        sDouble = s << 1
        rEven = (r.start):(sDouble):(r.stop-s)
        rOdd = (r.start+s):(sDouble):(r.stop)
        rrEven = (rResult.start):(rResult.start+nHalf-1)
        rrOdd = (rResult.start+nHalf):(rResult.stop)
        subTwiddleBasis = twiddleBasis * twiddleBasis
        _myfftTask!(x, rEven, result, rrEven, nHalf, subTwiddleBasis)
        _myfftTask!(x, rOdd, result, rrOdd, nHalf, subTwiddleBasis)

        # core calculation
        twiddle = twiddleBasis
        i = rrEven.start
        j = rrOdd.start
        xEven = result[i]
        xOdd = result[j]
        result[i] = xEven + xOdd
        result[j] = xEven - xOdd

        while i < rrEven.stop
            i = i + 1
            j = j + 1
            xEven = result[i]
            xOdd = twiddle * result[j]
            result[i] = xEven + xOdd
            result[j] = xEven - xOdd
            twiddle = twiddle * twiddleBasis
        end
    elseif n == 2
        xEven = x[r.start]
        xOdd = x[r.stop]
        result[rResult.start] = xEven + xOdd
        result[rResult.stop] = xEven - xOdd
    else
        result[rResult.start] = x[r.start]
    end

    nothing
end

function myfft{F<:Real}(x::AbstractArray{Complex{F}, 1})
    n = size(x)[1]

    if n == 0
        throw(ArgumentError)
    elseif n & (n-1) == 0
        r = 1:n
        twiddleBasis = cis(convert(F, -2) * pi / n)
        result = Array(Complex{F}, n)
        _myfftTask!(x, r, result, r, n, twiddleBasis)
    else
        result = _myBluestein(x, -n)
    end

    result
end

function myfftClearCache()
    _myBluesteinCache = Dict{(Type, Integer), AbstractArray}()
end

function myifft{F<:Real}(x::AbstractArray{Complex{F}, 1})
    n = size(x)[1]

    if n == 0
        throw(ArgumentError)
    elseif n & (n-1) == 0
        r = 1:n
        twiddleBasis = cis(convert(F, 2) * pi / n)
        result = Array(Complex{F}, n)
        _myfftTask!(x, r, result, r, n, twiddleBasis)
    else
        result = _myBluestein(x, n)
    end

    result/n
end

function myirealfft{F<:Real}(x::AbstractArray{Complex{F}, 1})
    n = size(x)[1]

    if n == 0
        throw(ArgumentError)
    elseif n % 2 != 0
        return real(ifft(x))
    end

    nHalf = n >> 1
    z = Array(Complex{F}, nHalf)
    twiddleBasis = cis(convert(F, 2) * pi / n)
    twiddle = twiddleBasis
    p = x[1]
    q = x[nHalf+1]
    xEven = 0.5 * (p + q)
    xOdd = 0.5 * (p - q)
    z[1] = xEven + 1im * xOdd

    for i = 1:nHalf-1
        p = x[i+1]
        q = x[nHalf+i+1]
        xEven = 0.5 * (p + q)
        xOdd = 0.5 * (p - q) * twiddle
        z[i+1] = xEven + 1im * xOdd
        twiddle = twiddle * twiddleBasis
    end

    z_ifft = myifft(z)
    result = Array(F, n)
    result[1:2:n-1] = real(z_ifft)
    result[2:2:n] = imag(z_ifft)
    result
end

function myrealfft{F<:Real}(x::AbstractArray{F, 1})
    n = size(x)[1]

    if n == 0
        throw(ArgumentError)
    elseif n % 2 != 0
        return fft(x + 0.0im)
    end

    nHalf = n >> 1
    z = x[1:2:n-1] + 1im * x[2:2:n]
    z_fft = myfft(z)
    result = Array(Complex{F}, n)
    twiddleBasis = cis(convert(F, -2) * pi / n)
    twiddle = -1im * twiddleBasis
    p = z_fft[1]
    q = conj(z_fft[1])
    result[1] = 0.5 * (p + q - 1im * (p - q))
    result[nHalf+1] = 0.5 * (p + q  + 1im * (p - q))
    i = 1
    j = nHalf - 1

    while i < nHalf
        p = z_fft[i+1]
        q = conj(z_fft[j+1])
        result[i+1] = 0.5 * (p + q  + twiddle * (p - q))
        result[n-i+1] = conj(result[i+1])
        twiddle = twiddle * twiddleBasis
        i = i + 1
        j = j - 1
    end

    result
end

end # module


text:
#!/usr/bin/env julia

using Base.Threads

function _picalc(numsteps)
    slice = 1.0 / numsteps

    sums = zeros(Float64, nthreads())
    n = cld(numsteps, nthreads())

    Threads.@threads for i in 1:nthreads()
        sum_thread = 0.0
        @simd for j in (1 + (i - 1) * n):min(numsteps, i * n)
            x = (j - 0.5) * slice
            sum_thread += 4.0 / (1.0 + x ^ 2)
        end
        sums[threadid()] = sum_thread
    end

    return sum(sums) * slice
end

function picalc(numsteps)

    println("Calculating PI using:")
    println("  ", numsteps, " slices")
    println("  ", nthreads(), " thread(s)")

    start = time()
    mypi = _picalc(numsteps)
    elapsed = time() - start

    println("Obtained value of PI: ", mypi)
    println("Time taken: ", elapsed, " seconds")

end

numsteps = if length(ARGS) > 0
    parse(Int, ARGS[1])
else
    1_000_000_000
end

# Warm up kernel
_picalc(10)

# Run the full example
picalc(numsteps)


text:
#=
Extremely simple method from paper:

https://journals.aps.org/pre/abstract/10.1103/PhysRevE.47.2401

where a delay embedding is done, and the corrected timeseries is in fact just
an average of nearest neigbhors in the embedded space
=#
using Neighborhood
using DelayEmbeddings
export ExtremelySimpleNL

"""
    ExtremelySimpleNL(k::Int, ℓ::Int, w::Int, τ::Int, ε::Real) <: Decomposition
Quite literally the "extremely simple nonlinear noise-reduction method"[^Schreiber1993].
It decomposes `s` into the **sum** `x + r` with `x` being the "noiseless" timeseries.
This is the average position of neighbors in the delay embedded space.

This method works well if your timeseries is composed by the
addition of a structured component (which follows deterministic and
stationary dynamics which the embedding should approximate) and some noise.

The cited paper has some info on choosing optimal `ε`.

Arguments:
* `k` amount of past delay
* `ℓ` amount of forward delay
* `w` Theiler window
* `τ` delay time
* `ε` radius of the neighborhood in the embedded space

[^Schreiber1993]: [Schreiber, (1993) Extremely simple nonlinear noise-reduction method. Physical Review E, 47(4)](https://doi.org/10.1103/PhysRevE.47.2401)
"""
struct ExtremelySimpleNL <: Decomposition
    k::Int
    ℓ::Int
    w::Int
    τ::Int
    ε::Float64
end

function decompose(t, s, method::ExtremelySimpleNL)
    k, ℓ, w, τ, ε = getproperty.(Ref(method), (:k, :ℓ, :w, :τ, :ε))
    delays = [i*τ for i in -k:1:ℓ]
    embedding = GeneralizedEmbedding(Tuple(delays))
    X = genembed(s, Tuple(delays))
    tree = searchstructure(KDTree, X.data, Chebyshev())
    c = copy(s) # corrected timeseries
    theiler = Theiler(w)

    vec_of_idxs = bulkisearch(tree, X.data, WithinRange(ε), theiler)
    zc = 0

    for (i, n) in enumerate(τrange(s, embedding))
        𝒰n = vec_of_idxs[i]
        K = length(𝒰n)
        if K == 0
            zc += 1
            continue
        end
        c[n] = (1/K) * sum(s[j+k] for j in 𝒰n)
    end
    zc > 0 && @warn "$(zc) points were not corrected because no neighbors where found."
    return c, s .- c
end


text:
# least squares penalty

export LeastSquares

### ABSTRACT TYPE

abstract type LeastSquares <: ProximableFunction end

is_smooth(f::LeastSquares) = true
is_quadratic(f::LeastSquares) = true

fun_name(f::LeastSquares) = "Least squares penalty"

### CONSTRUCTORS

"""
**Least squares penalty**

    LeastSquares(A, b, λ=1.0; iterative=false)

For a matrix `A`, a vector `b` and a scalar `λ`, returns the function
```math
f(x) = \\tfrac{\\lambda}{2}\\|Ax - b\\|^2.
```
By default, a direct method (based on Cholesky factorization) is used to evaluate `prox!`.
If `iterative=true`, then `prox!` is evaluated approximately using an iterative method instead.
"""
function LeastSquares(A::M, b::V, lam::R=R(1); iterative=false) where {R <: Real, RC <: RealOrComplex{R}, V <: AbstractArray{RC}, M}
    if iterative == false
        LeastSquaresDirect(A, b, lam)
    else
        LeastSquaresIterative(A, b, lam)
    end
end

infer_shape_of_x(A, ::AbstractVector) = (size(A, 2), )
infer_shape_of_x(A, b::AbstractMatrix) = (size(A, 2), size(b, 2))

### INCLUDE CONCRETE TYPES

include("leastSquaresDirect.jl")
include("leastSquaresIterative.jl")


text:
using GridInterpolations
"""
	dimensions of state and action spaces
"""
const state_dim = 5
const action_dim = 2

"""
	state discretization
	s = [ρ, θ, ϕ, v_own, v_int]
"""
const rho_ind   = 1
const theta_ind = 2
const phi_ind   = 3
const v_own_ind = 4
const v_int_ind = 5

const num_states = rho_dim * theta_dim * phi_dim * v_dim * v_dim + 1

states = zeros(state_dim, num_states)
istate = 1
for iv_int in 1:v_dim
	for iv_own in 1:v_dim
		for iphi in 1:phi_dim
			for itheta in 1:theta_dim
				for irho in 1:rho_dim
					states[rho_ind, istate] = discrete_rho[irho]
					states[theta_ind, istate] = discrete_theta[itheta]
					states[phi_ind, istate] = discrete_phi[iphi]
					states[v_own_ind, istate] = discrete_v[iv_own]
					states[v_int_ind, istate] = discrete_v[iv_int]
					istate += 1
				end
			end
		end
	end
end

const terminal_state_var = 1e5
const terminal_state = terminal_state_var * ones(state_dim)
states[:, end] = terminal_state

const grid_states = RectangleGrid(sort(unique(vec(states[rho_ind,   1:end - 1]))),
								  sort(unique(vec(states[theta_ind, 1:end - 1]))),
                      			  sort(unique(vec(states[phi_ind,   1:end - 1]))),
                      			  sort(unique(vec(states[v_own_ind, 1:end - 1]))),
                      			  sort(unique(vec(states[v_int_ind, 1:end - 1]))))

"""
	action space (individual action)
	non-cooperative, the ownship has no knowledge about the action of the intruder
"""
const advisories = [:SR, :WR, :KEEP, :WL, :SL, :COC]
const advisory_to_ind_dict = Dict(:SR => 1, :WR => 2, :KEEP => 3, :WL => 4,
								:SL => 5, :COC => 6)
const advisory_to_action_dict = Dict(:SR => -10.0, :WR => -5.0, :KEEP => 0.0,
									:WL => 5.0, :SL => 10.0, :COC => -1.0)
const num_actions = length(advisories)

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
@testset "Kalman filter system test" begin
    start_pt = 19
    start_vel = 2
    start_acc = 1
    σ_acc_noise = 0.0001
    σ_meas_noise = 0.25
    Δt = 0.1

    function measure(state, Δt, σ_meas_noise, σ_acc_noise)
        s, v, a = state[1], state[2], state[3]
        a_next = a + randn() * σ_acc_noise
        v_next = v + a * Δt
        s_next = a * Δt^2 / 2 + v * Δt + s
        s_next + randn() * σ_meas_noise, (s_next, v_next, a_next)
    end

    # State space matrices of discretized white noise acceleration model
    F = [1 Δt Δt^2/2; 0 1 Δt; 0 0 1]
    H = [1, 0, 0]'
    Q = [Δt^2/2; Δt; 1] * [Δt^2/2 Δt 1] * σ_acc_noise^2
    R = σ_meas_noise^2

    # Initialization
    maxiter = 20000
    counter = 1
    x_init = [0.0, 0.0, 0.0]
    P_init = [2.5 0.25 0.1; 0.25 2.5 0.2; 0.1 0.2 2.5]
    ỹ_over_time = Vector{Float64}(undef, maxiter)
    S_over_time = Vector{Float64}(undef, maxiter)
    states = (start_pt, start_vel, start_acc)
    s_over_time = Vector{Float64}(undef, maxiter)
    z_over_time = Vector{Float64}(undef, maxiter)
    s̃_over_time = Vector{Float64}(undef, maxiter)
    # run Kalman Filter
    tu = time_update(x_init, P_init, F, Q)
    for i = 1:maxiter
        measurement, states = measure(states, Δt, σ_meas_noise, σ_acc_noise)
        mu = measurement_update(get_state(tu), get_covariance(tu), measurement, H, R)
        tu = time_update(get_state(mu), get_covariance(mu), F, Q)

        #s̃_over_time[counter] = state(mu)[1]
        #s_over_time[counter] = states[1]
        #z_over_time[counter] = measurement
        ỹ_over_time[counter] = get_innovation(mu)
        S_over_time[counter] = get_innovation_covariance(mu)
        counter += 1
    end
    #using PyPlot
    #figure()
    #plot(s_over_time)
    #plot(z_over_time)
    #plot(s̃_over_time)
    #plot(ỹ_over_time)
    #plot(s_over_time .- s̃_over_time)
    #legend(("Ground truth", "Messung", "Schätzung", "Innovation", "Error"))

    # Statistical consistency testing
    @test sigma_bound_test(ỹ_over_time[4:end], S_over_time[4:end]) == true
    @test two_sigma_bound_test(ỹ_over_time[4:end], S_over_time[4:end]) == true

    window_start = 4
    window_length = 400
    window = window_start:window_start + window_length - 1
    dof = length(window) * size(ỹ_over_time[window_start], 1)
    nis_over_time_sys = map((x, σ²) -> calc_nis(x, σ²), ỹ_over_time[window], S_over_time[window])
    result_nis_test = nis_test(nis_over_time_sys, dof)
    @test result_nis_test == true
end