text:
using StochasticDiffEq, DiffEqNoiseProcess

@test_broken begin
  # Need to internally broadcast!
  oup = OrnsteinUhlenbeckProcess!(1.0, 0.0, 1.0, 0.0, [0.0], [0.0])
  noise = NoiseWrapper(oup)
  prob = SDEProblem((u,p,t) -> -u, (u,p,t) -> [1.0], [0.0], (0.0, 1.0), noise=noise)
  sol = solve(prob,SRIW1())
  oup = OrnsteinUhlenbeckProcess!(1.0, 0.0, 1.0, 0.0, [0.0], [0.0])
  noise = NoiseWrapper(oup)
  prob = SDEProblem((u,p,t) -> -u, (u,p,t) -> [1.0], [0.0], (0.0, 1.0), noise=noise)
  sol = solve(prob,SOSRI())
  oup = OrnsteinUhlenbeckProcess!(1.0, 0.0, 1.0, 0.0, [0.0], [0.0])
  noise = NoiseWrapper(oup)
  prob = SDEProblem((u,p,t) -> -u, (u,p,t) -> [1.0], [0.0], (0.0, 1.0), noise=noise)
  @test_broken sol = solve(prob,SRA1())
  oup = OrnsteinUhlenbeckProcess!(1.0, 0.0, 1.0, 0.0, [0.0], [0.0])
  noise = NoiseWrapper(oup)
  prob = SDEProblem((u,p,t) -> -u, (u,p,t) -> [1.0], [0.0], (0.0, 1.0), noise=noise)
  @test_broken sol = solve(prob,SRA2())
  oup = OrnsteinUhlenbeckProcess!(1.0, 0.0, 1.0, 0.0, [0.0], [0.0])
  noise = NoiseWrapper(oup)
  prob = SDEProblem((u,p,t) -> -u, (u,p,t) -> [1.0], [0.0], (0.0, 1.0), noise=noise)
  @test_broken sol = solve(prob,SRA3())

  println("Scalar g")
  A = [-1.0 0.0; 0.0 -0.5]
  u0 = [1.0, 1.0]; tspan = (0.0,1.0)
  _f = (u,p,t) -> t*(A*u)
  _g = (u,p,t) -> 1.0
  prob = SDEProblem(SDEFunction(_f, _g), _g, u0, tspan)
  integrator = init(prob, SKenCarp(); adaptive=false, dt=0.01)
  step!(integrator)

  println("Vector g")
  _g = (u,p,t) -> [1.0, 1.0]
  prob = SDEProblem(SDEFunction(_f, _g), _g, u0, tspan)
  println("Implicit EM")
  integrator = init(prob, ImplicitEM(); adaptive=false, dt=0.01)
  step!(integrator)
  println("SKenCarp")
  integrator = init(prob, SKenCarp(); adaptive=false, dt=0.01)
  step!(integrator)
  solve(prob, SOSRI(); adaptive=false, dt=0.01)
  solve(prob, SOSRA(); adaptive=false, dt=0.01)
  solve(prob, EM(); adaptive=false, dt=0.01)
  solve(prob, RKMil(); adaptive=false, dt=0.01)
  solve(prob, SOSRI2(); adaptive=false, dt=0.01)
  solve(prob, SOSRA2(); adaptive=false, dt=0.01)
end


text:
# This file was generated, do not modify it. # hide
const dim = 16
k = 15
db = MatrixDatabase(randn(Float32, dim, 10^5))
dist = SqL2Distance()
G = SearchGraph(; dist, db)
opt = OptimizeParameters(kind=MinRecall(0.9))
itime = @elapsed (index!(G); optimize!(G, opt))
nothing # hide

ext: jl

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
"""
$(DocStringExtensions.README)
"""
module DataDrivenDiffEq

using DocStringExtensions
using LinearAlgebra
using DiffEqBase
using ModelingToolkit
using QuadGK
using Statistics
using DSP
using FiniteDifferences, DataInterpolations
using Compat
using DocStringExtensions

abstract type AbstractKoopmanOperator <: Function end;

include("./optimizers/Optimize.jl")
using .Optimize

export set_threshold!, set_threshold
export STRRidge, ADMM, SR3

export ADM

include("./basis.jl")
export Basis
export dynamics
export free_parameters

include("./koopman/algorithms.jl")
export DMDPINV, DMDSVD, TOTALDMD

include("./koopman/koopman.jl")
export eigen, eigvals, eigvecs
export modes, frequencies
export is_discrete, is_continuous
export operator, generator
export inputmap, outputmap, updatable, isstable

include("./koopman/linearkoopman.jl")
export LinearKoopman, update!

include("./koopman/nonlinearkoopman.jl")
export NonlinearKoopman, reduce_basis

include("./koopman/exact_dmd.jl")
export DMD, gDMD

include("./koopman/dmdc.jl")
export DMDc, gDMDc

include("./koopman/extended_dmd.jl")
export EDMD, gEDMD

include("./sindy/results.jl")
export SparseIdentificationResult
export print_equations
export get_coefficients, get_error, get_sparsity, get_aicc

include("./sindy/sindy.jl")
export SINDy
export sparse_regression, sparse_regression!

function SInDy(Y, X, basis; opt = STRRidge(), kwargs...)
    @warn("SInDy has been deprecated. Use SINDy to recover the same functionality.")
    SINDy(Y, X, basis, opt; kwargs...)
end

function ISInDy(Y, X, basis; opt = ADM(), kwargs...)
    @warn("ISInDy has been deprecated. Use ISINDy to recover the same functionality.")
    ISINDy(Y, X, basis, opt; kwargs...)
end

export SInDy, ISInDy

include("./sindy/isindy.jl")
export ISINDy

include("./system_conversions.jl")

include("./utils.jl")
export AIC, AICC, BIC
export optimal_shrinkage, optimal_shrinkage!
export savitzky_golay
export burst_sampling, subsample

include("./basis_generators.jl")
export chebyshev_basis, monomial_basis, polynomial_basis
export sin_basis, cos_basis, fourier_basis

end # module


text:
function triangulate(polygon::Vector{Vector{Point{2, Float64}}})
    lengths = map(x-> UInt32(length(x)), polygon)
    len = UInt32(length(lengths))
    array = ccall(
        (:u32_triangulate_f64, earcut),
        Tuple{Ptr{GLTriangle}, Cint},
        (Ptr{Ptr{Float64}}, Ptr{UInt32}, UInt32),
        polygon, lengths, len
    )
    unsafe_wrap(Vector{GLTriangle}, array[1], array[2])
end

function triangulate(polygon::Vector{Vector{Point{2, Float32}}})
    lengths = map(x-> UInt32(length(x)), polygon)
    len = UInt32(length(lengths))
    array = ccall(
        (:u32_triangulate_f32, earcut),
        Tuple{Ptr{GLTriangle}, Cint},
        (Ptr{Ptr{Float32}}, Ptr{UInt32}, UInt32),
        polygon, lengths, len
    )
    unsafe_wrap(Vector{GLTriangle}, array[1], array[2])
end

function triangulate(polygon::Vector{Vector{Point{2, Int64}}})
    lengths = map(x-> UInt32(length(x)), polygon)
    len = UInt32(length(lengths))
    array = ccall(
        (:u32_triangulate_i64, earcut),
        Tuple{Ptr{GLTriangle}, Cint},
        (Ptr{Ptr{Int64}}, Ptr{UInt32}, UInt32),
        polygon, lengths, len
    )
    unsafe_wrap(Vector{GLTriangle}, array[1], array[2])
end

function triangulate(polygon::Vector{Vector{Point{2, Int32}}})
    lengths = map(x-> UInt32(length(x)), polygon)
    len = UInt32(length(lengths))
    array = ccall(
        (:u32_triangulate_i32, earcut),
        Tuple{Ptr{GLTriangle}, Cint},
        (Ptr{Ptr{Int32}}, Ptr{UInt32}, UInt32),
        polygon, lengths, len
    )
    unsafe_wrap(Vector{GLTriangle}, array[1], array[2])
end


text:
export IAAFT

include("../utils/powerspectrum.jl")
include("../utils/interpolation.jl")
export IAAFT

"""
    IAAFT(M = 100, tol = 1e-6, W = 75)

An iteratively adjusted amplitude-adjusted-fourier-transform surrogate[^SchreiberSchmitz1996].

IAAFT surrogates have the same linear correlation, or periodogram, and also
preserves the amplitude distribution of the original data, but are improved relative
to AAFT through iterative adjustment (which runs for a maximum of `M` steps).
During the iterative adjustment, the periodograms of the original signal and the
surrogate are coarse-grained and the powers are averaged over `W` equal-width
frequency bins. The iteration procedure ends when the relative deviation
between the periodograms is less than `tol` (or when `M` is reached).

IAAFT, just as AAFT, can be used to test the null hypothesis that the data
come from a monotonic nonlinear transformation of a linear Gaussian process.

[^SchreiberSchmitz1996]: T. Schreiber; A. Schmitz (1996). "Improved Surrogate Data for Nonlinearity Tests". [Phys. Rev. Lett. 77 (4)](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.77.635)
"""
struct IAAFT <: Surrogate
    M::Int
    tol::Real
    W::Int

    function IAAFT(;M::Int = 100, tol::Real = 1e-6, W::Int = 75)
        new(M, tol, W)
    end
end

Base.show(io::IO, x::IAAFT) = print(io, "IAAFT(M = $(x.M), tol = $(x.tol), W = $(x.W))")

function interpolated_spectrum(spectrum, n)
    intp_spectrum = zeros(n)
    interpolated_spectrum!(intp_spectrum, spectrum, n)
end

function interpolated_spectrum!(intp_spectrum, spectrum, n)
    t = 1:length(spectrum)
    r = getrange(t, n)
    iₓ = itp(spectrum)
    intp_spectrum .= iₓ.(r)
    return intp_spectrum
end

function surrogenerator(x, method::IAAFT, rng = Random.default_rng())
    m = mean(x)
    x_sorted = sort(x)
    forward = plan_rfft(x)
    inverse = plan_irfft(forward * x, length(x))
    𝓕 = forward * x
    r = abs.(𝓕)
    ϕ = abs.(𝓕)
    ix = zeros(Int, length(x))

    # Periodograms
    𝓕p = prepare_spectrum(x, forward)
    xpower = zeros(length(𝓕)); powerspectrum!(𝓕p, xpower, x, forward)
    spower = copy(xpower)

    # Binned periodograms
    xpowerᵦ = interpolated_spectrum(xpower, method.W)
    spowerᵦ = interpolated_spectrum(spower, method.W)

    init = (
        forward = forward,
        inverse = inverse,
        𝓕 = 𝓕,
        𝓕p = 𝓕p,
        r = r,
        ϕ = ϕ,
        m = m,
        x_sorted = x_sorted,
        xpower = xpower,
        spower = spower,
        xpowerᵦ = xpowerᵦ,
        spowerᵦ = spowerᵦ,
        ix = ix,
    )

    return SurrogateGenerator(method, x, similar(x), init, rng)
end
using LinearAlgebra
function (sg::SurrogateGenerator{<:IAAFT})()
    init_fields = (:forward, :inverse, :𝓕, :𝓕p, :r, :ϕ, :m, :x_sorted, :xpower, :spower, :xpowerᵦ, :spowerᵦ, :ix)
    forward, inverse, 𝓕, 𝓕p, r, ϕ, m, x_sorted, xpower, spower, xpowerᵦ, spowerᵦ, ix = getfield.(Ref(sg.init), init_fields)

    x, s, rng = sg.x, sg.s, sg.rng
    M, W = sg.method.M, sg.method.W
    tol = sg.method.tol

    # Surrogate starts out as a random permutation of `x`
    n = length(x)
    s .= x[sample(rng, 1:n, n)]

    sum_old, sum_new = 0.0, 0.0
    iter = 1
    while iter <= M
        mul!(𝓕, forward, s)
        ϕ .= angle.(𝓕)
        𝓕 .= r .* exp.(ϕ .* 1im)

        # TODO: Unfortunately, we can't simply do ldiv! here to avoid allocations.
        # But, although FFTW does not yet have irfft!, it will probably have.
        # See https://github.com/JuliaMath/FFTW.jl/pull/222.
        # Once that PR is merged, we should replace the following line with
        # the in-place version.
        ######################################################################
        s .= inverse * 𝓕
        sortperm!(ix, s)
        s[ix] .= x_sorted

        powerspectrum!(𝓕p, spower, s, forward)
        interpolated_spectrum!(spowerᵦ, spower, W)
        if iter == 1
            sum_old = sum((xpowerᵦ .- xpowerᵦ) .^ 2) / sum(xpowerᵦ .^ 2)
        else
            sum_new = sum((xpowerᵦ .- spowerᵦ) .^ 2) / sum(xpowerᵦ .^ 2)
            if abs(sum_old - sum_new) < tol
                iter = M + 1
            else
                sum_old = sum_new
            end
        end

        iter += 1
    end

    return s
end


text:
# one sided version of tictactoe (  | |  ), when any player has two adjatcent poles it wins
# (so basicly black can't win, and middle cell is win for white)

import AlphaZero.GI
using StaticArrays
# using Flux: onehot

const BOARD_SIDE = 5
const NUM_POSITIONS = BOARD_SIDE

const Player = Bool
const WHITE = true
const BLACK = false

