text:
function parameterization(inputs)
"""
Defines model parameters
"""

    # Parameters for beliefs about parameters
    # From /estimate_damage_parameters/table_1.do
    coeff_mean = .00556486
    coeff_var =  .00001433
    exp_mean = 1.882038
    exp_var = .4505272^2
    damage_shock_sd = 1.498293

    if !inputs[:robust]
        theta = 0.0
    end

    if inputs[:perfect_info]
        coeff_var = 1e-12
        exp_var = 1e-12
    end

    # Initialize model parameter dictionary
    mparams = Dict{Symbol,Any}( :theta_robust => inputs[:theta],
        :exp_mean => exp_mean, :exp_var => exp_var,
        :cscert => 3, :b2 => 0.0028388,
        :b3 => 2., :B0 => 1.1, :gB => -0.01, :a0 => 1.17,
        :a1 => 2., :a2 => 2.8, :gPsi => -0.005, :delta => 0.015, :beta => 1/(1.015),
        :L0 => 6514, :dL => 0.035, :Linf => 8600, :kappa => 0.3,
        :A0_N => 0.02722, :gA_N => 0.0092, :dA => 0.001,
        :sigma0 => 0.13418, :gsigma0 => -0.00730, :dsigma => 0.003,
        :deltak => 0.1, :EF0 => -0.06, :EF1 => 0.30,
        :lambda_0 => 0.315, :forcing_per_2xco2 => 3.8,
        :gtc_per_ppm => 2.16, :nqnodes => inputs[:nqnodes],
        :new_bounds_weight_up => inputs[:new_bounds_weight_up],
        :new_bounds_weight_down => inputs[:new_bounds_weight_down],
        :increase_adapt_iter => inputs[:increase_adapt_iter],
        :Mstart => 454.901, :Tstart => 0.7307,
        :increase_dims => inputs[:increase_dims], :increase_dims_adapt => inputs[:increase_dims_adapt],
        :increase_level => inputs[:increase_level], :increase_level_adapt => inputs[:increase_level_adapt],
        :initial_approximation_level => inputs[:initial_approximation_level],
        :approximation_level => inputs[:initial_approximation_level], :rho => -1., :deltak => 0.1,
        :ccr=> 0.0016, :time_horizon => inputs[:time_horizon], :damages => inputs[:damages],
        :sim_length => inputs[:sim_length], :num_runs => inputs[:num_runs], :root_dir => pwd(),
        :num_runs_adapt => inputs[:num_runs_adapt], :sim_length_adapt => inputs[:sim_length_adapt],
        :num_adapts => inputs[:num_adapts], :current_adapt => 0, :scaling_factor => inputs[:scaling_factor],
        :expand_domain => inputs[:expand_domain])

    # Initialize model options dictionary
    opts = Dict{Symbol,Bool}(:perfect_info => inputs[:perfect_info], :robust => inputs[:robust],
        :learn_e => inputs[:learn_e], :sim => false, :first_it => true,
        :adaptive_grid => inputs[:adapt], :adapted => false, :load => inputs[:load])

    # Parameters for misspecified runs: pure damage shocks and weitzman coefficients and exponents
    mparams[:shock_loc] = -0.59
    mparams[:shock_scale] = 1.18
    mparams[:weitz_c1] = 20.46
    mparams[:weitz_c2] = 6.081
    mparams[:weitz_e1] = 2
    mparams[:weitz_e2] = 6.754
    mparams[:dietz_c1] = 0.00284
    mparams[:dietz_c2] = 8.19e-5
    mparams[:dietz_e1] = 2
    mparams[:dietz_e2] = 6.754
    mparams[:extreme_c1] = 0.00284
    mparams[:extreme_c2] = 5.85e-4
    mparams[:extreme_e1] = 2
    mparams[:extreme_e2] = 6.754

    # Location and scale parameter for the joint d1-omega distribution
    mparams[:coeff_mean] = coeff_mean
    mparams[:coeff_loc] = log.(coeff_mean/sqrt(1 + coeff_var/coeff_mean^2)) .- log.(damage_shock_sd^2 + 1)/2.
    mparams[:coeff_scale] = log.(1 + coeff_var/coeff_mean^2) .+ log.(damage_shock_sd^2 + 1)

    # Initial production technology
    mparams[:A0] = mparams[:A0_N]^(1/(1-mparams[:kappa]))

    # Initial growth rate of production technology
    mparams[:gA0] = mparams[:gA_N]/(1-mparams[:kappa])

    # Initial capital
    mparams[:kstart] = 137/mparams[:A0]/mparams[:L0]

    # Initial fractional net output
    mparams[:lossstart] = 1/(1 + exp.(mparams[:coeff_loc])*mparams[:Tstart]^mparams[:exp_mean])

    # Quadrature scheme for joint d1-omega distribution
    de, dw = qnwlogn(mparams[:nqnodes], mparams[:coeff_loc], mparams[:coeff_scale])
    mparams[:coeff] = de
    mparams[:weights] = dw

    # Initial belief conditions
    mparams[:expstart] = mparams[:exp_mean]
    mparams[:varstart] = mparams[:exp_var]

    # Build domain
    if opts[:perfect_info]

        # Set exponent to mean
        mparams[:exp_var] = 1e-12
        mparams[:exp_mean] = 2.0

        # Restrict approximation domain
        domain_loc = mparams[:exp_mean]
        domain_scale = mparams[:exp_var]*1.01
        loc_adjust = .0001
        scale_adjust = .00001

    elseif opts[:learn_e]

        # Restrict approximation domain
        domain_loc = mparams[:exp_mean]
        domain_scale = mparams[:exp_var]*1.01
        loc_adjust = 1.25
        scale_adjust = domain_scale

    elseif ~opts[:perfect_info]

        # Restrict approximation domain
        domain_loc = mparams[:exp_mean]
        domain_scale = mparams[:exp_var]*1.01
        loc_adjust = .0001
        scale_adjust = .00001

    end

    sc = mparams[:expand_domain]

    # Set bounds on approximation domain
    # States in order: capital, cumulative emissions, prior_loc parameter, scale parameter, fractional net output
    if opts[:learn_e]
        mparams[:smin] = [0.00, mparams[:Mstart]-.01,  domain_loc-loc_adjust, domain_scale-scale_adjust, 0.4]
        mparams[:smax] = [6.00, 5000,                  domain_loc+loc_adjust, domain_scale,              1.0]
    else
        mparams[:smin] = [0.65/sc, (mparams[:Mstart]-.01)/sc, 0.3/sc/sc]
        mparams[:smax] = [4.30*sc, 3500*sc,                 1.0]
        mparams[:approximation_level] = mparams[:approximation_level][[1 2 5]]
        mparams[:increase_dims] = mparams[:increase_dims][[1 2 5]]
    end

    # Build approximation space
    space = fundefn(:cheb, mparams[:approximation_level], mparams[:smin], mparams[:smax])
    states = Array{Float64}(undef, prod(mparams[:approximation_level]), length(mparams[:approximation_level]))
    states, bounds = funnode(space)
    space[:num_nodes] = size(states, 1)
    space[:c_cont] = zeros(space[:num_nodes])
    Binv = funbase(space, states)\Matrix{Float64}(I, space[:num_nodes], space[:num_nodes])

    println("Approximating with $(space[:num_nodes]) nodes.")
    if opts[:robust]
        println("Penalty parameter: $(mparams[:theta_robust])")
    end
    println("Domain upper bound: $(mparams[:smax]).")
    println("Domain lower bound: $(mparams[:smin]).")

    ## Compute exogenously evolving variables
    t = 0:mparams[:time_horizon]

    # Non-industrial emissions
    exog_states = Dict{Symbol,Array{Float64}}(:Bs => mparams[:B0]*exp.(mparams[:gB]*t))

    # Labor
    exog_states[:Lt] = mparams[:L0] .+ (mparams[:Linf]-mparams[:L0])*(1 .- exp.(-mparams[:dL]*t))

    # Labor growth
    exog_states[:gL] = mparams[:dL]./((mparams[:Linf]/(mparams[:Linf] .- mparams[:L0]))*exp.(mparams[:dL]*t) .- 1)

    # Technological growth rate
    exog_states[:gA] = mparams[:gA0]*exp.(-mparams[:dA]*t)

    # Technology
    exog_states[:At] = mparams[:A0]*exp.( (mparams[:gA0]/mparams[:dA])*(1 .- exp.(-mparams[:dA]*t)))

    # Growth rate for exogenous growth in labor and productivity
    exog_states[:effkgrowth] = exp.(-(exog_states[:gA] .+ exog_states[:gL]));

    # Effective discount factor
    exog_states[:betaeff] = exp.((-mparams[:delta] .+ exog_states[:gA]*mparams[:rho] .+ exog_states[:gL]))

    # Exogenous decarbonization of production
    exog_states[:sigma] = mparams[:sigma0].*exp.((mparams[:gsigma0]./mparams[:dsigma])*(1 .- exp.(-mparams[:dsigma].*t)))
    exog_states[:emint] = exog_states[:sigma].*exog_states[:At].*exog_states[:Lt]

    # Abatement cost function coefficient
    exog_states[:Psi] = (mparams[:a0].*exog_states[:sigma]./mparams[:a2]).*(1 .- (1 .- exp.(mparams[:gPsi].*t))./mparams[:a1])

    # External forcing
    exog_states[:EF] = mparams[:EF0] .+ 0.01*(mparams[:EF1] .- mparams[:EF0])*min.(t, 100)

    # External forcing in next year
    exog_states[:EF_next] = mparams[:EF0] .+ 0.01*(mparams[:EF1] .- mparams[:EF0])*min.(t .+ 1, 100)

    return states, mparams, exog_states, opts, space, Binv

end

max_stars_repo_name: irudik/optimal-climate-policy-aej

max_issues_repo_name: irudik/optimal-climate-policy-aej

max_forks_repo_name: irudik/optimal-climate-policy-aej


text:
module TF

import ..d; import ..dxasy

sigmoid{F <: AbstractFloat}(x::F) = one(F)/(one(F) + exp(-x))
dasysigmoid{F <: AbstractFloat}(y::F) = y * (one(F) - y)
#explicitly reinterpret the dasysigmoid.
dsigmoid{F <: AbstractFloat}(x::F) = dasysigmoid(sigmoid(x))

id{F <: AbstractFloat}(x::F) = x

softplus{F <: AbstractFloat}(x::F) = log(one(F) + exp(x))
dasysoftplus{F <: AbstractFloat}(y::F) = (exp(y) - one(F)) / exp(y)

#broadcast functions over arrays.
sigmoid{F <: AbstractFloat}(a::Array{F})      = sigmoid.(a)
dasysigmoid{F <: AbstractFloat}(a::Array{F})  = dasysigmoid.(a)
dsigmoid{F <: AbstractFloat}(a::Array{F})     = dsigmoid.(a)
id{F <: AbstractFloat}(a::Array{F})           = a
softplus{F <: AbstractFloat}(a::Array{F})     = softplus.(a)
dasysoftplus{F <: AbstractFloat}(a::Array{F}) = dasysoftplus.(a)

#non-simple transfer functions.
function softmax{F <: AbstractFloat}(a::Array{F})
  #go over the entire array and add up the values.
  totalexpsum = zero(F)
  for idx = 1:length(a)
    totalexpsum += exp(a[idx])
  end
  exp.(a) / totalexpsum
end
dsoftmax{F <: AbstractFloat}(a::Array{F}) = dasysigmoid.(softmax(a))

dxasy(::typeof(sigmoid)) = dasysigmoid
dxasy(::typeof(id)) = one
dxasy(::typeof(softplus)) = dasysoftplus
dxasy(::typeof(softmax)) = dasysigmoid

d(::typeof(sigmoid)) = dsigmoid
d(::typeof(id)) = one
d(::typeof(softplus)) = sigmoid
d(::typeof(softmax)) = dsoftmax

end

ext: jl

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
function prefix_sum(arr::Vector{T})where T<:Number
    pre=[]
    preans=zero(T)
    for i in arr
        preans+=i
        push!(pre,preans)
    end
    return pre
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text:
# SPDX-License-Identifier: MIT

"""
    Base.union(f::AbstractFMS{T}...)

Take the union of the provided FSMs.
"""
function Base.union(fsm1::AbstractFSM{T}, fsm2::AbstractFSM{T}) where T
    allstates = union(collect(states(fsm1)), collect(states(fsm2)))
    newfsm = VectorFSM{T}()

    smap1 = Dict()
    for state in states(fsm1)
        smap1[state] = addstate!(newfsm, state.label,
                                initweight = state.initweight,
                                finalweight = state.finalweight)
    end

    smap2 = Dict()
    for state in states(fsm2)
        smap2[state] = addstate!(newfsm, state.label,
                                initweight = state.initweight,
                                finalweight = state.finalweight)
    end

    for src in states(fsm1)
        for arc in arcs(fsm1, src)
            addarc!(newfsm, smap1[src], smap1[arc.dest], arc.weight)
        end
    end

    for src in states(fsm2)
        for arc in arcs(fsm2, src)
            addarc!(newfsm, smap2[src], smap2[arc.dest], arc.weight)
        end
    end

    newfsm
end
Base.union(f::AbstractFSM{T}, o::AbstractFSM{T}...) where T =
    foldl(union, o, init = f)

