text:
### Baseline policy for the occlusion scenario
# move forward
# stop
# check for N seconds
# if free go, else wait

@with_kw mutable struct WaitAndGo
    env::CrosswalkEnv = CrosswalkEnv()
    initialstate::Vehicle = initial_ego(env, MersenneTwister(1))
    max_acc::Float64 = 2.
    max_dec::Float64 = -4.0
    Δt::Float64 = 0.1
    N0::Int64 = 3# Number of steps to wait
    threshold::Float64 = 2.5
    dist_delta::Float64 = 0.
    verbose::Bool = false

    # states
    reaching::Bool = true
    wait::Bool = false
    go::Bool = false
    N::Int64 = 0.
end

function action(policy::WaitAndGo, b::Vector{Vehicle})
    verbose = policy.verbose
    crosswalk_x = policy.env.params.roadway_length/2 - policy.env.params.crosswalk_width - policy.dist_delta
    v_max = policy.env.params.speed_limit
    ego = b[1].state
    @assert b[1].id == 1

    init_dec = -policy.initialstate.state.v^2/(2*(crosswalk_x - policy.initialstate.state.posG.x))
    @assert init_dec > policy.max_dec # negative numbers

    acc = 0.

    if policy.reaching #TODO use a smarter controller
        # if ego.v*cos(ego.posG.θ) < 0.
        #     if verbose; println("adjusting") end
        #     acc = -ego.v*cos(ego.posG.θ)/policy.Δt
        if ego.v ≈ 0.
            acc = 0.
        else
            acc = init_dec
        end
        # acc = init_dec
        if verbose; println("reaching phase, acc = $acc") end
        if !policy.wait && (ego.v ≈ 0. && ego.posG.x > policy.initialstate.state.posG.x) # change state
            policy.reaching = false
            policy.wait = true
        end
        return acc
    end

    if !policy.go
        if policy.N == policy.N0
            if verbose; println("GO!") end
            policy.go = true
            policy.wait = false
        else
            # measure ttc
            ttc_min = Inf
            for ped in b[2:end]
                ttc = (-env.params.lane_width/2 - ped.state.posG.y)/ped.state.v
                if verbose; println("pedestrian $(ped.id): ttc $(ttc) ; v $(ped.state.v)") end
                if 0 < ttc < ttc_min
                    ttc_min = ttc
                end
            end
            if verbose; println("minimum ttc: $ttc_min ; $(length(b)-1) pedestrians detected") end
            # check
            if ttc_min > policy.threshold
                if verbose; println("keep checking $(policy.N)") end
                # policy.wait = false
                policy.N += 1
            else
                if verbose; println("reset") end
                policy.N = 0
            end
        end
    end

    if policy.go
        acc = policy.max_acc
        return acc
    end
    return 0.
end

mutable struct WaitAndGoUpdater
    policy::WaitAndGo
    problem::OCPOMDP
end

"""
Trivial belief updater for the baseline policy, just return the observation
"""
function POMDPs.update(up::WaitAndGoUpdater, bold::Vector{Vehicle}, a::Float64, o::Vector{Vehicle})
    return o
end

"""
    reset_policy!(policy::waitAndGo, ego::Vehicle)
Reset the state of the FSM, must be run before each simulation. Change in place the fields of the
policy
"""
function reset_policy!(policy::WaitAndGo, ego::Vehicle)
    policy.initialstate = ego
    policy.reaching = true
    policy.wait = false
    policy.go = false
    policy.N = 0
end

const WaitAndGoModel = CrosswalkDriver{Float64,Vector{Vehicle},Vector{Vehicle},WaitAndGo,WaitAndGoUpdater,SimpleSensor}

function reset_model!(model::WaitAndGoModel,
                             ego::Vehicle)
        model.tick = 0
        model.a = 0.
        model.b = Vehicle[]
        model.o = Vehicle[]
        model.ego = ego
        reset_policy!(model.policy, ego)
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
module ITF1788

using IntervalArithmetic, LinearAlgebra

export parse_block, parse_command, generate

# Write your package code here.

include("generate.jl")
include("parse.jl")
include("functions.jl")
end

ext: jl


text:
using BandedMatrices

function sband(n::Integer, a::Vector{T}) where T
    k = length(a)
    SH = T<:Complex ? Hermitian : Symmetric
    SH(BandedMatrix([-i => fill(a[i+1], n-i) for i = 0:min(k,n)-1]...), :L)
end

function example_1(n, T=Float64)
    tt = sband(n, T[4, 1])
    ss = sband(n, T[2, 1])
    tt, ss
end

det0(A, B, s) = det(Matrix(A - s * B))
function det1(A, B, s)
    h = sqrt(eps(max(abs(s), 1.0)))
    ( det0(A, B, s + h) - det0(A, B, s - h) ) / (2*h)
end
function det2(A, B, s)
    h = eps((max(abs(s), 1)))^(1/3)
    ( det0(A, B, s + h) + det0(A, B, s - h) -2*det0(A, B, s) ) / h^2
end

@testset "basics n=$n s=$s" for n = [1, 2, 3, 10], s = [1.0, 0.5]
    A, B = sband(n, [0.1]), sband(n, [1.0])
    κ, η, ζ = detderivates(A, B, s)
    @test κ == n
    @test η ≈ -n / (0.1 - s)
    @test ζ ≈ n / (0.1 - s)^2
end

@testset "small real $n" for n = 1:4
    A, B = sband(n, [6.0, -4, 1]), sband(n, [2.0, 1])
    s = 0.001
    d0 = det0(A, B, s)
    κ, η, ζ = detderivates(A, B, s)
    @test κ == 0
    @test isapprox(η, det1(A, B, s) / d0; rtol=1e-3)
    @test isapprox(η^2 - ζ, det2(A, B, s) / d0; atol=1e-5, rtol=1e-3)
end

@testset "small complex $n" for n = 1:4
    A, B = sband(n, [6.0, -4, 1+0im]), sband(n, [2.0, 1im])
    s = 0.001
    d0 = det0(A, B, s)
    κ, η, ζ = detderivates(A, B, s)
    @test η isa Real
    @test ζ isa Real
    @test κ == 0
    @test isapprox(η, det1(A, B, s) / d0; rtol=1e-3)
    @test isapprox(η^2 - ζ, det2(A, B, s) / d0; atol=1e-5, rtol=1e-3)
end

kappa(x) = x < 0.0 ? 0 : 2
kappa1(x) = x < 1.0 ? 0 : 3
kappa2(x) = x < 0.2 ? 0 : x < 1 ? 1 : x < 2 ? 3 : typemax(Int)

MatrixAlgebra.detderivates!(f::Function, x) = f(x), 0.0, 0.0

@testset "bisection bounds 0-2" begin
    lb, ub = eigbounds(1, 2, kappa, 1.0)
    @test  all(isapprox.(lb, 0.0; atol=eps()))
    @test  all(isapprox.(ub, 0.0; atol=eps()))
end
@testset "bisection bounds 1-3" begin
    lb, ub = eigbounds(1, 3, kappa1, 1.0)
    @test  all(isapprox.(lb, 1.0; atol=eps()))
    @test  all(isapprox.(ub, 1.0; atol=eps()))
end
@testset "bisection bounds 2-4" begin
    lb, ub = eigbounds(1, 4, kappa2, 1.0)
    @test all(isapprox.(lb, [0, 1, 1, 1]))
    @test all(isapprox.(ub, [0.5, 1, 1, 3]))
end
@testset "bisection bounds 2-5" begin
    lb, ub = eigbounds(1, 5, kappa2, 1.0)
    @test all(isapprox.(lb, [0, 1, 1, 2, 2]))
    @test all(isapprox.(ub, [0.5, 1, 1, 2, 2]))
end


text:
using GMC_NS, Distributions

prior=Uniform(0,1)
priors=[prior,prior]

box=[0. 1.
     0. 1.]

gmc=GMC_DEFAULTS
gmc[1]=180
gmc[2]=eps()

e=Eggbox_Ensemble("Eggbox_test", 1000, priors, box, gmc...)

uds=Vector{Vector{Function}}([[evidence_display],[liwi_display],[convergence_display]])
lds=Vector{Vector{Function}}([[ensemble_display]])

converge_ensemble!(e,backup=(true,100),upper_displays=uds,lower_displays=lds, disp_rot_its=1000, converge_factor=1e-6)


text:
using NCDatasets
using Formatting
using ArgParse
using JSON
include("./HOOM_z_res.jl")

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--resolution"
            help = "Resolution keywords. Currently accept: POP2, Standard."
            arg_type = String
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

zs = getHOOMResolution(parsed["resolution"])

Nz = length(zs) - 1
N_zs = length(zs)

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "Nz", Nz)
    defDim(ds, "zs", N_zs)

    defVar(ds, "zs", zs, ("zs",))
    defVar(ds, "mid_zs", (zs[1:end-1] + zs[2:end]) / 2.0, ("Nz",))

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
#=
julia --sysimage sys_plots.so
=#

# Trying out GR, StatsPlots, Plots
# Need:
# - Heatmaps
# - KDE
# - Histograms
# - Grouped line plots
# - Grouped scatter plots

import Pkg; Pkg.activate(".")
@time begin
  using Plots; gr()
end

@time begin
  x, y = randn(100), randn(100)
  p1 = plot(x, y, label="p1");
  p2 = plot(x, y, seriestype=:scatter, label="p2");
  p3 = plot(x, seriestype=:hist, label="p3");
  p4 = plot(y, seriestype=:hist, label="p4", color="red");
  p = plot(p1, p2, p3, p4, layout=(2, 2));
  savefig(p, "test.pdf")
end

@time begin
  a = plot(randn(6, 6), seriestype=:heatmap);
  b = plot(randn(6, 6), seriestype=:heatmap);
  pp = plot(a, b, layout=(2, 1));
  pp
end

max_stars_repo_name: luiarthur/TuringBnpBenchmarks

max_issues_repo_name: luiarthur/TuringBnpBenchmarks

max_forks_repo_name: luiarthur/TuringBnpBenchmarks


text:
# to do the initialiteration
"""
let's think what I need
* well, maybe I can use to get the singular values
* maybe I can get the relative entopy by the singular values
* do the initial iteration on the trotter direction
* do I need to do it on the spatial direction too?
"""

"""
	getSingularValues
Return the singular values

----
notes:
maybe it should be done in the way, take the tensorT
and get the singular values
(when a simulator is passed) any simulator??
	1. get the tensorT
	2. pass it to a function.

the unit to get the singular values when a tensor is passed.
	1. form MM+ : this depends on the which type of the simulator...
		can I deal with it?
	2. do the svd
	3. the middle vector.

or it could be done from on the way..
I mean anyway in the process they get some singular values.. so.
But I guess it should be like 'printSingularValues'...
Besides, when should I do? well, the first answer would be
when I form the Core tensor.
But, no, forget it, it goes down too deep...

Maybe I can recycle some functions I created before..
in the case of Quantum2dFractal: calculateCoreTensors(tensorT, tensorTtilde, dimM)
Actually more properly: getTensorUy(dimM, "x", tensorT, mirror=mirror)
	that's true that it depends on the direction... well, well, well

	calculateCoreTensors <- add arg getsvd
	getTensorUy
	getTensorUz
	getProgectorFromMMd: return tensorU, lambdaVector - DONE
"""
function getSingularValues(simulator::Simulator)

end

function getSingularValues{T}(tensorT::Array{T,N})
end


text:
using Base.Iterators
using Particles: instantiate
using Distributions, ConjugatePriors

@testset "Labeled observations" begin
    ys = rand(MixtureModel([Normal(-2.), Normal(2.)]), 20)

    labels = ifelse.(ys .> 0, 1, 2)

    ylab = Labeled.(ys, labels)

    prior = NormalInverseChisq(0., 1., 2., 3.)

    p = InfiniteParticle(Component.([prior, prior]),
                         nothing,
                         0,
                         1.,
                         Component(prior),
                         NStatePrior([1., 1.]))

    @testset "putatives with labeled" begin
        lab1, lab2 = instantiate.(putatives(p, ys[1]))
        @test lab1.components == putatives(p, Labeled(ys[1], 1))[1].components
        @test lab2.components == putatives(p, Labeled(ys[1], 2))[1].components

        lab1m, lab2m = instantiate.(putatives(p, Labeled(ys[1])))
        @test lab1m.components == lab1.components
        @test lab2m.components == lab2.components
    end

    @testset "filter equivalence with missing" begin
        Random.seed!(1)
        fil1 = reduce(fit!, ys, init=FearnheadParticles([p], 10))

        Random.seed!(1)
        fil2 = reduce(fit!, Labeled.(ys, missing), init=FearnheadParticles([p], 10))

        @test getfield.(fil1.particles, :components) ==
            getfield.(fil2.particles, :components)
    end

    @testset "filtering with labels" begin
        fil = reduce(fit!, ylab, init=FearnheadParticles([p], 10))
        # only one particle:
        @test length(fil.particles) == 1
        pp = first(fil.particles)
        @test nobs(pp) == length(ylab)
        @test nobs.(pp.components) == [sum(labels.==1), sum(labels.==2)]

        # same results if missing labels are before or after:
        some_unlab = Labeled.([-1., 0, 1], missing)
        fil_before = reduce(fit!, append!(copy(some_unlab), ylab),
                            init=FearnheadParticles([p], 10))
        fil_after = reduce(fit!, prepend!(copy(some_unlab), ylab),
                           init=FearnheadParticles([p], 10))

        # test equality of component parameters
        for (p1,p2) in zip(fil_before.particles, fil_after.particles)
            for (c1,c2) in zip(p1.components, p2.components)
                @test all(params(c1) .≈ params(c2))
            end
        end

    end
end


text:
#####
##### Interface and implementation for split-apply-combine.
#####

export by, RepeatRow, ignoreindex, aggregator

"""
RepeatRow(row)

A row repeated as many times as needed. Can be `merge`d to a `FunctionalTable`, or
instantiated with `FunctionalTable(len, repeat_row)`.
"""
struct RepeatRow{T <: NamedTuple}
    row::T
end

function FunctionalTable(len::Int, repeat_row::RepeatRow,
                         ordering::OrderingRule = TrustOrdering())
    FunctionalTable(TrustLength(len), map(v -> RepeatValue(v, len), repeat_row.row),
                    ordering)
end

Base.merge(row::RepeatRow, ft::FunctionalTable; kwargs...) =
    merge(FunctionalTable(length(ft), row), ft; kwargs...)