const Cell = Union{Nothing, Player}
const Board = SVector{NUM_POSITIONS, Cell}
const INITIAL_BOARD = Board(repeat([nothing], NUM_POSITIONS))
const INITIAL_STATE = (board=INITIAL_BOARD, curplayer=BLACK) # BLACK starts

struct GameSpec <: GI.AbstractGameSpec end

mutable struct GameEnv <: GI.AbstractGameEnv
  board :: Board
  curplayer :: Player
end

GI.init(::GameSpec, state=INITIAL_STATE) = GameEnv(state.board, state.curplayer)

GI.spec(::GameEnv) = GameSpec()

GI.two_players(::GameSpec) = true

function GI.set_state!(g::GameEnv, state)
  g.board = state.board
  g.curplayer = state.curplayer
end

const ALIGNMENTS = [i-1:i for i in 2:BOARD_SIDE]

function has_won(g::GameEnv, player)
  any(ALIGNMENTS) do al
    all(al) do pos
      g.board[pos] == player
    end
  end
end

const ACTIONS = collect(1:NUM_POSITIONS)

GI.actions(::GameSpec) = ACTIONS

GI.actions_mask(g::GameEnv) = map(isnothing, g.board)

GI.current_state(g::GameEnv) = (board=g.board, curplayer=g.curplayer)

GI.white_playing(g::GameEnv) = g.curplayer

function terminal_white_reward(g::GameEnv)
  has_won(g, WHITE) && return 1.
  has_won(g, BLACK) && return -1.
  isempty(GI.available_actions(g)) && return 0.
  return nothing
end

GI.game_terminated(g::GameEnv) = !isnothing(terminal_white_reward(g))

function GI.white_reward(g::GameEnv)
  z = terminal_white_reward(g)
  return isnothing(z) ? 0. : z
end

function GI.play!(g::GameEnv, pos)
  g.board = setindex(g.board, g.curplayer, pos)
  g.curplayer = !g.curplayer
end

function flip_colors(board)
  flip(cell) = isnothing(cell) ? nothing : !cell
  # Inference fails when using `map`
  return @SVector Cell[flip(board[i]) for i in 1:NUM_POSITIONS]
end

function GI.vectorize_state(::GameSpec, state)
  board = state.curplayer == WHITE ? state.board : flip_colors(state.board)
  return Float32[
    board[i] == c
    for i in 1:BOARD_SIDE,
        c in [nothing, WHITE, BLACK]]
end

GI.heuristic_value(::GameEnv) = 0.


text:
function ask_oracle(oracle::PoolOracle, query_id::Int)
    return OneClassActiveLearning.ask_oracle(oracle, [query_id])[1]
end

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama


text:
# Bogumił Kamiński, 2019-03-25

using Serialization

d = Dict{Int, String}()

for i in 1:10
    d[i] = string(i)
end

if VERSION >= v"1.1"
    serialize("d.bin", d)
    d2 = deserialize("d.bin")
else
    with("d.bin", "w") do f
        serialize(f, d)
    end
    f = open("d.bin")
    d2 =deserialize("d.bin")
    close(f)
end

max_stars_repo_path: 8. IO/ioexample_1.jl

max_issues_repo_path: 8. IO/ioexample_1.jl

max_forks_repo_path: 8. IO/ioexample_1.jl


text:
#!/usr bin/env julia

#function to generate alphabet list for fasta naming purposes

function letter_list(letter_count::Int64)

    #set alphabet
    alpha_array = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", "N", "O", "P",
    "Q", "R", "S", "T", "U", "V", "W", "X", "Y", "Z"]

    alpha_array_n = length(alpha_array)

    #initialize array to return_array
    return_array = String[]

    #set iteration counter
    x = 1
    j = 1
    for i in 1:letter_count
        if j > 26
            j = 1
            x = x + 1
        end
        value = repeat(alpha_array[j], x)
        #println(string("i, j, x =", i, j, x))
        #println(string("value is: ", value))
        push!(return_array, value)
        j = j + 1
    end

    return return_array
end


text:
@testset "seek_magic_bytes" begin

    the_magic = BTCParser.to_byte_tuple(BTCParser.MAGIC)
    # const MAGIC = 0xf9be_b4d9
    bcio = BTCParser.BCIterator()
    BTCParser.seek_magic_bytes!(bcio.io, the_magic)
    @test position(bcio.io) == 0x4
    BTCParser.seek_magic_bytes!(bcio.io, the_magic)
    @test position(bcio.io) == 0x129
    BTCParser.seek_magic_bytes!(bcio.io, the_magic)
    @test position(bcio.io) == 0x208
    close(bcio)

end

@testset "seek_magic_bytes_boyer_moore" begin

    the_magic = BTCParser.to_byte_tuple(BTCParser.MAGIC)
    # const MAGIC = 0xf9be_b4d9
    bcio = BTCParser.BCIterator()
    BTCParser.seek_magic_bytes_boyer_moore!(bcio.io, the_magic)
    @test position(bcio.io) == 0x4
    BTCParser.seek_magic_bytes_boyer_moore!(bcio.io, the_magic)
    @test position(bcio.io) == 0x129
    BTCParser.seek_magic_bytes_boyer_moore!(bcio.io, the_magic)
    @test position(bcio.io) == 0x208
    close(bcio)

end

@testset "to_byte_tuple" begin

    bt = BTCParser.to_byte_tuple(0x1234_5678)
    @test length(bt) == 4
    @test bt[4] == 0x12
    @test bt[3] == 0x34
    @test bt[2] == 0x56
    @test bt[1] == 0x78

end


text:
using NeuroIO, MAT, BinDeps, Base.Test

# Unfortunately, these test files come from Black Rock Microsystems and
# are assumed not to be redistributable
testdir = joinpath(dirname(@__FILE__), "NSS Data")
if !isdir(testdir)
    parentdir = dirname(@__FILE__)
    bundlezip = joinpath(parentdir, "NSS Data.zip")
    run(download_cmd("http://cl.ly/1B3y1C0U380G/download/NSS%20Data.zip", bundlezip))
    cd(parentdir) do
        run(`unzip $bundlezip`)
    end
end

if !isfile(joinpath(testdir, "nev.mat"))
    using MATLAB
    npmkdir = joinpath(dirname(@__FILE__), "NPMK")
    if !isdir(npmkdir)
        cd(dirname(@__FILE__)) do
            run(`git clone https://github.com/BlackrockMicrosystems/NPMK.git`)
        end
    end
    @eval mat"""
    addpath(genpath($(joinpath(npmkdir, "NPMK"))))
    cd($testdir);
    NEV = openNEV('./SampleNSSData.nev');
    save('nev.mat', 'NEV');
    NS5 = openNSx('./SampleNSSData.ns5');
    save('ns5.mat', 'NS5');
    """
end

nev = readnev(open(joinpath(testdir, "SampleNSSData.nev")))
nev_mat = matread(joinpath(testdir, "nev.mat"))["NEV"]
@test nev.header.bytes_in_headers == nev_mat["MetaTags"]["HeaderOffset"]
@test nev.header.bytes_in_data_packets == nev_mat["MetaTags"]["PacketBytes"]
@test nev.header.time_resolution_of_time_stamps == nev_mat["MetaTags"]["TimeRes"]
@test nev.header.time_resolution_of_samples == nev_mat["MetaTags"]["SampleRes"]
@test nev.header.application_to_create_file == nev_mat["MetaTags"]["Application"]

info = nev_mat["ElectrodesInfo"]
spikes = nev_mat["Data"]["Spikes"]
for i = 1:length(info["ElectrodeID"])
    id = info["ElectrodeID"][i]
    ch = nev.spike_channels[id]
    @test ch.physical_connector == info["ConnectorBank"][i][1]
    @test ch.connector_pin == info["ConnectorPin"][i]
    # What is this?
    # @test ch.digitization_factor == info["DigitalFactor"][i]
    @test ch.energy_threshold == info["EnergyThreshold"][i]
    @test ch.high_threshold == info["HighThreshold"][i]
    @test ch.low_threshold == info["LowThreshold"][i]
    @test ch.number_of_sorted_units == info["Units"][i]
    @test ch.label == rstrip(string(info["ElectrodeLabel"][i]...), '\0')
    @test ch.high_pass.cutoff == info["HighFreqCorner"][i]
    @test ch.high_pass.order == info["HighFreqOrder"][i]
    @test ch.high_pass.filter_type == info["HighFilterType"][i]
    @test ch.low_pass.cutoff == info["LowFreqCorner"][i]
    @test ch.low_pass.order == info["LowFreqOrder"][i]
    @test ch.low_pass.filter_type == info["LowFilterType"][i]

    spike_index = find(spikes["Electrode"] .== id)
    @test NeuroIO.times(ch) == spikes["TimeStamp"][spike_index]/nev.header.time_resolution_of_time_stamps
    @test ch.unit_numbers == spikes["Unit"][spike_index]
    @test NeuroIO.data(ch) == spikes["Waveform"][:, spike_index]
end

ns5 = readnsx(open(joinpath(testdir, "SampleNSSData.ns5")))
ns5_mat = matread(joinpath(testdir, "ns5.mat"))["NS5"]
@test ns5.header.label == rstrip(string(ns5_mat["MetaTags"]["SamplingLabel"]...), '\0')
@test ns5.header.period == 30000/ns5_mat["MetaTags"]["SamplingFreq"]
@test ns5.header.time_resolution_of_time_stamps == ns5_mat["MetaTags"]["TimeRes"]

info = ns5_mat["ElectrodesInfo"]
for i = 1:length(info["ElectrodeID"])
    id = info["ElectrodeID"][i]
    ch = ns5.continuous_channels[id]
    @test ch.label == rstrip(string(info["Label"][i]...), '\0')
    @test ch.physical_connector == info["ConnectorBank"][i][1]
    @test ch.connector_pin == info["ConnectorPin"][i]
    @test ch.min_digital_value == info["MinDigiValue"][i]
    @test ch.max_digital_value == info["MaxDigiValue"][i]
    @test ch.min_analog_value == info["MinAnalogValue"][i]
    @test ch.max_analog_value == info["MaxAnalogValue"][i]
    @test ch.units == rstrip(string(info["AnalogUnits"][i]...), '\0')
    @test ch.high_pass.cutoff == info["HighFreqCorner"][i]
    @test ch.high_pass.order == info["HighFreqOrder"][i]
    @test ch.high_pass.filter_type == info["HighFilterType"][i]
    @test ch.low_pass.cutoff == info["LowFreqCorner"][i]
    @test ch.low_pass.order == info["LowFreqOrder"][i]
    @test ch.low_pass.filter_type == info["LowFilterType"][i]

    @test NeuroIO.data(ch) == vec(ns5_mat["Data"][i, :])
end


text:
"""
Struct for SubordinatedBrownianMotion

		subordinatedBrownianMotion=SubordinatedBrownianMotion(sigma::num,drift::num1,underlying::abstrUnderlying)

Where:\n
		sigma       =	Volatility.
		drift       = 	drift.
		underlying  = 	underlying.
"""
mutable struct SubordinatedBrownianMotion{num <: Number, num1 <: Number, Distr <: Distribution{Univariate, Continuous}, numtype <: Number} <: AbstractMonteCarloProcess{numtype}
    sigma::num
    drift::num1
    subordinator_::Distr
    function SubordinatedBrownianMotion(sigma::num, drift::num1, dist::Distr) where {num <: Number, num1 <: Number, Distr <: Distribution{Univariate, Continuous}}
        if sigma <= 0.0
            error("volatility must be positive")
        else
            zero_typed = zero(num) + zero(num1)
            return new{num, num1, Distr, typeof(zero_typed)}(sigma, drift, dist)
        end
    end
end

export SubordinatedBrownianMotion;

function simulate!(X, mcProcess::SubordinatedBrownianMotion, mcBaseData::MonteCarloConfiguration{type1, type2, type3, SerialMode, type5}, T::numb) where {numb <: Number, type1 <: Integer, type2 <: Integer, type3 <: StandardMC, type5 <: Random.AbstractRNG}
    Nsim = mcBaseData.Nsim
    Nstep = mcBaseData.Nstep
    drift = mcProcess.drift
    sigma = mcProcess.sigma

    @assert T > 0.0

    type_sub = typeof(rand(mcBaseData.rng, mcProcess.subordinator_))
    isDualZero = drift * zero(type_sub) * 0.0
    @views X[:, 1] .= isDualZero

    for i = 1:Nstep
        # SUBORDINATED BROWNIAN MOTION (dt_s=time change)
        for j = 1:Nsim
            dt_s = rand(mcBaseData.rng, mcProcess.subordinator_)
            @views X[j, i+1] = X[j, i] + drift * dt_s + sigma * sqrt(dt_s) * randn(mcBaseData.rng)
        end
    end

    #return X;
end

function simulate!(X, mcProcess::SubordinatedBrownianMotion, mcBaseData::MonteCarloConfiguration{type1, type2, type3, SerialMode, type5}, T::numb) where {numb <: Number, type1 <: Integer, type2 <: Integer, type3 <: AntitheticMC, type5 <: Random.AbstractRNG}
    Nsim = mcBaseData.Nsim
    Nstep = mcBaseData.Nstep
    drift = mcProcess.drift
    sigma = mcProcess.sigma

    @assert T > 0.0

    type_sub = typeof(quantile(mcProcess.subordinator_, 0.5))
    isDualZero = drift * sigma * zero(type_sub) * 0.0
    @views X[:, 1] .= isDualZero
    Nsim_2 = div(Nsim, 2)

    for j = 1:Nstep
        # SUBORDINATED BROWNIAN MOTION (dt_s=time change)
        for i = 1:Nsim_2
            dt_s = rand(mcBaseData.rng, mcProcess.subordinator_)
            Z = randn(mcBaseData.rng)
            DW = sigma * sqrt(dt_s) * Z
            muDt = drift * dt_s
            @views X[2*i-1, j+1] = X[2*i-1, j] + muDt + DW
            @views X[2*i, j+1] = X[2*i, j] + muDt - DW
        end
    end

    nothing