"""
    renormalize(fsm::AbstractFSM{T}) where T<:Semiring

Ensure the that the weights of all the outgoing arcs leaving a
state sum up to `one(T)`.
"""
function renormalize(fsm::AbstractFSM{T}) where T<:Semiring
    total_initweight = zero(T)
    for s in filter(isinit, states(fsm)) total_initweight += s.initweight end

    totals = Dict()
    for s in states(fsm)
        v = s.finalweight
        for a in arcs(fsm, s)
            v += a.weight
        end
        totals[s] = v
    end

    newfsm = VectorFSM{T}()
    smap = Dict()
    for s in states(fsm)
        iw = total_initweight == zero(T) ? zero(T) : s.initweight / total_initweight
        fw = totals[s] == zero(T) ? zero(T) : s.finalweight / totals[s]
        smap[s] = addstate!(newfsm, s.label, initweight = iw, finalweight = fw)
    end

    for s in states(fsm)
        for a in arcs(fsm, s)
            addarc!(newfsm, smap[s], smap[a.dest], a.weight / totals[s])
        end
    end

    newfsm
end

function next_states(fsm, set_states)
    nextstates = []
    for s in set_states
        for arc in arcs(fsm, s)
            push!(nextstates, (arc.dest, arc.weight))
        end
    end
    nextstates
end

function nextlabels(T, statelist; init = false)
    labels = Dict()
    for (s, w) in statelist
        set_states, iw, fw, sw = get(labels, s.label, (Set(), zero(T), zero(T), zero(T)))
        push!(set_states, s)
        labels[s.label] = (set_states, iw+s.initweight, fw+s.finalweight, sw+w)
    end

    retval = []
    for (label, (set_states, iw, fw, sw)) in labels
        if ! init
            push!(retval, (label, set_states, zero(T), fw, sw))
        else
            push!(retval, (label, set_states, iw, fw, sw))
        end
    end
    retval
end

"""
    determinize(fsm::AbstractFSM)

Determinize the FSM w.r.t. the state labels. A FSM is deterministic
if for any state there is no next states with the same label.
"""
function determinize(fsm::AbstractFSM{T}; renormalize_output = true) where T
    dfsm = VectorFSM{T}()
    smap = Dict()
    newarcs = Dict()
    visited = Set()

    init_ws = [(s, one(T)) for s in filter(isinit, collect(states(fsm)))]
    queue = nextlabels(T, init_ws; init = true)
    while ! isempty(queue)
        label, set_states, iw, fw, w = popfirst!(queue)

        nl= nextlabels(T, next_states(fsm, set_states); init = false)
        for (next_label, next_set_states, next_iw, next_fw, next_w) in nl
            newarcs[(set_states, next_set_states)] = next_w
            if next_set_states ∉ visited
                val = (next_label, next_set_states, next_iw, next_fw, next_w)
                push!(queue, val)
                push!(visited, next_set_states)
            end
        end

        if set_states ∉ keys(smap)
            smap[set_states] = addstate!(dfsm, label, initweight = iw, finalweight = fw)
        end
    end

    for ((src, dest), weight) in newarcs
        addarc!(dfsm, smap[src], smap[dest], weight)
    end

    renormalize_output ? dfsm |> renormalize : dfsm
end
unnorm_determinize(fsm) = determinize(fsm, renormalize_output = false)

function previous_states(r_arc_map, set_states)
    nextstates = []
    for s in set_states
        for (src, dest, weight) in r_arc_map[s]
            push!(nextstates, (src, weight))
        end
    end
    nextstates
end

function reverse_determinize(fsm::AbstractFSM{T}) where T
    r_arc_map = Dict()
    for s in states(fsm)
        for a in arcs(fsm, s)
            arclist = get(r_arc_map, a.dest, [])
            push!(arclist, (s, a.dest, a.weight))
            r_arc_map[a.dest] = arclist
        end
    end

    dfsm = VectorFSM{T}()
    smap = Dict()
    newarcs = Dict()
    visited = Set()

    init_ws = [(s, one(T)) for s in filter(isfinal, collect(states(fsm)))]
    queue, _ = nextlabels(T, init_ws; init = false)
    while ! isempty(queue)
        label, set_states, iw, fw, w = popfirst!(queue)

        nl, total = nextlabels(T, previous_states(r_arc_map, set_states);
                               init = true, finalweight = fw)
        for (next_label, next_set_states, next_iw, next_fw, next_w) in nl
            newarcs[(set_states, next_set_states)] = next_w
            if next_set_states ∉ visited
                val = (next_label, next_set_states, next_iw, next_fw, next_w)
                push!(queue, val)
                push!(visited, next_set_states)
            end
        end

        if set_states ∉ keys(smap)
            #nfw = (total ≠ zero(T) ) ? fw / total : fw
            nfw = fw
            smap[set_states] = addstate!(dfsm, label, initweight = iw,
                                         finalweight = nfw)
        end
    end

    for ((dest, src), weight) in newarcs
        addarc!(dfsm, smap[src], smap[dest], weight)
    end

    dfsm
end

"""
    transpose(fsm::AbstractFSM)

Reverse the direction of the arcs and, for each state, inverse the
initial and final weight.
"""
function Base.transpose(fsm::AbstractFSM{T}) where T
    newfsm = VectorFSM{T}()
    smap = Dict()
    for s in states(fsm)
        ns = addstate!(newfsm, s.label, initweight = s.finalweight,
                       finalweight = s.initweight)
        smap[s] = ns
    end

    for src in states(fsm)
        for arc in arcs(fsm, src)
            addarc!(newfsm, smap[arc.dest], smap[src], arc.weight)
        end
    end

    newfsm
end

"""
    minimize(fsm::AbstractFSM)

Return a minimal equivalent fsm.
"""
minimize(fsm::AbstractFSM{T}) where T =
    (renormalize ∘ transpose ∘ unnorm_determinize ∘
     transpose ∘ unnorm_determinize)(fsm)

"""
    label_closure!(closure, fsm, state, label;  [weight=1], [visited=[]])

Find label closure from `state` in `fsm`.
"""
function label_closure!(
        closure::Vector, fsm::AbstractFSM{T}, state::State, label;
        weight::T=one(T), visited::Vector{State} = State[],
        matchfn = (x,y) -> x != y
) where T <: Semiring

    if state in visited
        return closure
    end
	push!(visited, state)

    for l in arcs(fsm, state)
        if matchfn(l.dest.label, label)
            push!(closure, (l.dest, l.weight * weight))
        else
            label_closure!(closure, fsm, l.dest, label; weight=l.weight * weight,
                           visited=visited, matchfn = matchfn)
        end
    end
    return closure
end

"""
	remove_label(fsm, label)

Removes all states from`fsm` with label `label`.
"""
function remove_label(fsm::AbstractFSM{T}, label) where T <: Semiring
    nfsm = VectorFSM{T}()
    label_states = []

    iw, fw = Dict(), Dict()
    for s in states(fsm)
        if s.label != label
            iw[s] = s.initweight

            closure = label_closure!([], fsm, s, label)
            unique!(closure)
            nfw = zero(T)
            for (ns, w) in closure
                nfw += w*ns.finalweight
            end
            fw[s] = s.finalweight + nfw
        end
    end

    label_closures = Dict{State, Vector}()
    for s in states(fsm)
        if s.label == label
            closure = label_closure!([], fsm, s, label)
            label_closures[s] = unique!(closure)

            for (ns, nw) in label_closures[s]
                iw[ns] = iw[ns] + s.initweight*nw
            end
        end
    end

    smap = Dict{State, State}()
    for s in keys(iw)
        smap[s] = addstate!(nfsm, s.label; initweight = iw[s], finalweight = fw[s])
    end

    for s in states(fsm)
        for l in arcs(fsm, s)
            if s.label != label && l.dest.label != label
                addarc!(nfsm, smap[s], smap[l.dest], l.weight)
            elseif s.label != label
                for (ns, w) in label_closures[l.dest]
                    addarc!(nfsm, smap[s], smap[ns], l.weight * w)
                end
            end
        end
    end
    return nfsm
end

max_stars_repo_name: BUTSpeechFIT/HiddenMarkovModel

max_issues_repo_name: BUTSpeechFIT/HiddenMarkovModel

max_forks_repo_name: BUTSpeechFIT/ProbabilisticFSMs


text:
"""
    Check feasibility of Matpower solution
"""

include(joinpath(ROOT,"..","src_PowSysMod", "PowSysMod_body.jl"))
include("read_QCQP.jl")

instances = [name[1:end-4] for name in readdir("instances")]

for instance_name in instances
  # Building complex and real problems
  filepath = abspath(joinpath("instances", instance_name*".dat"))
  pb = read_QCQP(filepath)
  pb_R = pb_cplx2real(pb)

  # Building complex and real solutions
  filepath = abspath(joinpath("solutions", instance_name*".sol"))
  pt_sol_C, pt_sol_R = build_QCQP_sol(filepath)

  # Compute and build minSlack and objective for both cplx and real problems
  minSlack, minCstrName = get_minslack(pb, pt_sol_C)
  val = get_objective(pb, pt_sol_C)
  @printf("%19s min slack: %.3e at %10s\t objective value : %.10e + im %.10e\n", instance_name, minSlack, minCstrName[1], real(val), imag(val))

  minSlack_R, minCstrName_R = get_minslack(pb_R, pt_sol_R)
  val_R = get_objective(pb_R, pt_sol_R)
  @printf("%19s min slack: %.3e at %10s\t objective value : %.10e + im %.10e\n", "C2R "*instance_name, minSlack_R, minCstrName_R[1], real(val_R), imag(val_R))
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
# Autogenerated wrapper script for GLFW_jll for powerpc64le-linux-gnu
export libglfw