Base.merge(ft::FunctionalTable, row::RepeatRow; kwargs...) =
    merge(ft, FunctionalTable(length(ft), row); kwargs...)

"""
$(SIGNATURES)

Prepend the `index` as repeated columns to `f(index, tables...)`.
"""
fuse(f, index::NamedTuple, tables...) =
    merge(RepeatRow(index), FunctionalTable(f(index, tables...)))

"""
$(TYPEDEF)

Implements [`by`](@ref) using an iterator.

# Internals

Each rows from the underlying `FunctionalTable` is split into `index` and `rest`.

Iterator state is

1. `nothing` when the rows of the underlying FunctionalTable have been exhausted,

2. `index`, `rest`, `itrstate` for the next block, where `index` and `rest` are the first
(mismatching) row that has *not* been pushed to the buffers.
"""
struct SplitTable{K <: NamedTuple, B <: NamedTuple, O <: TableOrdering, T <: FunctionalTable}
    block_buffers::B
    block_ordering::O
    ft::T
    function SplitTable{K}(block_buffers::B, block_ordering::O, ft::T) where {K, B, O, T}
        new{K, B, O, T}(block_buffers, block_ordering, ft)
    end
end

"""
$(SIGNATURES)

Helper function to set up a `SplitTable`. *Internal*.
"""
function split_table(ft::FunctionalTable{C}, splitkeys::Keys) where {C}
    @argcheck is_prefix(splitkeys, orderkey.(ordering(ft)))
    block_ordering = ()     # FIXME calculate this
    T = NamedTuple{fieldnames(C)}(map(eltype, fieldtypes(C)))
    index, rest = split_namedtuple(NamedTuple{splitkeys}, T)
    block_buffers = map(V -> Vector{V}(), rest)
    K = NamedTuple{keys(index), Tuple{values(index)...}}
    SplitTable{K}(block_buffers, block_ordering, ft)
end

Base.IteratorSize(::Type{<:SplitTable}) = Base.SizeUnknown()

Base.IteratorEltype(::Type{<:SplitTable}) = Base.HasEltype()

Base.eltype(::Type{<: SplitTable{K,B,O}}) where {K,B,O} = Tuple{K, FunctionalTable{B,O}}

ordering(st::SplitTable{K}) where K = mask_ordering(ordering(st.ft), fieldnames(K))

function Base.iterate(g::SplitTable{K}) where K
    row, itrstate = @ifsomething iterate(g.ft)
    index, rest = split_namedtuple(K, row)
    iterate(g, (index, rest, itrstate))
end

function Base.iterate(g::SplitTable{K}, state) where K
    @unpack ft, block_buffers, block_ordering = g
    map(b -> resize!(b, 0), block_buffers) # FIXME use foreach?
    index, rest, itrstate = @ifsomething state
    len = 0
    _block() = (index, FunctionalTable(TrustLength(len),
                                       # FIXME not entirely sure we need to copy
                                       # if that is emphasized in the semantics
                                       map(copy, block_buffers),
                                       TrustOrdering(block_ordering)))
    while true
        map(push!, block_buffers, rest) # FIXME use foreach?
        len += 1
        y = iterate(ft, itrstate)
        y ≡ nothing && return _block(), nothing # done
        row, itrstate = y
        row_index, row_rest = split_namedtuple(K, row)
        if index == row_index
             rest = row_rest
        else
           return _block(), (row_index, row_rest, itrstate)
        end
    end
end

####
#### by and its implementation
####

"""
$(SIGNATURES)

An iterator that groups rows of tables by the columns `splitkeys`, returning
`(index::NamedTupe, table::FunctionalTable)` for each contiguous block of the index keys.

The function has a convenience form `by(ft, splitkeys...; ...)`.
"""
function by(ft::FunctionalTable, splitkeys::Keys)
    # TODO by could be very clever here by only sorting the subgroup which is unsorted,
    # and giving views of the underlying vectors instead of creating new tables
    sorted_ft = is_prefix(splitkeys, orderkey.(ordering(ft))) ? ft :
        sort(ft, split_compatible_ordering(ordering(ft), splitkeys))
    split_table(sorted_ft, splitkeys)
end

by(ft::FunctionalTable, splitkeys::Symbol...; kwargs...) = by(ft, splitkeys; kwargs...)

"""
$(SIGNATURES)

Map a table split with [`by`](@ref) using `f`.

Specifically, `f(index, table)` receives the split index (a `NamedTuple`) and a
`FunctionalTable`.

It is supposed to return an *iterable* that returns rows (can be a `FunctionalTable`). These
will be prepended with the corresponding index, and collected into a `FunctionalTable` with
`cfg`.

When `f` returns just a single row (eg aggregation), wrap by `Ref` to create a
single-element iterable.
"""
function Base.map(f, st::SplitTable; cfg = SINKCONFIG)
    # FIXME: idea: custom ordering override? would that make sense?
    FunctionalTable(Iterators.flatten(imap(args -> fuse(f, args...), st)),
                    TrustOrdering(ordering(st)); cfg = cfg)
end

"""
$(SIGNATURES)

Wrap a function returning a closure that ignores the first argument (the index in
[`by`](@ref) mappings).
"""
@inline ignoreindex(f) = (_, args...) -> f(args...)

"""
$(SIGNATURES)

Wrap a function so that it maps columns of a `FunctionalTable` to a table with a single row,
columwise, ignoring the index. Returns a closure.

# Example

```julia
map(aggregator(mean), by(ft, :col))
```

will calculate means after grouping by `:col`.
"""
@inline aggregator(f) = ignoreindex(ft -> Ref(map(f, columns(ft))))


text:
import MixedModels: MixedModel
import Tables: ColumnTable
import DataFrames: DataFrame
import RCall: rcopy,
              RClass,
              rcopytype,
              reval,
              S4Sxp,
              sexp,
              protect,
              unprotect,
              sexpclass,
              @rput,
              @rget

function sexp(::Type{RClass{:merMod}}, x::MixedModel{T}) where T
    throw(ArgumentError("You must use a tuple with the data -- (MixedModel, DataFrame) -- not just a model"))
end

sexpclass(x::MixedModel{T}) where T = RClass{:merMod}

# we could in theory support the other ordering by re-ordering and passing on
# but that introduces additional maintainance work and encourages messy style
function sexp(::Type{RClass{:merMod}}, x::Tuple{Union{DataFrame, ColumnTable}, MixedModel})
    throw(ArgumentError("The order in your Tuple is reversed. It should be (Model, Data)"))
end

sexpclass(x::Tuple{Union{DataFrame, ColumnTable}, MixedModel}) where T = RClass{:merMod}


text:
using LinearAlgebraicRepresentation, ViewerGL
Lar = LinearAlgebraicRepresentation
GL = ViewerGL

filename = "/home/dario/Plasm.jl/test/svg/Lar.svg"
V,EV = Lar.svg2lar(filename)

GL.VIEW([
	GL.GLLines(V,EV,GL.COLORS[1]),
	GL.GLFrame2
]);

function pointsinout(V,EV, n=10000)
	point_in = []
	point_out = []
	point_on = []
	classify = Lar.pointInPolygonClassification(V,EV)
	for k=1:n
		queryPoint = rand(1,2)
		inOut = classify(queryPoint)
		# println("k = $k, queryPoint = $queryPoint, inOut = $inOut")
		if inOut=="p_in"
			push!(point_in, queryPoint);
		elseif inOut=="p_out"
			push!(point_out, queryPoint);
		elseif inOut=="p_on"
			push!(point_on, queryPoint);
		end
	end
	#GL.GLPoints(queryPoint,GL.COLORS[2])
	return point_in,point_out,point_on
end

points_in, points_out, points_on = pointsinout(V,EV);
pointsin = [vcat(points_in...) zeros(length(points_in),1)]
pointsout = [vcat(points_out...) zeros(length(points_out),1)]

polygon = [GL.GLLines(V,EV,GL.COLORS[1])];
in_mesh = [GL.GLPoints(pointsin, GL.COLORS[2])]
out_mesh = [GL.GLPoints(pointsout, GL.COLORS[3])]

result = cat([polygon,in_mesh,out_mesh])
GL.VIEW(result);


text:
@testset "kron_fast" begin
    p, q, n, m = 3, 4, 5, 6
    A = randn(p, q)
    B = randn(q*n, m)
    # Do not use I(n) to allow compatibility with Julia v1.0
    eye = Matrix{Float64}(I, n, n)
    @test kron(eye, A) * B ≈ kron_fastl(A, B)
    @test kron(A, eye) * B ≈ kron_fastr(A, B)
end

@testset "getscore" begin
    T, K, N = 5, 3, 2
    X = randn(T, K)
    resid = randn(T, N)
    @test getscore(X, resid) == kron(X, ones(1, N)).*repeat(resid, outer=(1,K))
    resid = randn(T)
    @test getscore(X, resid) == X.*resid
end

@testset "hamilton_filter" begin
    df = exampledata(:hp)
    pa = 100.0.*log.(df.payrolls_a)
    c, t, e = hamilton_filter(pa, 8, 4)
    @test length(c) == length(pa)
    @test all(isnan(c[i]) for i in 1:11)
    # Compare results with Matlab code from replication files
    # Need to replace `b = inv(X'*X)*X'*y(p+h:T,1)` in `regcyc.m` with `b = X\y(p+h:T,1)`
    @test c[12] ≈ -7.828544973785029
    @test c[200] ≈ 1.462770379340782
    @test all(~e[i] for i in 1:11)
    @test sum(e) == length(c) - 11

    c1, t1 = hamilton_filter(pa, :q)
    @test isequal(c1, c)
    @test isequal(t1, t)
    c1, t1 = hamilton_filter(pa, 24, 12)
    c2, t2 = hamilton_filter(pa, :m)
    @test isequal(c2, c1)
    @test isequal(t2, t1)
    c1, t1 = hamilton_filter(pa, 2, 1)
    c2, t2 = hamilton_filter(pa, :y)
    @test isequal(c2, c1)
    @test isequal(t2, t1)

    T = length(pa)
    pa1 = copy(pa)
    pa1[2] = NaN
    c1, t1, e1 = hamilton_filter(pa1, :q)
    @test isfinite.(c1) == e1
    @test isfinite.(t1) == e1
    @test e1 == ((1:T).>13)

    c1, t1, e1 = hamilton_filter(pa1, :q, subset=1:T.<=100)
    @test isfinite.(c1) == e1
    @test isfinite.(t1) == e1
    @test e1 == ((1:T).>13).&((1:T).<=100)

    @test_throws ArgumentError hamilton_filter(pa, 0, 1)
    @test_throws ArgumentError hamilton_filter(pa, 1, 0)
    @test_throws ArgumentError hamilton_filter(pa[1:24], 8, 4)
    @test_throws ArgumentError hamilton_filter(pa, :n)
end

@testset "Cum" begin
    df = exampledata(:gk)
    T = size(df, 1)
    c = Cum(:ff4_tc)
    @test _geto(c) == :ff4_tc
    v = getcolumn(df, c)
    @test v === Cum(df.ff4_tc, nothing)
    @test size(v, 1) == size(df.ff4_tc, 1)
    @test length(v) == length(df.ff4_tc)
    @test isequal(vec(Cum(df.ff4_tc), nothing, :x, 0, Float64), coalesce.(df.ff4_tc, NaN))

    r = copy(df.ff4_tc)
    r[1:end-2] += df.ff4_tc[2:end-1] + df.ff4_tc[3:end]
    r[end-1:end] .= NaN
    v = vec(Cum(df.ff4_tc), nothing, :x, 2, Float64)
    @test v[.~isnan.(v)] ≈ r[.~isnan.(v)]
    v = vec(Cum(df.ff4_tc), (1:T.<=200).|(1:T.>220), :x, 2, Float64)
    @test v[.~isnan.(v)] ≈ r[.~isnan.(v)]
    @test isnan.(v[191:230]) == ((1:40).>8).&((1:40).<=30)
    r = copy(df.ff4_tc)
    r[3:end] += df.ff4_tc[1:end-2] + df.ff4_tc[2:end-1]
    r[1:2] .= NaN
    v = vec(Cum(df.ff4_tc), nothing, :y, 2, Float64)
    @test v[.~isnan.(v)] ≈ r[.~isnan.(v)]
    ss = (1:T.<=200).|(1:T.>220)
    v = vec(Cum(df.ff4_tc), ss, :y, 2, Float64)
    @test v[.~isnan.(v)] ≈ r[.~isnan.(v)]
    @test isnan.(v[191:230]) == ((1:40).>10).&((1:40).<=32)

    df.s = df.ff4_tc.>0
    c1 = Cum(:ff, :s)
    @test _geto(c1) == :ff
    v1 = getcolumn(df, c1)
    @test v1 === Cum(df.ff, df.s)
    @test size(v1, 1) == size(df.ff, 1)
    @test length(v1) == length(df.ff)

    @test isequal(vec(v1, nothing, :x, 0, Float64), coalesce.(df.ff.*df.s, NaN))
    @test isequal(vec(v1, nothing, :y, 0, Float64), coalesce.(df.ff, NaN))
    @test isequal(view(vec(v1, ss, :x, 0, Float64), ss), coalesce.(view(df.ff.*df.s, ss), NaN))
    @test isequal(view(vec(v1, ss, :y, 0, Float64), ss), coalesce.(view(df.ff, ss), NaN))

    @test _toint(df, c) === Cum(16)
    @test _toint(df, c1) === Cum(5, 20)
    @test _toname(df, Cum(16)) === c
    @test _toname(df, Cum(5, 20)) === c1
    @test sprint(show, c) == "Cum(ff4_tc)"
    @test sprint(show, c1) == "Cum(ff, s)"
end


text:
using Revise, Random, LinearAlgebra, Optim, StatsBase,
  Statistics, Distributions, COSMO, BenchmarkTools, Plots

Random.seed!(11)

#compute the minvol portfolio given a return target, expected returns z, and Σ
function minvolw(μ, z, Σ, Σinv = Σ\I)
  K= length(z)
  some1s = ones(K)

  A = some1s' * Σinv * some1s
  B = some1s' * Σinv * z
  C = z' * Σinv * z

  λ =(μ * B - C) /(B^2 - A*C)
  γ = (μ * A - B) /(A*C -B^2)

  w = λ * Σinv * some1s .+ γ * Σinv * z

  return w
end