end


text:
#TODO: Adds continuous constraints
# adds regions and literals to prec and eff
function build_funnel!(hybrid_action, params)
    name = hybrid_action.name
    if name == :pick

    elseif name == :place

    elseif name == :move

    end
end

function build_action(act, domain, vs)
    if isdiscrete(act)
        a = Discrete_Action(act.name)
        a.pos_prec, a.neg_prec = get_preconditions(domain, act, vs)
        a.pos_eff, a.neg_eff = get_effects(domain, act, vs)
        a.params = vs
    else
        a = Funnel(act.name)
        build_funnel!(a, vs)
    end
    return a

end

function get_all_actions(domain, problem)
    actions = []
    obs = collect(PDDL.get_objects(problem))
    for act in values(PDDL.get_actions(domain))
        vars = act.args
        for vs in collect(permutations(obs, length(vars)))
            a = build_action(act, domain, vs)
            push!(actions, a)
        end
    end
    return actions
end

function get_preconditions(domain, act, vars)
    action = PDDL.get_action(domain, act.name)
    vardict = Dict()
    [vardict[action.args[i]] = vars[i] for i=1:length(action.args)]
    pos, neg = [], []
    for pre in action.precond.args
        if pre.name == :not
            vs = []
            [push!(vs, vardict[arg]) for arg in pre.args[1].args]
            prop = fill_proposition(pre.args[1], vs)
            push!(neg, prop)
        else
            vs = []
            [push!(vs, vardict[arg]) for arg in pre.args]
            prop = fill_proposition(pre, vs)
            push!(pos, prop)
        end
    end
    return pos, neg
end

function get_effects(domain, act, vars)
    action = PDDL.get_action(domain, act.name)
    vardict = Dict()
    [vardict[action.args[i]] = vars[i] for i=1:length(action.args)]
    pos, neg = [], []
    for eff in action.effect.args
        if eff.name == :not
            vs = []
            [push!(vs, vardict[arg]) for arg in eff.args[1].args]
            prop = fill_proposition(eff.args[1], vs)
            push!(neg, prop)
        else
            vs = []
            [push!(vs, vardict[arg]) for arg in eff.args]
            prop = fill_proposition(eff, vs)
            push!(pos, prop)
        end
    end
    return pos, neg
end

function fill_proposition(proposition, objs)
    if isempty(objs) objs=Term[] end
    prop = Compound(Symbol(proposition.name), objs)
    return prop
end

#=

#TODO: Augment with continuous constraint checking
function goal_reached!(domain, problem, graph)
    goal_set = goalstate(domain, problem).facts
    goal_set = [goal_set...]
    index = graph.num_levels - 1
    props = graph.props[index]
    μprops = graph.μprops[index]
    goal_found = false
    if issubset(goal_set, props)
        goal_found = true
        for goal_pair in collect(permutations(goal_set, 2))
            if goal_pair in μprops
                goal_found = false
                break
            end
        end
    elseif index > 0 && graph.props[index-1] == props
        graph.leveled = true
    end
    return goal_found
end

#TODO: Add continuous mutex checking
function is_mutex_acts(act_pair, μprops)
    a = act_pair[1]
    b = act_pair[2]

    if !isempty(intersect(a.neg_eff, union(b.pos_prec, b.pos_eff))) ||
       !isempty(intersect(b.neg_eff, union(a.pos_prec, a.pos_eff)))
        return true
    end

    if !isempty(μprops)
        for μ in μprops
            p = μ[1]
            q = μ[2]
            if (p in a.pos_prec && q in b.pos_prec)
                return true
            end
        end
    end
    return false
end

#TODO: Add continuous mutex checking
function is_mutex_props(prop_pair, action_list, μacts)
    p = prop_pair[1]
    q = prop_pair[2]

    for a in action_list
        if p in a.pos_eff && q in a.pos_eff
            return false
        end
    end

    actions_with_p = Set()
    for a in action_list
        if p in a.pos_eff
            push!(actions_with_p, a)
        end
    end

    actions_with_q = Set()
    for a in action_list
        if q in a.pos_eff
            push!(actions_with_q, a)
        end
    end

    μall = true
    for p_action in actions_with_p
        for q_action in actions_with_q
            if p_action == q_action return false end
            if !([p_action, q_action] in μacts)
                μall = false
                break
            end
        end
        if !μall break end
    end
    return μall
end

#TODO: Add continuous prop expansion
function expand!(domain, problem, graph)
    level = graph.num_levels
    #As
    action_list = []
    for action in get_all_actions(domain, problem)
        if action_is_applicable(action, graph.props[level-1], graph.μprops[level-1])
            push!(action_list, action)
        end
    end
    for prop in graph.props[level-1]
        push!(action_list, HybridNoOp(prop))
    end
    graph.acts[level] = action_list

    #Ps
    proposition_list = Set()
    for action in action_list
        for eff in action.pos_eff
            push!(proposition_list, eff)
        end
    end
    graph.props[level] = collect(proposition_list)
    proposition_list = collect(proposition_list)

    #μA
    action_μ_list = []
    for act_pair in collect(permutations(action_list, 2))
        if is_mutex_acts(act_pair, graph.μprops[level-1])
            push!(action_μ_list, act_pair)
        end
    end
    graph.μacts[level] = action_μ_list

    #μP
    proposition_μ_list = []
    for prop_pair in collect(permutations(proposition_list, 2))
        if is_mutex_props(prop_pair, action_list,action_μ_list)
            if !(prop_pair in proposition_μ_list)
                swapped = [prop_pair[2], prop_pair[1]]
                if !(swapped in proposition_μ_list)
                    push!(proposition_μ_list, prop_pair)
                end
            end
        end
    end
    graph.μprops[level] = proposition_μ_list

    graph.num_levels = level + 1
    if graph.props[level-1] == proposition_list
        graph.leveled = true
    end

    return graph
end

#TODO: add continuous init state
function get_init_propositions(domain, problem)
    initprops=[]
    inits = collect(initstate(domain, problem).facts)
    for init in inits
        objs = init.args
        push!(initprops, fill_proposition(init, objs))
    end
    return initprops

end

#TODO: add continuous goal state
function get_goal_propositions(domain, problem)
    goalprops=[]
    goals = collect(goalstate(domain, problem).facts)
    for goal in goals
        objs = goal.args
        push!(goalprops, fill_proposition(goal, objs))
    end
    return goalprops

end

function create_graph(domain, problem; max_levels=10)
    graph = Graph()
    graph.num_levels = 1
    graph.props[0] =  get_init_propositions(domain, problem)

    for _ in 1:max_levels
        expand!(domain, problem, graph)
        if goal_reached!(domain, problem, graph) break end
        if graph.leveled break end
    end
    graph.initprops = get_init_propositions(domain, problem)
    graph.goalprops = get_goal_propositions(domain, problem)
    return graph
end
=#


text:
cd(@__DIR__)
using Revise, DifferentialEvolutionMCMC, Random, Distributions
using LinearAlgebra
Random.seed!(50514)

# number of variables
n_μ = 30
# number of observations per variable
n_d = 100
# random μ parameters
μs = rand(Normal(0.0, 1.0), n_μ)
# data
data = rand(MvNormal(μs, 1.0 * I), n_d)

# function for initial values
function sample_prior()
    μ = rand(Normal(0, 1), n_μ)
    σ = rand(truncated(Cauchy(0, 1), 0, Inf))
    return as_union([μ,σ])
end

# returns prior log likelihood
function prior_loglike(μ, σ)
    LL = 0.0
    LL += sum(logpdf.(Normal(0, 1), μ))
    LL += logpdf(truncated(Cauchy(0, 1), 0, Inf), σ)
    return LL
end

# likelihood function
function loglike(data, μs, σ)
    return sum(logpdf(MvNormal(μs, σ^2 * I), data))
end

# upper and lower bounds of parameters
bounds = ((-Inf,Inf),(0.0,Inf))
# parameter names
names = (:μ,:σ)

# model object
model = DEModel(;
    sample_prior,
    prior_loglike,
    loglike,
    data,
    names
)

# DEMCMC sampler
de = DE(;
    sample_prior,
    bounds,
    sample = resample,
    burnin = 5000,
    n_initial = (n_μ + 1) * 4,
    Np = 3,
    n_groups = 1,
    θsnooker = 0.1,
)
# sample from the posterior distribution
n_iter = 50_000
chains = sample(model, de, MCMCThreads(), n_iter, progress=true)


text:
module REPLTests
    include("repl.jl")
end
module REPLCompletionsTest
    include("replcompletions.jl")
end
module LineEditTest
    include("lineedit.jl")
end

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
"""
```julia
    coldstartdelta( G::AbstractMatrix{Float64}, g::AbstractVector{Float64}; H::AbstractMatrix{Float64} = zeros( Float64, (1,1) ), L::Float64 = eigmax(H) )
```

Compute the Δ parameter needed when computing the upper iteration bound for a cold-started
fast gradient method.

Either the Hessian `H` or its maximum eigenvalue `L` must be passed as a parameter to the function.
"""
function coldstartdelta( G::AbstractMatrix{Float64}, g::AbstractVector{Float64}; H::AbstractMatrix{Float64} = zeros( Float64, (1,1) ), L::Float64 = eigmax(H) )
    (n, m) = size( G )

    model = OSQP.Model()

    options = Dict( :verbose => false )

    OSQP.setup!( model; P = sparse( triu( 1.0*I(m) ) ),
                        q = zeros( Float64, (m) ),
                        A = sparse( G ),
                        u = Vector( g ),
                        options... )

    results = OSQP.solve!( model )

    r = norm( results.x )
    Δ  = L/2 * r^2

    return Δ
end

"""
```julia
    upperiterationbound( ϵ::Float64, G::AbstractMatrix{Float64}, g::AbstractVector{Float64}; H::AbstractMatrix{Float64} = zeros( Float64, (1,1) ), L::Float64 = eigmax(H), κ::Float64 = cond(H) )
```

Compute the maximum number of iterations needed for the fast gradient method to reach `ϵ` sub-optimality level when projecting
into the constraint set ``G x ≦ g``.

To compute the iteration bound, either the Hessian `H` must be specified, or its largest eigenvalue `L`
and condition number `κ` must be supplied instead.
"""
function upperiterationbound( ϵ::Float64, G::AbstractMatrix{Float64}, g::AbstractVector{Float64}; H::AbstractMatrix{Float64} = zeros( Float64, (1,1) ), L::Float64 = eigmax(H), κ::Float64 = cond(H) )
    Δ = coldstartdelta( G, g, L = L )

    a₁ = ceil( ( log(ϵ) - log(Δ) ) / ( log( 1- sqrt(1/κ) ) ) )
    a₂ = ceil( 2 * sqrt(Δ/ϵ) - 2 )

    return max( 0, min( a₁, a₂ ) )
end


text:
using Setfield
using Setfield: IndexLens, PropertyLens, ComposedLens

struct Lens!{L <:Lens} <: Lens
    pure::L
end

Setfield.get(o, l::Lens!) = Setfield.get(o, l.pure)
function Setfield.set(o, l::Lens!{<: ComposedLens}, val)
    o_inner = get(o, l.pure.outer)
    set(o_inner, Lens!(l.pure.inner), val)
end
function Setfield.set(o, l::Lens!{PropertyLens{prop}}, val) where {prop}
    setproperty!(o, prop, val)
    o
end
function Setfield.set(o, l::Lens!{<:IndexLens}, val) where {prop}
    o[l.pure.indices...] = val
    o
end

using Setfield: setmacro, lensmacro

macro myset!(ex)
    setmacro(Lens!, ex)
end

macro mylens!(ex)
    lensmacro(Lens!, ex)
end

mutable struct M
    a
    b
end

o = M(1,2)
l = Lens!(@lens _.b)
set(o, l, 20)
@test o.b == 20

l = Lens!(@lens _.foo[1])
o = (foo=[1,2,3], bar=:bar)
set(o, l, 100)
@test o == (foo=[100,2,3], bar=:bar)

using Setfield: setmacro, lensmacro

macro myset!(ex)
    setmacro(Lens!, ex)
end

macro mylens!(ex)
    lensmacro(Lens!, ex)
end

o = M(1,2)
@myset! o.a = 2
@myset! o.b += 98
@test o.a == :hi
@test o.b == 100

deep = [[[[1]]]]
@myset! deep[1][1][1][1] = 2
@test deep[1][1][1][1] === 2

l = @mylens! _.foo[1]
o = (foo=[1,2,3], bar=:bar)
set(o, l, 100)
@test o == (foo=[100,2,3], bar=:bar)

max_stars_repo_name: Cortalak/Pissaro

max_issues_repo_name: Cortalak/Pissaro

max_forks_repo_name: Cortalak/Pissaro


text:
__precompile__(true)