using Libglvnd_jll
using Xorg_libXcursor_jll
using Xorg_libXi_jll
using Xorg_libXinerama_jll
using Xorg_libXrandr_jll
JLLWrappers.@generate_wrapper_header("GLFW")
JLLWrappers.@declare_library_product(libglfw, "libglfw.so.3")
function __init__()
    JLLWrappers.@generate_init_header(Libglvnd_jll, Xorg_libXcursor_jll, Xorg_libXi_jll, Xorg_libXinerama_jll, Xorg_libXrandr_jll)
    JLLWrappers.@init_library_product(
        libglfw,
        "lib/libglfw.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import PredictMD
import Test

Test.@test PredictMD.is_make_examples() isa Bool
Test.@test PredictMD.is_make_docs() isa Bool
Test.@test PredictMD.is_deploy_docs() isa Bool
Test.@test PredictMD.is_docs_or_examples() == PredictMD.is_make_examples() ||
    PredictMD.is_make_docs() ||
    PredictMD.is_deploy_docs()


text:
"""
    ComtradeData

`ComtradeData` holds data and configuration information from a set of COMTRADE
files. Fields include:

* `filename::String` -- Base file name (without extension).
* `cfg::ComtradeCfg` -- Configuration information.
* `dat::DataFrame` -- Table of results with each channel as a column.
* `hdr::String` -- Contents of the "hdr" file.
* `inf::String` -- Contents of the "inf" file.
"""

using CSV

mutable struct ComtradeData
    filename::String
    cfg::ComtradeCfg
    dat::DataFrame
    hdr::String
    inf::String
end

function addcol!(d::DataFrame, x)
    d[!,Symbol(string("x", ncol(d) + 1))]=x
end

"""
    read_comtrade(basename)

Read contents of a set of COMTRADE files with base file name `basename` (no
extension).

Returns a `COMTRADE.ComtradeData` with the following fields:

* `filename::String` -- Base file name (without extension).
* `cfg::ComtradeCfg` -- Configuration information.
* `dat::DataFrame` -- Table of results with each channel as a column.
* `hdr::String` -- Contents of the "hdr" file.
* `inf::String` -- Contents of the "inf" file.

See also `COMTRADE.ComtradeCfg` for a description of the `cfg` field.
"""
function read_comtrade(basename)
    cfg = read_cfg("$basename.cfg")
    if cfg.ft == "ASCII"
        d = CSV.read("$basename.dat", delim = ',', header = false)
        #for i in 1:cfg.nA+2   # convert DataArray columns to Arrays
        #    d[i] = d[i]
        #end
        for i in (1 + 2 + cfg.nA):(cfg.nD + 2 + cfg.nA)
            d[!,i] = map(Bool, d[!,i])
        end
    elseif cfg.ft == "BINARY"
        x = reinterpret(Int16, read(open("$basename.dat")))
        n = cfg.nA + div(cfg.nD, 16) + 2 + 2
        x = reshape(x[1:(n * cfg.npts)], (n, cfg.npts))
        global d = DataFrame()
        addcol!(d, reshape(reinterpret(UInt32, x[1:2, :]), (cfg.npts,)))
        #addcol!(d, reinterpret(UInt32, x[1:2, :]))

        addcol!(d, reshape(reinterpret(UInt32, x[3:4, :]), (cfg.npts,)))
        #addcol!(d, reinterpret(UInt32, x[3:4, :]))

        for i = 1:cfg.nA
            addcol!(d, reshape(x[i+4, :], (cfg.npts,)))
        end
        for i = 1:cfg.nD
            addcol!(d, Bool[(x >>> (mod1(i, 16) - 1)) & 0x01 for x in x[cfg.nA + 4 + cld(i, 16), :]])
            # addcol!(d, Bool[(x >>> (16 - mod1(i, 16))) & 0x01 for x in x[cfg.nA + 4 + cld(i, 16), :]])
        end
    else
        error("Unsupported data type $(x.ft)")
    end
    newnames = DataFrames.make_unique(vcat(:n, :time, [Symbol(x) for x in vcat(cfg.A[!,:ch_id], cfg.D[!,:ch_id])]),makeunique=true)
    rename!(d, newnames)
    # analog multipliers
    for i in 1:cfg.nA
        d[!,i+2] = cfg.A[i,:a] * d[!,i+2] .+ cfg.A[i,:b]
    end
    # analog multipliers
    d[!,2] = 1e-6 * cfg.timemult * d[!,2]
    hdr = isfile("$basename.hdr") ? read("$basename.hdr",String) : ""
    inf = isfile("$basename.inf") ? read("$basename.inf",String) : ""
    hdr= ""
    inf= ""
    return ComtradeData(basename, cfg, d, hdr, inf)
end

"""
    S"some name"

Returns a symbol. Equivalent to `symbol("some name")`. This helps when COMTRADE
data columns use nonstandard names.
"""
macro S_str(arg)
    Meta.quot(Symbol(arg))
end


text:
# Utils to lift LinearMaps to LinearMaps acting on larger vector spaces
import LinearMaps
import LinearAlgebra

struct LiftedMap{T,TI,TJ} <: LinearMap{T}
    A::LinearMap{T}
    I::TI
    J::TJ
end

function LinearAlgebra.mul!(y::AbstractVector, L::LiftedMap, x::AbstractVector, α::Number, β::Number)

    yI = view(y, L.I)
    xJ = view(x, L.J)
    AIJ = L.A
    LinearAlgebra.mul!(yI, AIJ, xJ, α, β)
end


text:
struct RestrictedMeasure{F,M} <: AbstractMeasure
    f::F
    base::M
end

@inline function logdensity_def(d::RestrictedMeasure, x)
    d.f(x) || return -Inf
    return 0.0
end

function density_def(d::RestrictedMeasure, x)
    d.f(x) || return 0.0
    return 1.0
end

basemeasure(μ::RestrictedMeasure) = μ.base

basemeasure_type(::Type{RestrictedMeasure{F,M}}) where {F,M} = M

insupport(μ::RestrictedMeasure, x) = μ.f(x) && insupport(μ.base, x)

restrict(f, μ) = RestrictedMeasure(f, μ)


text:
using DH32ParallelUtils
using Test

# write your own tests here

### split_int ###
ix_start, ix_end = DH32ParallelUtils.split_int(2000, 3)
@test ix_start == [1,667,1334]
@test ix_end == [666,1333,2000]

ix_start, ix_end = DH32ParallelUtils.split_int(2000, 2)
@test ix_start == [1,1000]
@test ix_end == [999,2000]

### split_int_lengths ###
Ls = DH32ParallelUtils.split_int_lengths(2000, 3)
@test sum(Ls) == 2000
@test length(Ls) == 3

Ls = DH32ParallelUtils.split_int_lengths(2000, 2)
@test sum(Ls) == 2000
@test length(Ls) == 2


text:
# function bytes_to_unicode()
#   cs = vcat(collect(('!'):('~')), collect(('¡'):('¬')), collect(('®'):('ÿ')))
#   bs = map(Int, cs)
#   cs = bs[:]
#   n = 0
#   for b in 0:2^8 - 1
#     if !(b in bs)
#       push!(bs, b)
#       push!(cs, 2^8 + n)
#       n += 1
#     end
#   end
#   cs = map(Char, cs)
#   return Dict(zip(bs, cs))
# end

"the codemap used by openai gpt2"
default_codemap() = CodeMap([(0:32, 256:288), (127:160, 289:322), 173=>323])

"the tokenizer used by openai gpt2"
function gpt2_tokenizer(text)
  pattern = r"'s|'t|'re|'ve|'m|'ll|'d| ?\p{L}+| ?\p{N}+| ?[^\s\p{L}\p{N}]+|\s+(?!\S)|\s+"
  return map(x->intern(x.match), eachmatch(pattern, text))
end

"simply the built-in split function for the origin tokenize method in subword-nmt"
whitespace_tokenize(str::AbstractString) = split(str)


text:
###################################################
## geo/projectors/nearestnode.jl
## implementation of network projector with nn projection
###################################################

"""
    `NearestNode` : Projects trips to the nearest nodes, give the associated prediction
"""
mutable struct NearestNode <: NetworkProjector
    # compulsory attributes
    network::Network
    trips::GeoData

    "projection KD-tree"
    tree::KDTree
    "precomputation (trip->node projections)"
    proj::Vector{Tuple{Int,Int}}

    function NearestNode(n::Network)
        obj = new()
        obj.network = n
        obj.trips = GeoData()
        obj.proj  = Tuple{Int,Int}[]
        # Constructing tree
        dataPos = Array{Float64}(undef, (2,length(n.nodes)))
        for (i,node) in enumerate(n.nodes)
           dataPos[1,i] = node.x
           dataPos[2,i] = node.y
        end
        obj.tree = KDTree(dataPos)
        return obj
    end
end

"""
    `NearestNode` extended constructor: also preloads the data
"""
function NearestNode(n::Network, trips::GeoData)
    nn = NearestNode(n)
    preloadData!(nn, trips)
    return nn
end

"""
    `preloadData!`: project all trips onto their nearest node
"""
function preloadData!(nn::NearestNode, trips::GeoData)
    nTrips = length(trips)
    println("Projecting $nTrips geo-trips onto network-nodes...")
    #initializing containers
    nn.trips = trips
    nn.proj = Array{Tuple{Int,Int}}(undef, nTrips)
    for (i,t) in enumerate(trips)
        if i%100_000 == 0
            @printf("\r%.2f%% trips projected     ",100*i/nTrips)
        end
        pX, pY = toENU(t.pLon, t.pLat, nn.network)
        idP = knn(nn.tree,[pX,pY],1)[1][1]
        dX, dY = toENU(t.dLon, t.dLat, nn.network)
        idD = knn(nn.tree,[dX,dY],1)[1][1]
        nn.proj[i] = (idP,idD)
    end
    println("\r100.00% trips projected     ")
    return nn
end

"""
    `getNetworkTrips`: constructs list of NetworkTrip objects
"""
function getNetworkTrips(nn::NearestNode, tIds::AbstractArray{Int64,1})
    if length(nn.trips) < length(tIds)
        error("Trips not loaded yet")
    end
    #Store result in dictionary : (orig,dest)->NetworkTrip
    netTrips = Dict{Tuple{Int,Int},NetworkTrip}()
    for id in tIds
        t = nn.trips[id]
        od = nn.proj[id]
        #if same origin/destination
        if od[1] == od[2]
            continue
        end
        if haskey(netTrips, od)
            nt = netTrips[od]
            time = nt.time ^ (nt.weight/(nt.weight + 1.)) * t.time^(1/(nt.weight + 1.))
            netTrips[od] = NetworkTrip(nt.orig, nt.dest, time, nt.weight+1., false)
        else
            netTrips[od] = NetworkTrip((od[1], od[1], 0.), (od[2], od[2], 0.), t.time, 1., false)
        end
    end
    return collect(values(netTrips))
end

"""
    `getTripTiming`: returns a trip timing estimation from a NetworkTimings result
"""
function getTripTiming(nn::NearestNode, timings::NetworkTimings, tId::Int)
    o,d = nn.proj[tId]
    return timings.pathTimes[o,d]
end
function getTripTiming(nn::NearestNode, timings::NetworkTimings,
            pLon::Float32, pLat::Float32, dLon::Float32, dLat::Float32)
    pX, pY = toENU(pLon, pLat, nn.network)
    o = knn(nn.tree,[pX,pY],1)[1][1]
    dX, dY = toENU(dLon, dLat, nn.network)
    d = knn(nn.tree,[dX,dY],1)[1][1]
    return timings.pathTimes[o,d]
end


text:
# Test functionality with continuous features.
@testset "Evaporative Cooling ReliefF - Continuous Features" begin
    data = rand(1000, 10)
    for idx1 = 1:size(data, 2) - 1
        for idx2 = idx1+1:size(data, 2)
            target = Int64.(data[:, idx1] .> data[:, idx2])
            rank = Relief.ecrelieff(data, target, f_type="continuous")
            @test Set([rank[idx1], rank[idx2]]) == Set([1, 2])
        end
    end
end

# Test functionality with discrete features.
@testset "Evaporative Cooling ReliefF - Discrete Features" begin
    data = rand([0, 1, 2, 3], 1000, 10)
    for idx1 = 1:size(data, 2) - 1
        for idx2 = idx1+1:size(data, 2)
            target = Int64.(data[:, idx1] .> data[:, idx2])
            rank = Relief.ecrelieff(data, target, f_type="discrete")
            @test Set([rank[idx1], rank[idx2]]) == Set([1, 2])
        end
    end
end


text:
"""
	Implements all needed parts for constructing DeepSVDD model for detecting anomalies.
		(Depp One-Class Classification -> http://proceedings.mlr.press/v80/ruff18a/ruff18a.pdf)
	Code is inspired by original Pytorch implementation https://github.com/lukasruff/Deep-SVDD-PyTorch.

"""

mutable struct DeepSVDD
	encoder
	decoder
	objective
	c
	R
	nu
end

Flux.@functor DeepSVDD

function (svdd::DeepSVDD)(x)
	z = svdd.encoder(x)
	return z
end

function AE(svdd::DeepSVDD, x)
	z = svdd.encoder(x)
	x̂ = svdd.decoder(z)
	return x̂
end

mse_score(x̂,x;dims) = vec(Flux.sum((x̂ .- x).^2, dims=dims))

function anomaly_score(svdd::DeepSVDD, x)
	z = svdd(x)
	dist = mse_score(z, svdd.c, dims=1)
	#vec(Flux.sum((z .- svdd.c).^2, dims=1))
	score = (svdd.objective == "soft-boundary") ? dist .- svdd.R.^2 : dist
	return score
end

function init_c(svdd::DeepSVDD, data, eps=0.1; batchsize=50, to_testmode::Bool=true)
	(to_testmode == true) ? Flux.testmode!(svdd) : nothing
	c_ = deepcopy(svdd.c)
	for X in Flux.Data.DataLoader(data, batchsize=batchsize)
		X = getobs(X)|>gpu
		ŷ = svdd(X)
		c_ += Flux.sum(ŷ, dims=2)
	end
	c_ = c_ |> cpu
	c_ /= size(data, 4)
	c_[(abs.(c_) .< eps) .& (c_ .< 0)] .= -eps
	c_[(abs.(c_) .< eps) .& (c_ .> 0)] .= eps
	(to_testmode == true) ? Flux.testmode!(svdd, false) : nothing
	return c_
end

function get_radius(svdd::DeepSVDD, X)
	output = svdd(X)
	dist = vec(Flux.sum((output .- svdd.c).^2, dims=1))
	return Statistics.quantile(sqrt.(dist), 1-svdd.nu)
end

# LOSSES

# autoencoder loss function
ae_loss(ŷ, y; dims=(1,2,3)) = Flux.mean(Flux.sum((ŷ .- y).^2, dims=dims))

# objective / loss function for Soft-boundary "SVM"
function sb_loss(ŷ, c, R, nu)
	dist = Flux.sum((ŷ .- c).^2, dims=1)
	scores = dist .- R^2 #
	loss = R^2 + (1 / nu) * Flux.mean(max.(0, scores))
end

# objective / loss function for One-Class "SVM"
oc_loss(ŷ, c, R, nu) = Flux.mean(Flux.sum((ŷ .- c).^2, dims=1))

# anomaly score on GPU
function anomaly_score_gpu(svdd::DeepSVDD, real_input; batch_size=64, to_testmode::Bool=true)
	real_input = Flux.Data.DataLoader(real_input, batchsize=batch_size)
	(to_testmode == true) ? Flux.testmode!(svdd) : nothing
	svdd = svdd |> gpu
	scores = Array{Float32}([])
	for x in real_input
		x = x |> gpu
		score = anomaly_score(svdd, x)
		scores = cat(scores, score |> cpu, dims=1)
	end
	(to_testmode == true) ? Flux.testmode!(svdd, false) : nothing
	return scores
end

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

function conv_svdd_constructor(
	;idim=(2,2,1),
	zdim::Int=1,
	activation="relu",
	kernelsizes=(1,1),
	channels=(1,1),
	scalings=(1,1),
	init_seed=nothing,
	objective="soft-boundary",
	batchnorm=false,
	nu::AbstractFloat = 0.1f0,
	R::AbstractFloat = 0f0,
	c = nothing,
	kwargs...
)
	# if seed is given, set it
	(init_seed !== nothing) ? Random.seed!(init_seed) : nothing
	# encoder
	encoder = Chain(
		conv_encoder(idim, zdim, kernelsizes, channels, scalings; activation=activation, batchnorm=batchnorm)...
	)
	# decoder
	decoder = Chain(
		conv_decoder(idim, zdim, reverse(kernelsizes), reverse(channels), reverse(scalings),
			activation=activation, vec_output=false, vec_output_dim=nothing, batchnorm=batchnorm)...,
			x->sigmoid.(x) # our dataset images are scaled from 0-1
	)
	# reset seed
	(init_seed !== nothing) ? Random.seed!() : nothing

	c_ = (c !== nothing) ? c : zeros(Float32, zdim)

	model = DeepSVDD(encoder, decoder, objective, c_, R, nu)
end

function svdd_constructor(
	;idim::Int=1,
	zdim::Int=1,
	hdim::Int=64,
	activation = "relu",
	nlayers::Int=3,
	init_seed=nothing,
	objective="soft-boundary",
	nu::AbstractFloat = 0.1f0,
	R::AbstractFloat = 0f0,
	c = nothing,
	kwargs...
	)
	(init_seed !== nothing) ? Random.seed!(init_seed) : nothing

	encoder = build_mlp(idim, hdim, zdim, nlayers, activation=activation, lastlayer="linear")

	decoder = build_mlp(zdim, hdim, idim, nlayers, activation=activation, lastlayer="linear")

	(init_seed !== nothing) ? Random.seed!() : nothing

	c_ = (c !== nothing) ? c : zeros(Float32, zdim)

	return DeepSVDD(encoder, decoder, objective, c_, R, nu)
end

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

# fit Autoencoder
function fit_ae(svdd::DeepSVDD, optim, data, params)
	train_loader, val_loader = prepare_dataloaders(data, batch_size=params.batch_size, iters=params.ae_iters)
	history = MVHistory()
	total_iters = length(train_loader)
	progress = Progress(total_iters)
	# early stopping
	best_svdd = deepcopy(svdd)
	patience = params.patience
	val_batches = length(val_loader)
	best_val_loss = Inf

	if ndims(first(train_loader)) == 2
		dims = 1
	elseif ndims(first(train_loader)) == 4
		dims = (1,2,3)
	else
		error("unknown data type (no 2D or 4D tensor)")
	end

	svdd = svdd |> gpu
	ps_ae = Flux.params(svdd.encoder, svdd.decoder)
	for (iter, X) in enumerate(train_loader)
		X = getobs(X)|>gpu
		loss_ae, back = Flux.pullback(ps_ae) do
			ae_loss(AE(svdd, X), X, dims=dims)
		end
		grad_ae = back(1f0)
		Flux.Optimise.update!(optim, ps_ae, grad_ae)

		push!(history, :loss_ae, iter, loss_ae)

		next!(progress; showvalues=[
			(:iters, "$(iter)/$(total_iters)"),
			(:ae_loss, loss_ae)
			])

		if mod(iter, params.ae_check_every) == 0
			total_val_loss = 0
			Flux.testmode!(svdd)
			for X_val in val_loader
				X_val = X_val |> gpu
				total_val_loss += ae_loss(AE(svdd, X), X)
			end
			Flux.testmode!(svdd, false)
			push!(history, :validation_ae_loss, iter, total_val_loss/val_batches)
			if total_val_loss < best_val_loss
				best_val_loss = total_val_loss
				patience = params.patience
				best_svdd = deepcopy(svdd)
			else
				patience -= 1
				if patience == 0
					@info "Stopped training of Autoencoder after $(iter) iterations"
					break
				end
			end
		end
	end
	return best_svdd, history
end

function StatsBase.fit!(svdd::DeepSVDD, data, params)
	# prepare batches & loaders
	train_loader, val_loader = prepare_dataloaders(data, params)
	svdd_history = MVHistory()
	# prepare for early stopping
	best_val_loss = 1e10

	patience = params.patience
	val_batches = length(val_loader)

	#model to gpu
	svdd = svdd |> gpu

	# ADAMW(η = 0.001, β = (0.9, 0.999), decay = 0) = Optimiser(ADAM(η, β), WeightDecay(decay))
	opt_ae = haskey(params, :decay) ? ADAMW(params.lr_ae, (0.9, 0.999), params.decay) : ADAM(params.lr_ae)
	opt_svdd = haskey(params, :decay) ? ADAMW(params.lr_svdd, (0.9, 0.999), params.decay) : ADAM(params.lr_svdd)

	# loss_function for SVDD
	loss = (svdd.objective == "soft-boundary") ? sb_loss : oc_loss

	# Pretraining of AE
	svdd, ae_history = fit_ae(svdd, opt_ae, data, params)
	best_svdd = deepcopy(svdd)

	# init center of hypershpere
	svdd.c = init_c(svdd, data[1][1]) |> gpu

	# train "SVDD"
	progress = Progress(params.iters)
	ps_svdd = Flux.params(svdd.encoder)
	for (iter, X) in enumerate(train_loader)
		X = getobs(X)|>gpu
		loss_svdd, back = Flux.pullback(ps_svdd) do
			loss(svdd(X), svdd.c, svdd.R, svdd.nu)
		end
		grad_svdd = back(1f0)
		Flux.Optimise.update!(opt_svdd, ps_svdd, grad_svdd)

		# after warm up
		svdd.R = ((iter >= 0.05*params.iters) & (svdd.objective == "soft-boundary")) ? get_radius(svdd, X) : svdd.R

		push!(svdd_history, :loss_svdd, iter, loss_svdd)
		next!(progress; showvalues=[
			(:iters, "$(iter)/$(params.iters)"),
			(:svdd_loss, loss_svdd)
			])

		if mod(iter, params.check_every) == 0
			total_val_loss = 0
			Flux.testmode!(svdd)
			for X_val in val_loader
				X_val = X_val |> gpu
				total_val_loss += loss(svdd(X), svdd.c, svdd.R, svdd.nu)
			end
			Flux.testmode!(svdd, false)
			push!(svdd_history, :validation_svdd_loss, iter, total_val_loss/val_batches)
			if total_val_loss < best_val_loss
				best_val_loss = total_val_loss
				patience = params.patience
				best_svdd = deepcopy(svdd)
			else
				patience -= 1
				if patience == 0
					@info "Stopped training after $(iter) iterations"
					global iters = iter - params.check_every*params.patience
					break
				end
			end
		end
		if iter == params.iters
			global iters = params.iters
		end
	end

	return (svdd_history, ae_history), best_svdd, sum(map(p->length(p), Flux.params(svdd))), iters
end


text:
## context

## A context stores objects where named values are looked up.
type Context
    view ## of what? a Dict, Module, CompositeKind, DataFrame
    parent ## a context or nothing
    _cache::Dict
end
function Context(view, parent=nothing)
    d = Dict()
    d["."] = view
    Context(view, parent, d)
end

## used by renderTokens
function ctx_push(ctx::Context, view)
    Context(view, ctx) ## add context as a parent
end

## Lookup value by key in the context
function lookup(ctx::Context, key)
    if haskey(ctx._cache, key)
        value = ctx._cache[key]
    else
        context = ctx
        value = nothing
        while value == nothing && context != nothing
            ## does name have a .?
            if ismatch(r"\.", key)
                ## do something with "."
                error("Not implemented. Can use Composite Kinds in the view.")
            else
                ## strip leading, trailing whitespace in key
                value = lookup_in_view(context.view, stripWhitepace(key))
            end

            context = context.parent

        end

        ## cache
        ctx._cache[key] = value
    end

    if is(value, Function)
        value = value()
    end

    return(value)
end

## Lookup value in an object by key
## This of course varies based on the view.
## we special case dataframes here, so that we don't have to assume package is loaded
function lookup_in_view(view, key)
    if is_dataframe(view)
        if ismatch(r":", key)  key = key[2:end] end
        key = Symbol(key)
        out = nothing
        if haskey(view, key)
            out = view[1, key] ## first element only
        end
        return out
    else
        _lookup_in_view(view, key)
    end
end

function _lookup_in_view(view::Dict, key)

    ## is it a symbol?
    if ismatch(r"^:", key)
        key = Symbol(key[2:end])
    end

    out = nothing
    if haskey(view, key)
        out = view[key]
    end

    out
end

function _lookup_in_view(view::Module, key)

    hasmatch = false
    re = Regex("^$key\$")
    for i in names(view, true)
        if ismatch(re, string(i))
            hasmatch = true
            break
        end
    end

    out = nothing
    if hasmatch
        out = getfield(view, Symbol(key))  ## view.key
    end
    out

end

## Default is likely not great, but we use CompositeKind
function _lookup_in_view(view, key)
    nms = fieldnames(view)
    re = Regex(key)
    has_match = false
    for i in nms
        if ismatch(Regex(key), string(i))
            has_match=true
            break
        end
    end

    out = nothing
    if has_match
        out = getfield(view, Symbol(key))  ## view.key
    end

    out
end


text:
#=
# 4. The Kutta condition

For surfaces that contain convex edges, the vortex sheet strength assumes a singular behavior in the vicinity of these edges, with a strength that depends on the interior angle of the edge: sharper edges have more singular behavior. In the discrete representation of the surface, edges are only approximately represented by the sudden disruptions of positions in clusters of adjacent points. The behavior in this discrete form is not quite singular, but the solution for the vortex sheet strength $\mathfrak{f}$ nonetheless exhibits a large and rapid change of amplitude.
=#

# Let's illustrate this with the example of a flat plate. First create the grid.
#md # ```@setup 4.-The-Kutta-condition
#md # using GridPotentialFlow
#md # using Plots
#md # ```
#!md using GridPotentialFlow
#!md using Plots
Δx = 0.01
Lx = 2.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx);