#same as before but acquires the solution numerically
function nminvolw(μ, z, Σ)
  K= length(z)
  some1s = ones(K)

  m = COSMO.Model()
  c1 = COSMO.Constraint(ones(1,K), -1.0, COSMO.ZeroSet)
  c2 = COSMO.Constraint(z', -μ, COSMO.ZeroSet)
  COSMO.assemble!(m, Σ, zeros(K), [c1; c2])
  res = COSMO.optimize!(m)

  return res.x

  #return w
end

function gengraphs(;N=10^5, K=15, Δ=0.001)
  Erm = 0.08 #expected market reutn
  rf = 0.02 #risk free rate
  σm = 0.1 #market standard deviation

  ivol = rand(K)*0.4 .+ 0.1 #stock ivols
  beta = rand(Normal(1,0.2), K) #stock betas

  rm = rand(Normal(Erm, σm), N) #the market returns
  noise = reduce(hcat, (σe->rand(Normal(0, σe), N)).(ivol))
  r = rm .* beta' .+ noise

  z = mean(r, dims=1) |> vec
  Σ = cov(r)

  #sanity check
  Σinv = Σ\I
  w = minvolw(0.1, z, Σ, Σinv)
  nw = nminvolw(0.1, z, Σ)
  @assert w ≈ nw

  minvolσ(μ) = minvolw(μ, z, Σ, Σinv) |> (w)->(w'*Σ*w)^0.5
  nminvolσ(μ) = nminvolw(μ, z, Σ) |> (w)->(w'*Σ*w)^0.5

  #compute the minimum variance frontier
  μs = 0.01:Δ:0.2 |> collect
  σs = similar(μs)
  Threads.@threads for i ∈ 1:length(μs)
    σs[i] = minvolσ(μs[i])
  end

  #compute global min var portfolio
  some1s = ones(K)
  A = some1s' * Σinv * some1s
  B = some1s' * Σinv * z

  wg = Σinv * some1s / A
  σg = (wg'*Σ*wg)^(0.5)
  rg = z' * wg

  #compute the tangency portfolio
  wt = Σinv * (z .- some1s .* rf) ./ (B-rf*A)
  rt = z'*wt
  σt = (wt'*Σ*wt)^(0.5)

  pth = "C:\\Users\\Clinton\\Dropbox\\AAtawork\\ChernovInvestments\\Investments\\2020\\Exams\\final"
  p = plot(σs, μs,
    yaxis = ("E(r)", (0.0,0.16), :none),
    yticks=(0.0:0.02:0.16),
    xaxis=("vol", (0.0,0.4),),
    xticks=(0.0:0.05:0.4),
    title="Minimum Variance Frontier",
    legend=:none,
    linestyle=:dot,
    linecolor=:black)

  pans = p |> deepcopy

  plot!(pans, [0.0; σt; σs[μs .> rt]], [rf; rt; μs[μs .> rt]],
    linestyle=:dash, seriestype=:path, linecolor=:black, linewidth=4,
    annotations=(σt*1.5, rt*1.2,
    Plots.text("min var frontier\n (rf asset, no borrowing)", :left, 7,)))
  plot!(pans, [0.0; σt; σt*10], [rf; rt; (rt-rf)*10],
    linestyle=:solid, seriestype=:path, linecolor=:black, linewidth=2,
    annotations=(σt*1.1, rt*1.6,
    Plots.text("min varfrontier\n (rf asset, borrowing)", :left, 7,)))

  #scatter!(pans, [0.0], [rf], markersize=6)
  scatter!(pans, [σt], [rt], markersize=6,
    annotations=(σt+0.005, rt-0.0025, Plots.text("tangency portfolio", :left, 7,)))
  scatter!(pans, [σg], [rg],markersize=6,
    annotations=(σg+0.005, rg-0.0025, Plots.text("global min var portfolio", :left, 7,)))
  savefig(p, "$pth\\problem3.pdf")
  savefig(pans, "$pth\\problemans3.pdf")

end

@time gengraphs()

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
@testset "Projection math" begin
    @test eltype(AbstractPlotting.rotationmatrix_x(1)) == Float64
    @test eltype(AbstractPlotting.rotationmatrix_x(1f0)) == Float32
end


text:
doc"""
    rand!(a::nf_elem, U::AbstractArray) -> nf_elem
> Inplace, set the coefficients of $a$ to random elements in $U$.
> $a$ is returned.
"""
function rand!(a::nf_elem, U::AbstractArray)
  for i=1:degree(parent(a))
    _num_setcoeff!(a, i-1, rand(U))
  end
  return a
end

doc"""
    rand(K::AnticNumberField, U::AbstractArray) -> nf_elem
> Find an element in $K$ where the coefficients are selceted at random in $U$.
"""
function rand(K::AnticNumberField, U::AbstractArray)
  a = K()
  return rand!(a, U)
end

doc"""
    rand!(A::Generic.Mat{nf_elem}, U::AbstractArray) -> Generic.Mat{nf_elem}
> Inplace, replace each element in $A$ by an element where the coefficients are
> sected at random in $U$.
> Returns $A$.
"""
function rand!(A::Generic.Mat{nf_elem}, U::AbstractArray)
  for i=1:rows(A)
    for j=1:cols(A)
      rand!(A[i,j], U)
    end
  end
  return A
end

doc"""
    rand(A::Generic.MatSpace{nf_elem}, U::AbstractArray) -> Generic.Mat{nf_elem}
> Create a random matrix in $A$ where the coefficients are selected from $U$.
"""
function rand(A::Generic.MatSpace{nf_elem}, U::AbstractArray)
  return rand!(A(), U)
end

doc"""
    modular_lift(ap::Array{fq_nmod_mat, 1}, me::modular_env) -> Array
> Given an array of matrices as computed by \code{modular_proj},
> compute a global pre-image using some efficient CRT.
"""
function modular_lift(ap::Array{fq_nmod_mat, 1}, me::modular_env)
  A = zero_matrix(me.K, rows(ap[1]), cols(ap[1]))
  for i=1:rows(A)
    for j=1:cols(A)
      A[i,j] = modular_lift([ap[k][i,j] for k=1:length(ap)], me)
    end
  end
  return A
end

doc"""
    mod!(A::Generic.Mat{nf_elem}, m::fmpz)
> Inplace: reduce all entries of $A$ modulo $m$, into the positive residue system.
"""
function mod!(A::Generic.Mat{nf_elem}, m::fmpz)
  for i=1:rows(A)
    for j=1:cols(A)
      mod!(A[i,j], m)
    end
  end
end

doc"""
    mod_sym!!(A::Generic.Mat{nf_elem}, m::fmpz)
> Inplace: reduce all entries of $A$ modulo $m$, into the symmetric residue system.
"""
function mod_sym!(A::Generic.Mat{nf_elem}, m::fmpz)
  for i=1:rows(A)
    for j=1:cols(A)
      mod_sym!(A[i,j], m)
    end
  end
end

doc"""
    rational_reconstruction(A::Generic.Mat{nf_elem}, M::fmpz) -> Bool, Generic.Mat{nf_elem}
> Apply \code{rational_reconstruction} to each entry of $M$.
"""
function rational_reconstruction(A::Generic.Mat{nf_elem}, M::fmpz)
  B = similar(A)
  for i=1:rows(A)
    for j=1:cols(A)
      fl, B[i,j] = rational_reconstruction(A[i,j], M)
      if !fl
        return fl, B
      end
    end
  end
  return true, B
end

doc"""
    divexact!(A::Generic.Mat{nf_elem}, p::fmpz)
> Inplace: divide each entry by $p$.
"""
function divexact!(A::Generic.Mat{nf_elem}, p::fmpz)
  for i=1:rows(A)
    for j=1:cols(A)
      A[i,j] = A[i,j]//p
    end
  end
end

#TODO/ To experiment:
# - vector reconstruction ala Storjohan
# - reconstruction with algebraic denominators
# - vector reconstruction with algebraic denominators
# - Talk to Bill: fq_nmod_mat is missing in Nemo, the inv is dreadfully slow
# - extend to non-unique solutions
# - make Aip*D mult faster, A*y as well?
#
function Nemo.solve_dixon(A::Generic.Mat{nf_elem}, B::Generic.Mat{nf_elem})
  p = next_prime(p_start)
  K = base_ring(A)

  me = modular_init(K, p)
  ap = modular_proj(A, me)
  @time ap = [inv(x) for x= ap]
  Aip = modular_lift(ap, me)
  sol = 0*B
  D = B
  pp = fmpz(1)
  last_SOL = false
  while true
    y = Aip*D
    mod!(y, fmpz(p))
    sol += y*pp
    pp *= p
    fl, SOL = rational_reconstruction(sol, pp)
#    t = A*sol-B
#    mod!(t, pp)
#    @assert iszero(t)
    if fl
      if last_SOL == SOL && A*SOL == B
        return SOL
      else
        last_SOL = SOL
      end
    end
    D = D - A*y
    divexact!(D, fmpz(p))
#    if nbits(pp) > 10000 # a safety device to avoid infinite loops
#      error("not work")
#    end
  end
end


text:
@testset "precision" begin
    @test precision(Double64) == 113
    @test precision(Double32) == 53
    @test precision(Double32) == 24
end

@testset "maxintfloat" begin
    @test maxintfloat(Double64) >   maxintfloat(Double64) - 1
    @test maxintfloat(Double64) == (maxintfloat(Double64) - 1) + 1
    @test maxintfloat(Double32) >   maxintfloat(Double32) - 1
    @test maxintfloat(Double32) == (maxintfloat(Double32) - 1) + 1
    @test maxintfloat(Double16) >   maxintfloat(Double16) - 1
    @test maxintfloat(Double16) == (maxintfloat(Double16) - 1) + 1
end


text:
module DTCInitium

using Sockets
using AbstractDAQs

export Initium, SD1, SD2, SD3, SD5, PC4, CA2, AD0, AD2, socket
export readresponse, readresponse!
export scannerlist, SD1cmd, daqparams, setparams, SD2cmd
export addscanners
export daqaddinput, daqconfig, daqacquire, daqconfig, daqconfigdev
export daqstart, daqread, daqstop, daqzero
export daqchannels, numchannels, samplesread, isreading
export dtcsetstbl!, setfastdaq!
export devname, devtype
export DAQTask

include("ports.jl")

mutable struct DTCChannels
    "Number of channels"
    nchans::Int
    "List of port ranges"
    plst::Vector{PortRange}
    "Channels availables"
    channels::Vector{Int}
    "Names of each channel (pressure port)"
    channames::Vector{String}
    "Maps a channel name to the index of the pressure port"
    chanidx::Dict{String,Int}
end

"""
`DTCChannels()`

Stores information on available ports.

The ports available in the DTC Initium can be specified as range or
a sequence of ranges. Check command SD3 in the user's manual.

"""
DTCChannels() = DTCChannels(0, PortRange[], Int[], String[], Dict{String,Int}())

mutable struct Initium <: AbstractPressureScanner
    "IP address of the device"
    ipaddr::IPv4
    "TCP/IP port, 8400"
    port::Int32
    "Device name"
    devname::String
    "Socket used to communicate with Initium"
    sock::TCPSocket
    "Cluster/Rack/slot"
    crs::String
    "Attached scanners"
    scanners::Vector{Tuple{Int,Int,Int}}
    "Active data acquisition setup table"
    stbl::Int
    "Intermittency of temp-sets"
    actx::Int
    "DAQ Task handler - stores binary data"
    task::DAQTask
    "Buffer to store acquired data"
    buffer::CircMatBuffer{UInt8}
    "Data acquisition parameters"
    params::Dict{Symbol,Int}
    "Channels configured"
    chans::DTCChannels
    "Device configuration"
    conf::DAQConfig
    "Devices with `stbl` configured"
    stbldev::Dict{Int, Initium}
    "Have channels been added to device?"
    haschans::Bool
    "Has data acquisition been configured?"
    isconfigured::Bool
    usethread::Bool
    unit::Int
end

"""
`Initium(devname, ip; stbl=1, crs="111")`
`Initium(devname, scanners...; ip=ip stbl=1, crs="111",
    npp=64, lrn=1, bufsize, addallports)`
`Initium(devname, ip, stbl)`

Create a device that communicates with the DTC Initium.

The first constructor just establishes the communication with the
DTC Initium but has no configuration.

A more typical use is to add scanners, initialize buffers, and add all
possible pressure ports. The second constructor above does just that.

The DTC Initium can hold 5 different configurations. Usually a single
configuration is used. In some situations, more than one configuration
can be useful and a "sub-device" can be created for another DA setup table
(stbl) using the third constructor above.

## Description of the paramteters

 * `devname` String containing the device name. The device name is used to reference individual devices.
 * `ip` String with ip address of the DTC Initium
 * `stbl` Integer 1-5 that selects the DA Setup table to be used
 * `crs` Not used, exists for compatibility with othe Pressure Systems devices
 * `npp` Default number of pressure ports per scanner. Used as a default values
 * `lrn` Logical range number (see manual for further information)
 * `bufsize` Length of buffer that can store pressure data
 * `addallports` Boolean that specifies if all available pressure ports should be added to the device
 * `scanners` Scanners connected to the system. See method [`addscanners`](@ref) for detailed information on the format used. But if integers or integer ranges are used, parameters `npp` and `lrn` are used as default values.
 * `unit` Index of the unit to be used for pressure. See manual. 3 is for Pascal.
## Examples

Basic low level example:

```julia-repl
julia> dev = Initium("press", "192.168.129.7", stbl=1)
DTC Initium
    IP: 192.168.129.7
    stbl: 1
    Scanners:

julia> addscanners(dev, 1:4)

julia> print(dev)
DTC Initium
    IP: 192.168.129.7
    stbl: 1
    Scanners:
        1: npp=64, lrn=1
        2: npp=64, lrn=1
        3: npp=64, lrn=1
        4: npp=64, lrn=1

julia> numchannels(dev)
0

julia> daqaddinput(dev, 101:464)
0

julia> numchannels(dev)
256
```

```julia-repl
julia> # Let's do the same thing more conveniently...

julia> dev = Initium("press", 1:4, stbl=1, addallports=true)
DTC Initium
    IP: 192.168.129.7
    stbl: 1
    Scanners:
        1: npp=64, lrn=1
        2: npp=64, lrn=1
        3: npp=64, lrn=1
        4: npp=64, lrn=1

julia> numchannels(dev)
256
```

```julia-repl

julia> # Let's do the same thing more conveniently...

julia> dev1 = Initium("otherconfig", dev, 2)
DTC Initium
    IP: 192.168.129.7
    stbl: 2
    Scanners:
        1: npp=64, lrn=1
        2: npp=64, lrn=1
        3: npp=64, lrn=1
        4: npp=64, lrn=1

julia> numchannels(dev1)
0

julia> daqaddinput(dev1, 101:116, 201:216)
0

julia> numchannels(dev1)
32
```

Now we will create another device for anothe
"""
function Initium(devname::String, ip::String; stbl=1, crs="111", usethread=true)

    if !(1 ≤ stbl ≤ 5)
        error("DA Setup table (`stbl`) should be 1-5!")
    end

    ip1 = IPv4(ip)
    port = 8400
    sock = opensock(ip1, port)

    try
        tsk = DAQTask()

        ipars = Dict{String,Int}("stbl"=>stbl, "actx"=>1)
        fpars = Dict{String,Float64}()
        spars = Dict{String,String}()
        conf = DAQConfig(ipars, fpars, spars ;devname=devname,
                         model="DTCInitium", sn="", tag="", ip=ip)

        scanners = Tuple{Int,Int,Int}[]
        buffer = CircMatBuffer{UInt8}()
        params = Dict{Symbol,Int}()
        chans = DTCChannels()
        stbldev = Dict{Int,Initium}()

        dev = Initium(ip1, port, devname, sock, crs, scanners, stbl, 1, tsk,
                      buffer, params, chans, conf, stbldev,
                      false, false, usethread, 3)
        dev.stbldev[stbl] = dev
        return dev
    catch e
        isopen(sock) && close(sock)
        throw(e)
    end

end

function Initium(devname::String, scanners...; stbl=1,
                 ip="192.168.129.7", crs="111", npp=64, lrn=1,
                 bufsize=65_000, addallports=true, usethread=true, unit=3)

    if !(1 ≤ stbl ≤ 5)
        error("DA Setup table (`stbl`) should be 1-5!")
    end

    sock = TCPSocket()

    try
        dev = Initium(devname, ip; crs=crs, stbl=stbl,
                      usethread=usethread)
        sock = socket(dev)
        addscanners(dev, scanners...; npp=npp, lrn=lrn, unit=unit)
        # Allocate buffer
        nchans = availablechans(dev)
        w = 24 + nchans*4  # Maximum number of bytes per frame
        resize!(dev.buffer, w, bufsize)

        # Default data acquisition parameters
        SD2(dev, stbl=stbl)
        updateconf!(dev)

        if addallports
            # Add all possible pressure ports
            addallpressports(dev)
        end

        return dev
    catch e
        isopen(sock) && close(sock)
        throw(e)
    end

end

function Initium(devname::String, dev::Initium, stbl::Int)
    if !(1 ≤ stbl ≤ 5)
        error("DA Setup table (`stbl`) should be 1-5!")
    end

    if haskey(dev.stbldev, stbl)
        # This subdevice has already been created!
        # Just return it!
        return dev.stbldev[stbl]
    end
    conf = DAQConfig(devname=devname, ip=daqdevip(dev),
                     model = daqdevmodel(dev), sn=daqdevserialnum(dev),
                     tag=daqdevtag(dev))
    conf.ipars["stbl"] = stbl
    conf.ipars["actx"] = dev.actx

    newdev = Initium(dev.ipaddr, dev.port, devname, dev.sock, dev.crs,
                     dev.scanners, stbl, dev.actx, dev.task,
                     dev.buffer, Dict{Symbol,Int}(), DTCChannels(),
                     conf, dev.stbldev, false, false, dev.usethread, dev.unit)
    dev.stbldev[stbl] = newdev
    # Get a default configuration
    SD2(newdev, stbl=stbl)
    return newdev
end

function Base.show(io::IO, dev::Initium)
    println(io, "DTC Initium")
    println(io, "    Dev Name: $(devname(dev))")
    println(io, "    IP: $(string(dev.ipaddr))")
    println(io, "    stbl: $(dev.stbl)")
    println(io, "    Scanners:")
    for s in dev.scanners
        snum = s[1]
        npp = s[2]
        lrn = s[3]
        println(io, "        $snum: npp=$npp, lrn=$lrn")
    end

end

"""
`updateconf!(dev)`

The configuration of a `DTCInitium` device is stored in the `params` field of
an `Initium` object and other fields (such as `actx` or `stbl`).
This function updates the `conf::DAQConfig` field so that
the actual configuration can be stored in HDF5 files.
"""
function updateconf!(dev)
    stbl = dev.stbl

    p1 = dev.params
    p = dev.conf
    ipars = p.ipars
    ipars["stbl"] = stbl
    ipars["nfr"] = p1[:nfr]
    ipars["nms"] = p1[:nms]
    ipars["msd"] = p1[:msd]
    ipars["trm"] = p1[:trm]
    ipars["scm"] = p1[:scm]
    ipars["ocf"] = p1[:ocf]
    ipars["actx"] = dev.actx
    return

end

"""
`setfastdaq!(dev, actx=1)`

When measuring pressure with the DTC Initium, individual temperature of the
pressure sensors are used to correct the pressure. In the default configuration,
each time a pressure is read, so is the temperature. While this can be more accurate
it has a cost: it takes time to read each temperature.

To achieve a higher sampling rate, the temperature of each pressure port doesn't have
to be read for every pressure sample. If `actx=0`, it is read only at the beginning.
When `actx > 1`, this parameter will acquire the temperature intermittently
(once every  `actx` pressure samples are read).
"""
function setfastdaq!(dev::Initium, actx=1)

    SD5(dev, actx)
    # This is a global parameter! Should update ALL subdevices!
    for (stbl, xdev) in dev.stbldev
        xdev.actx = actx
        xdev.conf.ipars["actx"] = actx
    end
    return
end

"""
`addpressports(dev, plst; names="P")`

Add pressure port defined by `plst`. See [`PortRange`](@ref) and
[`portlist`](@ref) to see how `plst` is defined.

To give specific names to the pressure ports, use argument `names`:
 * String or symbol: The string will be preppended to the channel number
 * Vector: specific names to each pressure channel.

"""
function addpressports(dev, plst::AbstractVector{PortRange}; names="P")

    chans = defscanlist(scanners(dev), plst)
    # Check if there are ports in plst previously added:
    if length(intersect(chans, dev.chans.channels)) > 0
        throw(ArgumentError("Some ports have already been added!"))
    end

    SD3(dev, dev.stbl, plst)

    # Get channels names and numbers
    if isa(names, AbstractString) || isa(names, Symbol)
        chs = string.(names, chans)
    elseif isa(names, AbstractVector)
        length(chans) != length(names) && throw(ArgumentError("If `names` is a vector it should have the length of the number of channels"))
        chs = string.(names)
    else
        throw(ArgumentError("`names` should be either a string, a symbol or a vector"))
    end

    # Add the new channels to dev.chans
    nch = length(chans)
    dev.chans.nchans = length(chans)
    dev.chans.plst =  plst
    dev.chans.channels = chans
    dev.chans.channames = chs
    for i in 1:nch
        chname = chs[i]
        dev.chans.chanidx[chname] = i
    end

    dev.haschans = true

end

"""
`addallpressports(dev)`

Use every available pressure port when acquiring data. This is a convenience function
that can be used instead of [`daqaddinput`](@ref)  if every available pressure port
should be acquired.
"""
function addallpressports(dev)

    plst = PortRange[]

    for (s,n,l) in dev.scanners
        p1 = s*100 + 1
        p2 = s*100 + n
        push!(plst, PortRange(p1,p2,true))
    end
    addpressports(dev, plst; names="P")
end

"Total number of available channels in the scanners"
availablechans(dev::Initium) = availablechans(dev.scanners)

Base.open(dev::Initium) = dev.sock = opensock(ipaddr(dev), portnum(dev))
Base.close(dev::Initium) = close(dev.sock)

"""
`addscanners(dev::Initium, lst...; npp=64, lrn=1)`

Add scanners to the DTC Initium device. Before any channels can be added,
the scanners should be added. This function calls [`scannerlist`](@ref)
 to create a list of scanners attached to the DTC Initium.

 * `dev::Initium` Device to which the scanners should be added
 * `lst...` Sequence of scanners to be added (see details and [`scannerlist`](@ref))
 * `npp` Default number of pressure ports per scanner
 * `lrn` Default logical range number.

Each scanner attached to the DTC Initium is specified by the slot it is connected to
(1-8), the number of pressure ports and the logical range number.

Often, similar scanners (same number of pressure ports) are attached and the default
values of number of pressure ports can be specified by argument `npp`. The same goes for
argument `lrn`

## Examples

```julia-repl
julia> using DTCInitium

julia> dev = Initium("press", "192.168.129.7", stbl=1)
DTC Initium
    IP: 192.168.129.7
    stbl: 1
    Scanners:

julia> addscanners(dev, 1:2, (3,64,1), (4,64,1))

julia> print(dev)
DTC Initium
    IP: 192.168.129.7
    stbl: 1
    Scanners:
        1: npp=64, lrn=1
        2: npp=64, lrn=1
        3: npp=64, lrn=1
        4: npp=64, lrn=1

julia> close(dev)

julia> dev = Initium("press", "192.168.129.7", stbl=1)
DTC Initium
    IP: 192.168.129.7
    stbl: 1
    Scanners:

julia> addscanners(dev, 1:4)

julia> print(dev)
DTC Initium
    IP: 192.168.129.7
    stbl: 1
    Scanners:
        1: npp=64, lrn=1
        2: npp=64, lrn=1
        3: npp=64, lrn=1
        4: npp=64, lrn=1
```

"""
function addscanners(dev::Initium, lst...; npp=64, lrn=1, unit=3)

    scnlst = scannerlist(lst...; npp=npp, lrn=lrn)

    !isopen(socket(dev)) && open(dev)
    try
        SD1(dev, scnlst)
        dev.scanners = scnlst
        # Set buffer width size
        nchans = availablechans(scnlst)
        w = 24 + nchans*4  # Maximum number of bytes per frame
        dev.buffer.width = w
        if unit < 1 || unit > 12
            unit = 3 # Pa
        end

        # Set the default unit to Pascal (3)

        for lrn in unique([s[3] for s in scnlst])
            PC4(dev, unit, 0, lrn=lrn)
        end
        dev.unit=unit

    catch e
        if isa(e, DTCInitiumError)
            close(dev.sock)
            throw(e)
        else
            throw(e)
        end
    end

    return
end

"Return crs parameter"
getcrs(dev::Initium) = dev.crs

"Return list of scanners attached to the DTC Initium"
scanners(dev::Initium) = dev.scanners

"Return the socket used to communicate with the DTC Initium"
socket(dev::Initium) = dev.sock

"Return configuration parameters"
daqparam(dev::Initium) = dev.params

"Return IPv4 address"
ipaddr(dev::Initium) = dev.ipaddr

"Return TCP/IP port number used by the DTC Initium"
portnum(dev::Initium) = dev.port

include("errorcodes.jl")
include("packets.jl")
include("buildcmd.jl")
include("commands.jl")
include("daq.jl")

end


text:
#!/usr/bin/env julia
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Base: thispatch, thisminor, nextpatch, nextminor
import LinearAlgebra: checksquare
import UUIDs
import Pkg
using Pkg.Types
using Pkg.Types: uuid_package, uuid_registry, uuid5
import Pkg.Pkg2.Reqs: Reqs, Requirement
import Pkg.Pkg2.Pkg2Types: VersionInterval, VersionSet

## Loading data into various data structures ##

struct Require
    versions::VersionSet
    systems::Vector{Symbol}
end

struct Version
    sha1::String
    requires::Dict{String,Require}
end

struct Package
    uuid::UUID
    url::String
    versions::Dict{VersionNumber,Version}
end

Require(versions::VersionSet) = Require(versions, Symbol[])
Require(version::VersionInterval) = Require(VersionSet([version]), Symbol[])
Version(sha1::AbstractString) = Version(String(sha1), Dict{String,Require}())

function load_requires(path::String)
    requires = Dict{String,Require}()
    requires["julia"] = Require(VersionInterval())
    isfile(path) || return requires
    for r in filter!(r->r isa Requirement, Reqs.read(path))
        new = haskey(requires, r.package)
        versions, systems = VersionSet(r.versions.intervals), r.system
        if haskey(requires, r.package)
            versions = versions ∩ requires[r.package].versions
            systems  = systems  ∪ requires[r.package].systems
        end
        requires[r.package] = Require(versions, Symbol.(systems))
    end
    return requires
end

function load_versions(dir::String)
    versions = Dict{VersionNumber,Version}()
    isdir(dir) || return versions
    for ver in readdir(dir)
        path = joinpath(dir, ver)
        sha1 = joinpath(path, "sha1")
        isfile(sha1) || continue
        requires = load_requires(joinpath(path, "requires"))
        versions[VersionNumber(ver)] = Version(readchomp(sha1), requires)
    end
    return versions
end

const JULIA_VERSIONS = Dict(
    v"0.1.1"  => Version("bddb70374fdbb2256782398c13430fa4d8b55d0b"),
    v"0.1.2"  => Version("1187040d027210cc466f0b6a6d54118fd692cf2d"),
    v"0.2.0"  => Version("05c6461b55d9a66f05ead24926f5ee062b920d6b"),
    v"0.2.1"  => Version("e44b5939057d87c1e854077108a1a6d66203f4fa"),
    v"0.3.0"  => Version("768187890c2709bf2ff06818f40e1cdc79bd44b0"),
    v"0.3.1"  => Version("c03f413bbdb46c00033f4eaad402995cfe3b7be5"),
    v"0.3.2"  => Version("8227746b95146c2921f83d2ae5f37ecd146592d8"),
    v"0.3.3"  => Version("b24213b893ace6ca601722f6cdaac15f2b034b7b"),
    v"0.3.4"  => Version("33920266908e950936ef3f8503810a2d47333741"),
    v"0.3.5"  => Version("a05f87b79ad62beb033817fdfdefa270c9557aaf"),
    v"0.3.6"  => Version("0c24dca65c031820b91721139f0291068086955c"),
    v"0.3.7"  => Version("cb9bcae93a32b42cec02585c387396ff11836aed"),
    v"0.3.8"  => Version("79599ada444fcc63cb08aed64b4ff6a415bb4d29"),
    v"0.3.9"  => Version("31efe690bea1fd8f4e44692e205fb72d34f50ad3"),
    v"0.3.10" => Version("c8ceeefcc1dc25953a644622a895a3adcbc80dad"),
    v"0.3.11" => Version("483dbf5279efe29f94cf35d42de0ab0a673fb34e"),
    v"0.3.12" => Version("80aa77986ec84067f4696439f10fc77c93703bc7"),
    v"0.4.0"  => Version("0ff703b40afddf9b705bd6a06d3a59cb4c089ea5"),
    v"0.4.1"  => Version("cbe1bee3a8f27ba4f349556857bf4615ee5fa68f"),
    v"0.4.2"  => Version("bb73f3489d837e3339fce2c1aab283d3b2e97a4c"),
    v"0.4.3"  => Version("a2f713dea5ac6320d8dcf2835ac4a37ea751af05"),
    v"0.4.4"  => Version("ae683afcc4cc824dcf68251569a6cf8760e4fb6f"),
    v"0.4.5"  => Version("2ac304dfba75fad148d4070ef4f8a2e400c305bb"),
    v"0.4.6"  => Version("2e358ce975029ec97aba5994c17d4a2169c3b085"),
    v"0.4.7"  => Version("ae26b25d43317d7dd3ca05f60b70677aab9c0e08"),
    v"0.5.0"  => Version("3c9d75391c72d7c32eea75ff187ce77b2d5effc8"),
    v"0.5.1"  => Version("6445c82d0060dbe82b88436f0f4371a4ee64d918"),
    v"0.5.2"  => Version("f4c6c9d4bbbd9587d84494e314f692c15ff1f9c0"),
    v"0.6.0"  => Version("903644385b91ed8d95e5e3a5716c089dd1f1b08a"),
    v"0.6.1"  => Version("0d7248e2ff65bd6886ba3f003bf5aeab929edab5"),
    v"0.6.2"  => Version("d386e40c17d43b79fc89d3e579fc04547241787c"),
    v"0.7.0"  => Version("a4cb80f3edcf8cea00bd9660e3b65f544f41462f"),
    v"1.0.0"  => Version("5d4eaca0c9fa3d555c79dbacdccb9169fdf64b65"),
    v"1.0.1"  => Version("0d713926f85dfa3e4e0962215b909b8e47e94f48"),
    v"1.0.2"  => Version("d789231e9985537686052db9b2314c0d51656308"),
    v"1.1.0"  => Version("6ef1d76c50a39c0ce68b4e42948a1499c1551415"), # dummy commit, not actual 1.1
    v"1.2.0"  => Version("217d330296debe0567bb07addabf66b00602e325"), # dummy commit, not actual 1.2
)

function load_packages(dir::String)
    pkgs = Dict{String,Package}()
    for pkg in readdir(dir)
        path = joinpath(dir, pkg)
        url = joinpath(path, "url")
        versions = joinpath(path, "versions")
        isfile(url) || continue
        pkgs[pkg] = Package(
            uuid5(uuid_package, pkg),
            readchomp(url),
            load_versions(versions),
        )
    end
    pkgs["julia"] = Package(
        uuid5(uuid_package, "julia"),
        "https://github.com/JuliaLang/julia.git",
        JULIA_VERSIONS,
    )
    return pkgs
end

@eval julia_versions() = $(sort!(collect(keys(JULIA_VERSIONS))))
julia_versions(f::Function) = filter(f, julia_versions())
julia_versions(vi::VersionInterval) = julia_versions(v->v in vi)

macro clean(ex) :(x = $(esc(ex)); $(esc(:clean)) &= x; x) end

function prune!(pkgs::AbstractDict{String,Package})
    # remove unsatisfiable versions
    while true
        clean = true
        filter!(pkgs) do (pkg, p)
            filter!(p.versions) do (ver, v)
                @clean ver == thispatch(ver) > v"0.0.0" &&
                all(v.requires) do kv
                    req, r = kv
                    haskey(pkgs, req) &&
                    any(w->w in r.versions, keys(pkgs[req].versions))
                end
            end
            @clean !isempty(p.versions)
        end
        clean && break
    end
    return pkgs
end

## Load package data ##

const pkgs = load_packages(Pkg.Pkg2.dir("METADATA"))

# delete packages whose repos that no longer exist:
for pkg in [
    "CardinalDicts"
    "CreateMacrosFrom"
    "GSDicts"
    "S3Dicts"
    "ChainMap"
    "NanoTimes"
    "OnlineMoments"
    "ProjectiveDictionaryPairLearning"
    "React"
    "PackageGenerator"
    "Sparrow"
    "KeyedTables"
    "Arduino"
    "ControlCore"
    "DynamicalBilliardsPlotting"
    "GLUT"
    "GetC"
]

    delete!(pkgs, pkg)
end

# cap julia version for probably-0.7-incompatible packages

const passing = readlines(joinpath(@__DIR__, "passing.txt"))
for pkg in passing
    haskey(pkgs, pkg) || continue
    p = pkgs[pkg]
    ver = maximum(keys(p.versions))
    v = p.versions[ver]
    for dep in keys(v.requires)
        dep == "julia" && continue
        dep in passing || push!(passing, dep)
    end
end
sort!(passing, by=lowercase)

const meta_dir = Pkg.Pkg2.dir("METADATA")
const time_map = Dict{Tuple{String,VersionNumber},Int}()
let t = 0
    for line in eachline(`git -C $meta_dir log --format=%ct --name-only`)
        if (m = match(r"^(\d+)$", line)) !== nothing
            t = parse(Int, line)
        elseif (m = match(r"^([^/]+)/versions/(\d+\.\d+\.\d+)/requires$", line)) !== nothing
            pkg = String(m.captures[1])
            ver = VersionNumber(m.captures[2])
            haskey(time_map, (pkg, ver)) && continue
            time_map[pkg,ver] = t
        end
    end
end

const all_vers = julia_versions()
const old_vers = julia_versions(v -> v < v"0.7")
const jul_14 = 1531526400 # Jul 14, 2018
const oct_10 = 1539129600 # Oct 10, 2018

function cap_compat!(pkg::String, ver::VersionNumber, reqs::Dict{String,Require})
    jvers = reqs["julia"].versions
    ivals = jvers.intervals
    isempty(ivals) && return
    t = get(time_map, (pkg, ver), 0)
    if pkg in passing && ver ≥ maximum(keys(pkgs[pkg].versions)) || !isempty(ivals[end]) &&
        (ivals[end].upper < v"∞" || !any(v->v in ivals[end] && v < v"0.7", all_vers))
        # in the "passing list" from pkgeval and maxiumal version => leave alone
        # has final interval with explicit upper bound => leave alone
        # or interval only containing 0.7+ versions => 1.0 compatible
        return # no change
    elseif (v"0.7" in jvers || v"1.0" in jvers) && t ≥ oct_10
        # recently tagged and claims to support 0.7/1.0 => face value
        return # no change
    elseif pkg != "Compat" && !haskey(reqs, "Compat") && any(v in jvers for v in old_vers)
        # supports an older julia & doesn't use Compat => 0.7 incompatible
        # fall through
    elseif (v"0.7" in jvers || v"1.0" in jvers) && t ≥ jul_14
        # claims to support 0.7+ & tagged after date cutoff => 1.0 compatible
        return # no change
    end
    # cap supported Julia versions at 0.7
    ivals[end] = VersionInterval(ivals[end].lower, v"0.7+")
    return
end

for (pkg, p) in pkgs
    (pkg == "julia" || pkg == "Example") && continue
    for (ver, v) in p.versions
        cap_compat!(pkg, ver, v.requires)
    end
end

# prune versions that can't be satisfied
prune!(pkgs)


text:
####################################################
gaugexy2mat(gauges::Vector{VisClaw.Gauge}) = permutedims(hcat(getfield.(gauges, :loc)...), [2,1])

####################################################
"""
    gmtgaugelocation(gauge::VisClaw.Gauge; kwargs...)
    gmtgaugelocation!(gauge::VisClaw.Gauge; kwargs...)
    gmtgaugelocation(gauges::Vector{VisClaw.Gauge}; kwargs...)
    gmtgaugelocation!(gauges::Vector{VisClaw.Gauge}; kwargs...)
"""
gmtgaugelocation(gauge::VisClaw.Gauge; kwargs...) = GMT.scatter([gauge.loc[1] gauge.loc[2]]; kwargs...)
####################################################
"""
$(@doc gmtgaugelocation)
"""
gmtgaugelocation!(gauge::VisClaw.Gauge; kwargs...) = GMT.scatter!([gauge.loc[1] gauge.loc[2]]; kwargs...)
####################################################
gmtgaugelocation(gauges::Vector{VisClaw.Gauge}; kwargs...) = GMT.scatter(gaugexy2mat(gauges); kwargs...)
gmtgaugelocation!(gauges::Vector{VisClaw.Gauge}; kwargs...) = GMT.scatter!(gaugexy2mat(gauges); kwargs...)

####################################################
"""
    gmtgaugeannotation!(gauge::VisClaw.Gauge; kwargs...)
"""
gmtgaugeannotation!(gauge::VisClaw.Gauge, annot::AbstractString=gauge.label; R="", offset=(0.0,0.0), kwargs...) =
GMT.text!(GMT.text_record([gauge.loc[1]+offset[1] gauge.loc[2]+offset[2]], annot); R=R, kwargs...)
####################################################

####################################################
"""
    gmtgaugewaveform(gauge::VisClaw.Gauge; kwargs...)
    gmtgaugewaveform!(gauge::VisClaw.Gauge; kwargs...)
"""
gmtgaugewaveform(gauge::VisClaw.Gauge; kwargs...) = GMT.plot(gauge.time, gauge.eta; kwargs...)
####################################################
"""
$(@doc gmtgaugewaveform)
"""
gmtgaugewaveform!(gauge::VisClaw.Gauge; kwargs...) = GMT.plot!(gauge.time, gauge.eta; kwargs...)
####################################################
"""
    gmtgaugevelocity(gauge::VisClaw.Gauge; kwargs...)
    gmtgaugevelocity!(gauge::VisClaw.Gauge; kwargs...)
"""
gmtgaugevelocity(gauge::VisClaw.Gauge; kwargs...) = GMT.plot(gauge.time, sqrt.(gauge.u.^2 + gauge.v.^2); kwargs...)
####################################################
"""
$(@doc gmtgaugevelocity)
"""
gmtgaugevelocity!(gauge::VisClaw.Gauge; kwargs...) = GMT.plot!(gauge.time, sqrt.(gauge.u.^2 + gauge.v.^2); kwargs...)
####################################################


text:
# dqmc.jl called with arguments: whatever.in.xml
using Dates, LinearAlgebra, Pkg
start_time = now()
println("\nStarted: ", Dates.format(start_time, "d.u yyyy HH:MM"))
println("Hostname: ", gethostname())

# support for symbol link to app/dqmc.jl in root of repo
__FILE__ = joinpath(dirname(@__FILE__), splitpath(dirname(@__FILE__))[end] == "app" ? "" : "app")

function is_dqmc_env_activated()
    project_file = Base.active_project()
    project = Pkg.Types.read_project(project_file)
    return !isnothing(project.name) && lowercase(project.name) == "qmc"
end

if !is_dqmc_env_activated()
  println("Activating DQMC environment.")
  haskey(ENV, "JULIA_DQMC") || error("DQMC environment not loaded and JULIA_DQMC env variable not set!")
  Pkg.activate(ENV["JULIA_DQMC"])
  is_dqmc_env_activated() || error("Invalid JULIA_DQMC env variable value.")
end

try
  nthreads = parse(Int, ENV["OMP_NUM_THREADS"])
  LinearAlgebra.BLAS.set_num_threads(nthreads)
  ENV["MKL_NUM_THREADS"] = nthreads
  ENV["JULIA_NUM_THREADS"] = 1
  println("Using $(nthreads) threads for BLAS, OMP, and MKL.")
catch err
  # Single core simulation
  LinearAlgebra.BLAS.set_num_threads(1)
  ENV["OMP_NUM_THREADS"] = 1
  ENV["MKL_NUM_THREADS"] = 1
  ENV["JULIA_NUM_THREADS"] = 1
  println("Using single thread for BLAS, OMP, and MKL.")
end

# -------------------------------------------------------
#             Input/Output file preparation
# -------------------------------------------------------
if length(ARGS) == 1
  # ARGS = ["whatever.in.xml"]
  input_xml = ARGS[1]
  output_file = input_xml[1:first(findfirst(".in.xml", input_xml))-1]*".out.h5.running"
elseif length(ARGS) == 2
  # Backward compatibility
  # ARGS = ["sdwO3_L_4_B_2_dt_0.1_2", 1]
  prefix = convert(String, ARGS[1])
  idx = 1
  try global idx = parse(Int, ARGS[2]); catch end # SLURM_ARRAY_TASK_ID
  output_file = prefix * ".task" * string(idx) * ".out.h5.running"

  println("Prefix is ", prefix, " and idx is ", idx)
  input_xml = prefix * ".task" * string(idx) * ".in.xml"
else
  input_xml = "dqmc.in.xml"
  output_file = input_xml[1:first(findfirst(".in.xml", input_xml))-1]*".out.h5.running"
end
@assert isfile(input_xml)

# hdf5 write test/ dump git commit
using Git
branch = Git.branch(dir=dirname(@__FILE__)).string[1:end-1]
if branch != "master"
  println("!!!Not on branch master but \"$(branch)\"!!!")
  flush(stdout)
end

# TIMING parameter "hack"
using LightXML #, Iterators
include(joinpath(__FILE__, "../src/tools/xml_parameters.jl"))
params = xml2dict(input_xml, false)
haskey(params, "TIMING") && (parse(Bool, lowercase(params["TIMING"])) == true) && (global const TIMING = true)

# -------------------------------------------------------
#                    Includes
# -------------------------------------------------------
include(joinpath(__FILE__, "../src/dqmc_framework.jl"))
using JLD

# -------------------------------------------------------
#      Parse input parameters + check for resuming
# -------------------------------------------------------
p = Params()
p.output_file = output_file
xml2parameters!(p, input_xml)
set_walltimelimit!(p, start_time)

# mv old .out.h5 to .out.h5.running (and then try to resume below)
# (isfile(output_file[1:end-8]) && !isfile(output_file)) && mv(output_file[1:end-8], output_file)
isfile(output_file[1:end-8]) && begin println("Nothing to do here. Found .out.h5 file."); exit(); end

# check if there is a resumable running file
if isfile(output_file)
  println("\nFound old output file. Let's see if we can resume.")
  alreadydone = false
  resumable = false
  try
    jldopen(output_file) do f
      if HDF5.has(f.plain, "resume") && HDF5.has(f.plain, "obs/configurations/count") && read(f["GIT_BRANCH_DQMC"]) == branch
        nconfs = read(f["obs/configurations/count"])
        pmeasurements = read(f["params/measurements"])
        pwrite_every_nth = read(f["params/write_every_nth"])
        measurements = nconfs * pwrite_every_nth

        (measurements >= pmeasurements) && (measurements >= p.measurements) && (global alreadydone = true)
        (nconfs > 0) && (global resumable = true)
      end
    end
  catch
  end

  if alreadydone
    mv(output_file, output_file[1:end-8])
    println("Nothing to do here. There are already enough configurations present.")
    exit();
  else
    if resumable
      println("Will resume.")
      p.resume = true
    else
      println("Couldn't resume.")
    end
  end
end

if !p.resume

  # Start with a prethermalized conf?
  if isfile(p.output_file)
    try
      jldopen(output_file) do f
        if HDF5.has(f.plain, "thermal_init")
          println("Using thermal_init/conf as starting configuration.")
          global start_conf = read(f["thermal_init/conf"])

          if HDF5.has(f.plain, "thermal_init/prethermalized")
            p.prethermalized = read(f["thermal_init/prethermalized"])
            println("Using thermal_init/prethermalized. Will do $(p.thermalization - p.prethermalized) ud-sweeps.")
          end

          if HDF5.has(f.plain, "thermal_init/box")
            box = read(f, "thermal_init/box")
            p.box = box
            println("Using thermal_init/box.")
          end

          if HDF5.has(f.plain, "thermal_init/box_global")
            box_global = read(f, "thermal_init/box_global")
            p.box_global = box_global
            println("Using thermal_init/box_global.")
          end
        end
      end
    catch
    end
  end

  # overwrite (potential) running file
  jldopen(output_file, isfile(p.output_file) ? "r+" : "w") do f
    HDF5.has(f.plain, "GIT_COMMIT_DQMC") && o_delete(f.plain, "GIT_COMMIT_DQMC")
    HDF5.has(f.plain, "GIT_BRANCH_DQMC") && o_delete(f.plain, "GIT_BRANCH_DQMC")
    HDF5.has(f.plain, "RUNNING") && o_delete(f.plain, "RUNNING")
    f["GIT_COMMIT_DQMC"] = Git.head(dir=dirname(@__FILE__)).string
    f["GIT_BRANCH_DQMC"] = branch
    f["RUNNING"] = 1

    # we will dump the parameters below
    HDF5.has(f.plain, "params") && o_delete(f.plain, "params")
    HDF5.has(f.plain, "obs") && o_delete(f.plain, "obs")
  end
end

if !p.resume
  parameters2hdf5(p, p.output_file)
else
  println()
  println("RESUMING MODE -----------------")
  h5open(output_file, "r+") do f
    !HDF5.has(f, "RUNNING") || HDF5.o_delete(f, "RUNNING")
    f["RUNNING"] = 1
    !HDF5.has(f, "RESUME") || HDF5.o_delete(f, "RESUME")
    f["RESUME"] = 1
  end

  confs = loadobs_frommemory(output_file, "obs/configurations")
  prevcount = length(confs)
  println("Found $(prevcount) configurations.")
  println()
  lastmeasurements = h5read(output_file, "params/measurements")
  global prevmeasurements = prevcount * p.write_every_nth
  if prevmeasurements < lastmeasurements
    println("Finishing last run (i.e. overriding p.measurements)")
    p.measurements = lastmeasurements - prevmeasurements
  elseif prevmeasurements < p.measurements
    println("Extension run.")
    p.measurements = p.measurements - prevmeasurements
  end
  @show prevmeasurements
  @show p.measurements
  global lastconf = confs[end]
end

# -------------------------------------------------------
#                      Simulation
# -------------------------------------------------------
mc = DQMC(p)

println()
println("HoppingEltype = ", heltype(mc))
println("GreensEltype = ", geltype(mc))
println()
@printf("It took %.2f minutes to prepare everything. \n", (now() - start_time).value/1000/60)

if !mc.p.resume
  if @isdefined start_conf
    init!(mc, start_conf)
  else
    init!(mc)
  end
  run!(mc)
else
  resume!(mc, lastconf, prevmeasurements)
end

h5open(output_file, "r+") do f
  HDF5.o_delete(f, "RUNNING")
  f["RUNNING"] = 0
end
mv(output_file, output_file[1:end-8], force=true) # .out.h5.running -> .out.h5

end_time = now()
println("Ended: ", Dates.format(end_time, "d.u yyyy HH:MM"))
@printf("Duration: %.2f minutes", (end_time - start_time).value/1000/60)

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
module aoc_01
    using Test
    export read_input, part1, part2
    export test_part1, test_part2, test
    export input, demo

    function read_input(filename = "input.txt")
        lines = readlines(filename)
        # for line in eachline(filename)
        [parse(Int,line) for line in lines]
    end

    const input = read_input("input.txt")
    const demo = read_input("demo.txt")

    function part1(input=input)
        for i in eachindex(input), j in eachindex(input[i+1:end])
            if input[i] + input[j] == 2020
                return input[i]*input[j]
            end
        end
        error("Could not find two numbers which add up to 2020")
    end

    function test_part1()
    end

    function part2()
        for i in eachindex(input),
            j in eachindex(input[i+1:end]),
            k in eachindex(input[i+2:end])
            if input[i] + input[j] + input[k] == 2020
                return input[i]*input[j]*input[k]
            end
        end

    end

    function test_part2()
    end
end


text:
module ApplyNoise

using Random
using ImageCore
using ImageCore.MappedArrays
using ImageCore: NumberLike, GenericGrayImage, GenericImage

import ..NoiseAPI: AbstractImageNoise, apply_noise, apply_noise!
import Base: show

# avoid InexactError for Bool array input by promoting it to float
apply_noise(img::AbstractArray{T},
            n::AbstractImageNoise,
            args...; kargs...) where T<:Union{Bool, Gray{Bool}} =
    apply_noise(of_eltype(floattype(eltype(img)), img), n, args...; kargs...)

include("AdditiveWhiteGaussianNoise.jl")

export
    apply_noise, apply_noise!,

    AdditiveWhiteGaussianNoise

end # end module


text:
type BanditStats
  regret::Float64
  arm_counts::AbstractVector{Float64}
end
+(a::BanditStats, b::BanditStats) = BanditStats(a.regret + b.regret, a.arm_counts + b.arm_counts)
/(a::BanditStats, n::Integer) = BanditStats(a.regret / n, a.arm_counts / n)

function simulate(bandit::Bandit, agent::Agent, periods::Integer)
  total_regret = 0
  while agent.period < periods
    action = choose(agent)
    reward = pull(bandit, action)
    total_regret += regret(bandit, action)
    observe(agent, reward)
  end
  BanditStats(total_regret, agent.action_attempts)
end

function aggregate_simulate(bandit::Bandit, agent::Agent, periods::Integer, n::Integer)
  qs = @parallel (+) for i=1:n
    simulate(bandit, agent, periods)
  end
  return qs / n
end


text:
# using Revise
using Profile
import CSV
import DataFrames
using DecisionTrees

include(joinpath("..", "test\\runtests.jl"))

datafile = string("data\\freMTPL2\\freMTPL2.csv")
@assert isfile(datafile);
@time fullData = CSV.read(datafile, rows_for_type_detect=100000, allowmissing=:none, categorical=false);

areasSorted = sort(unique(fullData[:Area]))
AreaInteger = map(x->findall((in)([x]), areasSorted)[1], fullData[:Area])
fullData[:AreaInteger] = AreaInteger

# correct for unreasonable observations
for i = 1:size(fullData, 1)
    fullData[:ClaimNb][i] = min(4, fullData[:ClaimNb][i])
    fullData[:Exposure][i] = min(1.0, fullData[:Exposure][i])
end

# set independent variables
selected_explanatory_vars = ["Area","AreaInteger","VehPower","VehAge","DrivAge","BonusMalus","VehBrand","VehGas","Density","Region"]

dtmtable, sett, dfprepped = prepare_dataframe_for_dtm!(fullData, trnvalcol="trnTest", numcol="ClaimNb", denomcol="Exposure", weightcol="Exposure", independent_vars=selected_explanatory_vars);

updateSettingsMod!(sett,minWeight=-0.03,model_type="boosted_tree",iterations=40,learningRate=0.025,subsampling_features_prop=.7,calculatePoissonError=true)

resultingFiles, resM = dtm(dtmtable, sett) # 31s on notebook (x260 with profiling disabled but Revise on!)

Profile.clear()
@profile dtm(dtmtable, sett)

thisIO = open("c:\\temp\\prof.txt", "w")
Profile.print(IOContext(thisIO, :displaysize => (24, 500)),format=:flat)
# Profile.print(IOContext(s, :displaysize => (24, 500)))
close(thisIO);

# old stuff?
li, lidict = Profile.retrieve()

profile_loc = "R:\\temp\\profile.jlprof"
dd = splitdir(profile_loc)[1]
if isdir(dd)
    @warn("This is currently failing with JLD2. I do no want to install JLD as it causes some compatibility issues (for now)")
#    JLD2.@save "mi.jlddd" li lidict
else
    @warn("Folder $(dd) not found!")
end

# consider
Profile.print(noisefloor=2.0,sortedby=:count)

# open("C:\\temp\\logs\\profile.txt", "w") do s
#    Profile.print(IOContext(s, :displaysize => (24, 500)))
# end

#=
#to open and read the data
using HDF5, JLD2, ProfileView
JLD2.@load "/tmp/profdata.jld"
ProfileView.view(li, lidict=lidict) =#


text:
struct MatchRule{T} end

# default transform is to do nothing
transform(fn::Function, value) = value
transform(fn::Function, vector::Vector) = [transform(fn,n) for n in nodes]

fndefault(node,children,label) = Node(node.name,node.value,node.first,node.last,children,node.ruleType)

function transform(fn::Function, node::Node)
  transformedchildren = [transform(fn, child) for child in node.children]

  if hasmethod(fn, (Node, Any, MatchRule{Symbol(node.name)}))
    label = MatchRule{Symbol(node.name)}()
  elseif hasmethod(fn, (Node, Any, MatchRule{:default}))
    label = MatchRule{:default}()
  else
    return fndefault(node,transformedchildren,nothing)
  end

  return fn(node, transformedchildren, label)
end


text:
module Test13

import Automa
import Automa.RegExp: @re_str
const re = Automa.RegExp
using Test

@testset "Test13" begin
    abra = re"abra"
    ca = re"(ca)+"
    ca.actions[:enter] = [:ca_enter]
    ca.actions[:exit] = [:ca_exit]
    dabra = re"dabra"
    machine = Automa.compile(re.cat(abra, ca, dabra))
    ctx = Automa.CodeGenContext(generator=:inline, clean=true)
    @eval function validate(data)
        logger = Symbol[]
        $(Automa.generate_init_code(ctx, machine))
        p_end = p_eof = sizeof(data)
        $(Automa.generate_exec_code(ctx, machine, :debug))
        return logger, cs == 0 ? :ok : cs < 0 ? :error : :incomplete
    end
    @test validate(b"a") == ([], :incomplete)
    @test validate(b"abrac") == ([:ca_enter], :incomplete)
    @test validate(b"abraca") == ([:ca_enter], :incomplete)
    @test validate(b"abracad") == ([:ca_enter, :ca_exit], :incomplete)
    @test validate(b"abracadabra") == ([:ca_enter, :ca_exit], :ok)
    @test validate(b"abracacadabra") == ([:ca_enter, :ca_exit], :ok)
    @test validate(b"abrad") == ([], :error)
end

end


text:
#

using LinearAlgebra
using Gadfly
Gadfly.set_default_plot_size(35cm,20cm)

p1 = 0.7 + 0.7*im

b1 = 0.4 + 0.2*im
b2 = 0.3 + 0.9*im
b3 = 0.8 + 0.1*im

r1 = 100*norm(b1-p1)
r2 = 100*norm(b2-p1)
r3 = 100*norm(b3-p1)

N = 50000

# random range samples
Cc1 = (rand(Complex{Float64},N) .- b1)*100
Cc2 = (rand(Complex{Float64},N) .- b2)*100
Cc3 = (rand(Complex{Float64},N) .- b3)*100

A = zeros(100,100)

kr1 = 1 ./ (abs.(r1 .- norm.(Cc1)) )
kr2 = 1 ./ (abs.(r2 .- norm.(Cc2)) )
kr3 = 1 ./ (abs.(r3 .- norm.(Cc3)) )

kr1 .= exp.(kr1) .- 1
kr2 .= exp.(kr2) .- 1
kr3 .= exp.(kr3) .- 1

mask1 = rand(N) .< kr1
mask2 = rand(N) .< kr2
mask3 = rand(N) .< kr3

X1 = Cc1[mask1] .|> real .|> x->round(Int,x)
Y1 = Cc1[mask1] .|> imag .|> x->round(Int,x)
X1 .+= 100*real(b1)
Y1 .+= 100*imag(b1)
X1[X1 .== 0] .= 1
Y1[Y1 .== 0] .= 1
for i in 1:sum(mask1)
  A[X1[i],Y1[i]] += 1
end

X2 = Cc2[mask2] .|> real .|> x->round(Int,x)
Y2 = Cc2[mask2] .|> imag .|> x->round(Int,x)
X2 .+= 100*real(b2)
Y2 .+= 100*imag(b2)
X2[X2 .== 0] .= 1
Y2[Y2 .== 0] .= 1
for i in 1:sum(mask2)
  A[X2[i],Y2[i]] += 1
end

X3 = Cc3[mask3] .|> real .|> x->round(Int,x)
Y3 = Cc3[mask3] .|> imag .|> x->round(Int,x)
X3 .+= 100*real(b3)
Y3 .+= 100*imag(b3)
X3[X3 .== 0] .= 1
Y3[Y3 .== 0] .= 1
for i in 1:sum(mask3)
  A[X3[i],Y3[i]] += 1
end

pl = spy(A)

#


text:
@testset "BoxWedgeTail tests" begin

using DiffEqNoiseProcess, Test, Random
import Distributions
using Cubature
using Statistics
using DiffEqBase

# Test generation of boxes
W = BoxWedgeTail(0.0,zeros(2), box_grouping = :Columns)
# number of boxes with column-wise packing
@test length(W.boxes.boxes) == 2894
@test sum(W.boxes.probability) ≈ 0.911857 rtol=1e-5
@test DiffEqNoiseProcess.entropy(W.boxes.probability) ≈ 10.13 rtol=1e-5

W = BoxWedgeTail(0.0,zeros(2), box_grouping = :none)
# number of boxes of smallest size
@test length(W.boxes.boxes) == 119519
@test minimum(2*W.Δr*W.Δa*W.Δz .== W.boxes.probability)
@test sum(W.boxes.probability) ≈ 0.911857 rtol=1e-5
@test DiffEqNoiseProcess.entropy(W.boxes.probability) ≈ Distributions.entropy(W.boxes.dist,2) rtol=1e-5

W = BoxWedgeTail(0.0,zeros(2), box_grouping = :MinEntropy)
# number of boxes with min-entropy packing
@test length(W.boxes.boxes) == 2975
@test sum(W.boxes.probability) ≈ 0.911857 rtol=1e-5
@test DiffEqNoiseProcess.entropy(W.boxes.probability) ≈ 7.14 rtol=1e-3

# test generation of wedges

(val,err) = hcubature(x-> W.jpdf(x[1],x[2]), [zero(W.rM),zero(W.aM)],  [W.rM,W.aM])
@test 2*val ≈ 0.99732 rtol=1e-5

# number of columns
@test (W.aM/W.Δa) * (W.rM/W.Δr) == 4096
# probability for wedges
@test 2*val - sum(W.boxes.probability) ≈ 0.0855 rtol=1e-3
# on average less than half of the samples have to be rejected
@test sum(W.wedges.probability) <= 2*(2*val - sum(W.boxes.probability))
@test length(W.wedges.boxes) == 4096

# without sqeezing method
W = BoxWedgeTail(0.0,zeros(2), box_grouping = :MinEntropy, sqeezing=false)
# probability for wedges
@test 2*val - sum(W.boxes.probability) ≈ 0.0855 rtol=1e-3
# on average less than half of the samples have to be rejected
@test sum(W.wedges.probability) <= 2*(2*val - sum(W.boxes.probability))
@test length(W.wedges.boxes) == 4096

# check probabilities for tail approximation
# region 1
val1, _ = hcubature(x->W.jpdf(x[1],x[2]), [W.rM,zero(W.aM)], [12*one(W.rM),W.aM])
@test 2*val1 ≈ 2.98*1e-4 rtol=1e-3
@test W.tails.tail1.p ≈ 2*val1 rtol=1e-10
# region 2
val2, _ = hcubature(x->W.jpdf(x[1],x[2]), [W.rM,W.aM], [8*one(W.rM),8*one(W.aM)])
@test 2*val2 ≈ 3.65*1e-5 rtol=1e-3
@test W.tails.tail2.p ≈ 2*val2 rtol=1e-10
# region 3
val3, _ = hcubature(x->W.jpdf(x[1],x[2]), [2*one(W.rM),W.aM], [W.rM,8*one(W.aM)])
@test 2*val3 ≈ 1.80*1e-3 rtol=1e-2
@test W.tails.tail3.p ≈ 2*val3 rtol=1e-10
# region 4
val4, _ = hcubature(x->W.jpdf(x[1],x[2]), [zero(W.rM),W.aM], [one(W.rM)/2,6*one(W.aM)])
@test 2*val4 ≈ 1.61*1e-6 rtol=1e-2
@test W.tails.tail4.p ≈ 2*val4 rtol=1e-10
# region 5
val5, _ = hcubature(x->W.jpdf(x[1],x[2]), [one(W.rM)/2,W.aM], [one(W.rM),6*one(W.aM)])
@test 2*val5 ≈ 1.96*1e-5 rtol=1e-2
@test W.tails.tail5.p ≈ 2*val5 rtol=1e-10
# region 6
val6, _ = hcubature(x->W.jpdf(x[1],x[2]), [one(W.rM),W.aM], [3*one(W.rM)/2,6*one(W.aM)])
@test 2*val6 ≈ 1.20*1e-4 rtol=1e-2
@test W.tails.tail6.p ≈ 2*val6 rtol=1e-10
# region 7
val7, _ = hcubature(x->W.jpdf(x[1],x[2]), [3*one(W.rM)/2,W.aM], [2*one(W.rM),6*one(W.aM)])
@test 2*val7 ≈ 3.86*1e-4 rtol=1e-3
@test W.tails.tail7.p ≈ 2*val7 rtol=1e-10
# region 8
val8, _ = hcubature(x->W.jpdf(x[1],x[2]), [one(W.rM),6*one(W.aM)], [2*one(W.rM),8*one(W.aM)])
@test 2*val8 ≈ 6.57*1e-6 rtol=1e-3
@test W.tails.tail8.p ≈ 2*val8 rtol=1e-10
# region 9
val9, _ = hcubature(x->W.jpdf(x[1],x[2]), [2*one(W.rM),8*one(W.aM)], [5*one(W.rM),10*one(W.aM)])
@test 2*val9 ≈ 4.15*1e-6 rtol=1e-3
@test W.tails.tail9.p ≈ 2*val9 rtol=1e-10

# test remainder
@test 1.0-(2*(val+val1+val2+val3+val4+val5+val6+val7+val8+val9)) ≈ 3.81e-7 atol=1e-3

@test W.distBWT.p[1] ≈ 0.911857 rtol=1e-5
@test W.distBWT.p[2] ≈ 0.0855 rtol=1e-3
@test W.distBWT.p[3] ≈ 0.0027 rtol=1e-2
# tests if sampling works (samples from expected domain)

# boxes
#W = BoxWedgeTail(0.0,zeros(2), box_grouping = :MinEntropy)
samples = [DiffEqNoiseProcess.sample_box(W, W.boxes) for i=1:1000]
@test minimum(getindex.(samples, 1)) > 0
@test minimum(getindex.(samples, 2)) > 0
@test maximum(getindex.(samples, 1)) < W.rM
@test maximum(getindex.(samples, 2)) < W.aM

# using KernelDensity, StatsPlots, Plots
# r = getindex.(samples,1)
# a = getindex.(samples,2)
# BKDE = kde((r, a), boundary=((0,4),(0,4)))
# @test one(eltype(r)) ≈ sum(BKDE.density)*BKDE.x[2]*BKDE.y[2] rtol=1e-3
# pl1 = heatmap(BKDE)
# pl2 = marginalhist(r, a, fc = :plasma,  xlabel="r",  ylabel="a", bins=(64,64))
# pl3 = plot(0:W.Δr:4,0:W.Δa:4,W.jpdf,st =:contourf, xlabel="r", ylabel="a", title="f(r,a)")
# pl = plot(pl1, pl2, pl3, layout = (1, 3), legend = false, size=(800,250))
# savefig(pl, "boxes.png")

# wedges
samples = [DiffEqNoiseProcess.sample_wedge(W, W.wedges) for i=1:1_000]
@test minimum(getindex.(samples, 1)) > 0
@test minimum(getindex.(samples, 2)) > 0
@test maximum(getindex.(samples, 1)) < W.rM
@test maximum(getindex.(samples, 2)) < W.aM

# r = getindex.(samples,1)
# a = getindex.(samples,2)
# BKDE = kde((r, a), boundary=((0,4),(0,4)))
# @test one(eltype(r)) ≈ sum(BKDE.density)*BKDE.x[2]*BKDE.y[2] rtol=1e-3
# pl1 = heatmap(BKDE)
# pl2 = marginalhist(r, a, fc = :plasma,  xlabel="r",  ylabel="a", bins=(64,64))
# pl3 = plot(0:W.Δr:4,0:W.Δa:4,W.jpdf,st =:contourf, xlabel="r", ylabel="a", title="f(r,a)")
# pl = plot(pl1, pl2, pl3, layout = (1, 3), legend = false, size=(800,250))
# savefig(pl, "wedges.png")

# tails
# samples = [DiffEqNoiseProcess.sample_tail(W, W.tails) for i=1:100_000]
# r = getindex.(samples,1)
# a = getindex.(samples,2)
# BKDE = kde((r, a), boundary=((0,12),(0,10)))
# @test one(eltype(r)) ≈ sum(BKDE.density)*BKDE.x[2]*BKDE.y[2] rtol=1e-3
# pl1 = heatmap(BKDE)
# pl2 = marginalhist(r, a, fc = :plasma,  xlabel="r",  ylabel="a", bins=(64,64))
# pl3 = plot(0:W.Δr:12,0:W.Δa:10,W.jpdf,st =:contourf, xlabel="r", ylabel="a", title="f(r,a)")
# pl = plot(pl1, pl2, pl3, layout = (1, 3), legend = false, size=(800,250))
# savefig(pl, "tails.png")

samples = [DiffEqNoiseProcess.sample_tail(W.rng, W.jpdf, W.tails.tail1) for i=1:1_000]
@test minimum(getindex.(samples, 1)) > 4
@test minimum(getindex.(samples, 2)) > 0
@test maximum(getindex.(samples, 1)) < 12
@test maximum(getindex.(samples, 2)) < 4

samples = [DiffEqNoiseProcess.sample_tail(W.rng, W.jpdf, W.tails.tail2) for i=1:1_000]
@test minimum(getindex.(samples, 1)) > 4
@test minimum(getindex.(samples, 2)) > 4
@test maximum(getindex.(samples, 1)) < 8
@test maximum(getindex.(samples, 2)) < 8

samples = [DiffEqNoiseProcess.sample_tail(W.rng, W.jpdf, W.tails.tail3) for i=1:1_000]
@test minimum(getindex.(samples, 1)) > 2
@test minimum(getindex.(samples, 2)) > 4
@test maximum(getindex.(samples, 1)) < 4
@test maximum(getindex.(samples, 2)) < 8

samples = [DiffEqNoiseProcess.sample_tail(W.rng, W.jpdf, W.tails.tail4) for i=1:1_000]
@test minimum(getindex.(samples, 1)) > 0
@test minimum(getindex.(samples, 2)) > 4
@test maximum(getindex.(samples, 1)) < 0.5
@test maximum(getindex.(samples, 2)) < 6

samples = [DiffEqNoiseProcess.sample_tail(W.rng, W.jpdf, W.tails.tail5) for i=1:1_000]
@test minimum(getindex.(samples, 1)) > 0.5
@test minimum(getindex.(samples, 2)) > 4
@test maximum(getindex.(samples, 1)) < 1
@test maximum(getindex.(samples, 2)) < 6

samples = [DiffEqNoiseProcess.sample_tail(W.rng, W.jpdf, W.tails.tail6) for i=1:1_000]
@test minimum(getindex.(samples, 1)) > 1
@test minimum(getindex.(samples, 2)) > 4
@test maximum(getindex.(samples, 1)) < 1.5
@test maximum(getindex.(samples, 2)) < 6

samples = [DiffEqNoiseProcess.sample_tail(W.rng, W.jpdf, W.tails.tail7) for i=1:1_000]
@test minimum(getindex.(samples, 1)) > 1.5
@test minimum(getindex.(samples, 2)) > 4
@test maximum(getindex.(samples, 1)) < 2
@test maximum(getindex.(samples, 2)) < 6

samples = [DiffEqNoiseProcess.sample_tail(W.rng, W.jpdf, W.tails.tail8) for i=1:1_000]
@test minimum(getindex.(samples, 1)) > 1
@test minimum(getindex.(samples, 2)) > 6
@test maximum(getindex.(samples, 1)) < 2
@test maximum(getindex.(samples, 2)) < 8

samples = [DiffEqNoiseProcess.sample_tail(W.rng, W.jpdf, W.tails.tail9) for i=1:1_000]
@test minimum(getindex.(samples, 1)) > 2
@test minimum(getindex.(samples, 2)) > 8
@test maximum(getindex.(samples, 1)) < 5
@test maximum(getindex.(samples, 2)) < 10

# compare samples with pdf and check some statistics

samples = [DiffEqNoiseProcess.sample_distribution(W,one(W.dt),nothing,nothing,nothing,W.rng) for i=1:100_000]

dW1 = getindex.(getindex.(samples,1),1)
dW2 = getindex.(getindex.(samples,1),2)
r = sqrt.(dW1.^2 + dW2.^2)
a = abs.(getindex.(getindex.(samples,2),1))

@test mean(dW1) ≈ zero(eltype(r)) atol=1e-2
@test mean(dW2) ≈ zero(eltype(r)) atol=1e-2

@test var(dW1) ≈ one(eltype(r)) rtol=1e-1
@test var(dW2) ≈ one(eltype(r)) rtol=1e-1

# BKDE = kde((r, a), boundary=((0,6),(0,6)))
# @test one(eltype(r)) ≈ sum(BKDE.density)*BKDE.x[2]*BKDE.y[2] rtol=1e-3
# pl1 = heatmap(BKDE)
# pl2 = marginalhist(r, a, fc = :plasma,  xlabel="r",  ylabel="a", bins=(64,64))
# pl3 = plot(0:0.1:6,0:0.1:6,W.jpdf,st =:contourf, xlabel="r", ylabel="a", title="f(r,a)")
# pl = plot(pl1, pl2, pl3, layout = (1, 3), legend = false, size=(800,250))

# simulate oop noise problem

dt = 0.1
calculate_step!(W,dt,nothing,nothing)

for i in 1:10
  accept_step!(W,dt,nothing,nothing)
end

prob = NoiseProblem(W,(0.0,1.0));
sol = solve(prob;dt=0.1)

# simulate iip noise problem
W = BoxWedgeTail!(0.0,zeros(2))

dt = 0.1
calculate_step!(W,dt,nothing,nothing)

for i in 1:10
  accept_step!(W,dt,nothing,nothing)
end

prob = NoiseProblem(W,(0.0,1.0));
sol = solve(prob;dt=0.1)

end


text:
#
# tlibs2 -- julia module
# @author Tobias Weber <tobias.weber@tum.de>, <tweber@ill.fr>
# @date 2017 -- 2018
# @license GPLv3, see <<LICENSE>> file
# @desc Forked on 7-Nov-2018 from my privately and TUM-PhD-developed <<tlibs>> project (https://github.com/t-weber/tlibs).
#
# ----------------------------------------------------------------------------
# tlibs
# Copyright (C) 2017-2021  Tobias WEBER (Institut Laue-Langevin (ILL),
#                          Grenoble, France).
# Copyright (C) 2015-2017  Tobias WEBER (Technische Universitaet Muenchen
#                          (TUM), Garching, Germany).
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, version 3 of the License.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
# ----------------------------------------------------------------------------
#

__precompile__()
module tl2

t_real = Float64
g_enable_debug = 1

#
# initialises tlibs2
#
function __init__()
	ccall((:load_tl2, :tl2_jl), Cvoid, (Cint,), g_enable_debug)
end

#
# loads instrument data files
#
function loadinstr(strFile::String) :: Array{Any, 1}
	scandata = ccall((:load_instr, :tl2_jl), Array{Any, 1}, (Cstring,), strFile)
	if length(scandata) == 0 return scandata end

	thedict = Dict{String,String}()
	for (key, val) in zip(scandata[3], scandata[4])
		thedict[key] = val
	end

	# header names, data matrix, property dict
	return [ scandata[1], scandata[2], thedict ]
end

#
# build the @cfunction call as string which requires literal arguments tuples
#
function __build_cfunc_string(fkt, num_args)
	ty = string(t_real)

	argtup = "("
	for iArg in 1:num_args
		argtup *= ty * ", "
	end
	argtup *= ")"

	str = "@cfunction($fkt" * ", " * ty * ", " * argtup * ")"
	#println("Function: " * str)
	return str
end

#
# function fitting
#
function fit(fkt, x, y, yerr; fixed = [], values = Dict(), errors = Dict())
	# find number of function arguments
	meth = methods(fkt).ms[1]
	num_args = length(meth.sig.parameters) - 1
	num_free_params = num_args - 1

	# get function argument names
	strMeth = repr(meth)
	strArgs = strMeth[first(something(findfirst("(", strMeth), 0:-1))+1 : first(something(findfirst(")", strMeth), 0:-1))-1]
	arrArgs = map(strip, split(strArgs, ","))
	arrArgs = map(String, arrArgs)
	arrParams = arrArgs[2 : length(arrArgs)]	# only free params

	# build values/errors array
	iParam = 0
	arrValues = Array{t_real, 1}(undef, num_free_params)
	arrErrors = Array{t_real, 1}(undef, num_free_params)

	for strArg in arrArgs
		# skip x parameter
		if(iParam == 0)
			iParam += 1
			continue
		end

		valArg = get(values, strArg, 0.)
		errArg = get(errors, strArg, valArg*0.1)

		arrValues[iParam] = valArg
		arrErrors[iParam] = errArg
		iParam += 1
	end

	# map to a C function pointer with num_args arguments
	#cfkt = cfunction(fkt, t_real, NTuple{num_args, t_real})
	cfuncstr = __build_cfunc_string(fkt, num_args)
	cfkt = eval(Meta.parse(cfuncstr))

	# call C function pointer
	bOk = ccall((:fit, :tl2_jl),
		# return type
		Cint,

		# arg types
		(Ptr{Cvoid}, Csize_t,
		Ptr{t_real}, Ptr{t_real}, Ptr{t_real}, Csize_t,
		Array{String, 1}, Array{String, 1},
		Ptr{t_real}, Ptr{t_real}),

		# args
		cfkt, num_free_params, x, y, yerr, length(x), arrArgs, fixed, arrValues, arrErrors)

	# build map of values & errors
	dictRet = Dict()

	for (strParam, valArg, valErr) in zip(arrParams, arrValues, arrErrors)
		strErr = strParam * "_err"
		dictRet[strParam] = valArg
		dictRet[strErr] = valErr
	end

	dictRet["<valid>"] = bOk
	dictRet["<func>"] = fkt
	dictRet["<args>"] = arrParams
	dictRet["<args_val>"] = arrValues
	dictRet["<args_err>"] = arrErrors
	return dictRet
end

#
# function minimisation
#
function minimise(fkt; fixed = [], values = Dict(), errors = Dict())
	# find number of function arguments
	meth = methods(fkt).ms[1]
	num_args = length(meth.sig.parameters) - 1

	# get function argument names
	strMeth = repr(meth)
	strArgs = strMeth[first(something(findfirst("(", strMeth), 0:-1)) + 1 : first(something(findfirst(")", strMeth), 0:-1)) - 1]
	arrArgs = map(strip, split(strArgs, ","))
	arrArgs = map(String, arrArgs)
	arrParams = arrArgs[1 : length(arrArgs)]

	# build values/errors array
	arrValues = Array{t_real, 1}(undef, num_args)
	arrErrors = Array{t_real, 1}(undef, num_args)

	iParam = 1
	for strArg in arrArgs
		valArg = get(values, strArg, 0.)
		errArg = get(errors, strArg, valArg*0.1)

		arrValues[iParam] = valArg
		arrErrors[iParam] = errArg
		iParam += 1
	end

	# map to a C function pointer with num_args arguments
	#cfkt = cfunction(fkt, t_real, NTuple{num_args, t_real})
	cfuncstr = __build_cfunc_string(fkt, num_args)
	cfkt = eval(Meta.parse(cfuncstr))

	# call C function pointer
	bOk = ccall((:minimise, :tl2_jl),
		# return type
		Cint,

		# arg types
		(Ptr{Cvoid}, Csize_t,
		Array{String, 1}, Array{String, 1},
		Ptr{t_real}, Ptr{t_real}),

		# args
		cfkt.ptr, num_args, arrArgs, fixed, arrValues, arrErrors)

	# build map of values & errors
	dictRet = Dict()

	for (strParam, valArg, valErr) in zip(arrParams, arrValues, arrErrors)
		strErr = strParam * "_err"
		dictRet[strParam] = valArg
		dictRet[strErr] = valErr
	end

	dictRet["<valid>"] = bOk
	dictRet["<func>"] = fkt
	dictRet["<args>"] = arrParams
	dictRet["<args_val>"] = arrValues
	dictRet["<args_err>"] = arrErrors
	return dictRet
end

#
# evaluate the fitted function using the fitted parameters
#
function eval_func(fitresult, xmin, xmax, count=128)
	fkt = fitresult["<func>"]

	xs = collect(range(xmin, stop=xmax, length=count))
	ys = fkt(xs, fitresult["<args_val>"]...)

	return [ Array(xs), Array(ys) ]
end

#
# Gauss model
#
function gauss_model_amp(x, x0, sigma, amp, offs)
	return amp.*exp.(-0.5 .* ((x.-x0)./sigma).^2.) .+ offs
end

#
# Lorentz model
#
function lorentz_model_amp(x, x0, hwhm, amp, offs)
	return amp .* hwhm.^ 2. ./ ((x.-x0).^ 2. .+ hwhm.^2.) .+ offs
end

#
# Poisson error
#
function poisson_err(x) :: t_real
	if(x == 0)
		return 1
	end

	return sqrt(x)
end

#
# Conversions
#
SIGMA2HWHM = sqrt(2. * log(2.))
SIGMA2FWHM = 2. * SIGMA2HWHM
HWHM2SIGMA = 1. / SIGMA2HWHM
FWHM2SIGMA = 1. / SIGMA2FWHM

end		# module tl2


text:
@enum Button_Names begin
    app_close_button
    dropdown_button_datei
    dropdown_button_anzeigen
    dropdown_button_navigation
    dropdown_button_tools
    dropdown_button_hilfe
    button_abbrechen
    button_aendern
    button_bearbeiten
    button_clipboard
    button_speichern
    button_zuruecksetzen
    menu_button_1
    menu_button_preferences
    menu_ueber
    small_button_1
    small_button_2
    small_button_3
    small_button_4
    small_button_5
    small_button_6
    small_button_7
    small_button_8
    small_button_9
    small_button_10
    small_button_11
    small_button_12
    small_button_13

end
@enum Window_Names begin
    main_window
    window_preferences
    window_ueber
    dropdown_window
    no_window
end
@enum On_click_listeners begin
    placeholder
end


text:
"""
Interface between `Iris` and `NonlinearEigenproblems`
"""
module NonlinearEigenproblemsInterface

export IrisLinSolverCreator

using NonlinearEigenproblems
using ..Common

################################################################################
# overload all the NonlinearEigenproblems solvers

import ..HelmholtzNEP
import ..MaxwellNEP

fnames = names(NonlinearEigenproblems.NEPSolver)
for fn ∈ fnames
    @eval begin
        """
            $($fn)(nep::HelmholtzNEP, args...; kwargs...) -> λ, s::ScalarField
        """
        function NonlinearEigenproblems.$(fn)(nep::HelmholtzNEP, args...; kwargs...)
            λ, v = $(fn)(nep.nep,args...;kwargs...)
            return λ, ScalarField(nep,v)
        end
		"""
            $($fn)(nep::MaxwellNEP, args...; kwargs...) -> λ, e::ElectricField
        """
        function NonlinearEigenproblems.$(fn)(nep::MaxwellNEP, args...; kwargs...)
            λ, v = $(fn)(nep.nep,args...;kwargs...)
            return λ, ElectricField(nep,v)
        end
    end
end

################################################################################
# DEFINE LINSOLVERCREATOR, per instructions in NonlinearEigenproblems docs

using LinearAlgebra
using SparseArrays

"""
	IrisLinSolverCreator{S}
"""
struct IrisLinSolverCreator{S} <: LinSolverCreator solver::S end

"""
	IrisLinSolverCreator([lupack])

`lupack` defaults to $DEFAULT_LUPACK.
For use in [`NonlinearEigenproblems`](https://github.com/nep-pack/NonlinearEigenproblems.jl).
"""
IrisLinSolverCreator(s::T) where T<:AbstractLUPACK = IrisLinSolverCreator{T}(s)
IrisLinSolverCreator(::Type{T}) where T<:AbstractLUPACK = IrisLinSolverCreator{T}(T())
IrisLinSolverCreator() where T<:AbstractLUPACK = IrisLinSolverCreator(DEFAULT_LUPACK)

struct IrisLinSolver{S,L,T} <: LinSolver
    LU::L
    x::Vector{T}
    Z::SparseMatrixCSC{T,Int}
end

function NonlinearEigenproblems.create_linsolver(ILS::IrisLinSolverCreator{S}, nep, λ) where S
    LU = lu(compute_Mder(nep,λ),S)
    x = Vector{ComplexF64}(undef,size(nep,1))
    return IrisLinSolver{S,typeof(LU),eltype(x)}(LU,x,spzeros(eltype(x),length(x),length(x)))
end

function NonlinearEigenproblems.lin_solve(solver::IrisLinSolver{S}, b::Vector; tol=eps()) where S
    ldiv!(solver.x, solver.LU, b)
    return solver.x
end

################################################################################
# maxwelleigen and helmholtzeigen, NEP

import ..AbstractNonlinearEigenproblem
import ..INDEX_OFFSET
import ..orthogonalize!

function iris_eigen_nonlineareigenproblems(
            nep::AbstractNonlinearEigenproblem,
            ω::Number,
            args...;
            verbose::Bool = false,
            lupack::AbstractLUPACK = DEFAULT_LUPACK,
            linsolvercreator::LinSolverCreator = IrisLinSolverCreator(lupack),
            kwargs...
            ) where S<:AbstractLUPACK

	# is background fields ωs and ψs are provided, saturate nep
    isempty(args) ? nothing : nep(ω,args...)
	# use iar_chebyshev as default
    λ::Vector{ComplexF64}, v::Matrix{ComplexF64} = iar_chebyshev(nep; σ=ω, logger=Int(verbose), linsolvercreator=linsolvercreator, kwargs...)
	# wrap v in Vector field, depending on Helmholtz vs Maxwell
	if typeof(nep) <: HelmholtzNEP
		ψ = ScalarField(nep,v)
	elseif typeof(nep) <: MaxwellNEP
		ψ = ElectricField(nep,v)
	end
	# Normalize according to (ψ₁,ψ₂)=δ₁₂, fixing size(ψ,1)÷2+INDEX_OFFSET to be real
    normalize!(nep.simulation, ψ, nep.nep.A[end-1], size(ψ,1)÷2+INDEX_OFFSET)
	# orthogonalize in case of degeneracy
	orthogonalize!(ψ,nep.simulation, λ, nep.nep.A[end-1], nep.ky, nep.kz)
    return λ,ψ
end

import ..DOC_NEP_H
import ..DOC_NEP_M
import ..DEFAULT_NONLINEAR_EIGENSOLVER
import ..helmholtzeigen
import ..maxwelleigen

if DEFAULT_NONLINEAR_EIGENSOLVER == :NonlinearEigenproblems
	@doc """
	$DOC_NEP_H
	`neigs` Number of eigenpairs (`6`) |
	`maxit` Maximum iterations (`30`) |
	`verbose` (`false`) |
	`v` Initial eigenvector approximation (`rand(size(nep,1),1)`)

	Advanced keywords:
	`logger` Integer print level (`verbose`) |
	`tolerance` Convergence tolerance (`eps()*10000`) |
	`lupack` LU solver (`$DEFAULT_LUPACK`) |
	`linsolvercreator` Define linear solver (`IrisLinSolverCreator(lupack)`) |
	`check_error_every` Check for convergence every this number of iterations (`1`) |
	`γ` Eigenvalue scaling in `λ=γs+Ω` (`1`) |
	`orthmethod` Orthogonalization method, see `IterativeSolvers` (`DGKS`) |
	`errmeasure` See `NonlinearEigenproblems` documentation |
	`compute_y0_method` How to find next vector in Krylov space, see `NonlinearEigenproblems` docs (`ComputeY0ChebAuto`) |
	a (`-1`) |
	b (`1`)
	""" ->
	helmholtzeigen(nep::HelmholtzNEP,args...;kwargs...) = iris_eigen_nonlineareigenproblems(nep,args...;kwargs...)

	@doc """
	$DOC_NEP_M
	`neigs` Number of eigenpairs (`6`) |
	`maxit` Maximum iterations (`30`) |
	`verbose` (`false`) |
	`v` Initial eigenvector approximation (`rand(size(nep,1),1)`)

	Advanced keywords:
	`logger` Integer print level (`verbose`) |
	`tolerance` Convergence tolerance (`eps()*10000`) |
	`lupack` LU solver (`$DEFAULT_LUPACK`) |
	`linsolvercreator` Define linear solver (`IrisLinSolverCreator(lupack)`) |
	`check_error_every` Check for convergence every this number of iterations (`1`) |
	`γ` Eigenvalue scaling in `λ=γs+Ω` (`1`) |
	`orthmethod` Orthogonalization method, see `IterativeSolvers` (`DGKS`) |
	`errmeasure` See `NonlinearEigenproblems` documentation |
	`compute_y0_method` How to find next vector in Krylov space, see `NonlinearEigenproblems` docs (`ComputeY0ChebAuto`) |
	a (`-1`) |
	b (`1`)
	""" ->
	maxwelleigen(nep::MaxwellNEP,args...;kwargs...) = iris_eigen_nonlineareigenproblems(nep,args...;kwargs...)
end

end # module

using .NonlinearEigenproblemsInterface