module Stochatto
    include("model/constant.jl")
    include("model/scale.jl")
    include("model/interval.jl")
    include("model/chord.jl")
    include("model/key.jl")
    include("model/beat.jl")
    include("compute/utility.jl")
    include("compute/operation.jl")
    include("compute/engine.jl")
    export
        # constants
        N_OCTAVES,
        NOTE_RANGE,
        RANGE_MIN,
        RANGE_MAX,
        MIDDLE_C,
        SHARP,
        # intervals
        Interval,
        ROOT,
        OCTAVE,
        MINOR_THIRD,
        MAJOR_THIRD,
        PERFECT_FOURTH,
        PERFECT_FIFTH,
        MINOR_SIXTH,
        MAJOR_SIXTH,
        MINOR_SEVENTH,
        MAJOR_SEVENTH,
        # chords
        Chord,
        MAJOR_CHORD,
        MINOR_CHORD,
        POWER_CHORD,
        # scales
        CHROMATIC,
        MAJOR,
        MINOR,
        NATURAL_MINOR,
        HARMONIC_MINOR,
        IONIAN,
        DORIAN,
        PHRYGIAN,
        LYDIAN,
        MIXOLYDIAN,
        AEOLIAN,
        LOCRIAN,
        SCALES,
        MODES,
        # key signatures
        Key,
        # rhythm
        TPQ,
        THIRTY_SECOND,
        SIXTEENTH,
        TRIPLET,
        EIGHTH,
        QUARTER,
        HALF,
        WHOLE,
        # utility
        acf_series,
        find_note,
        to_sharp_convention,
        notefreq,
        # generation
        Engine,
        generate_note,
        generate_notes,
        generate_rhythm,
        generate
end


text:
"""
Export each given currency symbol into the current namespace. The individual
unit exported will be a full unit of the currency specified, not the smallest
possible unit. For instance, `@usingcurrencies EUR` will export `EUR`, a
currency unit worth 1€, not a currency unit worth 0.01€.

    @usingcurrencies EUR, GBP, AUD
    7AUD  # 7.00 AUD

There is no sane unit for certain currencies like `XAU` or `XAG`, so this macro
does not work for those. Instead, define them manually:

    const XAU = Monetary(:XAU; precision=4)
"""
macro usingcurrencies(curs)
    if isa(curs, Symbol)
        curs = Expr(:tuple, curs)
    end
    @assert Meta.isexpr(curs, :tuple)

    quote
        $([:(const $cur = Monetary($(Expr(:quote, cur))))
            for cur in curs.args]...)
        nothing
    end |> esc
end


text:
"""
    UniGrid(m::Int)

where `m` is the number of points on each dimension.
Adaptive uniform grid based on [1].
The resulting inducing points are stored in the memory-efficient custom type
`UniformGrid`.

[1] Moreno-Muñoz, P., Artés-Rodríguez, A. & Álvarez, M. A. Continual Multi-task Gaussian Processes. (2019).
"""
struct UniGrid <: OnIPSA
    m::Int # Number of points per dimension
end

Base.show(io::IO, Z::UniGrid) = print(io, "Uniform grid with side length $(Z.m).")

function inducingpoints(
    ::AbstractRNG,
    alg::UniGrid,
    X::Union{AbstractVector{<:Real},AbstractVector{<:AbstractVector{<:Real}}};
    kwargs...,
)
    ndim = length(first(X)) # Take the dimensionality
    bounds = [extrema(x -> getindex(x, i), X) for i in 1:ndim]
    Z = map(bounds) do lims
        LinRange(lims..., alg.m)
    end
    proditer = Iterators.product(Z...)
    return UniformGrid(proditer)
end

function updateZ!(
    ::AbstractRNG, Z::AbstractVector, alg::UniGrid, X::AbstractVector; kwargs...
)
    ndim = length(X[1])
    new_bounds = [extrema(x -> getindex(x, i), X) for i in 1:ndim]
    Zn = map(Z.proditer.iterators, new_bounds) do old_Z, new_b
        x_start = min(old_Z.start, new_b[1])
        x_stop = max(old_Z.stop, new_b[2])
        return LinRange(x_start, x_stop, alg.m)
    end
    Z.proditer = Iterators.product(Zn...)
    return Z
end

function updateZ(
    rng::AbstractRNG, Z::AbstractVector, alg::UniGrid, X::AbstractVector; kwargs...
)
    Zn = deepcopy(Z)
    return updateZ!(rng, Zn, alg, X; kwargs...)
end

"""
    UniformGrid{T,Titer} <: AbstractVector{T}

A memory-efficient custom object representing a wrapper around a `Iterators.ProductIterator`.
Supports broadcasting and other relevant array methods, and avoids explicitly computing all points on the grid.
"""
mutable struct UniformGrid{T,Titer} <: AbstractVector{T}
    proditer::Titer
end

function UniformGrid(proditer::Iterators.ProductIterator{NTuple{N,S}}) where {N,S}
    T = Vector{eltype(first(proditer.iterators))}
    return UniformGrid{T,typeof(proditer)}(proditer)
end

import Base: getindex, broadcastable, eachindex, length, size, iterate, eltype, IndexStyle

_getelement(x::NTuple{1,<:Real}) = only(x)
_getelement(x::NTuple) = collect(x)

function Base.getindex(ug::UniformGrid, i::Integer)
    return _getelement(first(Iterators.drop(ug.proditer, i - 1)))
end

function Base.getindex(ug::UniformGrid, r::UnitRange)
    return [
        _getelement(x) for
        x in Iterators.take(Iterators.drop(ug.proditer, first(r) - 1), length(r))
    ]
end

Base.getindex(ug::UniformGrid, ::Colon) = Base.getindex(ug, 1:length(ug))

function Base.broadcastable(ug::UniformGrid)
    r = similar(1:1, eltype(ug.proditer), length(ug))
    return copyto!(r, ug.proditer)
end

Base.length(ug::UniformGrid) = prod(length, ug.proditer.iterators)
Base.size(ug::UniformGrid) = (prod(length, ug.proditer.iterators),)

# slows iteration, but makes it consistent with iterating over an AbstractVector
function Base.iterate(ug::UniformGrid, state...)
    r = Base.iterate(ug.proditer, state...)
    r === nothing && return nothing
    return (_getelement(r[1]), r[2])
end

import KernelFunctions: pairwise, pairwise!
function pairwise(d::PreMetric, x::UniformGrid)
    return KernelFunctions.Distances_pairwise(d, x.proditer)
end

function pairwise!(out::AbstractMatrix, d::PreMetric, x::UniformGrid)
    return KernelFunctions.Distances.pairwise!(out, d, x.proditer)
end

Base.eltype(ug::UniformGrid) = typeof(ug[1])

Base.IndexStyle(::UniformGrid) = IndexLinear()

## show needs more improvement, maybe
function Base.show(io::IO, ug::UniformGrid)
    return print(io, "$(length.(ug.proditer.iterators)) uniform grid")
end

function Base.show(io::IO, ::MIME"text/plain", ug::UniformGrid)
    println(io, "$(length.(ug.proditer.iterators)) uniform grid with edges")
    pr(iter) = println(io, [iter[1], iter[end]])
    return pr.(ug.proditer.iterators)
end


text:
function hdr_hash(S::SeisData, i::Int64)
  h = hash(zero(UInt64))
  for f in (:id, :name, :loc, :fs, :gain, :resp, :units)
    h = hash(getfield(S, f)[i], h)
  end
  return h
end

function track_hdr!(S::GphysData, hashes::Array{UInt64, 1}, fmt::String, fname::String, opts::String)
  to_track = Array{Int64, 1}(undef, 0)

  # Add new channels
  L = length(hashes)
  if S.n > L
    append!(hashes, zeros(UInt64, S.n-L))
  end

  # Check existing channels
  for i in 1:S.n
    if hdr_hash(S, i) != hashes[i]
      push!(to_track, i)
    end
  end

  # note new meta source in to_track
  note!(S, to_track, string( "+meta ¦ read_meta!(S, \"",
                            fmt, "\", \"",
                            fname, "\", ",
                            opts) )
  return nothing
end


text:
function test_fft_cosmx(g::OneDGrid)
    f1, f1h, f1hr, f1hr_mul, f1h_th, f1hr_th = create_testfuncs(g)
    isapprox(f1h, f1h_th, rtol=rtol_fft)
end

function test_rfft_cosmx(g::OneDGrid)
    f1, f1h, f1hr, f1hr_mul, f1h_th, f1hr_th = create_testfuncs(g)
    isapprox(f1hr, f1hr_th, rtol=rtol_fft)
end

function test_rfft_mul_cosmx(g::OneDGrid)
    f1, f1h, f1hr, f1hr_mul, f1h_th, f1hr_th = create_testfuncs(g)
    isapprox(f1hr_mul, f1hr_th, rtol=rtol_fft)
end

function test_fft_cosmxcosny(g::TwoDGrid)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul, f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    isapprox(f1h, f1h_th, rtol=rtol_fft)
end

function test_rfft_cosmxcosny(g::TwoDGrid)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul, f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    isapprox(f1hr, f1hr_th, rtol=rtol_fft)
end

function test_fft_sinmxny(g::TwoDGrid)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mulPk, f2hr_mul, f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    isapprox(f2h, f2h_th, rtol=rtol_fft)
end

function test_rfft_sinmxny(g::TwoDGrid)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul, f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    isapprox(f2hr, f2hr_th, rtol=rtol_fft)
end

function test_rfft_mul_cosmxcosny(g::TwoDGrid)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul, f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    isapprox(f1hr_mul, f1hr_th, rtol=rtol_fft)
end

function test_rfft_mul_sinmxny(g::TwoDGrid)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul, f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    isapprox(f2hr_mul, f2hr_th, rtol=rtol_fft)
end


text:
import GeometryTypes: type_immutable, with_immutable, make_immutable
import GeometryTypes: gjk0, ⊖, support_vector_max

@testset "gjk" begin

    @testset "gjk examples" begin
        @testset "two simplices" begin
            c1 = Simplex(Vec(-1.))
            c2 = Simplex(Vec(4.))
            @test gjk(c1,c2) ≈5
            @test min_euclidean(c1,c2) ≈5

            c1 = Simplex(Vec(-1.,0.,0.))
            c2 = Simplex(Vec(4.,0.,0.))
            @test gjk(c1,c2) ≈5
            @test min_euclidean(c1,c2) ≈5

            c1 = FlexibleConvexHull([Vec(0.,0.), Vec(0.,1.), Vec(1.,0.),Vec(1.,1.)])
            c2 = Simplex(Vec(4.,0.5))
            @test gjk(c1, c2) ≈ 3
            @test min_euclidean(c1,c2) ≈ 3

            pt1 = Vec(1.,2.,3.)
            pt2 = Vec(3.,4.,5.)
            @test gjk(pt1, pt2) ≈norm(pt1-pt2)
            @test min_euclidean(pt1, pt2) ≈norm(pt1-pt2)
        end

        @testset "gjk intersecting lines" begin
            c1 = Simplex(Vec(1.,1.), Vec(1., 2.))
            @test gjk(c1, c1) == 0.
            @test min_euclidean(c1,c1) == 0.

            c2 = Simplex(Vec(1.,1.), Vec(10., 2.))
            @test gjk(c1, c2) == 0.
            @test min_euclidean(c1,c2) == 0.

            c3 = Simplex(Vec(0., 1.), Vec(2.,2.))
            md = vertices(c1 ⊖ c3)
            @test md == [Vec(1.0,0.0),Vec(-1.0,-1.0),Vec(1.0,1.0),Vec(-1.0,0.0)]
            @test gjk0(FlexibleConvexHull(md)) == (Vec(0., 0.), 0.)
            @test gjk(c1, c3) == 0.
            @test min_euclidean(c1,c3) == 0.

        end

        @testset "Cube" begin
            c = HyperCube(Vec(0.5,0.5,0.5), 1.)
            @test min_euclidean(Vec(2.,2.,2.), c) ≈ gjk(Vec(2.,2.,2.), c) ≈ √(3/4)

            s = Simplex(Vec(1., 0.5, 0.5), Vec(1.,2.,3.))
            @test 0 <= min_euclidean(s, c) <= 1e-14
            @test 0 <= gjk(s, c) <= 1e-14

            s = Simplex(Vec(2.,2.,2.), Vec(2.,3.,2.), Vec(2.,2.,7.), Vec(3.,4.,5.))
            @test min_euclidean(s,c) ≈ gjk(s, c) ≈ √(3/4)
        end

    end

    @testset "support_vector_max" begin
        r = HyperRectangle(Vec(-0.5, -1.), Vec(1., 2.))
        @test support_vector_max(r, Vec(1.,0.)) == (Vec(0.5,-1.), 0.5)
        @test support_vector_max(r, Vec(2.,0.)) == (Vec(0.5,-1.), 1.)
        @test support_vector_max(r, Vec(-1.,0.)) == (Vec(-0.5,-1.), 0.5)
        @test support_vector_max(r, Vec(0., 1.)) == (Vec(-0.5,1.), 1.)
        @test support_vector_max(r, Vec(1., 1.)) == (Vec(0.5,1.), 1.5)
        @test support_vector_max(FlexibleConvexHull(r), Vec(1., 1.)) == (Vec(0.5,1.), 1.5)

        c1 = Simplex(Vec(1.,1.), Vec(1., 2.))
        c3 = Simplex(Vec(0., 1.), Vec(2.,2.))
        md = c1 ⊖ c3
        fh = FlexibleConvexHull(md)
        for v in [Vec(1.,0.), Vec(12.,-10.), Vec(0.,-1.), Vec(1.,1.)]
            v_md, s_md = support_vector_max(md, v)
            v_fh, s_fh = support_vector_max(fh, v)
            @test s_md ≈ s_fh
            @test v_md ≈ v_fh
        end
    end

    @testset "make immutable" begin
        T = Vec{2, Float64}
        n = 3
        S = Simplex{n, T}
        FS = FlexibleSimplex{T}
        fs = FS([Vec(0.0,0.0), Vec(0.0,1.0), Vec(1.0,0.0)])
        s = S((Vec(0.0,0.0), Vec(0.0,1.0), Vec(1.0,0.0)))
        @test (@inferred type_immutable(FS, Val{3})) == S
        @test type_immutable(fs) == S
        @test (@inferred type_immutable(fs, Val{3})) == S
        @test make_immutable(fs) == s
        @test (@inferred make_immutable(fs, Val{3})) == s
        @test with_immutable(identity, fs) == s
    end