# To create a potential flow model with flat plate, we use the `Plate` constructor from `RigidBodyTools.jl` and transform the body such that its angle of attack is `α`.
c = Lx/2 # chord length
α = -π/6 # angle of attack
plate = Plate(c,4*cellsize(g))
Tr = RigidTransform((0.0,0.0),α)
Tr(plate)
pfb = PotentialFlowBody(plate);
model = VortexModel(g,bodies=[pfb]);

# To simulate the plate in a uniform flow, we have to set the `U∞` field of the vortex model to our desired value.
setU∞(model,(1.0,0.0))
sol = solve(model);

# By plotting a large number streamlines, we see how the flow accelerates as it tries to make its way around the sharp edges of the flat plate. In the continuous case, the velocity field would be singular at the edges.
plot(sol.ψ,g,levels=60);
plot!(plate,linecolor=:black,linewidth=2,xlabel="x",ylabel="y")

# The continous vortex sheet strength shows similar singularities at the edges, and also in our discrete case, vortex sheet strength the indeed assumes an approximate singular behavior near those edges.
plot(sol.f,xlabel="body point index",ylabel="f",legend=false)

#=
If we seek to eliminate this behavior, we must first have some means of exposing it. In fact, for any discretized surface, the essence of this nearly-singular behavior lies in the previously defined vector $\mathfrak{f}_0$, and all other bound vortex sheets associated with the same surface share the same nearly-singular behavior. Thus, we will use a multiplicative decomposition of the vortex sheet strength:

$\mathfrak{f} = \mathfrak{f}_0 \circ \tilde{\mathfrak{f}}$

where $\circ$ is the Hadamard product. This decomposed form isolates the singular behavior into $\mathfrak{f}_0$, and $\tilde{\mathfrak{f}}$ is a relatively smoother vector of surface point data.

=#

# For our flat plate, we create $\mathfrak{f}_0$ in the same way as we did on the previous page.
ones = ScalarData(length(plate))
ones .= 1.0
f₀ = model.system.ibp.Sfact\ones
plot(f₀)

# Then we can plot the smooth $\tilde{\mathfrak{f}}$ as follows.
plot(sol.f./f₀,xlabel="body point index",ylabel="f̃",legend=false)

#=
## Using the Kutta condition in a steady-state problem

The Kutta condition corresponds to annihilating the nearly-singular behavior at a surface point. At such points, we will set the corresponding value of $\tilde{\mathfrak{f}}$ to zero. Suppose we wish to enforce the Kutta condition at an edge corresponding to surface point $k$. The condition is

$\mathfrak{e}_{k}^{T} \tilde{\mathfrak{f}} = 0.$

In the regularization operation on $\mathfrak{f}$, we can absorb $\mathfrak{f}_0$ into $\mathsf{R}$, first noting that the Hadamard product can alternatively be written with the help of a diagonal matrix,

$\mathfrak{f} = \mathfrak{f}_0 \circ \tilde{\mathfrak{f}} = \mathsf{D}{\mathfrak{f}_0} \tilde{\mathfrak{f}}.$

Then, we can define a re-scaled regularization operator,

$\mathsf{R}\mathfrak{f} = \mathsf{R} \mathsf{D}{\mathfrak{f}_0} \tilde{\mathfrak{f}} = \tilde{\mathsf{R}} \tilde{\mathfrak{f}}.$

We will first take the steady-state approach to enforce the Kutta condition: allow the bound circulation to be set appropriately, with the implicit understanding that there is a starting vortex of equal and opposite circulation at infinity that preserves the Kelvin circulation theorem. The Lagrange multiplier for this constraint will not be $\Gamma_b$, but $s_0$.

$\begin{bmatrix}
\mathsf{L} & \tilde{\mathsf{R}} & 0\\
\mathsf{E} &  0 & \mathfrak{1}\\
0 & \mathfrak{e}_{k}^{T} & 0\\
\end{bmatrix} \begin{pmatrix} \mathsf{s} \\ \tilde{\mathfrak{f}} \\ s_0 \end{pmatrix} =
\begin{pmatrix} -\mathsf{w} \\ \mathfrak{s}'_b \\ 0 \end{pmatrix}.$

=#

# In `GridPotentialFlow.jl`, we create a regularized potential flow model by specifying the index of the regularized edges using the `edges` keyword when creating a `PotentialFlowBody`. In this case, we want to regularize the trailing edge, which is the last point of `Plate`, so its index is `length(plate)`.
pfb = PotentialFlowBody(plate,edges=[length(plate)])
model = VortexModel(g,bodies=[pfb],U∞=(1.0,0.0))
sol = solve(model);

# When there are no vortices present in the model, the package solves the above saddle point system and the solution is the steady state solution. Inspection of the streamlines shows that flow now indeed leaves the trailing edge smoothly.
plot(sol.ψ,g);
plot!(plate,linecolor=:black,linewidth=2,xlabel="x",ylabel="y")

# By plotting the vortex sheet strength, we can see that by constraining the trailing edge point of $\tilde{\mathfrak{f}}$, its whole distribution is shifted upward compared to the unregularized case such that the last value equals zero.
plot(plot(sol.f,xlabel="body point index",ylabel="f",legend=false),plot(sol.f./f₀,xlabel="body point index",ylabel="f̃",legend=false),size=[800,300])

#=
## Using the Kutta condition to set a new vortex element

In the previous section, we used the Kutta condition to set the bound circulation but did not explicitly create a new vortex element. This vortex element was assumed to lie at infinity so that its effect was negligible except insofar as it left equal but opposite circulation about the body.

In this section, we will create a new vortex element in the vicinity of the edge at which we are applying the Kutta condition. We will thus seek to establish the strength of this new element and to do so in such a manner that the overall circulation of the flow is conserved. Once the element is created, it will be allowed to advect with the local fluid velocity (minus its own contribution to this velocity).

Let us assume that the new vortex element (which we label with the subscript 1) is introduced at some point in physical space, and that its immersion into the Cartesian grid is described by a grid vector $\mathsf{d}_1$ and that its strength (i.e., its circulation) is $\delta\Gamma_1$. Thus, the fluid vorticity after this new element's introduction can be written as

$\mathsf{w} + \delta\Gamma_1 \mathsf{d}_1.$