end


text:
type RacosClassification
    solution_space
    sample_region
    label
    positive_solution
    negative_solution
    x_positive
    uncertain_bit

    function RacosClassification(dim, positive, negative; ub=1)
        solution_space = dim
        sample_region = []
        label = []
        # solution
        positive_solution = positive
        negative_solution = negative
        x_positive = Nullable()
        uncertain_bit = ub

        regions = dim.regions
        for i in 1:dim.size
            if dim.types[i] == true
                temp = [convert(Float64, regions[i][1]), convert(Float64, regions[i][2])]
            else
                temp = [regions[i][1], regions[i][2]]
            end
            push!(sample_region, temp)
            push!(label, false)
        end
        return new(solution_space, sample_region, label, positive_solution,
          negative_solution, x_positive, uncertain_bit)
    end
end

function reset_classifier(classifier)
    regions = classfier.solution_space.regions
    for i in 1:classifier.solution_space.size
        classifier.sample_region[i][0] = regions[i][0]
        classifier.sample_region[i][1] = regions[i][1]
        classifier.label[i] = false
    end
    classifier.x_positive = Nullable()
end

# This algos always works, whether discrete or continuous, we always use this function.
function mixed_classification(classifier)
    classifier.x_positive = classifier.positive_solution[rand(rng,
      1:length(classifier.positive_solution))]
    len_negative = length(classifier.negative_solution)
    index_set = Array(1:classifier.solution_space.size)
    types = classifier.solution_space.types
    while len_negative > 0
        pos = rand(rng, 1:length(index_set))
        k = index_set[pos]
        x_pos_k = classifier.x_positive.x[k]
        # continuous
        if types[k] == true
            x_negative = classifier.negative_solution[rand(rng, 1:len_negative)]
            x_neg_k = x_negative.x[k]
            if x_pos_k < x_neg_k
                r = rand_uniform(rng, x_pos_k, x_neg_k)
                if r < classifier.sample_region[k][2]
                    classifier.sample_region[k][2] = r
                    i = 1
                    while i <= len_negative
                        if classifier.negative_solution[i].x[k] >= r
                            itemp = classifier.negative_solution[i]
                            classifier.negative_solution[i] = classifier.negative_solution[len_negative]
                            classifier.negative_solution[len_negative] = itemp
                            len_negative -= 1
                        else
                            i += 1
                        end
                    end
                end
            else
                r = rand_uniform(rng, x_neg_k, x_pos_k)
                if r > classifier.sample_region[k][1]
                    classifier.sample_region[k][1] = r
                    i = 1
                    while i <= len_negative
                        if classifier.negative_solution[i].x[k] <= r
                            itemp = classifier.negative_solution[i]
                            classifier.negative_solution[i] = classifier.negative_solution[len_negative]
                            classifier.negative_solution[len_negative] = itemp
                            len_negative -= 1
                        else
                            i += 1
                        end
                    end
                end
            end
            # discrete
        else
            delete = 0
            i = 1
            while i <= len_negative
                if classifier.negative_solution[i].x[k] != x_pos_k
                    delete += 1
                    itemp = classifier.negative_solution[i]
                    classifier.negative_solution[i] = classifier.negative_solution[len_negative]
                    classifier.negative_solution[len_negative] = itemp
                    len_negative -= 1
                else
                    i += 1
                end
            end
            if delete != 0
                splice!(index_set, pos)
            end
            if length(index_set) == 0
                push!(index_set, k)
            end
        end
    end
    set_uncertain_bit!(classifier, index_set)
end

function set_uncertain_bit!(classifier, iset)
    index_set = iset
    for i in 1:classifier.uncertain_bit
        pos = rand(rng, 1:length(index_set))
        index = index_set[pos]
        classifier.label[index] = true
        splice!(index_set, pos)
    end
end

function rand_sample(classifier)
    x = []
    for i in 1:classifier.solution_space.size
        if classifier.label[i] == true
            if classifier.solution_space.types[i] == true
                push!(x, rand_uniform(rng, classifier.sample_region[i][1], classifier.sample_region[i][2]))
            else
                push!(x, rand(rng, classifier.sample_region[i][1]:classifier.sample_region[i][2]))
            end
        else
            push!(x, classifier.x_positive.x[i])
        end
    end
    return x
end

# for dubugging
function print_neg(classifier)
    zoolog("------print neg------")
    for x in classifier.negative_solution
        sol_print(x)
    end
end

function print_pos(classifier)
    zoolog("------print pos------")
    for x in classifier.positive_solution
        sol_print(x)
    end
end

function print_sample_region(classifier)
    zoolog("------print sample region------")
    zoolog(classifier.sample_region)
end


text:
##############################################################################
# From client
const FailureHandlingKind = String
const FailureHandlingKinds = (Abort = "abort",
                              Transactional = "transactional",
                              TextOnlyTransactional = "textOnlyTransactional",
                              Undo = "undo")

const ResourceOperationKind = String
const ResourceOperationKinds = (Create = "create",
                                Rename = "rename",
                                Delete = "delete")

@dict_readable struct WorkspaceEditClientCapabilities <: Outbound
    documentChanges::Union{Bool,Missing}
    resourceOperations::Union{Vector{ResourceOperationKind},Missing}
    failureHandling::Union{FailureHandlingKind,Missing}
end

@dict_readable struct DidChangeConfigurationClientCapabilities <: Outbound
    dynamicRegistration::Union{Bool,Missing}
end

@dict_readable struct DidChangeWatchedFilesClientCapabilities <: Outbound
    dynamicRegistration::Union{Bool,Missing}
end

@dict_readable struct WorkspaceClientCapabilities <: Outbound
    applyEdit::Union{Bool,Missing}
    workspaceEdit::Union{WorkspaceEditClientCapabilities,Missing}
    didChangeConfiguration::Union{DidChangeConfigurationClientCapabilities,Missing}
    didChangeWatchedFiles::Union{DidChangeWatchedFilesClientCapabilities,Missing}
    symbol::Union{WorkspaceSymbolClientCapabilities,Missing}
    executeCommand::Union{ExecuteCommandClientCapabilities,Missing}
    workspaceFolders::Union{Bool,Missing}
    configuration::Union{Bool,Missing}
end

@dict_readable struct TextDocumentSyncClientCapabilities <: Outbound
    dynamicRegistration::Union{Bool,Missing}
    willSave::Union{Bool,Missing}
    willSaveWaitUntil::Union{Bool,Missing}
    didSave::Union{Bool,Missing}
end

@dict_readable struct TagClientCapabilities
    valueSet::Vector{DiagnosticTag}
end

@dict_readable struct PublishDiagnosticsClientCapabilities <: Outbound
    relatedInformation::Union{Bool,Missing}
    tagSupport::Union{TagClientCapabilities,Missing}
    versionSupport::Union{Bool,Missing}
end

@dict_readable struct TextDocumentClientCapabilities <: Outbound
    synchronization::Union{TextDocumentSyncClientCapabilities,Missing}
    completion::Union{CompletionClientCapabilities,Missing}
    hover::Union{HoverClientCapabilities,Missing}
    signatureHelp::Union{SignatureHelpClientCapabilities,Missing}
    declaration::Union{DeclarationClientCapabilities,Missing}
    definition::Union{DefinitionClientCapabilities,Missing}
    typeDefinition::Union{TypeDefinitionClientCapabilities,Missing}
    implementation::Union{ImplementationClientCapabilities,Missing}
    references::Union{ReferenceClientCapabilities,Missing}
    documentHighlight::Union{DocumentHighlightClientCapabilities,Missing}
    documentSymbol::Union{DocumentSymbolClientCapabilities,Missing}
    codeAction::Union{CodeActionClientCapabilities,Missing}
    codeLens::Union{CodeLensClientCapabilities,Missing}
    documentLink::Union{DocumentLinkClientCapabilities,Missing}
    colorProvider::Union{DocumentColorClientCapabilities,Missing}
    formatting::Union{DocumentFormattingClientCapabilities,Missing}
    rangeFormatting::Union{DocumentRangeFormattingClientCapabilities,Missing}
    onTypeFormatting::Union{DocumentOnTypeFormattingClientCapabilities,Missing}
    rename::Union{RenameClientCapabilities,Missing}
    publishDiagnostics::Union{PublishDiagnosticsClientCapabilities,Missing}
    foldingRange::Union{FoldingRangeClientCapabilities,Missing}
    selectionRange::Union{SelectionRangeClientCapabilities,Missing}
end

@dict_readable struct WindowClientCapabilities <: Outbound
    workDoneProgress::Union{Bool,Missing}
end

@dict_readable struct ClientCapabilities <: Outbound
    workspace::Union{WorkspaceClientCapabilities,Missing}
    textDocument::Union{TextDocumentClientCapabilities,Missing}
    window::Union{WindowClientCapabilities,Missing}
    experimental::Union{Any,Missing}
end
@dict_readable struct InfoParams <: Outbound
    name::String
    version::Union{String,Missing}
end

struct InitializeParams <: Outbound
    processId::Union{Int,Nothing}
    clientInfo::Union{InfoParams,Missing}
    rootPath::Union{DocumentUri,Nothing,Missing}
    rootUri::Union{DocumentUri,Nothing}
    initializationOptions::Union{Any,Missing}
    capabilities::ClientCapabilities
    trace::Union{String,Missing}
    workspaceFolders::Union{Vector{WorkspaceFolder},Nothing,Missing}
    workDoneToken::Union{Int,String,Missing} # ProgressToken
end

# Requires handwritten implementaiton to account for 3-part Unions
function InitializeParams(dict::Dict)
    InitializeParams(dict["processId"],
    haskey(dict, "clientInfo") ? InfoParams(dict["clientInfo"]) : missing,
    !haskey(dict, "rootPath") ? missing : dict["rootPath"] === nothing ? nothing : DocumentUri(dict["rootPath"]),
    dict["rootUri"] === nothing ? nothing : DocumentUri(dict["rootUri"]),
    get(dict, "initializationOptions", missing),
    ClientCapabilities(dict["capabilities"]),
    haskey(dict, "trace") ? String(dict["trace"]) : missing ,
    !haskey(dict, "workspaceFolders") ? missing : dict["workspaceFolders"] === nothing ? nothing : WorkspaceFolder.(dict["workspaceFolders"]),
    haskey(dict, "workDoneToken") ? dict["workDoneToken"] : missing)
end
##############################################################################

##############################################################################
# Server Response
struct SaveOptions <: Outbound
    includeText::Union{Bool,Missing}
end

struct ColorProviderOptions <: Outbound end

const TextDocumentSyncKind = Int
const TextDocumentSyncKinds = (None = 0,
                               Full = 1,
                               Incremental = 2)

struct TextDocumentSyncOptions <: Outbound
    openClose::Union{Bool,Missing}
    change::Union{TextDocumentSyncKind,Missing}
    willSave::Union{Bool,Missing}
    willSaveWaitUntil::Union{Bool,Missing}
    save::Union{SaveOptions,Missing}
end

struct WorkspaceFoldersOptions <: Outbound
    supported::Union{Bool,Missing}
    changeNotifications::Union{Bool,String,Missing}
end

struct WorkspaceOptions <: Outbound
    workspaceFolders::Union{WorkspaceFoldersOptions,Missing}
end

struct WorkspaceFoldersServerCapabilities <: Outbound
    supported::Union{Bool,Missing}
    changeNotifications::Union{String,Bool,Missing}
end

struct ServerCapabilities <: Outbound
    textDocumentSync::Union{TextDocumentSyncOptions,Int,Missing}
    completionProvider::Union{CompletionOptions,Missing}
    hoverProvider::Union{Bool,HoverOptions,Missing}
    signatureHelpProvider::Union{SignatureHelpOptions,Missing}
    declarationProvider::Union{Bool,DeclarationOptions,DeclarationRegistrationOptions,Missing}
    definitionProvider::Union{Bool,DefinitionOptions,Missing}
    typeDefinitionProvider::Union{Bool,TypeDefinitionOptions,TypeDefinitionRegistrationOptions,Missing}
    implementationProvider::Union{Bool,ImplementationOptions,ImplementationRegistrationOptions,Missing}
    referencesProvider::Union{Bool,ReferenceOptions,Missing}
    documentHighlightProvider::Union{Bool,DocumentHighlightOptions,Missing}
    documentSymbolProvider::Union{Bool,DocumentSymbolOptions,Missing}
    codeActionProvider::Union{Bool,CodeActionOptions,Missing}
    codeLensProvider::Union{CodeLensOptions,Missing}
    documentLinkProvider::Union{DocumentLinkOptions,Missing}
    colorProvider::Union{Bool,DocumentColorOptions,DocumentColorRegistrationOptions,Missing}
    documentFormattingProvider::Union{Bool,DocumentFormattingOptions,Missing}
    documentRangeFormattingProvider::Union{Bool,DocumentRangeFormattingOptions,Missing}
    documentOnTypeFormattingProvider::Union{DocumentOnTypeFormattingOptions,Missing}
    renameProvider::Union{Bool,RenameOptions,Missing}
    foldingRangeProvider::Union{Bool,FoldingRangeOptions,FoldingRangeRegistrationOptions,Missing}
    executeCommandProvider::Union{ExecuteCommandOptions,Missing}
    selectionRangeProvider::Union{Bool,SelectionRangeOptions,SelectionRangeRegistrationOptions,Missing}
    workspaceSymbolProvider::Union{Bool,Missing}
    workspace::Union{WorkspaceOptions,Missing}
    experimental::Union{Any,Missing}