The vector $\mathsf{d}_1$ is the discrete Dirac delta function, identical (or similar) to the function that constitutes the regularization $\mathsf{R}$ and interpolation $\mathsf{E}$ matrices.

The Kutta condition is still to be enforced. We also seek to ensure that the total circulation is zero. (We are assuming that the flow has started from rest.) Let us denote the circulation of the existing fluid vorticity $\mathsf{w}$ by

$\Gamma_{\mathsf{w}} = \Delta x \, \mathsf{1}^{T}\mathsf{w},$

where $\mathsf{1}$ is a grid vector of ones. Then, the circulation constraint is

$\mathfrak{1}^{T} \mathfrak{f} + \delta\Gamma_1 + \Gamma_{\mathsf{w}} = 0.$

The circulation of the bound vortex sheet $\mathfrak{f}$ can be re-written in terms of the smooth part of the sheet as $\mathfrak{1}^{T} \mathfrak{f} = \mathfrak{f}_0 \tilde{\mathfrak{f}}$.

With these two constraints, the overall saddle point system of equations is

$\begin{bmatrix}
\mathsf{L} & \tilde{\mathsf{R}} & 0 & \mathsf{d}_1\\
\mathsf{E} &  0 & \mathfrak{1} & 0\\
0 & \mathfrak{e}_{k}^{T} & 0 & 0\\
0 & \mathfrak{f}_{0}^{T} & 0 & 1
\end{bmatrix} \begin{pmatrix} \mathsf{s} \\ \tilde{\mathfrak{f}} \\ s_0 \\ \delta\Gamma_1 \end{pmatrix} =
\begin{pmatrix} -\mathsf{w} \\ \mathfrak{s}'_b \\ 0 \\ -\Gamma_{\mathsf{w}} \end{pmatrix}.$

=#

# To simulate this case for our flat plate, we need to insert a point vortex in the flow, preferably near the trailing edge. This situation corresponds to the flow right after impulsively starting a uniform flow around a flat plate and the point vortex now represents the starting vortex. We position the vortex behind and orthogonal the plate at a distance that is proportional to the convective length from the trailing edge. This position does not matter much here, but is of importance when we step in time and insert more point vortices every time step. Its strength can be set arbitrarily, because any vortices that are used for regularizing edges will be ignored when calculating $\mathsf{w}$ in the right-hand side for the saddle-point system.
Δt = 1e-2
vTE = Vortex(plate.x[end]+3*Δt*cos(α+π/2),plate.y[end]+3*Δt*sin(α+π/2),0.0);

# We include the point vortex when we create our model. Because the package detects that there is a vortex in the model, it solves the unsteady saddle point system and the solution will include the strength for the new point vortex.
pfb = PotentialFlowBody(plate,edges=[length(plate)])
model = VortexModel(g,bodies=[pfb],vortices=[vTE],U∞=(1.0,0.0))
sol = solve(model);
plot(sol.ψ,g);
plot!(plate,linecolor=:black,linewidth=2)
scatter!(model.vortices.x,model.vortices.y,color=:black,markersize=2,xlabel="x",ylabel="y")
# Note that the strength of the new point vortex is not automatically set when calling `solve` on the vortex model

# Because of the proximity of the point vortex to the flat plate, $\tilde{\mathfrak{f}}$ exhibits a quick variation in its value at the surface points that lie closest to the point vortex. The value at the trailing edge point itself is still constrained to zero.
plot(plot(sol.f,xlabel="body point index",ylabel="f",legend=false),plot(sol.f./f₀,xlabel="body point index",ylabel="f̃",legend=false),size=[800,300])

#=
## Applying more than one Kutta condition on a body

Suppose we wish to enforce the Kutta condition at two edges of the body---at points $k_{1}$ and $k_{2}$---instead of one. Each such point has a constraint,

$\mathfrak{e}_{k_j}^{T} \tilde{\mathfrak{f}} = 0,\quad j = 1,2.$

For two such constraints, we need two Lagrange multipliers: the strengths of two new vortices, $\delta\Gamma_1$ and $\delta\Gamma_2$, immersed into the grid with $\mathsf{d}_1$ and $\mathsf{d}_2$, respectively; and we still need the Lagrange multiplier $s_0$ to ensure that Kelvin's circulation theorem is also enforced. The system in the previous section is thus easily generalized to the following:

$\begin{bmatrix}
\mathsf{L} & \tilde{\mathsf{R}} & 0 & \mathsf{d}_1 & \mathsf{d}_2\\
\mathsf{E} &  0 & \mathfrak{1} & 0& 0\\
0 & \mathfrak{e}_{k_1}^{T} & 0 & 0 & 0\\
0 & \mathfrak{e}_{k_2}^{T} & 0 & 0 & 0\\
0 & \mathfrak{f}_{0}^{T} & 0 & 1 & 1
\end{bmatrix} \begin{pmatrix} \mathsf{s} \\ \tilde{\mathfrak{f}} \\ s_0 \\ \delta\Gamma_1 \\ \delta\Gamma_2 \end{pmatrix} =
\begin{pmatrix} -\mathsf{w} \\ \mathfrak{s}'_b \\ 0 \\ 0 \\ -\Gamma_{\mathsf{w}} \end{pmatrix}.$
=#

# We now simulate this case again for our flat plate problem. This solution corresponds to the flow right after impulsively starting a uniform flow around a flat plate, but unlike the previous case, the flow now also separates at the leading edge.

# We create two vortices, one near the leading edge and one near the trailing edge.
vLE = Vortex(plate.x[1]+3*Δt*plate.len*cos(plate.α+π/2),plate.y[1]+3*Δt*plate.len*sin(plate.α+π/2),0.0);
vTE = Vortex(plate.x[end]+3*Δt*cos(α+π/2),plate.y[end]+3*Δt*sin(α+π/2),0.0);

# The model now has two regularized edges, corresponding to the first body point and the last one.
pfb = PotentialFlowBody(plate,edges=[1,length(plate)])
model = VortexModel(g,bodies=[pfb],vortices=[vLE,vTE],U∞=(1.0,0.0))
sol = solve(model);
plot(sol.ψ,g);
plot!(plate,linecolor=:black,linewidth=2)
scatter!(model.vortices.x,model.vortices.y,markersize=2,xlabel="x",ylabel="y")

#

plot(plot(sol.f,xlabel="body point index",ylabel="f",legend=false),plot(sol.f./f₀,xlabel="body point index",ylabel="f̃",legend=false),size=[800,300])


text:
import LinearAlgebra
using LinearAlgebra
using SparseArrays
using Random
import StatsBase
using StatsBase

using jInv
export prepareSyntheticPointCloudData
export subv2ind,ind2subv
function subv2ind(shape, cindices)
    """Return linear indices given vector of cartesian indices.
    shape: d-tuple of dimensions.
    cindices: n-iterable of d-tuples, each containing cartesian indices.
    Returns: n-vector of linear indices.

    Based on:
    https://discourse.julialang.org/t/psa-replacement-of-ind2sub-sub2ind-in-julia-0-7/14666/8
    Similar to this matlab function:
    https://github.com/tminka/lightspeed/blob/master/subv2ind.m
    """
    lndx = LinearIndices(Dims(shape))
    n = length(cindices)
    out = Array{Int}(undef, n)
    for i = 1:n
        out[i] = lndx[cindices[i]...]
    end
    return out
end

function ind2subv(shape, indices)
    """Map linear indices to cartesian.
    shape: d-tuple with size of each dimension.
    indices: n-iterable with linear indices.
    Returns: n-vector of d-tuples with cartesian indices.

    Based on:
    https://discourse.julialang.org/t/psa-replacement-of-ind2sub-sub2ind-in-julia-0-7/14666/8
    Similar to this matlab function:
    https://github.com/probml/pmtk3/blob/master/matlabTools/util/ind2subv.m
    """
    n = length(indices)
    d = length(shape)
    cndx = CartesianIndices(Dims(shape))
    out = Array{Tuple}(undef, n)
    for i=1:n
        lndx = indices[i]
        out[i] = cndx[lndx]
    end
    return out
end

function loc2cs3D(loc1::Union{Int64,Array{Int64}},loc2::Union{Int64,Array{Int64}},loc3::Union{Int64,Array{Int64}},n::Array{Int64,1})
@inbounds cs = loc1 .+ (loc2.-1)*n[1] .+ (loc3.-1)*n[1]*n[2];
return cs;
end

function prepareSyntheticPointCloudData(P,Mesh::RegularMesh,npc::Int64,theta_phi_PC::Array{Float64,2} ,b::Array{Float64,2}, noiseTrans :: Float64, noiseAngle :: Float64 ,filename::String)
eps = 0.01171875*2.5 #0.0234375;
mid = (Mesh.domain[1:2:end] + Mesh.domain[2:2:end]) ./ 2.0;
#Add noise:
b = b + noiseTrans*ones(size(b,1),3)
println("noiseangle:",noiseAngle);
theta_phi_PC = theta_phi_PC +  noiseAngle*randn(size(theta_phi_PC));
b[1,:] .= 0.0;
theta_phi_PC[1,:] .= 0.0;
println("theta phi after noise added:",theta_phi_PC);