end

struct InitializeResult <: Outbound
    capabilities::ServerCapabilities
    serverInfo::Union{InfoParams,Missing}
end

##############################################################################
@dict_readable struct InitializedParams
end


text:
using EarthDataLab
using Test
using Random
using Statistics
using OnlineStats
@testset "OnlineStats" begin
c=Cube()

d = getCubeData(c,variable="air_temperature_2m",longitude=(30,31),latitude=(50,51),
              time=(Date("2002-01-01"),Date("2008-12-31")))

mytable = @CubeTable tair=d axes=(lat,)

d2=readcubedata(d)

@test mean(d2.data) ≈ oo.data[1]
@test mean(d2.data,dims=(1,3))[:] ≈ oo2.data

#Test KMeans
d2 = getCubeData(c,variable=["air_temperature_2m","gross_primary_productivity"],longitude=(30,31),latitude=(50,51),
              time=(Date("2002-01-01"),Date("2008-12-31")))

dm=readcubedata(d2)
xin=permutedims(dm.data,[4,1,2,3])
xin=reshape(xin,(2,length(xin) ÷ 2))
startVal=mean(xin,dims=2)[:].+rand(Float32,2,5)

x=mapCube(KMeans,d2,5,2,MDAxis=VariableAxis)

o2=KMeans(2,5)
o2.value[:]=startVal
fit!(o2,xin')
@test all(isapprox.(OnlineStats.value(o2),x.data))

#Test covariance Matrix
covmat,means = mapCube(CovMatrix,dm,MDAxis=VariableAxis)

@test all(isapprox.(covmat.data,cov(reshape(dm.data,length(dm.data) ÷ 2,2))))

using EarthDataLab
using DataStructures
c=Cube()
d2 = getCubeData(c,variable=["air_temperature_2m","gross_primary_productivity"],longitude=(30,31),latitude=(50,51),
              time=(Date("2002-01-01"),Date("2008-12-31")))
dm=readcubedata(d2)

randmask=CubeMem(dm.axes[1:2],rand(1:2,size(dm)[1:2]),zeros(UInt8,size(dm)[1:2]),
  Dict{String,Any}("labels"=>OrderedDict(1=>"one",2=>"two"),"name"=>"RandMask"))

#Test PCA
m1=randmask.data.==1
m2=randmask.data.==2
ii2=findall(m2)
for i=1:length(ii2)
  dm.data[ii2[i],:,:]=rand(322,2)
end

p=cubePCA(dm,by=[randmask],noutdims=2)

println((explained_variance(p)).data[1, 1])

@test all([0.984 0.45; 0.09 0.45] .< explained_variance(p).data[1,1])
rotation_matrix(p)
transformPCA(p,dm).data

Random.seed!(1)
d2=readcubedata(d)
mask=CubeMem(d2.axes,rand(1:10,size(d2.data)),zeros(UInt8,size(d2.data)),Dict("labels"=>Dict(i=>string(i) for i=1:10)))
mask2=CubeMem(d2.axes[1:2],rand(1:10,4,4),zeros(UInt8,4,4),Dict("labels"=>Dict(i=>string(i) for i=1:10)))

oogrouped = mapCube(Mean,d2,by=(mask,))
@test isa(oogrouped.axes[1],CategoricalAxis{String,:Label})

oogrouped2 = mapCube(Mean,d2,by=(mask2,))

for k=1:10
  know = findfirst(j->j==string(k),oogrouped2.axes[1].values)
  @test oogrouped.data[know] ≈ mean(d2.data[mask.data.==k])

  i=findall(mask2.data.==k)

  if length(i) > 0
    dhelp=Float32[]
    for jj in i
      append!(dhelp,d2.data[jj,:])
    end
    @test mean(dhelp) ≈ oogrouped2.data[know]
  end
end

end


text:
struct Weight <: NodeModifier
    w
end

n1 = Node(:a)
n2 = ModifiedNode(Node(:b), Weight(2.0))
e = Edge(n1, n2)

@testset "Node Ids" begin
    @test StenoGraphs.id(n1) == :a
    @test StenoGraphs.id(n2) == :b
end

@testset "Edge Ids" begin
    @test StenoGraphs.src_id(e) == :a
    @test StenoGraphs.dst_id(e) == :b
    @test StenoGraphs.id(e) == (:a, :b)
end

@testset "UndirectedEdge Ids" begin
    @test StenoGraphs.id(UndirectedEdge(Node(:a), Node(:b))) ==
    StenoGraphs.id(UndirectedEdge(Node(:b), Node(:a)))
end


text:
using LightGraphs
using MonteCarloX
using Random

include("utils.jl")

function test_sweep_random_element(;verbose = false)
    function test_random_element(list_prob::Vector{Float64})
        rng = MersenneTwister(1000)
        N = 10000000
        sum = 0.0
        for i in 1:N
            id = MonteCarloX.random_element(list_prob, rng)
            sum += id / N
        end
        return sum
    end
    function test_random_element(list_prob::ProbabilityWeights)
        rng = MersenneTwister(1000)
        N = 10000000
        sum = 0.0
        for i in 1:N
            id = StatsBase.sample(rng, list_prob)
            sum += id / N
        end
        return sum
    end
    function test_random_element(list_freq::FrequencyWeights)
        rng = MersenneTwister(1000)
        N = 10000000
        sum = 0.0
        for i in 1:N
            id = StatsBase.sample(rng, list_freq)
            sum += id / N
        end
        return sum
    end
    function test_binary_search(cum_prob::Vector{Float64})
        rng = MersenneTwister(1000)
        N = 10000000
        sum = 0.0
        for i in 1:N
            id = MonteCarloX.binary_search(cum_prob, rand(rng))
            sum += id / N
        end
        return sum
    end

    pass = true

    # list_prob = [0.1,0.1,0.4,0.2,0.2]
    list_prob = [0.4,0.2,0.2,0.1,0.1]
    list_freq = [4,2,2,1,1]
    cum_prob    = cumsum(list_prob)
    #############################################################################
    expected_id = 0
    for id in 1:length(list_prob)
        expected_id += id * list_prob[id]
    end
    if verbose
        println("expectation value id = $(expected_id)")
    end
    #############################################################################
    result_re = @time test_random_element(list_prob)
    result_pb = @time test_random_element(ProbabilityWeights(list_prob))
    result_fr = @time test_random_element(FrequencyWeights(list_freq))
    result_bs = @time test_binary_search(cum_prob)
    diff1 = abs(result_re - result_bs)
    diff2 = abs(result_pb - result_bs)
    diff3 = abs(result_fr - result_bs)
    if verbose
        println("result of loop over random element (deprecated) = $(result_re) vs loop over binary_search = $(result_bs) [abs diff = $(diff)]")
        println("result of loop over random element (prob) = $(result_pb) vs loop over binary_search = $(result_bs) [abs diff = $(diff)]")
        println("result of loop over random element (freq) = $(result_fr) vs loop over binary_search = $(result_bs) [abs diff = $(diff)]")
    end
    pass &= diff1 < 0.01
    pass &= diff2 < 0.01
    pass &= diff3 < 0.01
    #############################################################################
    return pass
end

# TODO: use this also for PT and adaptive overlap stories etc

""" Testing metropolis on unimodaL probability distribution"""
function test_unimodal_metropolis(;verbose = false)
    log_weight_unimodal(x, x0, var) = -(x - x0)^2 / 2 / var - log(sqrt(2 * pi * var))
    var = 1.0
    dx = 0.1
    list_x0 = [0.,1.,2.,3.]
    pass = true

    rng = MersenneTwister(1000)

    for x0 in list_x0
        log_weight_sampling(x) = log_weight_unimodal(x, x0, var)
        samples = Int(1e6)
        list_x = zeros(samples)
        x = x0
        # thermalization 100 updates
        for update in 1:samples + 100
            # Metropolis:
            x_new = step(x, dx, rng)
            if accept(rng, log_weight_sampling, x_new, x)
                x = x_new
            end
            if update > 100
                list_x[update - 100] = x
            end
        end
        # sort into binned histogram
        # bin size/ 2*sigma ~    #elements per bin (~100)/ samples*0.6/
        bin_x =    1e4 * sqrt(var) * 2 / 0.6 / samples
        P_meas = fit(Histogram, list_x, minimum(list_x):bin_x:maximum(list_x) + bin_x, closed = :left)
        P_meas = normalize(P_meas)
        kld = abs(kldivergence(P_meas, x->exp(log_weight_sampling(x))))
        if verbose
            println("result kldivergence P_meas to P_true = $(kld)")
        end
        pass &= kld < 0.1
        ###########################################################################
        est_x0 = sum(list_x) / length(list_x)
        diff_est_x0 = abs(est_x0 - x0)
        if verbose
            println("result estimate x0 = $(est_x0) vs x0 = $(x0) [diff = $(diff_est_x0)]")
        end
        pass &= diff_est_x0 < 0.1
    end
    return pass
end

""" Testing metropolis on unimodaL probability distribution"""
function test_2D_unimodal_metropolis(;verbose = false)
    log_weight_unimodal(x::Float64, y::Float64, x0::Float64, y0::Float64)::Float64 = -(x - x0)^2 / 2.0 - (y - y0)^2 / 2.0
    dx = dy = 0.1
    list_x0 = [0.0,1.0,2.0,3.0]
    y0 = 1.0
    pass = true

    rng = MersenneTwister(1000)

    for x0 in list_x0
        log_weight_sampling(x::Float64, y::Float64)::Float64 = log_weight_unimodal(x, y, x0, y0)
        samples = Int(1e6)
        list_x = zeros(samples)
        list_y = zeros(samples)
        x = x0
        y = y0
        # thermalization 100 updates
        for update in 1:samples + 100
            # Metropolis:
            x_new = x - dx + 2.0 * dx * rand(rng)
            y_new = y - dy + 2.0 * dy * rand(rng)
            if accept(rng, log_weight_sampling, (x_new, y_new), (x, y))
                x = x_new
                y = y_new
            end
            if update > 100
                list_x[update - 100] = x
                list_y[update - 100] = y
            end
        end
        # sort into binned histogram
        # bin size/ 2*sigma ~    #elements per bin (~100)/ samples*0.6/
        ###########################################################################
        est_x0 = sum(list_x) / length(list_x)
        diff_est_x0 = abs(est_x0 - x0)
        if verbose
            println("result estimate x0 = $(est_x0) vs x0 = $(x0) [diff = $(diff_est_x0)]")
        end
        pass &= diff_est_x0 < 0.1
    end
    return pass
end

# TODO: test speed difference
""" Testing sweep on unimodaL probability distribution"""
function test_unimodal_sweep(;verbose = false)
    log_weight_unimodal(x, x0, var) = -(x - x0)^2 / 2 / var - log(sqrt(2 * pi * var))
    var = 1.0

    rng = MersenneTwister(1000)

    pass = true
    list_x0 = [0.0,1.,2.,3.]
    for x0 in list_x0
        log_weight_sampling(x) = log_weight_unimodal(x, x0, var)
        samples = Int(1e5)
        list_x = zeros(samples)

        if verbose
            println("\nTesting sweep with multiple functions")
        end
        system = System(x0)
        stats1 = Stats(0, 0)
        stats2 = Stats(0, 0)
        # TODO: optimize this
        list_updates = [()->update(log_weight_sampling, system, 0.1, rng, stats1),
                        ()->update(log_weight_sampling, system, 0.2, rng, stats2),
                        ()->update(log_weight_sampling, system, 0.1, rng)]
        list_probabilities = ProbabilityWeights([0.2,0.4,0.4])

        # perfrom 100 updates as thermalization
        sweep(list_updates, list_probabilities, rng, number_updates = 100)

        stats2 = Stats(0, 0)
        for updates in 1:samples
            sweep(list_updates, list_probabilities, rng, number_updates = 10)
            list_x[updates] = system.x
        end

        ###########################################################################
        num_update1 = stats1.accept + stats1.reject
        num_update1_target = 0.2 * samples * 10 + 100
        rel_diff_num_update1 = abs(num_update1 - num_update1_target) / num_update1_target
        if verbose
            println("number calls to update1= $(num_update1) vs target = $(num_update1_target) [rel_diff = $(rel_diff_num_update1)]")
        end
        pass &= rel_diff_num_update1 < 0.01
        ###########################################################################
        num_update2 = stats2.accept + stats2.reject
        num_update2_target = 0.4 * samples * 10
        rel_diff_num_update2 = abs(num_update2 - num_update2_target) / num_update2_target
        if verbose
            println("number calls to update2= $(num_update2) vs target = $(num_update2_target) [rel_diff = $(rel_diff_num_update2)]")
        end
        pass &= rel_diff_num_update2 < 0.01
        ##########################################################################
        # sort into binned histogram
        # bin size/ 2*sigma ~    #elements per bin (~100)/ samples*0.6/
        bin_x =    1e4 * sqrt(var) * 2 / 0.6 / samples
        P_meas = fit(Histogram, list_x, minimum(list_x):bin_x:maximum(list_x) + bin_x, closed = :left)
        P_meas = normalize(P_meas)
        kld = abs(kldivergence(P_meas, x->exp(log_weight_sampling(x))))
        if verbose
            println("result kldivergence P_meas to P_true = $(kld)")
        end
        pass &= kld < 0.1
        ###########################################################################
        est_x0 = sum(list_x) / length(list_x)
        diff_est_x0 = abs(est_x0 - x0)
        if verbose
            println("result estimate x0 = $(est_x0) vs x0 = $(x0) [diff = $(diff_est_x0)]")
        end
        pass &= diff_est_x0 < 0.1

        if verbose
            println("\nTesting sweep with single function")
        end
        update1() = update(log_weight_sampling, system, 0.1, rng)
        for updates in 1:samples
            sweep(update1, rng, number_updates = 10)
            list_x[updates] = system.x
        end
        ##########################################################################
        # sort into binned histogram
        # bin size/ 2*sigma ~    #elements per bin (~100)/ samples*0.6/
        bin_x =    1e4 * sqrt(var) * 2 / 0.6 / samples
        P_meas = fit(Histogram, list_x, minimum(list_x):bin_x:maximum(list_x) + bin_x, closed = :left)
        P_meas = normalize(P_meas)
        kld = abs(kldivergence(P_meas, x->exp(log_weight_sampling(x))))
        if verbose
            println("result kldivergence P_meas to P_true = $(kld)")
        end
        pass &= kld < 0.1
        ###########################################################################
        est_x0 = sum(list_x) / length(list_x)
        diff_est_x0 = abs(est_x0 - x0)
        if verbose
            println("result estimate x0 = $(est_x0) vs x0 = $(x0) [diff = $(diff_est_x0)]")
        end
        pass &= diff_est_x0 < 0.1
    end
    return pass
end

""" Testing reweighting on unimodaL probability distribution"""
# TODO: write
function test_unimodal_reweighting()
    return pass
end

###############################################################################
###############################################################################
###############################################################################

mutable struct System
    x::Float64
end

mutable struct Stats
    accept::Float64
    reject::Float64
end

function update(log_weight_sampling, system::System, dx, rng::AbstractRNG)
    x_new = step(system.x, dx, rng)
    if accept(rng, log_weight_sampling, x_new, system.x)
        system.x = x_new
    end
end

# TODO: test second version with statistics could this be implemented in update w/o speed loss? -> could write a macro like @assert
function update(log_weight_sampling, system::System, dx, rng::AbstractRNG, stats::Stats)
    x_new = step(system.x, dx, rng)
    if accept(rng, log_weight_sampling, x_new, system.x)
        system.x = x_new
        stats.accept += 1
    else
        stats.reject += 1
    end
end

function step(x, dx, rng)
    return x - dx + 2 * dx * rand(rng)
end


text:
using Plots

N = 1000
dt = 0.001
x0 = 1.0

t = zeros(N)

x = zeros(N)
x[1] = 1.0

for i in 1:N
    t[i] = i * dt
end

for i in 1:N - 1
    x[i + 1] = (1 - dt) * x[i]
end

scatter(t, x)
savefig("Julia-file/output/result_array.png")


text:
using BlockDiagonals
using Documenter
using Test

@testset "BlockDiagonals" begin
    doctest(BlockDiagonals)
    include("blockdiagonal.jl")
    include("base_maths.jl")
    include("linalg.jl")
end  # tests


text:
# Lotka-Volterra model
using Logging

using StructuralIdentifiability

logger = Logging.SimpleLogger(stdout, Logging.Debug)
global_logger(logger)

ode = @ODEmodel(
    x1'(t) = a * x1(t) - b * x1(t) * x2(t) + u(t),
    x2'(t) = -c * x2(t) + d * x1(t) * x2(t),
    y(t) = x1(t)
)

@time println(assess_identifiability(ode))


text:
using DataDeps
addd
# WIP


text:
"""
    FeedbackNets

Implements deep networks with feedback operations from higher to lower layers.
Uses Flux as a backend.
"""
module FeedbackNets
using Reexport

include("Splitters.jl")
include("Mergers.jl")
include("FeedbackChains.jl")
include("FeedbackTrees.jl")
include("ModelFactory.jl")
@reexport using .Splitters
@reexport using .Mergers
@reexport using .FeedbackChains
@reexport using .FeedbackTrees
@reexport using .ModelFactory

end # module FeedbackNets


text:
#converts gradient field into heightmap using Fourier method
function convert_gradient(scheme::Frankot, pin::AbstractArray, qin::AbstractArray)
    p = Complex{Float64}.(pin)
    q = Complex{Float64}.(qin)
    wx, wy = setup_transform_values(last(size(p)),first(size(p)))
    fft!(p)
    fft!(q)
    Z = zeros(Complex{Float64}, size(p))
    Z = (-1im .* wx .* p .+ 1im .* wy .* q) ./ (2 .* π .* (wx.^2 .+ wy.^2 .+ eps()))
    ifft!(Z)
    Z = abs.(Z)
    return Z
end

#converts gradient field into heightmap using average of two integrals
function convert_gradient(scheme::Path, p::AbstractArray, q::AbstractArray)
    R, C = size(p)
    Z₁ = zeros(Float64, size(p))
    Z₂ = zeros(Float64, size(p))

    #path 1
    for i = 2:R
        Z₁[i,1] = Z₁[i-1,1] - q[i,1]
    end
    for i = 1:R
        for j = 2:C
            Z₁[i,j] = Z₁[i,j-1] + p[i,j]
        end
    end

    #path 2
    for i = 2:C
        Z₂[1,i] = Z₂[1,i-1] + p[1,i]
    end
    for i = 2:R
        for j = 1:C
            Z₂[i,j] = Z₂[i-1,j] - q[i,j]
        end
    end
    Z = (Z₁ .+ Z₂) / 2
    return Z
end

#converts gradient field into heightmap using integral culculated from average gradient at point
function convert_gradient(scheme::SplitPath, p::AbstractArray, q::AbstractArray)
    R, C = size(p)
    Z = zeros(Float64, size(p))
    for i = 2:R
        Z[i,1]=Z[i-1,1] + q[i,1]
        Z[1,i]=Z[1,i-1] - p[1,i]
    end
    for i = 2:R
        for j = i:R
            Z[i,j] = ((Z[i,j-1] + p[i,j]) + (Z[i-1,j] - q[i,j])) / 2
            Z[j,i] = ((Z[j,i-1] + p[j,i]) + (Z[j-1,i] - q[j,i])) / 2
        end
    end
    return Z
end

#Horn and Brooks
function convert_gradient(scheme::Horn, p::AbstractArray, q::AbstractArray, iter::Real = 10000, ϵ::Real = 1)
    Z = zeros(Float64, size(p))
    Zᵏ⁺¹ = zeros(Float64, size(p))
    h = zeros(Float64, size(p))
    v = zeros(Float64, size(p))
    hv = zeros(Float64, size(p))
    R, C = size(p)
    for i = 1:R
        Z[i,1] = p[i,1]
        Z[i,R] = p[i,R]
    end
    for i = 1:C
        Z[1,i] = q[1,i]
        Z[C,i] = q[C,i]
    end
    for i = 2:(R-1)
        for j = 2:(C-1)
            h[i,j] = (p[i,j+1] - p[i,j-1]) / 2
            v[i,j] = (q[i+1,j] - q[i-1,j]) / 2
            hv[i,j] = (ϵ / 4) * (h[i,j] - v[i,j])
        end
    end
    for k = 1:iter
        copyto!(Zᵏ⁺¹, Z)
        for i = 2:(R-1)
            for j = 2:(C-1)
                Zᵏ⁺¹[i,j] = (Z[i-1,j] + Z[i+1,j] + Z[i,j-1] + Z[i,j+1]) / 4
                Zᵏ⁺¹[i,j] = Zᵏ⁺¹[i,j] - hv[i,j]
            end
        end
        Z = Zᵏ⁺¹
    end
    return Z
end

#Durou and Courteille
function convert_gradient(scheme::Durou, p::AbstractArray, q::AbstractArray, iter::Real = 1000, ϵ::Real = 1)
    Z = zeros(Float64, size(p))
    Zᵏ⁺¹ = zeros(Float64, size(p))
    h = zeros(Float64, size(p))
    v = zeros(Float64, size(p))
    hv = zeros(Float64, size(p))
    R, C = size(p)
    for i = 1:(R-1)
        for j = 1:(C-1)
            h[i,j] = (p[i,j+1] + p[i,j]) / 2
            v[i,j] = (q[i+1,j] + q[i,j]) / 2
            hv[i,j] = (ϵ / 2) * (h[i,j] - v[i,j])
        end
    end
    for k = 1:iter
        copyto!(Zᵏ⁺¹, Z)
        for i = 1:(R-1)
            for j = 1:(C-1)
                Zᵏ⁺¹[i,j] = (Z[i+1,j] + Z[i,j+1]) / 2
                Zᵏ⁺¹[i,j] = Zᵏ⁺¹[i,j] - hv[i,j]
            end
        end
        Z = Zᵏ⁺¹
    end
    return Z
end


text:
#=
  Medical test.

  https://www.math.hmc.edu/funfacts/ffiles/30002.6.shtml
  """
  Suppose that you are worried that you might have a rare disease. You decide to get tested,
  and suppose that the testing methods for this disease are correct 99 percent of the time
  (in other words, if you have the disease, it shows that you do with 99 percent probability,
  and if you don't have the disease, it shows that you do not with 99 percent probability).
  Suppose this disease is actually quite rare, occurring randomly in the general population
  in only one of every 10,000 people.

  If your test results come back positive, what are your chances that you actually have the disease?

  Do you think it is approximately: (a) .99, (b) .90, (c) .10, or (d) .01?

  Surprisingly, the answer is (d), less than 1 percent chance that you have the disease!

  """

  Cf ~/cplint/medical_test.pl
     ~/blog/medical_test.blog
     ~/psi/medical_test.psi
     ~/webppl/medical_test.wppl

=#

# using Memoization
using Turing, StatsPlots, DataFrames
using ReverseDiff, Zygote, Tracker
using Printf
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)