Parray = Array{Array{Float64}}(undef,npc);
Normals = Array{Array{Float64}}(undef,npc);
global d = [];
for i=1:npc
	s = 0; e = 1;
	sz = round(Int64,size(P,1)/2)+1;
	ws  = Weights( vec([collect(LinRange(s,e,sz)) ; zeros(size(P,1)-sz,1)  ] ))
	if(i == 2)
		ws = Weights(reverse(vec([collect(LinRange(s,e,sz)) ; zeros(size(P,1)-sz,1)  ] )));
	end
	sampledPointIndices = sample(collect(1:size(P,1)), ws, round(Int64,size(P,1)/2),replace=false )
	# println("sampledPointIndices:",sampledPointIndices)
	P_curr = P[sampledPointIndices,:];
	Normals[i] = P_curr[:,4:6];
	cursubs = P_curr[:,1:3];
	writedlm(string("PC_",filename,"_",i,".txt"),cursubs);
	#println("size of cursubs:",size(cursubs))
	subsForward = cursubs + (eps .*Normals[i]);
	subsBackward = cursubs - (eps .*Normals[i]);
	curr_points = [subsBackward; cursubs; subsForward];
	R = getRotate3D(theta_phi_PC[i,1],theta_phi_PC[i,2]);
	curr_points = curr_points .- (mid')
	curr_points =  ( curr_points )*(R') .+ (mid')
	curr_points .+= b[i,:]';
	writedlm(string("PC_rotated_",filename,"_",i,".txt"),curr_points[round(Int64,size(curr_points,1)/3)+1:round(Int64,2*size(curr_points,1)/3),:]);

	println("size of P:",size(curr_points))
	Parray[i] = curr_points;
	vals = [ones(size(subsBackward,1));0.5.*ones(size(cursubs,1)); zeros(size(subsForward,1))];
	#vals = [ones(size(cursubsbwd,1)); zeros(size(cursubsfwd,1))];
	d = [d ;vals];
end

d = Array{Float64,1}(d);
#println("size of d:",size(d))

file = matopen(string(filename,"_data.mat"), "w");
write(file,"Data",d);
write(file,"P",Parray);
write(file,"Normals",Normals);
write(file,"b",b);
write(file,"theta_phi_PC",theta_phi_PC)
write(file,"npc",npc)
close(file);
return;
end

export readPointCloudDataFile
function readPointCloudDataFile(filename::String)
file = matopen(string(filename,"_data.mat"), "r");
Data = read(file,"Data");
P = read(file,"P");
Normals = read(file,"Normals");
b = read(file,"b");
theta_phi_PC = read(file,"theta_phi_PC");
npc = read(file,"npc")
close(file);
return Data,P,Normals,b,theta_phi_PC,npc;
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
using ForecastModels
using Test

@testset "ForecastModels.jl" begin
    # Write your tests here.
end


text:
using Base.Test, PyCall, Compat

roundtrip(T, x) = convert(T, PyObject(x))
roundtrip(x) = roundtrip(PyAny, x)
roundtripeq(T, x) = roundtrip(T, x) == x
roundtripeq(x) = roundtrip(x) == x

# test handling of type-tuple changes in Julia 0.4
import PyCall.pyany_toany
@test pyany_toany(Int) == Int
@test pyany_toany(PyAny) == Any
@test pyany_toany(@compat Tuple{Int,PyAny}) == @compat Tuple{Int,Any}
@test pyany_toany(@compat Tuple{Int,Tuple{PyAny,Int8}}) == @compat Tuple{Int,Tuple{Any,Int8}}
if VERSION >= v"0.4.0-dev+4319"
    @test pyany_toany(@compat Tuple{PyAny,Int,Vararg{PyAny}}) == @compat Tuple{Any,Int,Vararg{Any}}
end

@test roundtripeq(17)
@test roundtripeq(0x39)
@test roundtripeq(true) && roundtripeq(false)
@test roundtripeq(3.14159)
@test roundtripeq(1.3+4.5im)
@test roundtripeq(nothing)
@test roundtripeq("Hello world")
@test roundtripeq("Hëllö")
@test roundtripeq(Symbol, :Hello)
@test roundtripeq(C_NULL) && roundtripeq(convert(Ptr{Void}, 12345))
@test roundtripeq([1,3,4,5]) && roundtripeq([1,3.2,"hello",true])
@test roundtripeq([1 2 3;4 5 6]) && roundtripeq([1. 2 3;4 5 6])
@test roundtripeq((1,(3.2,"hello"),true)) && roundtripeq(())
@test roundtripeq(Int32)
@test roundtripeq(@compat Dict(1 => "hello", 2 => "goodbye")) && roundtripeq(Dict())
@test roundtripeq(Uint8[1,3,4,5])

@test pycall(PyObject(x -> x + 1), PyAny, 314158) == 314159
@test roundtrip(x -> x + 1)(314158) == 314159

testkw(x; y=0) = x + 2*y
@test pycall(PyObject(testkw), PyAny, 314157) == 314157
@test pycall(PyObject(testkw), PyAny, 314157, y=1) == 314159
@test roundtrip(testkw)(314157) == 314157
@test roundtrip(testkw)(314157, y=1) == 314159

if PyCall.npy_initialized
    @test PyArray(PyObject([1. 2 3;4 5 6])) == [1. 2 3;4 5 6]
end
@test PyVector(PyObject([1,3.2,"hello",true])) == [1,3.2,"hello",true]
@test PyDict(PyObject(@compat Dict(1 => "hello", 2 => "goodbye"))) == @compat Dict(1 => "hello", 2 => "goodbye")

let d = PyDict(@compat Dict(1 => "hello", "yes" => 34))
    @test get(d.o, 1) == "hello"
    set!(d.o, "yes", "goodbye")
    @test d["yes"] == "goodbye"
end

@test roundtripeq(Any[1 2 3; 4 5 6])
@test roundtripeq([])
@test convert(Array{PyAny,1}, PyObject(Any[1 2 3; 4 5 6])) == Any[Any[1,2,3],Any[4,5,6]]
if PyCall.npy_initialized
    @test roundtripeq(begin A = Array(Int); A[1] = 3; A; end)
end
@test convert(PyAny, PyObject(begin A = Array(Any); A[1] = 3; A; end)) == 3

array2py2arrayeq(x) = PyCall.py2array(Float64,PyCall.array2py(x)) == x
@test array2py2arrayeq(rand(3))
@test array2py2arrayeq(rand(3,4))
@test array2py2arrayeq(rand(3,4,5))

@test roundtripeq(2:10) && roundtripeq(10:-1:2)
@test roundtrip(2:2.0:10) == convert(Vector{Float64}, 2:2.0:10)

@pyimport math
@test_approx_eq math.sin(3) sin(3)

@test collect(PyObject([1,"hello",5])) == [1,"hello",5]

@test try @eval (@pyimport os.path) catch ex isa(ex, ArgumentError) end

@test PyObject("hello") == PyObject("hello")
@test PyObject("hello") != PyObject("hellö")
@test PyObject(hash) == PyObject(hash)
@test PyObject(hash) != PyObject(println)
@test hash(PyObject("hello")) == hash(PyObject("hello"))
@test hash(PyObject("hello")) != hash(PyObject("hellö"))
@test hash(PyObject("hello")) != hash("hellö")
@test hash(PyObject(hash)) == hash(PyObject(hash))
@test hash(PyObject(hash)) != hash(PyObject(println))
@test hash(PyObject(hash)) != hash(hash)

# issue #92:
let x = PyVector(PyAny[])
    pyeval("lambda x: x.append(\"bar\")")(x)
    @test x == ["bar"]
end

if pyversion >= v"2.7" && isdefined(PyCall, :PyDateTime_CAPI)
    # Dates is built-in in Julia 0.4
    if !isdefined(Base, :Dates)
        import Dates
    end
    @test roundtripeq(Dates.Date(2012,3,4))
    @test roundtripeq(Dates.DateTime(2012,3,4, 7,8,9,11))
    @test roundtripeq(Dates.Millisecond(typemax(Int32)))
    @test roundtripeq(Dates.Millisecond(typemin(Int32)))
    @test roundtripeq(Dates.Second, Dates.Second(typemax(Int32)))
    @test roundtripeq(Dates.Second, Dates.Second(typemin(Int32)))
    @test roundtripeq(Dates.Day, Dates.Day(999999999)) # max allowed day timedelta
    @test roundtripeq(Dates.Day, Dates.Day(-999999999)) # min allowed day timedelta
end

# fixme: is there any nontrivial mimewritable test we can do?
@test !mimewritable("text/html", PyObject(1))

# in Python 3, we need a specific encoding to write strings or bufferize them
# (http://stackoverflow.com/questions/5471158/typeerror-str-does-not-support-the-buffer-interface)
pyutf8(s::PyObject) = pycall(s["encode"], PyObject, "utf-8")
pyutf8(s::ByteString) = pyutf8(PyObject(s))

# IO (issue #107)
@test roundtripeq(STDOUT)
let buf = IOBuffer(false, true), obuf = PyObject(buf)
    @test !obuf[:isatty]()
    @test obuf[:writable]()
    @test !obuf[:readable]()
    @test obuf[:seekable]()
    obuf[:write](pyutf8("hello"))
    @test position(buf) == obuf[:tell]() == 5
    let p = obuf[:seek](-2, 1)
        @test p == position(buf) == 3
    end
    let p = obuf[:seek](0, 0)
        @test p == position(buf) == 0
    end
    @test takebuf_string(buf) == "hello"
    obuf[:writelines](["first\n", "second\n", "third"])
    @test takebuf_string(buf) == "first\nsecond\nthird"
    obuf[:write](convert(Vector{Uint8}, "möre stuff"))
    @test takebuf_string(buf) == "möre stuff"
    @test isopen(buf) == !obuf[:closed] == true
    obuf[:close]()
    @test isopen(buf) == !obuf[:closed] == false
end
let buf = IOBuffer("hello\nagain"), obuf = PyObject(buf)
    @test !obuf[:writable]()
    @test obuf[:readable]()
    @test obuf[:readlines]() == ["hello\n","again"]
end
let buf = IOBuffer("hello\nagain"), obuf = PyObject(buf)
    @test obuf[:readall]() == convert(Vector{Uint8}, "hello\nagain")
end
let buf = IOBuffer("hello\nagain"), obuf = PyTextIO(buf)
    @test obuf[:encoding] == "UTF-8"
    @test obuf[:readall]() == "hello\nagain"
end
let nm = tempname()
    open(nm, "w") do f
        @test roundtripeq(f)
        pf = PyObject(f)
        @test pf[:fileno]() == fd(f)
        @test pf[:writable]()
        @test !pf[:readable]()
        pf[:write](pyutf8(nm))
        pf[:flush]()
    end
    @test readall(nm) == nm
end

# issue #112
@test roundtripeq(Array, [1,2,3,4])
@test roundtripeq(Array{Int8}, [1,2,3,4])

# buffers
let b = PyCall.PyBuffer(pyutf8("test string"))
    @test ndims(b) == 1
    @test (length(b),) == (length("test string"),) == (size(b, 1),) == size(b)
    @test stride(b, 1) == 1
    @test PyCall.iscontiguous(b) == true
end


text:
#!/usr/bin/env julia6
using ReadGlobal

function coef(out,f1,f2,f3,m1,m2,m3,α)
    Threads.@threads for i in linearindices(out)
        @inbounds out[i] = α * muladd(f1[i],m1[i], muladd(f2[i],m2[i],f3[i]*m3[i])) / muladd(m1[i],m1[i], muladd(m2[i],m2[i],m3[i]*m3[i]))
    end
end

function main()
    f1 = readfield("flux1")
    f2 = readfield("flux2")
    f3 = readfield("flux3")
    N,F = getnfilter()
    m1 = readfield("$(F)drhodx_N$N")
    m2 = readfield("$(F)drhody_N$N")
    m3 = readfield("$(F)drhodz_N$N")

    nx,ny,nz,lx,ly,lz = getdimsize()
    Δ = parse(Float64,N)*2π*lx/nx
    S,_,_,_ = read_info("normD.info")

    coef(m1,f1,f2,f3,m1,m2,m3,1/(Δ*Δ*S))
    write("coef_1modelflux",m1)

    read!("2modelflux1",m1)
    read!("2modelflux2",m2)
    read!("2modelflux3",m3)

    coef(m1,f1,f2,f3,m1,m2,m3,S)
    write("coef_2modelflux",m1)

    read!("3modelflux1",m1)
    read!("3modelflux2",m2)
    read!("3modelflux3",m3)

    coef(m1,f1,f2,f3,m1,m2,m3,1/(Δ*Δ))
    write("coef_3modelflux",m1)

    read!("4modelflux1",m1)
    read!("4modelflux2",m2)
    read!("4modelflux3",m3)

    coef(m1,f1,f2,f3,m1,m2,m3,1/(Δ*Δ))
    write("coef_4modelflux",m1)

    read!("5modelflux1",m1)
    read!("5modelflux2",m2)
    read!("5modelflux3",m3)

    coef(m1,f1,f2,f3,m1,m2,m3,1/(Δ*Δ))
    write("coef_5modelflux",m1)

end

main()


text:
using Test
using NLsolve
using Distributions
# using Statistics

using IncrementalInference

@testset "Ensure lambda's work with NLsolve" begin

function f!(F, x, dummy)
    F[1] = (x[1]+3)*(x[2]^3-7)+18
    F[2] = sin(x[2]*exp(x[1])-1)
end

data = randn(2,3)
r = nlsolve( (res, x) -> f!(res, x, data), zeros(2), inplace=true )

@test norm(r.zero - [0, 1.0]) < 1e-10

end

# abstract functional defintion
mutable struct OneDimensionTest{T} <: FunctorPairwise
  Dx::T
  OneDimensionTest{T}(a::T) where T = new(a)
end
OneDimensionTest(a::T) where T = OneDimensionTest{T}(a)

# standardized parameter list used by IIF
function (Dp::OneDimensionTest)(res::Array{Float64},
                                userdata::FactorMetadata,
                                idx::Int,
                                meas::Tuple{Array{Float64,2}},
                                p1::Array{Float64,2},
                                p2::Array{Float64,2} )
  #
  res[1] = meas[1][1,idx] - (p2[1,idx] - p1[1,idx])
  nothing
end

@testset "minimalistic pattern on how NLsolve is used in IIF." begin

N = 10
global p1 = rand(1,N)
global p2 = rand(1,N)
global VARS = Array{Array{Float64,2},1}()
push!(VARS,p1)
push!(VARS,p2)

global odo = OneDimensionTest(Normal(10.0,1.0))
global meas = (reshape(rand(odo.Dx,N),1,N),)
global ccw = CommonConvWrapper(odo, VARS[1], 1, VARS, measurement=meas, varidx=2)

# Find definition of (::ccw)(res, x) here:
nlsolve(  ccw, [0.0], inplace=true )

end

#


text:
export LowStorageRungeKutta3N
export LS3NRK44Classic, LS3NRK33Heuns

"""
    LowStorageRungeKutta3N(f, RKA, RKB, RKC, RKW, Q; dt, t0 = 0)

This is a time stepping object for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

The constructor builds a low-storage Runge--Kutta scheme using 3N storage
based on the provided `RKA`, `RKB` and `RKC` coefficient arrays.
 `RKC` (vector of length the number of stages `ns`) set nodal points position;
 `RKA` and `RKB` (size: ns x 2) set weight for tendency and stage-state;
 `RKW` (unused) provides RK weight (last row in Butcher's tableau).

The 3-N storage formulation from Fyfe (1966) is applicable to any 4-stage,
fourth-order RK scheme. It is implemented here as:

```math
\\hspace{-20mm} for ~~ j ~~ in ~ [1:ns]: \\hspace{10mm}
  t_j = t^n + \\Delta t ~ rkC_j
```
```math
  dQ_j = dQ^*_j + f(Q_j,t_j)
```
```math
  Q_{j+1} = Q_{j} + \\Delta t \\{ rkB_{j,1} ~ dQ_j + rkB_{j,2} ~ dR_j \\}
```
```math
  dR_{j+1} = dR_j + rkA_{j+1,2} ~ dQ_j
```
```math
  dQ^*_{j+1} = rkA_{j+1,1} ~ dQ_j
```

The available concrete implementations are:

  - [`LS3NRK44Classic`](@ref)
  - [`LS3NRK33Heuns`](@ref)

### References

    @article{Fyfe1966,
       title = {Economical Evaluation of Runge-Kutta Formulae},
       author = {Fyfe, David J.},
       journal = {Mathematics of Computation},
       volume = {20},
       pages = {392--398},
       year = {1966}
    }
"""
mutable struct LowStorageRungeKutta3N{T, RT, AT, Nstages} <: AbstractODESolver
    "time step"
    dt::RT
    "time"
    t::RT
    "elapsed time steps"
    steps::Int
    "rhs function"
    rhs!
    "Storage for RHS during the `LowStorageRungeKutta3N` update"
    dQ::AT
    "Secondary Storage for RHS during the `LowStorageRungeKutta3N` update"
    dR::AT
    "low storage RK coefficient array A (rhs scaling)"
    RKA::Array{RT, 2}
    "low storage RK coefficient array B (rhs add in scaling)"
    RKB::Array{RT, 2}
    "low storage RK coefficient vector C (time scaling)"
    RKC::Array{RT, 1}
    "RK weight coefficient vector W (last row in Butcher's tableau)"
    RKW::Array{RT, 1}

    function LowStorageRungeKutta3N(
        rhs!,
        RKA,
        RKB,
        RKC,
        RKW,
        Q::AT;
        dt = 0,
        t0 = 0,
    ) where {AT <: AbstractArray}

        T = eltype(Q)
        RT = real(T)

        dQ = similar(Q)
        dR = similar(Q)
        fill!(dQ, 0)
        fill!(dR, 0)

        new{T, RT, AT, length(RKC)}(
            RT(dt),
            RT(t0),
            0,
            rhs!,
            dQ,
            dR,
            RKA,
            RKB,
            RKC,
            RKW,
        )
    end
end

"""
    dostep!(Q, lsrk3n::LowStorageRungeKutta3N, p, time::Real,
            [slow_δ, slow_rv_dQ, slow_scaling])

Use the 3N low storage Runge--Kutta method `lsrk3n` to step `Q` forward in time
from the current time `time` to final time `time + getdt(lsrk3n)`.

If the optional parameter `slow_δ !== nothing` then `slow_rv_dQ * slow_δ` is
added as an additional ODE right-hand side source. If the optional parameter
`slow_scaling !== nothing` then after the final stage update the scaling
`slow_rv_dQ *= slow_scaling` is performed.
"""
function dostep!(
    Q,
    lsrk3n::LowStorageRungeKutta3N,
    p,
    time,
    slow_δ = nothing,
    slow_rv_dQ = nothing,
    in_slow_scaling = nothing,
)
    dt = lsrk3n.dt

    RKA, RKB, RKC = lsrk3n.RKA, lsrk3n.RKB, lsrk3n.RKC
    rhs!, dQ, dR = lsrk3n.rhs!, lsrk3n.dQ, lsrk3n.dR

    rv_Q = realview(Q)
    rv_dQ = realview(dQ)
    rv_dR = realview(dR)
    groupsize = 256

    rv_dR .= -0
    for s in 1:length(RKC)
        rhs!(dQ, Q, p, time + RKC[s] * dt, increment = true)

        slow_scaling = nothing
        if s == length(RKC)
            slow_scaling = in_slow_scaling
        end
        # update solution and scale RHS
        event = Event(array_device(Q))
        event = update!(array_device(Q), groupsize)(
            rv_dQ,
            rv_dR,
            rv_Q,
            RKA[s % length(RKC) + 1, 1],
            RKA[s % length(RKC) + 1, 2],
            RKB[s, 1],
            RKB[s, 2],
            dt,
            slow_δ,
            slow_rv_dQ,
            slow_scaling;
            ndrange = length(rv_Q),
            dependencies = (event,),
        )
        wait(array_device(Q), event)
    end
end

@kernel function update!(
    dQ,
    dR,
    Q,
    rka1,
    rka2,
    rkb1,
    rkb2,
    dt,
    slow_δ,
    slow_dQ,
    slow_scaling,
)
    i = @index(Global, Linear)
    @inbounds begin
        if slow_δ !== nothing
            dQ[i] += slow_δ * slow_dQ[i]
        end
        Q[i] += rkb1 * dt * dQ[i] + rkb2 * dt * dR[i]
        dR[i] += rka2 * dQ[i]
        dQ[i] *= rka1
        if slow_scaling !== nothing
            slow_dQ[i] *= slow_scaling
        end
    end
end

"""
    LS3NRK44Classic(f, Q; dt, t0 = 0)

This function returns a [`LowStorageRungeKutta3N`](@ref) time stepping object
for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This uses the classic 4-stage, fourth-order Runge--Kutta scheme
in the low-storage implementation of Blum (1962)

### References
    @article {Blum1962,
       title = {A Modification of the Runge-Kutta Fourth-Order Method}
       author = {Blum, E. K.},
       journal = {Mathematics of Computation},
       volume = {16},
       pages = {176-187},
       year = {1962}
    }
"""
function LS3NRK44Classic(F, Q::AT; dt = 0, t0 = 0) where {AT <: AbstractArray}
    T = eltype(Q)
    RT = real(T)

    RKA = [
        RT(0) RT(0)
        RT(0) RT(1)
        RT(-1 // 2) RT(0)
        RT(2) RT(-6)
    ]

    RKB = [
        RT(1 // 2) RT(0)
        RT(1 // 2) RT(-1 // 2)
        RT(1) RT(0)
        RT(1 // 6) RT(1 // 6)
    ]

    RKC = [RT(0), RT(1 // 2), RT(1 // 2), RT(1)]

    RKW = [RT(1 // 6), RT(1 // 3), RT(1 // 3), RT(1 // 6)]

    LowStorageRungeKutta3N(F, RKA, RKB, RKC, RKW, Q; dt = dt, t0 = t0)
end

"""
    LS3NRK33Heuns(f, Q; dt, t0 = 0)

This function returns a [`LowStorageRungeKutta3N`](@ref) time stepping object
for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This method uses the 3-stage, third-order Heun's Runge--Kutta scheme.

### References
    @article {Heun1900,
       title = {Neue Methoden zur approximativen Integration der
       Differentialgleichungen einer unabh\"{a}ngigen Ver\"{a}nderlichen}
       author = {Heun, Karl},
       journal = {Z. Math. Phys},
       volume = {45},
       pages = {23--38},
       year = {1900}
    }
"""
function LS3NRK33Heuns(
    F,
    Q::AT;
    dt = nothing,
    t0 = 0,
) where {AT <: AbstractArray}
    T = eltype(Q)
    RT = real(T)

    RKA = [
        RT(0) RT(0)
        RT(0) RT(1)
        RT(-1) RT(1 // 3)
    ]

    RKB = [
        RT(1 // 3) RT(0)
        RT(2 // 3) RT(-1 // 3)
        RT(3 // 4) RT(1 // 4)
    ]

    RKC = [RT(0), RT(1 // 3), RT(2 // 3)]

    RKW = [RT(1 // 4), RT(0), RT(3 // 4)]

    LowStorageRungeKutta3N(F, RKA, RKB, RKC, RKW, Q; dt = dt, t0 = t0)
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
@testset "Hulls" begin
  @testset "GrahamScan" begin
    # basic test
    pset = PointSet(rand(P2, 100))
    chul = hull(pset, GrahamScan())
    @test all(pset .∈ Ref(chul))

    # duplicated points
    pset = PointSet([rand(P2, 100); rand(P2, 100)])
    chul = hull(pset, GrahamScan())
    @test all(pset .∈ Ref(chul))

    # corner cases
    pset = PointSet(P2[(0,0)])
    chul = hull(pset, GrahamScan())
    @test chul == P2(0,0)
    pset = PointSet(P2[(0,1),(1,0)])
    chul = hull(pset, GrahamScan())
    @test chul == Segment(P2(1,0), P2(0,1))
    pset = PointSet(P2[(1,0),(0,0),(0,1)])
    chul = hull(pset, GrahamScan())
    @test vertices(boundary(chul)) == P2[(0,0),(1,0),(0,1)]

    # original point set is already in hull
    pset  = PointSet(P2[(0,0),(1,0),(1,1),(0,1),(0.5,-1)])
    chul  = hull(pset, GrahamScan())
    verts = vertices(boundary(chul))
    @test verts == P2[(0.5,-1),(1,0),(1,1),(0,1),(0,0)]

    # random points in interior do not affect result
    p1 = P2[(0,0),(1,0),(1,1),(0,1),(0.5,-1)]
    p2 = P2[0.5.*(rand(), rand()) .+ 0.5 for _ in 1:10]
    pset  = PointSet([p1; p2])
    chul  = hull(pset, GrahamScan())
    verts = vertices(boundary(chul))
    @test verts == P2[(0.5,-1),(1,0),(1,1),(0,1),(0,0)]
  end
end


text:
__precompile__()

module TypeFuzz

export testtypeintersect, testcases

using Base: unwrap_unionall, rewrap_unionall

randtype(const_allowed=false) = randtype(Base.Random.GLOBAL_RNG, const_allowed)

function randtype(rng::AbstractRNG, const_allowed=false)
    i = rand(rng, 1:(const_allowed ? 5 : 4))
    if i == 1
        return Int64
    elseif i == 2
        return Float64
    elseif i == 3
        return Val{randtype(rng, true)}
    elseif i == 4
        num = rand(rng, 1:5)
        return Tuple{[randtype(rng, true) for _ in 1:num]...}
    else
        return rand(rng, 1:5)
    end
end

randsupertype(t, varprefix::String="T") = randsupertype(Base.Random.GLOBAL_RNG, t, varprefix)

function randsupertype(rng::AbstractRNG, t, varprefix::String="T")
    vars = Pair{TypeVar, Any}[]
    r = _randsupertype!(rng, t, vars, varprefix)
    for v in vars
        r = UnionAll(v[1], r)
    end
    if !(t <: r)
        @show t
        @show r
        @show vars
        error("randsupertype invalid")
    end
    return r
end

function _randsupertype!(rng::AbstractRNG, t::DataType, vars, varprefix)
    np = length(t.parameters)
    if rand(rng, 1:np+2) == 1
        return genvar(rng, t, vars, varprefix)
    else
        if np == 0
            return t
        else
            return t.name.wrapper{[_randsupertype!(rng, tp, vars, varprefix) for tp in t.parameters]...}
        end
    end
end

function _randsupertype!(rng::AbstractRNG, t::Any, vars, varprefix)
    if rand(rng, 1:2) == 1
        return genvar(rng, t, vars, varprefix)
    else
        return t
    end
end

function genvar(rng::AbstractRNG, t, vars, varprefix)
    tv = nothing
    for v in vars
        if t == v[2] && rand(rng, 1:2) == 1
            tv = v[1]
        end
    end
    if tv === nothing
        if length(vars) >= 8
            return t
        else
            tv = TypeVar(Symbol(varprefix, length(vars)+1))
            push!(vars, tv => t)
        end
    end
    return tv
end

istuple(x) = x isa DataType && x <: Tuple
isval(x) = x isa DataType && x <: Val

function tryminimize(f::Function, args...)
    if any(arg -> arg isa UnionAll, args)
        args´ = tryminimize(map(unwrap_unionall, args)...) do args´´...
            f(map(rewrap_unionall, args´´, args)...)
        end
        return map(rewrap_unionall, args´, args)
    end

    if any(arg -> arg isa DataType, args)
        # at least one DataType instance (i.e. not a TypeVar)
        if all(arg -> !(arg isa DataType) || length(arg.parameters) == 1, args)
            # try to unwrap single element DataTypes
            args´ = map(arg -> arg isa DataType ? arg.parameters[1] : arg, args)
            if f(args´...)
                return tryminimize(f, args´...)
            end
        end

        np = minimum(arg -> arg isa DataType ? length(arg.parameters) : typemax(Int), args)
        if np == 0
            return args
        end
        if all(arg -> arg isa TypeVar || istuple(arg), args) && any(istuple, args)
            # Tuples or TypeVars, but at least one Tuple
            for i in 1:np
                # try removing tuple elements
                args´ = map(arg -> arg isa TypeVar ? arg : Tuple{arg.parameters[[1:i-1; i+1:end]]...}, args)
                if f(args´...)
                    return tryminimize(f, args´...)
                end
            end
        end

        paramsmat = Base.typed_hcat(Any, (collect(args[j].parameters[1:np]) for j in 1:length(args) if args[j] isa DataType)...)
        for i in 1:np
            # minimize parameters
            paramsmat[i,:] = collect(tryminimize(paramsmat[i,:]...) do paramsi´...
                args´ = collect(args)
                j´ = 1
                for j in eachindex(args)
                    if args[j] isa DataType
                        args´[j] = args[j].name.wrapper{paramsmat[1:i-1,j´]...,  paramsi´[j´], paramsmat[i+1:end,j´]...}
                        j´ += 1
                    end
                end
                f(args´...)
            end)
        end
        args´ = collect(args)
        j´ = 1
        for j in eachindex(args)
            if args[j] isa DataType
                args´[j] = args[j].name.wrapper{paramsmat[:,j´]...}
                j´ += 1
            end
        end
        return args´
    end
    return args
end

mutable struct IntersectProc
    out::IO
    in::IO
    proc::IO
    function IntersectProc()
        (ioout, ioin, proc) = startiproc()
        new(ioout, ioin, proc)
    end
end

function startiproc()
    code_object = """
        while !eof(STDIN)
            eval(Main, deserialize(STDIN))
        end
        """
    readandwrite(`$(Base.julia_cmd()) --color=$(Base.have_color ? "yes" : "no")
                                      --eval $code_object`)
end

function Base.close(proc::IntersectProc)
    close(proc.in)
    wait(proc.proc)
end

function trytypeintersect(proc::IntersectProc, a, b)
    if !process_running(proc.proc)
        println("Subprocess not running, restarting...")
        (proc.out, proc.in, proc.proc) = startiproc()
    end
    serialize(proc.in, quote
        try
            r = typeintersect($a, $b)
            serialize(STDOUT, r)
        catch ex
            serialize(STDOUT, ex)
        end
    end)
    tim = Timer((::Timer) -> (println("Timeout!"); kill(proc.proc, Base.SIGKILL)), 15)
    if eof(proc.out)
        close(tim)
        println("Premature end of data - slave died?")
        error("oops")
    end
    res = deserialize(proc.out)
    close(tim)
    if res isa Exception
        throw(res)
    else
        return res
    end
end

testtypeintersect(n::Int=1) = testtypeintersect(Base.Random.GLOBAL_RNG, n)

function testtypeintersect(rng::AbstractRNG, n::Int=1)
    failures = Expr[]
    proc = IntersectProc()
    for i in 1:n
        _testtypeintersect!(proc, rng, failures)
    end
    close(proc)
    return failures
end

function _testtypeintersect!(proc::IntersectProc, rng::AbstractRNG, failures)
    x = randtype(rng)
    a = randsupertype(rng, x, "T")
    b = randsupertype(rng, x, "S")
    r1 = nothing
    r2 = nothing
    try
        r1 = trytypeintersect(proc, a, b)
    catch e
        a´, b´ = tryminimize(a, b) do a, b
            try
                trytypeintersect(proc, a, b)
            catch e
                return true
            end
            return false
        end
        push!(failures, :(typeintersect($a´, $b´)))
    end
    try
        r2 = trytypeintersect(proc, b, a)
    catch e
        a´, b´ = tryminimize(a, b) do a, b
            try
                trytypeintersect(proc, b, a)
            catch e
                return true
            end
            return false
        end
        push!(failures, :(typeintersect($b´, $a´)))
    end
    if r1 !== nothing && r2 !== nothing && r1 != r2
        a´, b´ = tryminimize(a, b) do a, b
            try
                trytypeintersect(proc, a,b) != trytypeintersect(proc, b, a)
            catch e
                return false
            end
        end
        push!(failures, Expr(:let, :(typeintersect(a, b) == typeintersect(b, a)), :(a = $(a´)), :(b = $(b´))))
    end
    if r1 !== nothing
        if !(x <: r1)
            a´, b´, x´ = tryminimize(a, b, x) do a, b, x
                try
                    x <: a && x <: b && !(x <: trytypeintersect(proc, a,b))
                catch e
                    return false
                end
            end
            push!(failures, :($(x´) <: typeintersect($(a´), $(b´))))
        end
        if !(r1 <: a)
            a´, b´ = tryminimize(a, b) do a, b
                try
                    !(trytypeintersect(proc, a,b) <: a)
                catch e
                    return false
                end
            end
            push!(failures, Expr(:let, :(typeintersect(a, $(b´)) <: a), :(a = $(a´))))
        end
        if !(r1 <: b)
            a´, b´ = tryminimize(a, b) do a, b
                try
                    !(trytypeintersect(proc, a,b) <: b)
                catch e
                    return false
                end
            end
            push!(failures, Expr(:let, :(typeintersect($(a´), b) <: b), :(b = $(b´))))
        end
    end
    if r2 !== nothing && r1 != r2
        if !(x <: r2)
            a´, b´, x´ = tryminimize(a, b, x) do a, b, x
                try
                    x <: a && x <: b && !(x <: trytypeintersect(proc, b, a))
                catch e
                    return false
                end
            end
            push!(failures, :($(x´) <: typeintersect($(b´), $(a´))))
        end
        if !(r2 <: a)
            a´, b´ = tryminimize(a, b) do a, b
                try
                    !(trytypeintersect(proc, b, a) <: a)
                catch e
                    return false
                end
            end
            push!(failures, Expr(:let, :(typeintersect($b´, a) <: a), :(a = $(a´))))
        end
        if !(r2 <: b)
            a´, b´ = tryminimize(a, b) do a, b
                try
                    !(trytypeintersect(proc, b, a) <: b)
                catch e
                    return false
                end
            end
            push!(failures, Expr(:let, :(typeintersect(b, $a´) <: b), :(b = $(b´))))
        end
    end
    return failures
end

if VERSION >= v"0.7.0-DEV.357"
    testcases(failures::Vector{Expr}) = Expr(:block, (Expr(:macrocall, Symbol("@test_broken"), nothing, failure) for failure in failures)...)
else
    testcases(failures::Vector{Expr}) = Expr(:block, (:(@test_broken $failure) for failure in failures)...)
end

end


text:
using IndexedTables
using TableTraits
using NamedTuples
using Base.Test

@testset "TableTraits" begin

source_nds = NDSparse(Columns(a=[1,2,3]), Columns(b=[1.,2.,3.], c=["A","B","C"]))

@test isiterable(source_nds) == true

target_array_nds = collect(getiterator(source_nds))

@test length(target_array_nds) == 3
@test target_array_nds[1] == @NT(a=1, b=1., c="A")
@test target_array_nds[2] == @NT(a=2, b=2., c="B")
@test target_array_nds[3] == @NT(a=3, b=3., c="C")

source_array = [@NT(a=1,b=1.,c="A"), @NT(a=2,b=2.,c="B"), @NT(a=3,b=3.,c="C")]

it1 = NDSparse(source_array)
@test length(it1) == 3
@test it1[1,1.].c == "A"
@test it1[2,2.].c == "B"
@test it1[3,3.].c == "C"

it2 = NDSparse(source_array, idxcols=[:a])
@test length(it2) == 3
@test it2[1] == @NT(b=1., c="A")
@test it2[2] == @NT(b=2., c="B")
@test it2[3] == @NT(b=3., c="C")

it3 = NDSparse(source_array, datacols=[:b, :c])
@test length(it3) == 3
@test it3[1] == @NT(b=1., c="A")
@test it3[2] == @NT(b=2., c="B")
@test it3[3] == @NT(b=3., c="C")

it4 = NDSparse([(1=>"A"), (2=>"B")])
@test length(it4) == 2
@test it4[1] == "A"
@test it4[2] == "B"

source_nt = table([1,2,3],[1.,2.,3.],["A","B","C"], names=[:a,:b,:c])

target_array_nt = collect(getiterator(source_nt))

@test length(target_array_nt) == 3
@test target_array_nt[1] == @NT(a=1, b=1., c="A")
@test target_array_nt[2] == @NT(a=2, b=2., c="B")
@test target_array_nt[3] == @NT(a=3, b=3., c="C")

it4 = table(source_array, copy=true)
@test length(it4) == 3
@test it4[1] == @NT(a=1,b=1.,c="A")
@test it4[2] == @NT(a=2,b=2.,c="B")
@test it4[3] == @NT(a=3,b=3.,c="C")

# Use a Set to really hit the iterator constructor
as_set = Set(source_array)
it4 = table(as_set, copy=true)
@test length(it4) == 3
@test contains(==, as_set, it4[1])
@test contains(==, as_set, it4[2])
@test contains(==, as_set, it4[3])

@testset "avoid too narrow dict value type" begin
    # The array constructor
    @test table([@NT(col1 = "A")], pkey = [1]) ==
            table(["A"], names = [:col1], pkey = [1])

    # The iterator constructor
    @test table(Iterators.repeated(@NT(col1 = "A"), 2), pkey = [1]) ==
            table(["A", "A"], names = [:col1], pkey = [1])
end

# Non NamedTuples iterators

@test table([(1, 2), (3, 4)]) == table([1, 3], [2, 4])
@test table([@NT(a=1) => @NT(b=2)]) == table([1], [2], names = [:a, :b], pkey = :a)

end


text:
using Zygote

function lkl_quantile(transformation, batch, percentage=(0.0, 1.0))
    score = logpdf(transformation, batch)
    n = nobs(batch)
    l, r = percentage
    if l <= 0.0
        lq = 1
    else
        lq = Int(ceil(l*n))
    end
    if r >= 1.0
        rq = n
    else
        rq = Int(ceil(r*n))
    end
    perms = sortperm(score)
    indxs = perms[lq:rq]
    return batch[:, indxs]
end
Zygote.@nograd lkl_quantile


text:
import Distributions

function random(a::Float64, b::Float64)::Float64
    return rand(Distributions.Uniform(a, b))
end


text:
export DeterministicDistributionModel,  DynamicDistributionModel, ExperienceSampleModel, TimeBasedSampleModel, PrioritizedSweepingSampleModel

include("deterministic_distribution_model.jl")
include("dynamic_distribution_model.jl")
include("experience_sample_model.jl")
include("time_based_sample_model.jl")
include("prioritized_sweeping_sample_model.jl")


text:
module TestConstraintBridge

using Test

using MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

struct DummyBridge <: MOI.Bridges.Constraint.AbstractBridge end

function test_AbstractBridge()
    model = MOI.Utilities.Model{Float64}()
    bridge = DummyBridge()
    attr = MOI.ConstraintPrimalStart()
    @test !MOI.supports(model, attr, typeof(bridge))
    @test_throws MOI.UnsupportedAttribute(attr) MOI.set(
        model,
        attr,
        bridge,
        1.0,
    )
    attr = MOI.ConstraintFunction()
    err = MOI.SetAttributeNotAllowed(attr)
    @test_throws err MOI.set(model, attr, bridge, MOI.EqualTo(1.0))
end

end  # module

TestConstraintBridge.runtests()


text:
function read_ba_instance(fn)
    fid = open(fn)
    lines = readlines(fid)
    close(fid)
    line=split(lines[1]," ")
    n = parse(Int,line[1])
    m = parse(Int,line[2])
    p = parse(Int,line[3])
    off = 2

    one_cam = zeros(Float64,11,1)
    line=split(lines[off]," ")
    for i in 1:11
        one_cam[i] = parse(Float64,line[i])
    end
    cams = repeat(one_cam,1,n)
    off += 1

    one_X = zeros(Float64,3,1)
    line=split(lines[off]," ")
    for i in 1:3
        one_X[i] = parse(Float64,line[i])
    end
    X = repeat(one_X,1,m)
    off += 1

    one_w = parse(Float64,lines[off])
    w = repeat([one_w],1,p)
    off += 1

    one_feat = zeros(Float64,2,1)
    line=split(lines[off]," ")
    for i in 1:2
        one_feat[i] = parse(Float64,line[i])
    end
    feats = repeat(one_feat,1,p)

    camIdx = 1
    ptIdx = 1
    obs = zeros(Int,2,p)
    for i in 1:p
        obs[1,i] = camIdx
        obs[2,i] = ptIdx
        camIdx = (camIdx%n) + 1
        ptIdx = (ptIdx%m) + 1
    end

    (cams,X,w,obs,feats)
end


text:
using GLVisualize, GeometryTypes, GLAbstraction
using Colors, Reactive, FileIO
if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0,1,50))
end

description = """
Animating some billboard particles in 2D.
"""

n = 30
const border = 50f0
function bounce_particles(pos_velo, _)
    positions, velocity = pos_velo
    dt = 0.1f0
    @inbounds for i=1:length(positions)
        pos,velo = positions[i], velocity[i]
        positions[i] = Point2f0(pos[1], pos[2] + velo*dt)
        if pos[2] <= border
            velocity[i] = abs(velo)
        else
            velocity[i] = velo - 9.8*dt
        end
    end
    positions, velocity
end
start_position = (rand(Point2f0, n)*700f0) + border
position_velocity = foldp(bounce_particles,
    (start_position, zeros(Float32, n)),
    timesignal
)
circle = HyperSphere(Point2f0(0), 40f0)
vis = visualize((circle, map(first, position_velocity)),
    image = loadasset("kitty.png"),
    stroke_width = 3f0,
    stroke_color = RGBA{Float32}(0.91,0.91,0.91,1),
    boundingbox = AABB(start_position)
)
_view(vis, window, camera=:orthographic_pixel)

if !isdefined(:runtests)
    renderloop(window)
end


text:
module TestHessian
    using Base.Test
    import FiniteDiff

    negative_sin(x) = -sin(x)
    negative_cos(x) = -cos(x)

    f1(x) = sin(x[1]) + cos(x[2])
    function h1!(out, x)
        out[1, 1] = negative_sin(x[1])
        out[1, 2] = 0.0
        out[2, 1] = 0.0
        out[2, 2] = negative_cos(x[2])
        return
    end

    f2(x) = 2.3 * sin(x[1]) + 3.1 * cos(x[2])
    function h2!(out, x)
        out[1, 1] = 2.3 * negative_sin(x[1])
        out[1, 2] = 0.0
        out[2, 1] = 0.0
        out[2, 2] = 3.1 * negative_cos(x[2])
        return
    end

    f3(x) = 0.1 * (x[1] - 1.0)^2 - 0.3 * (x[2] - 0.7)^2
    function h3!(out, x)
        out[1, 1] = 0.1 * 2.0
        out[1, 2] = 0.0
        out[2, 1] = 0.0
        out[2, 2] = -0.3 * 2.0
        return
    end

    f4(x) = x[1]^2 * x[2]^2
    function h4!(out, x)
        out[1, 1] = 2.0 * x[2]^2
        out[1, 2] = 4.0 * x[1] * x[2]
        out[2, 1] = 4.0 * x[1] * x[2]
        out[2, 2] = 2.0 * x[1]^2
        return
    end

    function check_error(y_true, y_approximate, err = 1 // 10_000)
        n = length(y_true)
        @test length(y_approximate) == n
        for i in 1:n
            # TODO: Decide how to handle very small inputs.
            if y_true[i] > eps(y_true[i])
                @test abs(y_true[i] - y_approximate[i]) / abs(y_true[i]) < err
            end
        end
    end

    function run_tests(n::Integer)
        funcs = (
            (f1, h1!),
            (f2, h2!),
            (f3, h3!),
            (f4, h4!),
        )

        x = Array(Float64, 2)
        output = Array(Float64, 2, 2)
        buffer = Array(Float64, 2)
        H = Array(Float64, 2, 2)

        n_dims = length(x)

        for (f, f′′!) in funcs
            for itr in 1:n
                for i in 1:n_dims
                    x[i] = itr
                end

                f′′!(H, x)

                FiniteDiff.hessian!(output, f, x, buffer)
                check_error(H, output)

                output2 = FiniteDiff.hessian(f, x)
                check_error(H, output2)

                tmp! = FiniteDiff.hessian(f, mutates = true)
                tmp!(output, x, buffer)
                check_error(H, output)

                tmp = FiniteDiff.hessian(f, mutates = false)
                output2 = tmp(x)
                check_error(H, output2)
            end
        end
    end

    @testset "hessian tests" begin
        run_tests(100)
    end
end


text:
import NonUniformRandomVariateGeneration: sampleSortedUniforms!,
  sampleCategorical, sampleCategoricalSorted!

# sample nans categorical r.v.s where Pr(X=i) ∝ p[i]
# required: length(Fs) = length(p), unifs[start+1:start+nans] are valid
@inline function _sampleCategoricalSorted!(vs::Vector{Int64}, p::Vector{Float64},
  nans::Int64, unifs::Vector{Float64}, Fs::Vector{Float64}, start::Int64,
  rng::RNG)
  sampleSortedUniforms!(unifs, start, nans, rng)
  cumsum!(Fs, p)
  @inbounds maxFs = Fs[length(p)]
  j = 1
  for i = 1:nans
    @inbounds while maxFs * unifs[start + i] > Fs[j]
      j += 1
    end
    @inbounds vs[start + i] = j
  end
end

@inline function _cresample!(smcio::SMCIO)
  smcio.internal.as[1] = 1
  _sampleCategoricalSorted!(smcio.internal.as, smcio.ws, smcio.N-1,
    smcio.internal.scratch1, smcio.internal.scratch2, 1, getRNG())
end

@inline function _resample!(smcio::SMCIO)
  sampleCategoricalSorted!(smcio.internal.as, smcio.ws,
    smcio.internal.scratch1, smcio.internal.scratch2, getRNG())
end