include("jl_utils.jl")

@model function medical_test(test, disease,prob_of_disease=1/10000.0,reliability=0.99)
    # Suppose this disease is actually quite rare, occurring randomly in the general population
    # in only one of every 10,000 people.
    # prob_of_disease = 1/10000.0
    # prob_of_disease = 1/1000.0

    # Probability that a person has a disease (and there's a reason to do a test)
    disease ~ flip(prob_of_disease)

    # The test is quite accurate: It shows correct result (test is positive if disease) in 99%.
    # However, in 1% of the case it shows incorrect result (positive even if there is no disease).
    # reliability = 0.99    # => disease:0.009803921568627446
    # reliability = 0.999  # A more reliable test => disease:0.09083469721767588
    # reliability = 0.9999 # An even more reliable test => disease:0.500000000000028

    test ~ disease ? flip(reliability) : flip(1-reliability)

end

function run_medical_test()
    for t in [missing,true,false],
        d in [missing,true,false],
        pd in [1/10000.0, 1/1000.0],
        r in [0.99, 0.999,0.9999]
        println("\n\ntest:$t disease:$d prob disease:$pd reliability:$r ")
        # model = medical_test(t, d,pd,r)
        model = medical_test(t, d)

        # chns = sample(model, Prior(), 40_000)
        chns = sample(model, MH(), 100_000)
        # display(chns)
        show_var_dist(chns,:test)
        show_var_dist(chns,:disease)

    end
end

@time run_medical_test()

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Base.Test

const testdatadir = dirname(@__FILE__)

REPO_URL = "https://github.com/yeesian/ArchGDALDatasets/blob/master/"

remotefiles = [
    "data/A.tif",
    "data/point.geojson",
    "data/utmsmall.tif",
    "gdalworkshop/world.tif",
    "ospy/data1/sites.dbf",
    "ospy/data1/sites.prj",
    "ospy/data1/sites.sbn",
    "ospy/data1/sites.sbx",
    "ospy/data1/sites.shp",
    "ospy/data1/sites.shx",
    "ospy/data2/ut_counties.dbf",
    "ospy/data2/ut_counties.prj",
    "ospy/data2/ut_counties.shp",
    "ospy/data2/ut_counties.shx",
    "ospy/data2/ut_counties.txt",
    "ospy/data3/cache_towns.dbf",
    "ospy/data3/cache_towns.sbn",
    "ospy/data3/cache_towns.sbx",
    "ospy/data3/cache_towns.shp",
    "ospy/data3/cache_towns.shx",
    "ospy/data3/sites.dbf",
    "ospy/data3/sites.sbn",
    "ospy/data3/sites.sbx",
    "ospy/data3/sites.shp",
    "ospy/data3/sites.shx",
    "ospy/data4/aster.img",
    "ospy/data4/aster.rrd",
    "ospy/data4/sites.dbf",
    "ospy/data4/sites.sbn",
    "ospy/data4/sites.sbx",
    "ospy/data4/sites.shp",
    "ospy/data4/sites.shx",
    "ospy/data5/aster.img",
    "ospy/data5/aster.rrd",
    "ospy/data5/doq1.img",
    "ospy/data5/doq1.rrd",
    "ospy/data5/doq2.img",
    "ospy/data5/doq2.rrd",
    "pyrasterio/example.tif",
    "pyrasterio/example2.tif",
    "pyrasterio/example3.tif",
    "pyrasterio/float_nan.tif",
    "pyrasterio/float.tif",
    "pyrasterio/RGB.byte.tif",
    "pyrasterio/shade.tif"
]

for f in remotefiles
    # create the directories if they don't exist
    currdir = dirname(f)
    isdir(currdir) || mkpath(currdir)
    # download the file
    currfile = joinpath(testdatadir, f)
    isfile(currfile) || download(REPO_URL*f*"?raw=true", currfile)
end

cd(dirname(@__FILE__)) do
    isdir("tmp") || mkpath("tmp")
    include("test_datastreams.jl")
    include("test_gdal_tutorials.jl")
    include("test_geometry.jl")
    include("test_types.jl")
    include("test_drivers.jl")
    include("test_feature.jl")
    include("test_featurelayer.jl")
    include("test_fielddefn.jl")
    include("test_styletable.jl")
    include("test_dataset.jl")
    include("test_rasterband.jl")
    include("test_rasterio.jl")
    include("test_rasterattrtable.jl")
    include("test_ospy_examples.jl")
    include("test_geos_operations.jl")
    include("test_cookbook_geometry.jl")
    include("test_cookbook_projection.jl")
end

for f in remotefiles
    currfile = joinpath(testdatadir, f)
    isfile(currfile) && rm(currfile)
end


text:
errval(::C.PyPtr) = C.PyNULL
errval(::T) where {T<:Number} = zero(T) - one(T)

iserrval(val) = val == errval(val)

iserrset() = C.PyErr_Occurred() != C.PyNULL
iserrset(val) = val == errval(val)

errcheck() = iserrset() ? pythrow() : nothing
errcheck(val) = iserrset(val) ? pythrow() : val

iserrset_ambig(val) = iserrset(val) && iserrset()

errcheck_ambig(val) = iserrset_ambig(val) ? pythrow() : val

errclear() = C.PyErr_Clear()

errmatches(t) = (@autopy t C.PyErr_ExceptionMatches(getptr(t_))) == 1

function errget()
    t = Ref(C.PyNULL)
    v = Ref(C.PyNULL)
    b = Ref(C.PyNULL)
    C.PyErr_Fetch(t, v, b)
    (pynew(t[]), pynew(v[]), pynew(b[]))
end

errset(t::Py) = GC.@preserve t C.PyErr_SetNone(getptr(t))
errset(t::Py, v::Py) = GC.@preserve t v C.PyErr_SetObject(getptr(t), getptr(v))
errset(t::Py, v::String) = GC.@preserve t C.PyErr_SetString(getptr(t), v)

function errnormalize!(t::Py, v::Py, b::Py)
    tptr = getptr(t)
    vptr = getptr(v)
    bptr = getptr(b)
    tref = Ref(tptr)
    vref = Ref(vptr)
    bref = Ref(bptr)
    C.PyErr_NormalizeException(tref, vref, bref)
    setptr!(t, tref[])
    setptr!(v, vref[])
    setptr!(b, bref[])
    (t, v, b)
end

"""
    PyException(x)

Wraps the Python exception `x` as a Julia `Exception`.
"""
mutable struct PyException <: Exception
    _t :: Py
    _v :: Py
    _b :: Py
    _isnormalized :: Bool
end
function PyException(v::Py=pybuiltins.None)
    if pyisnone(v)
        t = b = v
    elseif pyisinstance(v, pybuiltins.BaseException)
        t = pytype(v)
        b = pygetattr(v, "__traceback__", pybuiltins.None)
    else
        throw(ArgumentError("expecting a Python exception"))
    end
    PyException(t, v, b, true)
end
export PyException

ispy(x::PyException) = true
getpy(x::PyException) = x.v

pyconvert_rule_exception(::Type{R}, x::Py) where {R<:PyException} = PyException(Py(x))

function Base.show(io::IO, x::PyException)
    show(io, typeof(x))
    print(io, "(")
    show(io, x.v)
    print(io, ")")
end

function Base.getproperty(exc::PyException, k::Symbol)
    if k in (:t, :v, :b) && !exc._isnormalized
        errnormalize!(exc._t, exc._v, exc._b)
        pyisnull(exc._t) && setptr!(exc._t, incref(getptr(pybuiltins.None)))
        pyisnull(exc._v) && setptr!(exc._v, incref(getptr(pybuiltins.None)))
        pyisnull(exc._b) && setptr!(exc._b, incref(getptr(pybuiltins.None)))
        pyisnone(exc._v) || (exc._v.__traceback__ = exc._b)
        exc._isnormalized = true
    end
    k == :t ? exc._t : k == :v ? exc._v : k == :b ? exc._b : getfield(exc, k)
end

pythrow() = throw(PyException(errget()..., false))

file_to_pymodule(fname::String) = begin
    isfile(fname) || return nothing
    modules = pyimport("sys").modules
    for (n, m) in modules.items()
        if pyhasattr(m, "__file__")
            fname2 = pystr(String, m.__file__)
            if isfile(fname2) && realpath(fname) == realpath(fname2)
                return pystr(String, n)
            end
        end
    end
end

Base.showerror(io::IO, e::PyException) = _showerror(io, e, nothing, backtrace=false)

Base.showerror(io::IO, e::PyException, bt; backtrace=true) = _showerror(io, e, bt; backtrace=backtrace)

function _showerror(io::IO, e::PyException, bt; backtrace=true)
    print(io, "Python: ")

    if pyisnone(e.t)
        print(io, "mysterious error (no error was actually set)")
        return
    end

    if CONFIG.auto_sys_last_traceback
        try
            sys = pyimport("sys")
            sys.last_type = e.t
            sys.last_value = e.v
            sys.last_traceback = e.b
        catch err
            print(io, "<error while setting 'sys.last_traceback': $err>")
        end
    end

    if !pyisnull(pyJuliaError) && pyissubclass(e.t, pyJuliaError)
        # handle Julia exceptions specially
        try
            je, jb = pyconvert(Tuple{Any,Any}, e.v.args)
            print(io, "Julia: ")
            if je isa Exception
                showerror(io, je, jb, backtrace=backtrace && jb !== nothing)
            else
                print(io, je)
                backtrace && jb !== nothing && Base.show_backtrace(io, jb)
            end
        catch err
            println("<error while printing Julia exception inside Python exception: $err>")
        end
    else
        # print the type name
        try
            print(io, e.t.__name__)
        catch err
            print(io, "<error while printing type: $err>")
        end

        # print the error message
        if !pyisnone(e.v)
            print(io, ": ")
            try
                print(io, e.v)
            catch err
                print(io, "<error while printing value: $err>")
            end
        end
    end

    if backtrace
        # print the Python stacktrace
        println(io)
        printstyled(io, "Python stacktrace:")
        if pyisnone(e.b)
            printstyled(io, " none")
        else
            try
                fs = [(pystr(String, x.name), pystr(String, x.filename), pystr(String, x.lineno)) for x in pyimport("traceback").extract_tb(e.b)]
                if Base.VERSION < v"1.6.0-rc1"
                    for (i, (name, fname, lineno)) in enumerate(reverse(fs))
                        println(io)
                        printstyled(io, " [", i, "] ")
                        printstyled(io, name, bold = true)
                        printstyled(io, " at ")
                        printstyled(io, fname, ":", lineno, bold = true)
                    end
                else
                    mcdict = Dict{String, Symbol}()
                    mccyclyer = Iterators.Stateful(Iterators.cycle(Base.STACKTRACE_MODULECOLORS))
                    # skip a couple as a basic attempt to make the colours different from the Julia stacktrace
                    popfirst!(mccyclyer)
                    popfirst!(mccyclyer)
                    for (i, (name, fname, lineno)) in enumerate(reverse(fs))
                        println(io)
                        printstyled(io, " [", i, "] ")
                        printstyled(io, name, bold = true)
                        println(io)
                        printstyled(io, "   @ ", color = :light_black)
                        mod = file_to_pymodule(fname)
                        if mod !== nothing
                            # print the module, with colour determined by the top level name
                            tmod = first(split(mod, ".", limit=2))
                            color = get!(mcdict, tmod) do
                                popfirst!(mccyclyer)
                            end
                            printstyled(io, mod, " ", color = color)
                        end
                        if isfile(fname) && :stacktrace_contract_userdir in names(Base, all=true) && Base.stacktrace_contract_userdir()
                            if :replaceuserpath in names(Base, all=true)
                                fname = Base.replaceuserpath(fname)
                            elseif :contractuser in names(Base.Filesystem, all=true)
                                fname = Base.Filesystem.contractuser(fname)
                            end
                        end
                        printstyled(io, fname, ":", lineno, color = :light_black)
                    end
                end
            catch err
                print(io, "<error while printing stacktrace: $err>")
            end
        end

        # print the Julia stacktrace
        Base.show_backtrace(io, bt)
    end
end


text:
# This file was generated, do not modify it. # hide
y, X = unpack(df, ==(:Churn), !=(:customerID));
schema(X).names

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
function StatsBase.addcounts!{T,U}(cm::Dict{U,Int}, x::AbstractDataArray{T})
    for v in x
        cm[v] = get(cm, v, 0) + 1
    end
    return cm
end

function StatsBase.addcounts!{T,U,W}(cm::Dict{U,W}, x::AbstractDataArray{T}, wv::WeightVec{W})
    n = length(x)
    length(wv) == n || raise_dimerror()
    w = values(wv)
    z = zero(W)

    @inbounds for i in 1:n
        xi = x[i]
        wi = w[i]
        cm[xi] = get(cm, xi, z) + wi
    end
    return cm
end

function StatsBase.countmap{T}(x::AbstractDataArray{T})
    addcounts!(Dict{Union(T, NAtype), Int}(), x)
end

function StatsBase.countmap{T,W}(x::AbstractDataArray{T}, wv::WeightVec{W})
    addcounts!(Dict{Union(T, NAtype), W}(), x, wv)
end

function cut{S, T}(x::AbstractVector{S}, breaks::Vector{T})
    if !issorted(breaks)
        sort!(breaks)
    end
    min_x, max_x = minimum(x), maximum(x)
    if breaks[1] > min_x
        unshift!(breaks, min_x)
    end
    if breaks[end] < max_x
        push!(breaks, max_x)
    end
    refs = fill(zero(DEFAULT_POOLED_REF_TYPE), length(x))
    for i in 1:length(x)
        if x[i] == min_x
            refs[i] = 1
        else
            refs[i] = searchsortedfirst(breaks, x[i]) - 1
        end
    end
    n = length(breaks)
    from = map(x -> sprint(showcompact, x), breaks[1:(n - 1)])
    to = map(x -> sprint(showcompact, x), breaks[2:n])
    pool = Array(ASCIIString, n - 1)
    if breaks[1] == min_x
        pool[1] = string("[", from[1], ",", to[1], "]")
    else
        pool[1] = string("(", from[1], ",", to[1], "]")
    end
    for i in 2:(n - 1)
        pool[i] = string("(", from[i], ",", to[i], "]")
    end
    PooledDataArray(RefArray(refs), pool)
end

cut(x::AbstractVector, ngroups::Integer) = cut(x, quantile(x, [1 : ngroups - 1] / ngroups))

function rep{T <: Integer}(x::AbstractVector, lengths::AbstractVector{T})
    if length(x) != length(lengths)
        error("vector lengths must match")
    end
    res = similar(x, sum(lengths))
    i = 1
    for idx in 1:length(x)
        tmp = x[idx]
        for kdx in 1:lengths[idx]
            res[i] = tmp
            i += 1
        end
    end
    return res
end

function rep(x::AbstractVector, times::Integer = 1, each::Integer = 1)
    res = similar(x, each * times * length(x))
    i = 1
    for jdx in 1:times
        for idx in 1:length(x)
            tmp = x[idx]
            for kdx in 1:each
                res[i] = tmp
                i += 1
            end
        end
    end
    return res
end

function rep(x::AbstractVector; times::Integer = 1, each::Integer = 1)
    rep(x, times, each)
end

rep(x::Any, times::Integer) = fill(x, times)


text:
let
    mdp = SimpleGridWorld()

    pomdp = FullyObservablePOMDP(mdp)

    @test observations(pomdp) == states(pomdp)
    @test statetype(pomdp) == obstype(pomdp)

    @test observations(pomdp) == states(pomdp)
    @test statetype(pomdp) == obstype(pomdp)

    s_po = rand(MersenneTwister(1), initialstate(pomdp))
    s_mdp = rand(MersenneTwister(1), initialstate(mdp))
    @test s_po == s_mdp

    solver = ValueIterationSolver(max_iterations = 100)
    @test_skip begin
        mdp_policy = solve(solver, mdp)
        pomdp_policy = solve(solver, UnderlyingMDP(pomdp))
        mdp_policy.util == pomdp_policy.util
    end

    is = rand(MersenneTwister(3), initialstate(mdp))
    for (sp, o, r) in stepthrough(pomdp,
                               FunctionPolicy(o->:left),
                               PreviousObservationUpdater(),
                               is, is, "sp,o,r",
                               rng=MersenneTwister(2),
                               max_steps=10)
        @test sp == o
    end
end