text:
function read_asdf!(  S::GphysData,
                      hdf::String,
                      id::String,
                      s::TimeSpec,
                      t::TimeSpec,
                      msr::Bool,
                      v::Int64  )

  SX = SeisData() # for XML
  idr = isa(id, String) ? id_to_regex(id) : id
  (v > 2) && println("Reading IDs that match ", idr)

  if typeof(s) == String && typeof(t) == String
    d0 = s
    d1 = t
    ts = DateTime(s).instant.periods.value*1000 - dtconst
    te = DateTime(t).instant.periods.value*1000 - dtconst
  else
    (d0, d1) = parsetimewin(s, t)
    ts = DateTime(d0).instant.periods.value*1000 - dtconst
    te = DateTime(d1).instant.periods.value*1000 - dtconst
  end
  ts *= 1000
  te *= 1000
  Δ = 0
  fs = 0.0

  # this nesting is a mess
  netsta = netsta_to_regex(id)
  idr = id_to_regex(id)
  f = h5open(hdf, "r")
  W = f["Waveforms"]
  A = names(W)
  sort!(A)
  (v > 2) && println("Net.sta found: ", A)

  for i in A
    if occursin(netsta, i)
      w = W[i]
      N = names(w)
      sort!(N)
      (v > 2) && println("Traces found: ", N)
      for n in N
        if n == "StationXML"
          sxml = String(UInt8.(read(w[n])))
          read_station_xml!(SX, sxml, msr=msr, s=s, t=t, v=v)
        elseif occursin(idr, n)
          x = w[n]
          nx = length(x)
          t0 = read(x["starttime"])
          fs = read(x["sampling_rate"])

          # convert fs to sampling interval in ns
          Δ = round(Int64, 1.0e9/fs)
          t1 = t0 + (nx-1)*Δ
          (v > 2) && println("t0 = ", t0,"; t1 = ", t1, "; ts = ", ts, "; te = ", te)

          if (ts ≤ t1) && (te ≥ t0)
            xi = 0
            i0, i1, t2 = get_trace_bounds(ts, te, t0, t1, Δ, nx)
            ni = i1-i0+1
            cid = String(split(n, "_", limit=2, keepempty=true)[1])
            j = findid(cid, S.id)
            (v > 2) && println("cid = ", cid, " (found at index in S = ", j, ")")
            nX = div(te-ts, Δ)+1
            if j == 0
              T = eltype(x)
              push!(S, SeisChannel(id = cid,
                                   fs = fs,
                                   x = Array{T,1}(undef, nX)))
              j = S.n
              if has(x, "event_id")
                S.misc[j]["event_id"] = read(x["event_id"])
              end
            else
              ct = getindex(getfield(S, :t), j)
              nt = div(lastindex(ct), 2)
              L = lastindex(S.x[j])
              if nt > 0
                xi = getindex(ct, nt)
                check_for_gap!(S, j, div(t2, 1000), ni, v)
              end
              if xi + ni > L
                resize!(S.x[j], xi + max(ni, nX))
              end
              if S.fs[j] == 0.0
                S.fs[j] = fs
              end
            end

            if xi == 0
              ct = Array{Int64, 2}(undef, 2, 2)
              setindex!(ct, one(Int64), 1)
              setindex!(ct, ni, 2)
              setindex!(ct, div(t2, 1000), 3)
              setindex!(ct, zero(Int64), 4)
              setindex!(getfield(S, :t), ct, j)
            end

            si = xi+1
            ei = si+ni-1
            X = S.x[j]
            load_data!(X, x, i0:i1, si:ei)
          end
        end
      end
    end
  end

  # Ensure data source is logged accurately
  fill!(S.src, realpath(hdf))

  # merge in the XML that we read
  sxml_mergehdr!(S, SX, app=false, nofs=true, v=v)

  trunc_x!(S)

  # Done
  close(f)
  return S
end

function read_asdf(hdf::String, id::Union{String,Regex}, s::TimeSpec, t::TimeSpec, msr::Bool, v::Int64)
  S = SeisData()
  read_asdf!(S, hdf, id, s, t, msr, v)
  return S
end


text:
module TestContextManagers
include("../../../examples/contextmanagers.jl")
end  # module

ext: jl


text:
tests =
  [
    "common",
    "VariableStates",
    "BasicDiscUnvParameterState",
    "BasicDiscMuvParameterState",
    "BasicContUnvParameterState",
    "BasicContMuvParameterState",
    "VariableNStates",
    "ParameterNStates",
    "VariableIOStreams",
    "ParameterIOStreams",
    "BasicContUnvParameter",
    "BasicContMuvParameter",
    "GenericModel",
    "generators"
  ]

println("Running tests:")

for t in tests
  tfile = t*".jl"
  println("  * $(tfile) *")
  include(tfile)
end


text:
using Base.Test
using AutomotiveDrivingModels
using AutoViz
using NBInclude

nbinclude(Pkg.dir("AutoViz", "doc", "AutoViz.ipynb"))


text:
using Test
import BinaryTrees
bt = BinaryTrees


text:
# random interval
function rand(::Type{Interval}; N::Type{<:Real}=Float64,
              rng::AbstractRNG=GLOBAL_RNG)::Interval{N}

    x, y = randn(rng, N), randn(rng, N)
    return x < y ? Interval(x, y) : Interval(y, x)
end

# create a random interval for the given numeric type and random-number generator
@inline function _rand_interval(; N=Float64, rng::AbstractRNG=GLOBAL_RNG)
    x, y = randn(rng, N), randn(rng, N)
    return x < y ? Interval(x, y) : Interval(y, x)
end

"""
rand(::Type{IntervalMatrix}, m::Int=2, [n]::Int=m;
N=Float64, rng::AbstractRNG=GLOBAL_RNG)

Return a random interval matrix of the given size and numeric type.

### Input

- `IntervalMatrix` -- type, used for dispatch
- `m`              -- (optional, default: `2`) number of rows
- `n`              -- (optional, default: `m`) number of columns
- `rng`            -- (optional, default: `GLOBAL_RNG`) random-number generator

### Output

An interval matrix of size ``m × n`` whose coefficients are normally-distributed
intervals of type `N` with mean `0` and standard deviation `1`.

### Notes

If this function is called with only one argument, it creates a square matrix,
because the number of columns defaults to the number of rows.
"""
function rand(::Type{IntervalMatrix}, m::Int=2, n::Int=m;
    N = Float64, rng::AbstractRNG=GLOBAL_RNG)
    B = Matrix{Interval{N}}(undef, m, n)

    for j in 1:n
        for i in 1:m
            @inbounds B[i, j] = _rand_interval(N=N, rng=rng)
        end
    end

    return IntervalMatrix(B)
end

"""
sample(A::IntervalMatrix{T}; rng::AbstractRNG=GLOBAL_RNG) where {T}

Return a sample of the given random interval matrix.

### Input

- `A`   -- interval matrix
- `m`   -- (optional, default: `2`) number of rows
- `n`   -- (optional, default: `2`) number of columns
- `rng` -- (optional, default: `GLOBAL_RNG`) random-number generator

### Output

An interval matrix of size ``m × n`` whose coefficients are normally-distributed
intervals of type `N` with mean `0` and standard deviation `1`.
"""
function sample(A::IntervalMatrix{T}; rng::AbstractRNG=GLOBAL_RNG) where {T}
    m, n = size(A)
    B = Matrix{T}(undef, m, n)

    @inbounds for j in 1:n
        for i in 1:m
            itv = A[i, j]
            B[i, j] = (sup(itv) - inf(itv)) * rand(rng) + inf(itv)
        end
    end

    return B
end


text:
"""
	drawcircle!(img,i,j,r[,trange])
Modifies the image to superimpose a circle of specified center and radius. If `trange` is given, it's a 2-vector defining the range of angles to be used, measured ccw from "straight down" in the usual image visualization (i.e., vertically flipped).
"""
function drawcircle!(img,i,j,r,trange::AbstractVector=[-π,π];usecolor=RGB(1,0,0.8))
	t = π/2 .+ LinRange(trange[1],trange[2],5size(img,2))
	valid = i -> checkbounds(Bool,img,i)
	for r in 0.99r:1.01r
		idx = CartesianIndex.(zip(round.(Int,i.+r*sin.(t)),round.(Int,j.+r*cos.(t))))
		idx = filter(valid,idx)
		img[idx] .= usecolor
	end
	return img
end

"""
	makemovie(files,result,size;numframes=Inf)
Make a movie to show the detected cornea on each image named in the string vector `files`, using the detection results in `result` and the given `size` as [height,width]. If a value is given for `numframes`, frames will be skipped to bring the total frame count to be no greater than that value. Returns a vector of images.
"""
function makemovie(folder::AbstractVector{String},result,sz;numframes=Inf)
	imgstack = []
	circ = sz[1].*[ result.cenrow result.cencol result.radius ]
	N = length(folder)
	select = 1:max(1,ceil(Int,(N-1)/numframes)):N
	@showprogress for k in select
		img = imresize(RGB.(load(folder[k])),sz...)
		drawcircle!(img,circ[k,:]...)
		push!(imgstack,img)
	end
	return imgstack
end

max_stars_repo_name: tobydriscoll/cornea-detection

max_issues_repo_name: tobydriscoll/cornea-detection

max_forks_repo_name: tobydriscoll/cornea-detection


text:
using MLKernels

# Constructing Standard Kernels

GaussianKernel()
SquaredDistanceKernel()

ϕ = ExponentialClass(2.0, 1.0);
κ = SquaredDistanceKernel(1.0);
ψ = ϕ ∘ κ   # use \circ for '∘'

ismercer(κ)
isnegdef(κ)
ismercer(ψ)
isnegdef(ψ)

w = round(rand(5),3);
ARD(κ, w)

# Kernel Functions

x = rand(3); y = rand(3);
kernel(ψ, x[1], y[1])
kernel(ψ, x, y)

X = rand(5,3);
kernelmatrix(ψ, X)
kernelmatrix(ψ, X', true)

Y = rand(4,3);
kernelmatrix(ψ, X, Y)
kernelmatrix(ψ, X', Y', true)

kernelmatrix(ψ, [X; Y])[1:5, 6:9]

# Kernel Operations

κ1 = GaussianKernel();

2k1 + 3

κ2 = PolynomialKernel();

κ1 + κ2

κ1 * κ2

κ3 = SineSquaredKernel();

κ3^0.5

κ4 = 2ScalarProductKernel() + 3;

κ4^3

exp(κ4)

tanh(κ4)


text:
# # User Recipes

"""
    _process_userrecipes(plt, plotattributes, args)

Wrap input arguments in a `RecipeData' vector and recursively apply user recipes and type
recipes on the first element. Prepend the returned `RecipeData` vector. If an element with
empy `args` is returned pop it from the vector, finish up, and it to vector of `Dict`s with
processed series. When all arguments are processed return the series `Dict`.
"""
function _process_userrecipes!(plt, plotattributes, args)
    still_to_process = _recipedata_vector(plt, plotattributes, args)

    # For plotting recipes, we swap out the args and update the parameter dictionary.  We are keeping a stack of series that still need to be processed.
    #
    # On each pass through the loop, we pop one off and apply the recipe.
    # the recipe will return a list a Series objects. The ones that are
    # finished (no more args) get added to the `kw_list`, and the ones that are not
    # are placed on top of the stack and are then processed further.

    kw_list = KW[]
    while !isempty(still_to_process)
        # grab the first in line to be processed and either add it to the kw_list or
        # pass it through apply_recipe to generate a list of RecipeData objects
        # (data + attributes) for further processing.
        next_series = popfirst!(still_to_process)
        # recipedata should be of type RecipeData.
        # if it's not then the inputs must not have been fully processed by recipes
        if !(typeof(next_series) <: RecipeData)
            error("Inputs couldn't be processed... expected RecipeData but got: $next_series")
        end
        if isempty(next_series.args)
            _finish_userrecipe!(plt, kw_list, next_series)
        else
            rd_list =
                RecipesBase.apply_recipe(next_series.plotattributes, next_series.args...)
            warn_on_recipe_aliases!(plt, rd_list, :user, next_series.args...)
            prepend!(still_to_process, rd_list)
        end
    end

    # don't allow something else to handle it
    plotattributes[:smooth] = false
    kw_list
end

# TODO Move this to api.jl?

function _recipedata_vector(plt, plotattributes, args)
    still_to_process = RecipeData[]
    # the grouping mechanism is a recipe on a GroupBy object
    # we simply add the GroupBy object to the front of the args list to allow
    # the recipe to be applied
    if haskey(plotattributes, :group)
        args = (_extract_group_attributes(plotattributes[:group], args...), args...)
    end

    # if we were passed a vector/matrix of seriestypes and there's more than one row,
    # we want to duplicate the inputs, once for each seriestype row.
    if !isempty(args)
        append!(still_to_process, _expand_seriestype_array(plotattributes, args))
    end

    # remove subplot and axis args from plotattributes...
    # they will be passed through in the kw_list
    if !isempty(args)
        for (k, v) in plotattributes
            if is_subplot_attribute(plt, k) || is_axis_attribute(plt, k)
                reset_kw!(plotattributes, k)
            end
        end
    end

    still_to_process
end

function _expand_seriestype_array(plotattributes, args)
    sts = get(plotattributes, :seriestype, :path)
    if typeof(sts) <: AbstractArray
        reset_kw!(plotattributes, :seriestype)
        rd = Vector{RecipeData}(undef, size(sts, 1))
        for r in axes(sts, 1)
            dc = copy(plotattributes)
            dc[:seriestype] = sts[r:r, :]
            rd[r] = RecipeData(dc, args)
        end
        rd
    else
        RecipeData[RecipeData(copy(plotattributes), args)]
    end
end

function _finish_userrecipe!(plt, kw_list, recipedata)
    # when the arg tuple is empty, that means there's nothing left to recursively
    # process... finish up and add to the kw_list
    kw = recipedata.plotattributes
    preprocess_attributes!(plt, kw)
    # if there was a grouping, filter the data here
    _filter_input_data!(kw)
    process_userrecipe!(plt, kw_list, kw)
end

# --------------------------------
# Fallback user recipes
# --------------------------------

# These call `_apply_type_recipe` in type_recipe.jl and finally the `SliceIt` recipe in
# series.jl.

# handle "type recipes" by converting inputs, and then either re-calling or slicing
@recipe function f(x, y, z)
    wrap_surfaces!(plotattributes, x, y, z)
    did_replace = false
    newx = _apply_type_recipe(plotattributes, x, :x)
    x === newx || (did_replace = true)
    newy = _apply_type_recipe(plotattributes, y, :y)
    y === newy || (did_replace = true)
    newz = _apply_type_recipe(plotattributes, z, :z)
    z === newz || (did_replace = true)
    if did_replace
        newx, newy, newz
    else
        SliceIt, x, y, z
    end
end
@recipe function f(x, y)
    wrap_surfaces!(plotattributes, x, y)
    did_replace = false
    newx = _apply_type_recipe(plotattributes, x, :x)
    x === newx || (did_replace = true)
    newy = _apply_type_recipe(plotattributes, y, :y)
    y === newy || (did_replace = true)
    if did_replace
        newx, newy
    else
        SliceIt, x, y, nothing
    end
end
@recipe function f(y)
    wrap_surfaces!(plotattributes, y)
    newy = _apply_type_recipe(plotattributes, y, :y)
    if y !== newy
        newy
    else
        SliceIt, nothing, y, nothing
    end
end

# if there's more than 3 inputs, it can't be passed directly to SliceIt
# so we'll apply_type_recipe to all of them
@recipe function f(v1, v2, v3, v4, vrest...)
    did_replace = false
    newargs = map(
        v -> begin
            newv = _apply_type_recipe(plotattributes, v, :unknown)
            if newv !== v
                did_replace = true
            end
            newv
        end,
        (v1, v2, v3, v4, vrest...),
    )
    if !did_replace
        error("Couldn't process recipe args: $(map(typeof, (v1, v2, v3, v4, vrest...)))")
    end
    newargs
end

# helper function to ensure relevant attributes are wrapped by Surface
function wrap_surfaces!(plotattributes, args...) end
wrap_surfaces!(plotattributes, x::AMat, y::AMat, z::AMat) = wrap_surfaces!(plotattributes)
wrap_surfaces!(plotattributes, x::AVec, y::AVec, z::AMat) = wrap_surfaces!(plotattributes)
function wrap_surfaces!(plotattributes, x::AVec, y::AVec, z::Surface)
    wrap_surfaces!(plotattributes)
end
function wrap_surfaces!(plotattributes)
    if haskey(plotattributes, :fill_z)
        v = plotattributes[:fill_z]
        if !isa(v, Surface)
            plotattributes[:fill_z] = Surface(v)
        end
    end
end

# --------------------------------
# Special Cases
# --------------------------------

# --------------------------------
# 1 argument

@recipe function f(n::Integer)
    if is3d(plotattributes)
        SliceIt, n, n, n
    else
        SliceIt, n, n, nothing
    end
end

# return a surface if this is a 3d plot, otherwise let it be sliced up
@recipe function f(mat::AMat)
    if is3d(plotattributes)
        n, m = axes(mat)
        m, n, Surface(mat)
    else
        nothing, mat, nothing
    end
end

# if a matrix is wrapped by Formatted, do similar logic, but wrap data with Surface
@recipe function f(fmt::Formatted{<:AMat})
    if is3d(plotattributes)
        mat = fmt.data
        n, m = axes(mat)
        m, n, Formatted(Surface(mat), fmt.formatter)
    else
        nothing, fmt, nothing
    end
end

# assume this is a Volume, so construct one
@recipe function f(vol::AbstractArray{<:MaybeNumber, 3}, args...)
    seriestype := :volume
    SliceIt, nothing, Volume(vol, args...), nothing
end

# Dicts: each entry is a data point (x,y)=(key,value)
@recipe function f(d::AbstractDict)
    seriestype --> :line
    collect(keys(d)), collect(values(d))
end
# function without range... use the current range of the x-axis
@recipe function f(f::FuncOrFuncs{F}) where {F <: Function}
    plt = plotattributes[:plot_object]
    xmin, xmax = if haskey(plotattributes, :xlims)
        plotattributes[:xlims]
    else
        try
            get_axis_limits(plt, :x)
        catch
            xinv = inverse_scale_func(get(plotattributes, :xscale, :identity))
            xm = PlotUtils.tryrange(f, xinv.([-5, -1, 0, 0.01]))
            xm, PlotUtils.tryrange(f, filter(x -> x > xm, xinv.([5, 1, 0.99, 0, -0.01])))
        end
    end
    f, xmin, xmax
end

# --------------------------------
# 2 arguments

# if functions come first, just swap the order (not to be confused with parametric
# functions... as there would be more than one function passed in)
@recipe function f(f::FuncOrFuncs{F}, x) where {F <: Function}
    F2 = typeof(x)
    @assert !(F2 <: Function || (F2 <: AbstractArray && F2.parameters[1] <: Function))
    # otherwise we'd hit infinite recursion here
    x, f
end

# --------------------------------
# 3 arguments

# surface-like... function
@recipe function f(x::AVec, y::AVec, zf::Function)
    x, y, Surface(zf, x, y)  # TODO: replace with SurfaceFunction when supported
end

# surface-like... matrix grid
@recipe function f(x::AVec, y::AVec, z::AMat)
    if !is_surface(plotattributes)
        plotattributes[:seriestype] = :contour
    end
    x, y, Surface(z)
end

# parametric functions
# special handling... xmin/xmax with parametric function(s)
@recipe function f(f::Function, xmin::Number, xmax::Number)
    xscale, yscale = [get(plotattributes, sym, :identity) for sym in (:xscale, :yscale)]
    _scaled_adapted_grid(f, xscale, yscale, xmin, xmax)
end
@recipe function f(fs::AbstractArray{F}, xmin::Number, xmax::Number) where {F <: Function}
    xscale, yscale = [get(plotattributes, sym, :identity) for sym in (:xscale, :yscale)]
    unzip(_scaled_adapted_grid.(vec(fs), xscale, yscale, xmin, xmax))
end
@recipe f(
    fx::FuncOrFuncs{F},
    fy::FuncOrFuncs{G},
    u::AVec,
) where {F <: Function, G <: Function} = _map_funcs(fx, u), _map_funcs(fy, u)
@recipe f(
    fx::FuncOrFuncs{F},
    fy::FuncOrFuncs{G},
    umin::Number,
    umax::Number,
    n = 200,
) where {F <: Function, G <: Function} = fx, fy, range(umin, stop = umax, length = n)

function _scaled_adapted_grid(f, xscale, yscale, xmin, xmax)
    (xf, xinv), (yf, yinv) = ((scale_func(s), inverse_scale_func(s)) for s in (xscale, yscale))
    xs, ys = PlotUtils.adapted_grid(yf ∘ f ∘ xinv, xf.((xmin, xmax)))
    xinv.(xs), yinv.(ys)
end

# special handling... 3D parametric function(s)
@recipe function f(
    fx::FuncOrFuncs{F},
    fy::FuncOrFuncs{G},
    fz::FuncOrFuncs{H},
    u::AVec,
) where {F <: Function, G <: Function, H <: Function}
    _map_funcs(fx, u), _map_funcs(fy, u), _map_funcs(fz, u)
end
@recipe function f(
    fx::FuncOrFuncs{F},
    fy::FuncOrFuncs{G},
    fz::FuncOrFuncs{H},
    umin::Number,
    umax::Number,
    numPoints = 200,
) where {F <: Function, G <: Function, H <: Function}
    fx, fy, fz, range(umin, stop = umax, length = numPoints)
end

# list of tuples
@recipe f(v::AVec{<:Tuple}) = unzip(v)
@recipe f(tup::Tuple) = [tup]


text:
#specific multi-process tests of random generators, jobsegments and runners
@test nprocs() > 1 #only run this test with more than 1 process

#get the first worker process
where1 = workers()[1]
println("===============================================")
println("Multi-process tests with ",nprocs()," processes")
println("===============================================")
println("Local process: 1")
println("Remote process: ",where1)
println("========================================")
println("Multi-process tests of random generators")
println("========================================")
#but when the seed has been set in both, they should be the same
rng1loc = srand(345)
rng1rem = remotecall_fetch(srand,where1,345)
@test rng1loc.seed == rng1rem.seed
v1loc = rand()
v1rem = remotecall_fetch(rand,where1)
@test v1loc == v1rem

#each random generator should now be in the same state
v2loc = rand()
v2rem = remotecall_fetch(rand,where1)
@test v2loc == v2rem

#now bring them into different states
rng3loc = srand()
rng3rem = remotecall_fetch(srand,where1)
@test rng3loc.seed != rng3rem.seed
v3loc = rand()
v3rem = remotecall_fetch(rand,where1)
@test v3loc != v3rem

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

println("==============================================")
println("Multi-process tests of sampler state functions")
println("==============================================")
#create a sampler and sampler states
s1loc = sampler(:mh,:normal,0.1,2)
#create an indicatorstate
ind1loc = samplerstate(s1loc,1,Float64,Float64,false)
#create a local and remote auxiliary state
aux1loc = samplerstate(s1loc,3,Float64,Float64,true)
aux1rem = samplerstate(s1loc,3,Float64,Float64,true)

#set the local and remote random generator in the same state
rng4loc = srand(789)
rng4rem = remotecall_fetch(srand,where1,789)
@test rng4loc.seed == rng4rem.seed

#test that the propose! function gives the same results on local and remote
aux2loc = propose!(aux1loc)
aux2rem = remotecall_fetch(propose!,where1,aux1rem)
@test aux2loc.from.values == aux2rem.from.values
@test aux2loc.proposals.values == aux2rem.proposals.values
aux3loc = propose!(aux2loc)
aux3rem = remotecall_fetch(propose!,where1,aux2rem)
@test aux3loc.from.values == aux3rem.from.values
@test aux3loc.proposals.values == aux3rem.proposals.values

#set the from field of auxiliary states via indicator state
propose!(ind1loc)
aux4loc = prepareauxiliary!(ind1loc,aux3loc)
aux4rem = remotecall_fetch(prepareauxiliary!,where1,ind1loc,aux3rem)
@test aux4loc.from.values == ind1loc.proposals.values
@test aux4rem.from.values == ind1loc.proposals.values
@test aux4loc.proposals.values == aux4rem.proposals.values

#repeat propose!
rng5loc = srand(78932)
rng5rem = remotecall_fetch(srand,where1,78932)
aux5loc = propose!(aux4loc)
aux5rem = remotecall_fetch(propose!,where1,aux4rem)
@test aux5loc.from.values == aux5rem.from.values
@test aux5loc.proposals.values == aux5rem.proposals.values

#now set the local and remote random generators to different seeds
rng6loc = srand()
rng6rem = remotecall_fetch(srand,where1)
aux6loc = propose!(aux5loc)
aux6rem = remotecall_fetch(propose!,where1,aux5rem)
@test aux6loc.from.values == aux6rem.from.values
@test aux6loc.proposals.values != aux6rem.proposals.values

################################################################################
println("================================================")
println("Multi-process tests of remote segement functions")
println("================================================")
########Create all required objects
nprops1 = 5
niter1 = 6
nburnin1 = 4
ntunerperiod = 2
time1 = linspace(0.0,10.0,100)
params1 = parameters([:a,:b],[Distributions.Normal(3.0,0.2),Distributions.Normal(4.0,0.2)],[3.0,4.0])
nparams1 = length(params1)
data1 = data(:array,time1,cat(2,sin(3time1),cos(4time1)))
noise1 = noise(:gaussian,[0.01,0.01])
model1 = model(:target,params1,data1,noise1,(p,t)->cat(2,sin(p[1]*t),cos(p[2]*t)),time1;name="MultiProcessTestModel")
policy1 = policy(:mh,nprops1;jobsegments=:one,store=:all)
runner1 = runner(policy1,niter1,nprops1;numburnin =nburnin1)
@test GeneralizedMetropolisHastings._numjobsegments(policy1,nprops1) == 1
@test collect(GeneralizedMetropolisHastings._processnumbers(policy1,1)) == [where1]
tuner1 = tuner(:scale,ntunerperiod,0.5,:erf)
ind1loc = samplerstate(s1loc,1,Float64,Float64,false)
ind1rem = samplerstate(s1loc,1,Float64,Float64,false)
seg1loc = segment(policy1,model1,s1loc,nprops1)
seg1rem = remotesegments(policy1,model1,s1loc,nprops1)
indicator1loc,ts1loc,chain1loc = GeneralizedMetropolisHastings.createcommon(runner1,tuner1,nparams1,nprops1,nprops1)
indicator1rem,ts1rem,chain1rem = GeneralizedMetropolisHastings.createcommon(runner1,tuner1,nparams1,nprops1,nprops1)

#initialize indicator state and generate auxiliary variable
srand(383928) ; initialize!(runner1,model1,ind1loc,chain1loc,true)
srand(383928) ; initialize!(runner1,model1,ind1rem,chain1rem,true)

function testoneiteration(indloc,indrem,segloc,segrem,indicatorloc,indicatorrem,r,m,c,s1,s2,w,tl,tr)
    srand(s1) ; auxiliary!(r,m,indloc) ; indaccloc = acceptance!(indloc)
    srand(s1) ; auxiliary!(r,m,indrem) ; indaccrem = acceptance!(indrem)
    @test from(indloc) == from(indrem)
    @test proposals(indloc) == proposals(indrem)
    @test indaccloc == indaccrem
    srand(s2) ; accloc = iterate!(segloc,indloc)
    remotecall_wait(srand,w,s2) ; accrem = iterate!(segrem,indrem)
    sleep(1.0)
    @test isequal(accloc,accrem)
    transitionprobability!(indicatorloc,indaccloc,accloc) ; sampleindicator!(indicatorloc)
    transitionprobability!(indicatorrem,indaccrem,accrem) ; sampleindicator!(indicatorrem) ; store!(r,indrem,segrem,indicatorrem,c)
    @test indicatorloc.stationary == indicatorrem.stationary
    @test indicatorloc.samples == indicatorrem.samples
    ###Print
    println("Indicator state from: ",from(indrem).values)
    println("Indicator state proposals: ",proposals(indrem).values)
    println("Indicator state acceptance: ",indaccrem)
    println("Segment state acceptance: ",accrem)
    println("Stationary indicator: ",indicatorrem.stationary)
    println("Indicator samples: ",indicatorrem.samples)
    println("Remote segment acceptances: ",segrem.acceptances)
    println("Remote segment collected samples: ")
    dump(segrem.collectedsamples)
    println("Remote segment prop2collected: ")
    show(segrem.prop2collected)
    println()
    ###Prepare for next
    accepted!(tl,indicatorloc) ; indlocnew = preparenext!(r,model1,indloc,segrem,indicatorloc)
    accepted!(tr,indicatorrem) ; indremnew = preparenext!(r,model1,indrem,segrem,indicatorrem)
    @test from(indlocnew) == from(indremnew)
    @test proposals(indlocnew) == proposals(indremnew)
    println("Indicator state from new: ",from(indremnew).values)
    println("Indicator state proposals new: ",proposals(indremnew).values)
    indlocnew,indremnew
end

for i=1:10
    indloc,indrem = testoneiteration(ind1loc,ind1rem,seg1loc,seg1rem,indicator1loc,indicator1rem,runner1,model1,chain1rem,750-i,190+i,where1,ts1loc,ts1rem)
end

println("==========================")
println("Multi-process tests passed")
println("==========================")
nothing


text:
# Adapted from Knet's src/data.jl (author: Deniz Yuret)
using Random: AbstractRNG, shuffle!, GLOBAL_RNG

struct DataLoader{D,R<:AbstractRNG}
    data::D
    batchsize::Int
    nobs::Int
    partial::Bool
    imax::Int
    indices::Vector{Int}
    shuffle::Bool
    rng::R
end

"""
    Flux.DataLoader(data; batchsize=1, shuffle=false, partial=true, rng=GLOBAL_RNG)

An object that iterates over mini-batches of `data`,
each mini-batch containing `batchsize` observations
(except possibly the last one).

Takes as input a single data tensor, or a tuple (or a named tuple) of tensors.
The last dimension in each tensor is the observation dimension, i.e. the one
divided into mini-batches.

If `shuffle=true`, it shuffles the observations each time iterations are re-started.
If `partial=false` and the number of observations is not divisible by the batchsize,
then the last mini-batch is dropped.

The original data is preserved in the `data` field of the DataLoader.

# Examples
```jldoctest
julia> Xtrain = rand(10, 100);

julia> array_loader = Flux.DataLoader(Xtrain, batchsize=2);

julia> for x in array_loader
         @assert size(x) == (10, 2)
         # do something with x, 50 times
       end

julia> array_loader.data === Xtrain
true

julia> tuple_loader = Flux.DataLoader((Xtrain,), batchsize=2);  # similar, but yielding 1-element tuples

julia> for x in tuple_loader
         @assert x isa Tuple{Matrix}
         @assert size(x[1]) == (10, 2)
       end

julia> Ytrain = rand('a':'z', 100);  # now make a DataLoader yielding 2-element named tuples

julia> train_loader = Flux.DataLoader((data=Xtrain, label=Ytrain), batchsize=5, shuffle=true);

julia> for epoch in 1:100
         for (x, y) in train_loader  # access via tuple destructuring
           @assert size(x) == (10, 5)
           @assert size(y) == (5,)
           # loss += f(x, y) # etc, runs 100 * 20 times
         end
       end

julia> first(train_loader).label isa Vector{Char}  # access via property name
true

julia> first(train_loader).label == Ytrain[1:5]  # because of shuffle=true
false

julia> foreach(println∘summary, Flux.DataLoader(rand(Int8, 10, 64), batchsize=30))  # partial=false would omit last
10×30 Matrix{Int8}
10×30 Matrix{Int8}
10×4 Matrix{Int8}
```
"""
function DataLoader(data; batchsize=1, shuffle=false, partial=true, rng=GLOBAL_RNG)
    batchsize > 0 || throw(ArgumentError("Need positive batchsize"))

    n = _nobs(data)
    if n < batchsize
        @warn "Number of observations less than batchsize, decreasing the batchsize to $n"
        batchsize = n
    end
    imax = partial ? n : n - batchsize + 1
    DataLoader(data, batchsize, n, partial, imax, [1:n;], shuffle, rng)
end

@propagate_inbounds function Base.iterate(d::DataLoader, i=0)     # returns data in d.indices[i+1:i+batchsize]
    i >= d.imax && return nothing
    if d.shuffle && i == 0
        shuffle!(d.rng, d.indices)
    end
    nexti = min(i + d.batchsize, d.nobs)
    ids = d.indices[i+1:nexti]
    batch = _getobs(d.data, ids)
    return (batch, nexti)
end

function Base.length(d::DataLoader)
    n = d.nobs / d.batchsize
    d.partial ? ceil(Int,n) : floor(Int,n)
end

_nobs(data::AbstractArray) = size(data)[end]

function _nobs(data::Union{Tuple, NamedTuple})
    length(data) > 0 || throw(ArgumentError("Need at least one data input"))
    n = _nobs(data[1])
    for i in keys(data)
        ni = _nobs(data[i])
        n == ni || throw(DimensionMismatch("All data inputs should have the same number of observations, i.e. size in the last dimension. " *
            "But data[$(repr(first(keys(data))))] ($(summary(data[1]))) has $n, while data[$(repr(i))] ($(summary(data[i]))) has $ni."))
    end
    return n
end

_getobs(data::AbstractArray, i) = data[ntuple(i -> Colon(), Val(ndims(data) - 1))..., i]
_getobs(data::Union{Tuple, NamedTuple}, i) = map(Base.Fix2(_getobs, i), data)

Base.eltype(::DataLoader{D}) where D = D


text:
using TinyTimeModels
using FactCheck

srand(10)
rtol = 0.05

facts("Parameter estimation") do

    n, m = 10000, 5
    σ, η, β = 1, .1, randn(m)
    X = randn(n, m)

    context("Local level innovation and model error variances") do

        y = cumsum(η*randn(n)) + σ*randn(n)
        llfit = fit(y)

        @fact llfit.locallevel_innovation_std --> roughly(η, rtol=rtol)
        @fact llfit.residual_std --> roughly(σ, rtol=rtol)

    end

    context("Variances and regression coefficients") do

        y = X*β + cumsum(η*randn(n)) + σ*randn(n)
        llfit = fit(y, X)

        @fact llfit.regression_coefs --> roughly(β, rtol=rtol)
        @fact llfit.locallevel_innovation_std --> roughly(η, rtol=rtol)
        @fact llfit.residual_std --> roughly(σ, rtol=rtol)

    end

    context("Regression with parametrized covariate transformation") do
        function thresholdify(p::Vector, X::Matrix)
            above_p = X .>= p'
            X_above = X .- p'
            X_below = -X_above
            X_above[!above_p] = 0
            X_below[above_p] = 0
            return [X_above X_below]
        end #thresholdify

        β_new = 10randn(6)
        p = 0.1randn(1)
        X_new = thresholdify(p, X[:,1:3])

        y = X_new*β_new + cumsum(η*randn(n)) + σ*randn(n)
        llfit = fit(y, thresholdify, 0.1randn(1), X[:, 1:3])

        @fact llfit.regression_coefs --> roughly(β_new, rtol=rtol)
        @fact llfit.transform_params --> roughly(p, rtol=rtol)
        @fact llfit.locallevel_innovation_std --> roughly(η, rtol=rtol)
        @fact llfit.residual_std --> roughly(σ, rtol=rtol)

    end

end


text:
export CircularArrayBuffer, sample,
       SumTree,
       CircularTurnBuffer, CircularSARDBuffer, CircularSARDSBuffer, CircularSARDSABuffer,
       EpisodeTurnBuffer, EpisodeSARDBuffer, EpisodeSARDSBuffer, EpisodeSARDSABuffer

include("circular_array_buffer.jl")
include("sum_tree.jl")

include("circular_turn_buffer.jl")
include("episode_turn_buffer.jl")


text:
# ------------ Utilities ------------ #

function trace_retained(vcs::VectorizedCallSite,
                        s::K,
                        ks,
                        o_len::Int,
                        n_len::Int,
                        args::Vector) where K <: UnconstrainedSelection
    w_adj = -sum(map(vcs.trace.subrecords[n_len + 1 : end]) do cl
                     get_score(cl)
                 end)
    new = vcs.trace.subrecords[1 : n_len]
    new_ret = typeof(vcs.ret)(undef, n_len)
    for i in 1 : n_len
        if i in ks
            ss = get_sub(s, i)
            prev_cl = get_sub(vcs, i)
            ret, u_cl, u_w, rd, ds = regenerate(ss, prev_cl, UndefinedChange(), args[i]...)
            new_ret[i] = ret
            new[i] = u_cl
            w_adj += get_score(u_cl) - get_score(prev_cl)
        else
            new_ret[i] = get_ret(get_sub(vcs, i))
        end
    end
    return w_adj, new, new_ret
end

function trace_new(vcs::VectorizedCallSite,
                   s::K,
                   ks,
                   o_len::Int,
                   n_len::Int,
                   args::Vector) where K <: UnconstrainedSelection
    w_adj = 0.0
    new_ret = typeof(vcs.ret)(undef, n_len)
    new = vcs.trace.subrecords
    for i in o_len + 1 : n_len
        ss = get_sub(s, i)
        ret, cl = simulate(call, args[i]...)
        new_ret[i] = ret
        new[i] = cl
    end
    for i in 1 : o_len
        i in ks && begin
            ss = get_sub(s, i)
            prev_cl = get_sub(vcs, i)
            ret, u_cl, u_w, rd, d = regenerate(ss, prev_cl, UndefinedChange(), args[i]...)
            new_ret[i] = ret
            new[i] = u_cl
            w_adj += u_w
            continue
        end
    end
    return w_adj, new, new_ret
end

# ------------ Call sites ------------ #

@inline function (ctx::RegenerateContext{C, T})(c::typeof(plate),
                                                addr::Address,
                                                call::Function,
                                                args::Vector) where {C <: HierarchicalCallSite, T <: HierarchicalTrace}
    vcs = get_prev(ctx, addr)
    n_len, o_len = length(args), length(vcs.args)
    s = get_subselection(ctx, addr)
    _, ks = keyset(s, n_len)
    if n_len <= o_len
        w_adj, new, new_ret = trace_retained(vcs, s, ks, o_len, n_len, args)
    else
        w_adj, new, new_ret = trace_new(vcs, s, ks, o_len, n_len, args)
    end
    add_call!(ctx, addr, VectorizedCallSite{typeof(plate)}(VectorizedTrace(new), get_score(vcs) + w_adj, call, n_len, args, new_ret))
    increment!(ctx, w_adj)

    return new_ret
end

@inline function (ctx::RegenerateContext{C, T})(c::typeof(plate),
                                                addr::Address,
                                                call::Function,
                                                args::Vector) where {C <: VectorizedCallSite, T <: VectorizedTrace}
    vcs = ctx.prev
    n_len, o_len = length(args), length(vcs.args)
    s = get_subselection(ctx, addr)
    _, ks = keyset(s, n_len)
    if n_len <= o_len
        w_adj, new, new_ret = trace_retained(vcs, s, ks, o_len, n_len, args)
    else
        w_adj, new, new_ret = trace_new(vcs, s, ks, o_len, n_len, args)
    end

    # TODO: fix - allocate static vector.
    for n in new
        add_call!(ctx, n)
    end
    increment!(ctx, w_adj)

    return new_ret
end


text:
using AxisArrayVariables
using JuMP: Model, Variable
using AxisArrays: AxisArray, Axis
using Base.Test

@testset "axes as keywords" begin
    m = Model()
    @axis_variable(m, x[time=1:5])
    @test isa(x, AxisArray{Variable, 1, Array{Variable, 1},
                           Tuple{Axis{:time, UnitRange{Int}}}})
    @test size(x, Axis{:time}) == 5
    @test size(x) == (5,)

    m = Model()
    @axis_variable(m, x[time=1:5, side=[:left, :right]])
    @test isa(x, AxisArray{Variable, 2, Array{Variable, 2},
                           Tuple{Axis{:time, UnitRange{Int}},
                                 Axis{:side, Vector{Symbol}}}})
    @test size(x, Axis{:time}) == 5
    @test size(x, Axis{:side}) == 2
    @test size(x) == (5, 2)

    m = Model()
    @axis_variable(m, x[side=[:left, :right], time=1:5])
    @test isa(x, AxisArray{Variable, 2, Array{Variable, 2},
                           Tuple{Axis{:side, Vector{Symbol}},
                                 Axis{:time, UnitRange{Int}}}})
    @test size(x, Axis{:time}) == 5
    @test size(x, Axis{:side}) == 2
    @test size(x) == (2, 5)
end

@testset "axes as variables" begin
    time = Axis{:time}(1:5)
    side = Axis{:side}([:left, :right])

    m = Model()
    @axis_variable(m, x[time])
    @test isa(x, AxisArray{Variable, 1, Array{Variable, 1},
                           Tuple{Axis{:time, UnitRange{Int}}}})
    @test size(x, Axis{:time}) == 5
    @test size(x) == (5,)

    m = Model()
    @axis_variable(m, x[time, side])
    @test isa(x, AxisArray{Variable, 2, Array{Variable, 2},
                           Tuple{Axis{:time, UnitRange{Int}},
                                 Axis{:side, Vector{Symbol}}}})
    @test size(x, Axis{:time}) == 5
    @test size(x, Axis{:side}) == 2
    @test size(x) == (5, 2)

    m = Model()
    @axis_variable(m, x[side, time])
    @test isa(x, AxisArray{Variable, 2, Array{Variable, 2},
                           Tuple{Axis{:side, Vector{Symbol}},
                                 Axis{:time, UnitRange{Int}}}})
    @test size(x, Axis{:time}) == 5
    @test size(x, Axis{:side}) == 2
    @test size(x) == (2, 5)
end

@testset "mixing variables and keywords" begin
    time = Axis{:time}(1:5)

    m = Model()
    @axis_variable(m, x[time, side=[:left, :right]])
    @test isa(x, AxisArray{Variable, 2, Array{Variable, 2},
                           Tuple{Axis{:time, UnitRange{Int}},
                                 Axis{:side, Vector{Symbol}}}})
    @test size(x, Axis{:time}) == 5
    @test size(x, Axis{:side}) == 2
    @test size(x) == (5, 2)
end

@testset "upper and lower bounds" begin
    m = Model()
    @axis_variable(m, x[time=1:5] >= 0)
    @test all(JuMP.getlowerbound(x) .== 0)
    @test all(JuMP.getupperbound(x) .== Inf)

    m = Model()
    @axis_variable(m, x[time=1:5] <= 0)
    @test all(JuMP.getlowerbound(x) .== -Inf)
    @test all(JuMP.getupperbound(x) .== 0)

    m = Model()
    @axis_variable(m, 1 <= x[time=1:5] <= 3)
    @test all(JuMP.getlowerbound(x) .== 1)
    @test all(JuMP.getupperbound(x) .== 3)
end

@testset "variable types" begin
    m = Model()
    @axis_variable(m, x[time=1:5, side=[:left, :right]])
    @axis_variable(m, in_contact[time=1:5, side=[:left, :right]], Bin)
    @test all(JuMP.getcategory.(in_contact) .== :Bin)
    @test all(JuMP.getcategory.(x) .== :Cont)
end

function f()
    m = Model()
    @axis_variable(m, x[time=1:5, side=[:left, :right]])
    x
end

@testset "inference" begin
    @inferred(f())
end


text:
module TestAmbWarning
using PatternDispatch

@pattern f(x::Int, y) = 1
@pattern f(x, y::Int) = 2

@pattern g(1, y) = 1
@pattern g(x, 2) = 2

end


text:
using Test

tests = [
    # "edgelist",
    # "adjlist",
    # "bellman_test",
    # "inclist",
    # "graph",
    # "gmatrix",
    # "bfs",
    # "dfs",
    # "conn_comp",
    # "dijkstra",
    # "a_star_spath",
    # "a_star_implicit_test",
    # "a_star_implicit_constrained_test",
    "a_star_implicit_constrained_epsilon_test",
    # "mst",
    # "floyd",
    # "dot",
    # "dot2",
    # "cliques",
    # "random",
    # "generators",
    # "maximum_adjacency_visit",
    # "issue_related_tests",
    # "inclist_dict_delete"
    ]

testdir = joinpath(dirname(@__DIR__), "test")

for t in tests
    @testset "$(t)" begin
        tp = joinpath(testdir, "$(t).jl")
        include(tp)
    end
end


text:
struct Envelope
    dd_interval_::DataFrame
    dd_evals_::DataFrame
    dd_function_::DataFrame
    x_name_::Symbol
    output_for_envelope_::Symbol
    splines_::Dict{Symbol,Dict{Symbol,Schumaker}}
    chance_of_function_::Multinomial
    funcs_::Array{Symbol,1}
    function Envelope(dd_interval::DataFrame, dd_eval::DataFrame, dd_function::DataFrame, x_name::Symbol, output_for_envelope::Symbol)
         func_names = dd_function[:,:name]
         nams = Symbol.(names(dd_eval))
         val_names  = nams[(nams .!= :name) .& (nams .!= x_name)]
         splines = Dict{Symbol,Dict{Symbol,Schumaker}}()
         for var_name in val_names
             name_d = Dict{Symbol,Schumaker}()
             for func in func_names
                 data = dd_eval[dd_eval[:,:name] .== func,:]
                 if sum(ismissing.(data[:,var_name])) > 0 continue end
                 schum = Schumaker(data[:,x_name], data[:,var_name])
                 name_d[Symbol(func)] = schum
             end
             splines[Symbol(var_name)] = name_d
         end
         #mass_by_func = by(dd_interval, :func, mass = :length => sum)
         mass_by_func = combine(groupby(dd_interval, :func), [:length => (sum => :mass)]...)

         total_mass   = sum(mass_by_func[:,:mass])
         mass_by_func[!,:mass] .= mass_by_func[:,:mass] ./ total_mass
         chance_of_function_ = Multinomial(1, mass_by_func[:,:mass])
         return new(dd_interval, dd_eval, dd_function, x_name, output_for_envelope, splines, chance_of_function_, mass_by_func[:,:func])
    end

    function Envelope(dd_function::DataFrame, grid::Union{Array,StepRangeLen,UnitRange}, output_for_envelope::Symbol,  x_name::Symbol = :x;
                      do_all_at_ends::Bool = true, bracketing_parameter::Real = 0.8, max_interval_width::Real = 1E-08, additional_points_for_all_functions::AbstractArray = Array{typeof(grid[1]),1}([]))
        dd_eval = crawler(dd_function, grid, x_name, output_for_envelope; do_all_at_ends = do_all_at_ends, additional_points_for_all_functions = additional_points_for_all_functions)
        dd_interval, dd_eval = get_upper_envelope(dd_eval, dd_function, x_name, output_for_envelope, bracketing_parameter, max_interval_width; grid = sort(grid))
        return Envelope(dd_interval, dd_eval, dd_function, x_name, output_for_envelope)
    end
end

import SchumakerSpline.evaluate, SchumakerSpline.evaluate_integral

"""
    evaluate_integral(env::Envelope, lhs::Real, rhs::Real, var::Symbol = env.output_for_envelope_)
Evaluates the integral of an envelope between two spots.
"""
function evaluate_integral(env::Envelope, lhs::Real, rhs::Real, variable::Symbol = env.output_for_envelope_)
    first_func_index = max(1, searchsortedlast(env.dd_interval_[!,:interval_start], lhs))
    last_func_index  = searchsortedlast(env.dd_interval_[!,:interval_start], rhs)
    spl_start_name   = env.dd_interval_[first_func_index, :func]
    integral         = 0.0
    if first_func_index == last_func_index
        #  If all is in the same interval.
        integral  += evaluate_integral(env.splines_[variable][spl_start_name], lhs, rhs)
        return integral
    else
        # Start and end components
        integral  += evaluate_integral(env.splines_[variable][spl_start_name], lhs, env.dd_interval_[first_func_index,:interval_end])
        spl_end_name = env.dd_interval_[last_func_index, :func]
        integral  += evaluate_integral(env.splines_[variable][spl_end_name], env.dd_interval_[last_func_index,:interval_start], rhs)
    end
    if last_func_index - first_func_index == 1
        # If start and end components are all they are.
        return integral
    end
    other_intervals = (first_func_index+1):1:(last_func_index-1)
    for i in other_intervals
        spl_name = env.dd_interval_[i, :func]
        integral  += evaluate_integral(env.splines_[variable][spl_name], env.dd_interval_[i,:interval_start], env.dd_interval_[i,:interval_end])
    end
    return integral
end

function top_function_at_point(env::Envelope, x::Real)
    return env.dd_interval_[searchsortedlast(env.dd_interval_[!,:interval_start], x),:func]
end

"""
    evaluate(env::Envelope, x::Real, output_for_envelope::Symbol = env.output_for_envelope_; derivative::Integer = 0)
Evaluates an upper envelope at a particular coordinate. It can also estimate other inputs on the same intervals and derivatives.
"""
function evaluate(env::Envelope, x::Real, output_for_envelope::Symbol = env.output_for_envelope_; derivative::Integer = 0)
    top_spline_name = top_function_at_point(env, x)
    top_spline = env.splines_[output_for_envelope][top_spline_name]
    return evaluate(top_spline, x, derivative)
end
function (env::Envelope)(x::Real, derivative::Integer = 0)
    return evaluate(env, x, env.output_for_envelope_; derivative = derivative)
end

function sample_envelope_given_func(env::Envelope, func_name::Symbol, var::Symbol, twister::MersenneTwister; derivative::Integer = 0)
    eval_quantile = rand(twister,1)[1]
    x, y = sample_envelope_given_func(env, func_name, var, eval_quantile; derivative = 0)
    return x, y
end

function sample_envelope_given_func(env::Envelope, func_name::Symbol, var::Symbol, quantile::Real; derivative::Integer = 0)
    func_dd = env.dd_interval_[env.dd_interval_[!,:func] .== func_name,:]
    func_mass = sum(func_dd[!,:interval_end] .- func_dd[!,:interval_start])
    distance_to_point = func_mass * quantile
    final_x = NaN
    for i in 1:size(func_dd)[1]
        interval_width = func_dd[i,:interval_end] .- func_dd[i,:interval_start]
        if distance_to_point < interval_width
            final_x = func_dd[i,:interval_start] + distance_to_point
            break
        end
        distance_to_point = distance_to_point - interval_width
    end
    final_y = evaluate(env.splines_[var][func_name], final_x, derivative)
    return final_x, final_y
end

function sample_envelope(env::Envelope, var::Symbol, twister::MersenneTwister; derivative::Integer = 0)
    func_name = env.funcs_[findall(rand(twister, env.chance_of_function_) .== 1)[1]]
    x, y = sample_envelope_given_func(env, func_name, var, twister; derivative = 0)
    return func_name, x, y
end


text:
import Images
import TestImages
import ImageMetadata
import ImageAxes
import ImageDraw
import ImageView
import ImageTransformations
import ImageFiltering

ImageFromTestImages = TestImages.testimage("lighthouse")
ImageView.imshow(ImageFromTestImages)

ImageFromTestImagesGray = Images.Gray.(ImageFromTestImages)

shrinkedImage = Images.imresize(ImageFromTestImages, ratio = 0.5)
ImageView.imshow(shrinkedImage)

newSize = div.(size(ImageFromTestImages), 2)
blurValue = 0.1
sigma = map((o,n)->blurValue*o/n, size(ImageFromTestImages), newSize)
gaussFilter = ImageFiltering.KernelFactors.gaussian(sigma)
filtAndResizedImage = Images.imresize(Images.imfilter(ImageFromTestImages, gaussFilter, Images.NA()), newSize)
ImageView.imshow(filtAndResizedImage)


text:
function toh02()
#
#  generalization of K.C. Toh and L. N. Trefethen
#
    k = [(20:-1:0)...]
    p1 = exp10.(k)
    p2 = ones(21)
    for j = 20:-1:0, p2[21-j] = factorial(j); end
    p = p1 ./ p2;
    z = 2 .* (k .+0.5) ./ 19 + im .* sin.((2*pi) .* (k .+ 0.5) ./ 19)
    z = [1.08045842459084 + 0.92291979046775im;
        1.08045842459084 - 0.92291979046775im;
        0.57624370699834 + 0.97554880155012im;
        0.57624370699834 - 0.97554880155012im;
        0.23672740822527 + 0.94133582606706im;
        0.23672740822527 - 0.94133582606706im;
        -0.01684064074376 + 0.86388496107331im;
        -0.01684064074376 - 0.86388496107331im;
        -0.21255469312515 + 0.76040546396802im;
        -0.21255469312515 - 0.76040546396802im;
        -0.36406222857634 + 0.63986739119691im;
        -0.36406222857634 - 0.63986739119691im;
        -0.47903273599273 + 0.50773158889098im;
        -0.47903273599273 - 0.50773158889098im;
        -0.56209888019876 + 0.36770724207462im;
        -0.56209888019876 - 0.36770724207462im;
        -0.64273026408505 + 0.07449713941656im;
        -0.64273026408505 - 0.07449713941656im;
        -0.61611009709266 + 0.22255230874526im;
        -0.61611009709266 - 0.22255230874526im;
        ];
   z = [z ones(20)];
   p, PolyZeros(z)
end


text:
"""
Collection of functions which a practical to use when editing sourcs code, such as
* capitalizing words
* switch between snake case and camel case of an identifier
* toggle case
"""
module EditorUtils

export capitalize, length_of_longest, camel_case, snake_case, toggle, spaced_words, camel_split

"Capitalize strings such as \"HELLO\" to \"Hello\" or \"hello\" to \"Hello\""
function capitalize(s::AbstractString)
    if isempty(s)
        return s
    end
    uppercase(s[1:1]) * lowercase(s[2:end])
end

"Find the length of longest string in list. So for strings `[\"aa\", \"bbb\", \"c\"]` the result would be 3, since thel longest string \"bbb\" is 3 characters long"
function length_of_longest(strings::Vector{<: AbstractString})
    reduce(strings, init = 0) do len, s
        max(len, length(s))
    end
end

"Turns string \"hello_world\" into \"HelloWorld\""
camel_case(s::AbstractString) =  join(capitalize.(split(s, "_")))

const EOFChar = Char(0xC0) # Using illegal UTF8 as sentinel

"Split camel cased word into multiple strings. Turns \"HelloWorldWide\" into [\"hello\", \"world\", \"wide\"]"
function camel_split(s::AbstractString)
    words = String[]
    i = 1       # start position in string
    j = 0       # current position

    function peek()
        k = nextind(s, j)
        if k > lastindex(s)
            return EOFChar
        end
        s[k]
    end

    function scan(pred::Function)
        ch = peek()
        while ch != EOFChar && (pred(ch) || !isletter(ch))
            j += 1
            ch = peek()
        end
        ch
    end

    ch = peek()
    while ch != EOFChar
       ch = if islowercase(ch)
                scan(islowercase)
            else
                scan(isuppercase)
                scan(islowercase)
            end
       @debug "Added word '$(lowercase(s[i:j]))'" maxlog=6
       push!(words, lowercase(s[i:j]))
       i = j + 1
    end
    words
end

"""
    snake_case(string)
Turns string \"HelloWorld\" into \"hello_world\"
"""
snake_case(s::AbstractString) = join(camel_split(s), "_")

"""
    spaced_words(string)
Turns \"HELLO_THERE_FRIENDS\" into \"Hello There Friends\"
"""
spaced_words(s::AbstractString) = join(capitalize.(split(s, "_")), " ")

"""
    camel_case()
Convert string on clipboard in snake case into camel case. Result placed on clipboard.
"""
function camel_case()
    s = camel_case(clipboard())
    clipboard(s)
    s
end

"""
    snake_case()
Convert string on clipboard in camel case into snake case. Result placed on clipboard.
"""
function snake_case()
    s = snake_case(clipboard())
    clipboard(s)
    s
end

"""
    toggle_case(string)
Toggle between camel case and snake case
"""
function toggle_case(s::AbstractString)
    if isempty(search(s, "_")) || any(isuppercase, s)
        snake_case(s)
    else
        camel_case(s)
    end
end

"""
    toggle_case()
Toggle between camel case and snake case on string found on clipboard. Result placed
on clipboard.
"""
function toggle_case()
    s = toggle_case(clipboard())
    clipboard(s)
    s
end

end # module


text:
using BNP_WMReg_Joint

using BSON
import Clustering.hclust, Clustering.cutree
import Distances.pairwise, Distances.Euclidean
import BayesInference.embed, BayesInference.etr
using PDMats
using StatsBase
using Statistics
using Random
using Dates

using LinearAlgebra
BLAS.set_num_threads(1)

## Settings from shell script
const srnd = parse(Int64, ARGS[1])
const n = parse(Int64, ARGS[2])
const K = parse(Int64, ARGS[3])
const H = parse(Int64, ARGS[4])
const Σx_type = Symbol(ARGS[5])
const γ_type = Symbol(ARGS[6])
const γ_init = ARGS[7] # "all" or "none"
const SNR = parse(Float64, ARGS[8])
const whichdat = ARGS[9]

## MCMC Settings (used in article)
# n_burn = Int(300e3)
# n_keep = 5000
# n_adapt = 500
# n_thin = 100
# centerscale = false

## Alternate settings (for shorter runs)
n_burn = Int(200e3)
n_keep = 1000
n_adapt = 500
n_thin = 50
centerscale = false

const mesg = "$(whichdat)_n$(n)_K$(K)_SigX$(Σx_type)_VS$(γ_type)_gaminit$(γ_init)_snr$(Int(floor(SNR)))_$(srnd)"
Random.seed!(srnd)

dat = BSON.load("data/$(whichdat).bson")
y = deepcopy(dat[:y])
X = deepcopy(dat[:X])

y_fit = deepcopy(y[1:n])
X_fit = deepcopy(X[1:n, 1:K])

if centerscale
    X_fit = (X_fit .- mean(y_fit)) ./ std(y_fit) # time series
    y_fit = (y_fit .- mean(y_fit)) ./ std(y_fit)
end

## Customize the prior
range_y = maximum(y_fit) - minimum(y_fit)
center_y = mean(y_fit)
range_X = fill(deepcopy(range_y), K) # assuming time series application
center_X = fill(deepcopy(center_y), K)

prior = Prior_BNP_WMReg_Joint(K, H, range_y = range_y, center_y = center_y,
    range_X = range_X, center_X = center_X, snr=SNR, Σx_type=Σx_type, γ_type=γ_type)

if γ_type == :local
    prior.π_ξ = (0.1 .+ (0.5 .^ collect(1:K) .* 0.4/0.5 ))
end

init = init_state_BNP_WMReg_Joint(n, K, H, prior, random=1, Σx_type=Σx_type, γ_type=γ_type) # random: 0 for none, 1 for lower level params, 2 for all

if γ_init == "all"
    ## default is for all to initialize on
    all(init.γ) || throw("Reset default lag initialization to all on.")
elseif γ_init == "none"
    all(init.γ) || throw("Reset default lag initialization to all on.")
    init.γ = .!init.γ
else
    throw("Invalid initialization for lag selectors.")
end

if γ_type == :global
    init.π_γ = (0.1 .+ (0.5 .^ collect(1:K) .* 0.4/0.5 ))
end

## Clustering algorithm to initialize allocation
if γ_init == "all"
    all(init.γ) || throw("Reset default lag initialization to all on.")

    clust = hclust(pairwise(Euclidean(), transpose(hcat(y_fit, X_fit)), dims=2), linkage=:ward)
    init.S = cutree(clust, k=H)
    init.n_occup = length(unique(init.S))
    println(counts(init.S, 1:H))

elseif γ_init == "none"
    all(.!init.γ) || throw("Reset default lag initialization to all on.")

    clust = hclust(abs.(broadcast(-, y_fit, permutedims(y_fit))), linkage=:ward)
    init.S = cutree(clust, k=H)
    init.n_occup = length(unique(init.S))
    println(counts(init.S, 1:H))

else
    throw("Invalid initialization for lag selectors.")
end

## log file
logfilename = "postsim_progress/out_prog_$(mesg).txt"

## initialize model
model = Model_BNP_WMReg_Joint(y_fit, X_fit, H, prior, init, Σx_type=Σx_type, γ_type=γ_type)

monitor = Monitor_BNP_WMReg_Joint(true, # ηlω
                              false, # γ
                              false, # S
                              false) # G0

updatevars = Updatevars_BNP_WMReg_Joint(true,  # η
                                    true,  # lω
                                    false, # γ
                                    true,  # α
                                    true,  # S
                                    false)  # G0

@time sims, accptr = mcmc!(model, 1000,
    updatevars, monitor, logfilename, 1, 100)
println(model.state.iter)
println(accptr)

adapt!(model, n_iter_collectSS=n_adapt, n_iter_scale=n_adapt,
    adapt_thin=10, updatevars=updatevars,
    report_filename=logfilename, maxtries=10,
    accptr_bnds=[0.02, 0.15], tune_only=true)

if model.γ_type != :fixed
    updatevars.γ = true
end
monitor.γ = true # always on
updatevars.G0 = true

adapt!(model, n_iter_collectSS=n_adapt, n_iter_scale=n_adapt,
    adapt_thin=10, updatevars=updatevars,
    report_filename=logfilename, maxtries=10,
    accptr_bnds=[0.02, 0.15], tune_only=true)

## time to explore posterior
sims, accptr = mcmc!(model, Int(floor(n_burn/10)),
    updatevars, monitor, logfilename, 2, 30000)

report_file = open(logfilename, "a+")
if model.γ_type == :global
    write(logfilename, "\n\nCurrent variable selection rates: $([ mean( [ 1*sims[i][:γ][k] for i=1:length(sims) ] ) for k = 1:model.K ])\n\n")
    println([ mean( [ 1*sims[i][:γ][k] for i=1:length(sims) ] ) for k = 1:model.K ])
elseif model.γ_type == :local
    γ_occupied = counts(model.state.S, 1:model.H)'model.state.γ ./ float(model.n)
    write(logfilename, "Current variable selection of occupied clusters: $(γ_occupied)\n\n")
    println( mean( [ sims[i][:Scounts]'sims[i][:γ] ./ float(model.n) for i=1:length(sims) ] ) )
end
write(report_file, "Pre-burn-in complete. Now adapting with variable selection on.\n\n")
close(report_file)

for ii = 1:2
    ## adapt with variable selection on
    adapt!(model, n_iter_collectSS=n_adapt, n_iter_scale=n_adapt,
        adapt_thin=10, updatevars=updatevars,
        report_filename=logfilename, maxtries=10,
        accptr_bnds=[0.02, 0.20], tune_only=true)

    ## time to explore posterior
    global sims, accptr = mcmc!(model, Int(floor(n_burn/10)),
        updatevars, monitor, logfilename, 1, 10000)

    report_file = open(logfilename, "a+")
    if model.γ_type == :global
        write(logfilename, "\n\nCurrent variable selection rates: $([ mean( [ 1*sims[i][:γ][k] for i=1:length(sims) ] ) for k = 1:model.K ])\n\n")
        println([ mean( [ 1*sims[i][:γ][k] for i=1:length(sims) ] ) for k = 1:model.K ])
    elseif model.γ_type == :local
        γ_occupied = counts(model.state.S, 1:model.H)'model.state.γ ./ float(model.n)
        write(logfilename, "Current variable selection of occupied clusters: $(γ_occupied)\n\n")
        println( mean( [ sims[i][:Scounts]'sims[i][:γ] ./ float(model.n) for i=1:length(sims) ] ) )
    end
    write(report_file, "Adaptation $(ii) of 3 complete.\n\n")
    close(report_file)
end

## burn-in
timestartburn = Dates.now()
sims, accptr = mcmc!(model, Int(floor(n_burn/10)),
    updatevars, monitor, logfilename, 10, 50000)
etr(timestartburn, n_iter_timed=n_burn, n_keep=n_keep, thin=n_thin, outfilename=logfilename)

report_file = open(logfilename, "a+")
if model.γ_type == :global
    write(logfilename, "\n\nCurrent variable selection rates: $([ mean( [ 1*sims[i][:γ][k] for i=1:length(sims) ] ) for k = 1:model.K ])\n\n")
    println([ mean( [ 1*sims[i][:γ][k] for i=1:length(sims) ] ) for k = 1:model.K ])
elseif model.γ_type == :local
    γ_occupied = counts(model.state.S, 1:model.H)'model.state.γ ./ float(model.n)
    write(logfilename, "Current variable selection of occupied clusters: $(γ_occupied)\n\n")
    println( mean( [ sims[i][:Scounts]'sims[i][:γ] ./ float(model.n) for i=1:length(sims) ] ) )
end
write(report_file, "Burn-in complete. Beginning inference samples.\n\n")
close(report_file)

## inference runs
@time global sims, accptr = mcmc!(model, n_keep,
    updatevars, monitor, logfilename, n_thin, 50000)

## save samples
using BSON

if haskey(dat, :truth)

    bson("postsim/postsim_$(mesg).bson",
        model=deepcopy(model), sims=deepcopy(sims),
        init=deepcopy(init), mesg=deepcopy(mesg),
        truth=deepcopy(dat[:truth]) )

else

    bson("postsim/postsim_$(mesg).bson",
        model=deepcopy(model), sims=deepcopy(sims),
        init=deepcopy(init), mesg=deepcopy(mesg) )

end


text:
using Documenter, FixedPointSinCosApproximations

makedocs(;
    modules=[FixedPointSinCosApproximations],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/zsoerenm/FixedPointSinCosApproximations.jl/blob/{commit}{path}#L{line}",
    sitename="FixedPointSinCosApproximations.jl",
    authors="Soeren Zorn",
    assets=String[],
)

deploydocs(;
    repo="github.com/zsoerenm/FixedPointSinCosApproximations.jl",
)


text:
@testset "Delete Method Function Tests" begin

  @test isdefined(Julz, :delete_method) == true

end


text:
# julia
using FASTX
# main
print("collect proteins with length < 300 aa\nusage: julia fasta_collect_putative_effectome.jl input_fasta output_fasta")
reader = open(FASTA.Reader,ARGS[1])
open(ARGS[2],"a") do io
    for record in reader
        if length(FASTA.sequence(record)) < 300
            println(io,record)
        end
    end
    close(reader)
end


text:
################################################################################
#
#	ABSTRACT TYPE
#
#   AbstractEnergyManifold{hamiltonian}
#
#   FILE CONTAINS
#       - abstract type definition
#       - interface definition
#
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
abstract type AbstractEnergyManifold{
    H<:AbstractHamiltonian{L,UC,HB} where {L,UC,HB}
} end

# export the type
export AbstractEnergyManifold

################################################################################
#
#	INTERFACING / ACCESSING ENERGY MANIFOLDS
#	(functions have to be overwritten by concrete types)
#
################################################################################

# obtaining the hamiltonian
function hamiltonian(
            aem :: AEM
        ) :: H where {L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}, AEM <: AbstractEnergyManifold{H}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'hamiltonian' for energy manifold type " * string(typeof(aem)))
end

# export the function
export hamiltonian

# obtaining the k points of the manifold
function kpoints(
            aem :: AEM
        ) :: Vector{Vector{Float64}} where {L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}, AEM <: AbstractEnergyManifold{H}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'kpoints' for energy manifold type " * string(typeof(aem)))
end

# export the function
export kpoints

# obtaining the energy at which the manifold cuts
function energy(
            aem :: AEM
        ) :: Float64 where {L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}, AEM <: AbstractEnergyManifold{H}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'energy' for energy manifold type " * string(typeof(aem)))
end

# export the function
export energy

# obtaining the energy at which the manifold cuts
function numKPoints(
            aem :: AEM
        ) :: Float64 where {L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}, AEM <: AbstractEnergyManifold{H}}

    # return length of the point list
    return length(kpoints(aem))
end

# export the function
export numKPoints

# recalculate the energy manifold (k points)
function recalculate!(
            aem      :: AEM,
            n_points :: Integer
            ;
            resolution :: Union{Int64, Vector{Int64}} = 100
        ) where {L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}, AEM <: AbstractEnergyManifold{H}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'recalculate!' for energy manifold type " * string(typeof(aem)))
end

# export the function
export recalculate!


text:
# See also the file grid/subgrid.jl in BasisFunctions for the definition of
# SubGrid and IndexSubGrid.

"""
    struct MaskedGrid{G,M,I,T} <: SubGrid{T,1}

A MaskedGrid is a subgrid of another grid that is defined by a mask.
The mask is true or false for each point in the supergrid. The set of points
for which it is true make up the MaskedGrid.
"""
struct MaskedGrid{G,M,I,D<:Domain,T} <: SubGrid{T,1}
    supergrid   ::	G
    mask	    ::	M
    indices     ::  Vector{I}
    M           ::	Int				# Total number of points in the mask
    domain      ::  D

    MaskedGrid{G,M,I,D,T}(supergrid::AbstractGrid{T}, mask, indices, domain) where {G,M,I,D,T} =
        new(supergrid, mask, indices, sum(mask), domain)
end
# TODO: In MaskedGrid, perhaps we should not be storing pointers to the points of the underlying grid, but
# rather the points themselves. In that case we wouldn't need to specialize on the type of grid (parameter G can go).

function MaskedGrid(supergrid::AbstractGrid{T}, mask, indices, domain) where {T}
	@assert size(supergrid) == size(mask)

	MaskedGrid{typeof(supergrid),typeof(mask),eltype(indices),typeof(domain),T}(supergrid, mask, indices, domain)
end

# These are for the assignment to indices in the function below.
convert(::Type{NTuple{N,Int}},i::CartesianIndex{N}) where {N} = ntuple(k->i[k],N)

MaskedGrid(supergrid::AbstractGrid, domain::Domain) =
	MaskedGrid(supergrid, in.(supergrid, Ref(domain)), domain)

# MaskedGrid(maskedgrid::MaskedGrid, domain::Domain) =
#     MaskedGrid(supergrid(maskedgrid), mask(maskedgrid) .& in.(supergrid(maskedgrid), domain))

MaskedGrid(supergrid::AbstractGrid, mask, domain) =
    MaskedGrid(supergrid, mask, subindices(supergrid, mask), domain)

function subindices(supergrid, mask::BitArray)
    I = eltype(eachindex(supergrid))
    indices = Array{I}(undef, sum(mask))
    i = 1
    for m in eachindex(supergrid)
       if mask[m]
           indices[i] = m
           i += 1
       end
    end
    indices
end

size(g::MaskedGrid) = (g.M,)

mask(g::MaskedGrid) = g.mask

covering(g::MaskedGrid) = g.domain

subindices(g::MaskedGrid) = g.indices

similar_subgrid(g::MaskedGrid, g2::AbstractGrid) = MaskedGrid(g2, g.mask, g.indices)

# Check whether element grid[i] (of the underlying grid) is in the masked grid.
issubindex(i, g::MaskedGrid) = g.mask[i]

unsafe_grid_getindex(g::MaskedGrid, idx::Int) = @inbounds getindex(g.supergrid, g.indices[idx])

getindex(g::AbstractGrid, idx::BitArray) = MaskedGrid(g, idx, covering(g))

function subgrid(grid::MaskedGrid, domain::Domain)
    submask = in.(supergrid(grid), Ref(domain))
    MaskedGrid(supergrid(grid), submask .& mask(grid), domain)
    # points = grid.points[mask]
    # ScatteredGrid(points)
end


text:
using Test
using LinearAlgebra
using SparseArrays
using ProximalOperators

@testset "Epicompose (Gram-diagonal)" begin

n = 5

for R in [Float64] # TODO: enable Float32?
    for T in [R, Complex{R}]
        A = randn(T, n, n)
        Q, _ = qr(A)
        mu = R(2)
        L = mu*Q
        Lfs = [
            (L, NormL1(R(1))),
            (sparse(L), NormL1(R(1))),
        ]
        for (L, f) in Lfs
            g = Epicompose(L, f, mu)
            x = randn(T, n)
            prox_test(g, x, R(2))
        end
    end
end

end

@testset "Epicompose (Quadratic)" begin

n = 5
m = 3

for R in [Float64] # TODO: enable Float32?
    W = randn(R, n, n)
    Q = W' * W
    q = randn(R, n)
    L = randn(R, m, n)
    Lfs = [
        (L, Quadratic(Q, q)),
        (sparse(L), Quadratic(sparse(Q), q)),
    ]
    for (L, f) in Lfs
        g = Epicompose(L, f)
        x = randn(R, m)
        prox_test(g, x, R(2))
    end
end

end


text:
const obsnames = ["Energy", "Energy^2", "Specific Heat",
                  "|Magnetization x|", "Magnetization x^2", "Magnetization x^4",
                  "Susceptibility x", "Connected Susceptibility x", "Binder Ratio x",
                  "|Magnetization y|", "Magnetization y^2", "Magnetization y^4",
                  "Susceptibility y", "Connected Susceptibility y", "Binder Ratio y",
                  "|Magnetization|", "|Magnetization|^2", "|Magnetization|^4",
                  "Susceptibility", "Connected Susceptibility", "Binder Ratio",
                 ]
const Ns = [2,3,4]
const J = 1.0
const Ts = [0.3, 1.0, 3.0, 10.0]

"""
exact finite-T energy (per site) of XY model on `N` site fully connected network
"""
function integrate(J, N, Ts; Q::Integer=50)
    j = J/N
    Es = zeros(0)
    Ms = Vector{Float64}[]
    spins = zeros(Int, N)
    nst = Q^N
    for ist in 0:(nst-1)
        spins .= 1
        for i in 1:N
            spins[i] = mod(ist,Q) + 1
            ist = div(ist,Q)
        end
        E = 0.0
        M = zeros(2)
        nqs = [count(s->s==q, spins) for q in 1:Q]
        for q in 1:Q
            nq = nqs[q]
            E -= 0.5j*(nq*(nq-1))
            for q2 in (q+1):Q
                nq2 = nqs[q2]
                E -= j*nq*nq2*cospi(2*(q-q2)/Q)
            end
            M[1] += nq*cospi(2*q/Q)
            M[2] += nq*sinpi(2*q/Q)
        end
        M ./= N
        push!(Es, E)
        push!(Ms, M)
    end
    sorted = sortperm(Es, rev=true)
    E0 = Es[sorted[end]]

    mbetas = Ts .\ (-1.0)
    nT = length(Ts)
    res = Dict(n=>zeros(nT) for n in obsnames)
    res["Magnetization x"] = zeros(nT)
    res["Magnetization y"] = zeros(nT)
    Zs = zeros(nT)
    for i in sorted
        E = Es[i] - E0
        zs = exp.(mbetas.*E)
        Zs .+= zs
        E += E0
        res["Energy"] .+= zs.*(E/N)
        res["Energy^2"] .+= zs.*(E/N)^2

        M = Ms[i]
        for (d,c) in ((1,"x"),(2,"y"))
            m = M[d]
            res["Magnetization $(c)"] .+= zs.*m
            res["|Magnetization $(c)|"] .+= zs.*abs(m)
            res["Magnetization $(c)^2"] .+= zs.*m^2
            res["Magnetization $(c)^4"] .+= zs.*m^4
        end
        m2 = sum(abs2,M)
        res["|Magnetization|"] .+= zs.*sqrt(m2)
        res["|Magnetization|^2"] .+= zs.*m2
        res["|Magnetization|^4"] .+= zs.*m2^2
    end
    res["Magnetization x"] ./= Zs
    res["Magnetization y"] ./= Zs
    for n in obsnames
        res[n] ./= Zs
    end
    return res
end

betas = 1.0 ./ Ts

for N in Ns
    res = integrate(J,N,Ts)
    res["Specific Heat"] = N.*(res["Energy^2"] .- res["Energy"].^2).*betas.^2
    for c in ("x","y")
        res["Susceptibility $c"] = N.*(res["Magnetization $c^2"] .- res["Magnetization $c"].^2).*betas
        res["Connected Susceptibility $c"] = N.*(res["Magnetization $c^2"] .- res["|Magnetization $c|"].^2).*betas
        res["Binder Ratio $c"] = res["Magnetization $c^4"] ./ (res["Magnetization $c^2"].^2)
    end
    res["Susceptibility"] = N.*res["|Magnetization|^2"].*betas
    res["Connected Susceptibility"] = N.*(res["|Magnetization|^2"] .- res["|Magnetization|"].^2).*betas
    res["Binder Ratio"] = res["|Magnetization|^4"] ./ (res["|Magnetization|^2"].^2)
    open(@sprintf("J_%.1f__N_%d.dat", J, N), "w") do io
        for i in 1:length(Ts)
            @printf(io, "%.15f", Ts[i])
            for n in obsnames
                @printf(io, " %.15f", res[n][i])
            end
            println(io)
        end
    end
end


text:
using CalculiX
using Base.Test

# write your own tests here
@test 1 == 1
@test 1. == 1.


text:
using CSV
using Plots
using GR

#cd("D:\\test_data\\")
dat = CSV.read("tmp_table.txt", header=false)
#length(dat.Column14)
# using DataFrames
# nrow(dat)
# ncol(dat)
#a =indexin([12,147,240], dat.Column15)

# show the third column
# dat[!,5]
#dat[1:3, [1,5]]
#dat[!, [1,5]]

#last(dat, 3)
#first(dat, 3)
#size(dat)

x=1:22

gr()
plot!(x, dat.Column15, label="line")
scatter!(x, dat.Column15, label="points")

png("./test_a")


text:
module WriteToJSON

include("./../constraints/geometric/constants.jl")
using JuMP
using DelimitedFiles
using Dates
using JSON
using JSON3
using .Constants

export write_to_json_simple

function write_to_json_simple(m)

    for i in 1:length(boxes)-1
        print(boxes[i])
    end

    positions = []
    dimensions = []
    for i in 1:length(boxes)-1
        if (value((m[:lpx_length_of_box_parallel_to_x])[i]) == 1)
            push!(positions, [
                boxes[i][1]/2 + value(m[:flbx_front_left_bot_x][i]),
                boxes[i][3]/2 + value(m[:flbz_front_left_bot_z][i]),
                boxes[i][2]/2 + value(m[:flby_front_left_bot_y][i])
                ]
            )
            push!(dimensions,
                [
                    boxes[i][2],
                    boxes[i][3],
                    boxes[i][1]
                ]
            )
        else
            push!(positions, [
                boxes[i][2]/2 + value(m[:flbx_front_left_bot_x][i]),
                boxes[i][3]/2 + value(m[:flbz_front_left_bot_z][i]),
                boxes[i][1]/2 + value(m[:flby_front_left_bot_y][i])
                ]
            )
            push!(dimensions,
                [
                    boxes[i][1],
                    boxes[i][3],
                    boxes[i][2]
                ]
            )
        end
    end

    inside_json = []
    for i in 1:length(boxes)-1
        dim = Dict(
            "w" => dimensions[i][1],
            "h" => dimensions[i][2],
            "d" => dimensions[i][3],
        )
        pos = Dict(
            "x" => positions[i][1],
            "z" => positions[i][2],
            "y" => positions[i][3]

        )
        push!(inside_json, Dict(
                "dimensions" => dim,
                "position" => pos
            )
        )
    end

    open(resultjson, "w") do io
        JSON3.pretty(io, inside_json)
    end

end

end


text:
"""
An [`Embed`](@ref)'s thumbnail image information.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object-embed-thumbnail-structure).
"""
struct EmbedThumbnail
    url::Union{String, Missing}
    proxy_url::Union{String, Missing}
    height::Union{Int, Missing}
    width::Union{Int, Missing}
end
@boilerplate EmbedThumbnail :dict :docs :lower :merge

"""
An [`Embed`](@ref)'s video information.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object-embed-video-structure).
"""
struct EmbedVideo
    url::Union{String, Missing}
    height::Union{Int, Missing}
    width::Union{Int, Missing}
end
@boilerplate EmbedVideo :dict :docs :lower :merge

"""
An [`Embed`](@ref)'s image information.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object-embed-image-structure).
"""
struct EmbedImage
    url::Union{String, Missing}
    proxy_url::Union{String, Missing}
    height::Union{Int, Missing}
    width::Union{Int, Missing}
end
@boilerplate EmbedImage :dict :docs :lower :merge

"""
An [`Embed`](@ref)'s provider information.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object-embed-provider-structure).
"""
struct EmbedProvider
    name::Union{String, Missing}
    url::Union{String, Nothing, Missing}  # Not supposed to be nullable.
end
@boilerplate EmbedProvider :dict :docs :lower :merge

"""
An [`Embed`](@ref)'s author information.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object-embed-author-structure).
"""
struct EmbedAuthor
    name::Union{String, Missing}
    url::Union{String, Missing}
    icon_url::Union{String, Missing}
    proxy_icon_url::Union{String, Missing}
end
@boilerplate EmbedAuthor :dict :docs :lower :merge

"""
An [`Embed`](@ref)'s footer information.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object-embed-footer-structure).
"""
struct EmbedFooter
    text::String
    icon_url::Union{String, Missing}
    proxy_icon_url::Union{String, Missing}
end
@boilerplate EmbedFooter :dict :docs :lower :merge

"""
An [`Embed`](@ref) field.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object-embed-field-structure).
"""
struct EmbedField
    name::String
    value::String
    inline::Union{Bool, Missing}
end
@boilerplate EmbedField :dict :docs :lower :merge

"""
A [`Message`](@ref) embed.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object).
"""
struct Embed
    title::Union{String, Missing}
    type::Union{String, Missing}
    description::Union{String, Missing}
    url::Union{String, Missing}
    timestamp::Union{DateTime, Missing}
    color::Union{Int, Missing}
    footer::Union{EmbedFooter, Missing}
    image::Union{EmbedImage, Missing}
    thumbnail::Union{EmbedThumbnail, Missing}
    video::Union{EmbedVideo, Missing}
    provider::Union{EmbedProvider, Missing}
    author::Union{EmbedAuthor, Missing}
    fields::Union{Vector{EmbedField}, Missing}
end
@boilerplate Embed :dict :docs :lower :merge


text:
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

# Method (1) real
n, t=10, 100
X=genDataMatrix(n, t)
C=(X*X')/t
wC=whitening(Hermitian(C); simple=true)
# or, shortly
wC=whitening(ℍ(C); simple=true)

# Method (1) complex
Xc=genDataMatrix(ComplexF64, n, t)
Cc=(Xc*Xc')/t
wCc=whitening(Hermitian(Cc); simple=true)

# Method (2) real
wX=whitening(X; simple=true)
@test wC.F'*C*wC.F≈I
@test wX.F'*C*wX.F≈I
@test wX≈wC

# Method (2) complex
wXc=whitening(Xc; simple=true)
@test wCc.F'*Cc*wCc.F≈I
@test wXc.F'*Cc*wXc.F≈I
@test wXc≈wCc

# Method (3) real
k=10
Xset=[genDataMatrix(n, t) for i=1:k]

# whitening on the average covariance matrix
w=whitening(Xset)

# ... selecting subspace dimension allowing an explained variance = 0.5
w=whitening(Xset; eVar=0.5)

# ... averaging the covariance matrices using the logEuclidean metric
w=whitening(Xset; metric=logEuclidean, eVar=0.5)

# ... giving weights `w` to the covariance matrices
w=whitening(Xset; metric=logEuclidean, w=abs2.(randn(k)), eVar=0.5)

# ... subtracting the mean
w=whitening(Xset; meanX=nothing, metric=logEuclidean, w=abs2.(randn(k)), eVar=0.5)

# whitening on the average of the covariance matrices computed along dims 1
w=whitening(Xset; dims=1)

# explained variance
w.eVar

# name of the filter
w.name

using Plots
# plot regularized accumulated eigenvalues
plot(w.arev)

# plot the original covariance matrix and the whitened covariance matrix
 Cmax=maximum(abs.(C));
 h1 = heatmap(C, clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="C");
 D=wC.F'*C*wC.F;
 h2 = heatmap(D, clim=(0, 1), yflip=true, c=:amp, title="F'*C*F");
 📈=plot(h1, h2, size=(700, 300))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigWhitening.png")

# Method (3) complex
k=10
Xcset=[genDataMatrix(ComplexF64, n, t) for i=1:k]

# whitening on the average covariance matrix
wc=whitening(Xcset)

# ... selecting subspace dimension allowing an explained variance = 0.5
wc=whitening(Xcset; eVar=0.5)

# ... averaging the covariance matrices using the logEuclidean metric
wc=whitening(Xcset; metric=logEuclidean, eVar=0.5)

# ... giving weights `w` to the covariance matrices
wc=whitening(Xset; metric=logEuclidean, w=abs2.(randn(k)), eVar=0.5)

# ... subtracting the mean
wc=whitening(Xcset; meanX=nothing, metric=logEuclidean, w=abs2.(randn(k)), eVar=0.5)

# whitening on the average of the covariance matrices computed along dims 1
wc=whitening(Xcset; dims=1)

# explained variance
wc.eVar

# name of the filter
wc.name


text:
module Scratchy

module Utils

export clip, decision, accuracy

function clip(x::Real, low=1e-15, hi=1 - 1e-15)
    if x < low
        return low
    elseif x > hi
        return hi
    else
        return x
    end
end

function decision(logits::Array)
	# 0.5 decision threshold
	return Int.(logits .>= 0.5)
end

function accuracy(y::Array, y_hat::Array)
	# Compute accuracy
    return sum(y .== y_hat) / size(y, 1)
end

end  # Utils

module LayerOps

export Operation

struct Operation
    name::String       # The name of this layer operation
    forward::Function  # The function itself
    backward::Function # The derivative of this function
end

end  # LayerOps

module Activations

using ..LayerOps

export Sigmoid, ReLU, Tanh

# Sigmoid (Logistic)
function sigmoid(h::Array{Float64})
    return @. 1.0 / (1.0 + exp(-h))
end

function d_sigmoid(h::Array{Float64})
    forward_vals = sigmoid(h)
    return @. forward_vals * (1.0 - forward_vals)
end

Sigmoid = Operation("Sigmoid", sigmoid, d_sigmoid)

# ReLU
function relu(x)
    return replace(a -> a<0 ? 0 : a, x)
end

function d_relu(x)
    return replace(a -> a<0 ? 0 : 1, x)
end

ReLU = Operation("ReLU", relu, d_relu)

# Tanh
function tanh(x)
    ex = exp.(x)
    enx = exp.(-x)
    return @. (ex - enx) / (ex + enx)
end

function d_tanh(x)
    return @. 1 - (tanh(x)^2)
end

Tanh = Operation("Tanh", tanh, d_tanh)

end  # Activations

module Losses

using ..LayerOps
import ..Utils

export BinaryCrossEntropy

# Binary cross entropy
function binary_crossentropy(y::Array, o::Array{Float64})
    # y: the gold-standard labels 0 or 1
    # o: the real-valued output [0,1]
    n = size(o, 1)
    y_hat = Utils.clip.(o)
    ces = @. y * log(o) + (1 - y) * log(1 - o)
    return -(1 / n) * sum(ces)
end

function d_binary_crossentropy(y::Array, o::Array{Float64})
    # y: the gold-standard labels 0 or 1
    # o: the real-valued output [0,1]
    y_hat = Utils.clip.(o)
    return @. (o - y) / (o * (1 - o))
end

BinaryCrossEntropy = Operation("BinaryCrossEntropy", binary_crossentropy,
                               d_binary_crossentropy)

end  # Losses

module Layers

using ..LayerOps
import Random

export Dense, forward!

# Dense feed-forward layer
mutable struct Dense
    name::String            # What we'll call this layer.
    W::Array{Float64}       # The weight matrix
    b::Array{Float64}       # The bias term
    σ::Operation            # The activation function
    δ::Array{Float64}       # Used in backprop
    input::Array{Float64}   # The inputs from the previous layer
    z::Array{Float64}       # Weighted inputs XW + b
    output::Array{Float64}  # The output of this layer, i.e. the actvations

    # Constructor v1: When we have already initialized W and b
    function Dense(name::String, W::Array{Float64},
				   b::Array{Float64}, σ::Operation)
        if !isa(σ, Operation)
            error("Activation function '$σ' in layer '$name' is not an Operation")
        end
        if size(W, 2) != size(b, 2)
            error("Incompatible W and b shapes in layer '$name'")
        else
            # Since δ, input, and output are computed later,
			# we init empty arrays.
            new(name, W, b, σ, [], [], [], [])
        end
    end

    # Constructor v2: When we want to randomly initialize W and
	# b given the shape of W.
    function Dense(name::String, size::Tuple{Int,Int}, σ::Operation)
        W = Random.rand(Float64, size)
        b = Random.rand(Float64, (1, size[2]))
        # Since δ, input, and output are computed later, we init as empty arrays.
        new(name, W, b, σ, [], [], [], [])
    end
end

# Dense layer forward op
function forward!(l::Dense, X::Array)
    # It's always a good idea to double check that the shapes are compatible.
    if size(X, 2) != size(l.W, 1)
        error("Input shape $(size(X)) incompatible with layer '$(l.name)' shape $(size(l.W))")
    end
    l.input = X
    l.z = (X * l.W) .+ l.b  # Weighted inputs
    l.output = l.σ.forward(l.z)
    return l.output
end

end  # Layers

module Models

using ..LayerOps
import ..Layers

export Network, forward!, backward!, fit!, display

struct Network
    name::String
    layers::Array
    loss::Operation
    η::Real  # The learning rate

    # Constructor: Check that all layer shapes are compatible
    function Network(name::String, layers::Array, loss::Operation, η::Real)
        l1 = layers[1]
        for l2 in layers[2:end]
            l1_outdim = size(l1.W, 2)
            l2_indim = size(l2.W, 1)
            if l1_outdim != l2_indim
                error("""Output dimension $(size(l1.W, 2)) of layer
                      '$(l1.name)' incompatible with input dimension
                      $(size(l2.W, 1)) of layer '$(l2.name)'.""")
            end
        end
        new(name, layers, loss, η)
    end
end

# Network forward
function forward!(network::Network, X::Array)
    layer_input = X
    for l in network.layers
        layer_input = Layers.forward!(l, layer_input)
    end
    return layer_input
end

# Helper function for backward!
function compute_δ(network::Network, y::Array, i::Int)
    l = network.layers[i]
    dσ_dz = l.σ.backward(l.z)
    # The base case of the iteration: the last layer in the network.
    if i == length(network.layers)
        o = network.layers[i].output
        δ = network.loss.backward(y, o) .* dσ_dz
        return δ
    end
    l_next = network.layers[i+1]
    # l_next.δ should have already been computed because we're iterating
    # backwards through the network's layers in backward!() below.
    δ = (l_next.δ * transpose(l_next.W)) .* dσ_dz
    return δ
end

# Network backprop
function backward!(net::Network, y::Array)
    num_layers = length(net.layers)
    for i in num_layers:-1:1  # Iterate backwards through the network
        l = net.layers[i]
        l.δ = compute_δ(net, y, i)
        ∂C_∂W = transpose(l.δ) * l.input
        l.W = l.W .- (net.η * transpose(∂C_∂W))
    end
end

# Model estimation
function fit!(net::Network, X::Array, y::Array, epochs::Int)
    for i=1:epochs
        o = forward!(net, X)
        backward!(net, y)
        if i == 1 || i % 100 == 0
            loss = net.loss.forward(y, o)
            println("Loss ($i): $(loss)")
        end
    end
end

# Summarize a network's structure.
import Base.display
function display(nn::Network)
    println("Network: $(nn.name)")
    for layer in nn.layers
        println("$(layer.name) : $(size(layer.W)) : $(layer.σ.name)")
    end
    println("Loss: $(nn.loss.name)")
end

end  # Models

end  # Scratchy

max_stars_repo_name: jvasilakes/scratchy

max_issues_repo_name: jvasilakes/scratchy

max_forks_repo_name: jvasilakes/scratchy


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 4133e4a2-ed4d-11ea-11cd-0da48c09c2e8
using Pkg, DrWatson

# ╔═╡ 54bd3938-ed4d-11ea-29e4-c1c7020c1400
begin
	@quickactivate "SR2TuringPluto"
	using StatisticalRethinking
end

# ╔═╡ 0c5fe59e-ed4f-11ea-35a8-45c3aff83362
md"## Broadcasting.jl"

# ╔═╡ 54ccfc74-ed4d-11ea-366d-b738d54ff79a
p_grid = range(0, 1, length = 5)

# ╔═╡ 54d9bce8-ed4d-11ea-279e-3bc3bf465e28
prior = ones(5)

# ╔═╡ 54da91d4-ed4d-11ea-04f5-e9a4f8c659a3
likelihood = pdf.(Binomial.(9, p_grid), 6)

# ╔═╡ 54e45ff4-ed4d-11ea-0288-e3475c709007
posterior = likelihood .* prior

# ╔═╡ 54e8f906-ed4d-11ea-2fee-316c2df1d8f6
posterior ./= sum(posterior)

# ╔═╡ 54f1592c-ed4d-11ea-0210-0be7741d2604
md"##### Broadcasting is another important difference between Julia and R.

You don't have automatic
broadcasting (making functions work the same whether you run it on a single
number or on an array of numbers). However, fixing this is very easy,
you only have to annotate your function call with a dot:
  `f(single_number)` -> `f.(array_of_numbers)`.

That means the translation of the line `likelihood = pdf.(Binomial.(9, p_grid), 6)` goes something like:
1. First or every value in `p_grid` make a binomial distribution with that p value.
2. For every distribution then take the pdf of 6.

The same in the next line, this is elementwise multiplication.

For more infos read [this](https://julialang.org/blog/2017/01/moredots)."

# ╔═╡ e822d97e-ed4f-11ea-2c18-17ecfc81fdbd
md"## End of broadcasting.jl"

# ╔═╡ Cell order:
# ╟─0c5fe59e-ed4f-11ea-35a8-45c3aff83362
# ╠═4133e4a2-ed4d-11ea-11cd-0da48c09c2e8
# ╠═54bd3938-ed4d-11ea-29e4-c1c7020c1400
# ╠═54ccfc74-ed4d-11ea-366d-b738d54ff79a
# ╠═54d9bce8-ed4d-11ea-279e-3bc3bf465e28
# ╠═54da91d4-ed4d-11ea-04f5-e9a4f8c659a3
# ╠═54e45ff4-ed4d-11ea-0288-e3475c709007
# ╠═54e8f906-ed4d-11ea-2fee-316c2df1d8f6
# ╟─54f1592c-ed4d-11ea-0210-0be7741d2604
# ╟─e822d97e-ed4f-11ea-2c18-17ecfc81fdbd


text:
# ensure the requisite libraries are available

const noNemo = "Nemo.jl is not found:\n  Pkg.rm(\"Nemo\"); Pkg.add(\"Nemo\"); quit()\n  Pkg.rm(\"ArbFloats\");Pkg.add(\"ArbFloats\");"
const reNemo = "Nemo.jl is not as expected:\n  Pkg.rm(\"Nemo\"); Pkg.add(\"Nemo\"); quit()\n  Pkg.rm(\"ArbFloats\");Pkg.add(\"ArbFloats\");"

function package_directory(pkgName::String)
    pkgdir = Base.find_in_path(pkgName)
    nothing == pkgdir && throw(ErrorException(noNemo))
    return abspath(joinpath(split(pkgdir, pkgName)[1], pkgName))
end

function library_filepath(libsdir::String, filenames::Vector{String}, libname::String)
    libfile = filenames[ findfirst([startswith(x,libname) for x in filenames]) ]
    return joinpath( libsdir, libfile )
end

NemoLibsDir = abspath(joinpath( package_directory("Nemo"), "local/lib"))
libFiles = readdir(NemoLibsDir)

libarb   = library_filepath( NemoLibsDir, libFiles, "libarb"  )
libflint = library_filepath( NemoLibsDir, libFiles, "libflint")

isfile(libarb)   || throw(ErrorException(reNemo))
isfile(libflint) || throw(ErrorException(reNemo))

# prepare the libraries for use

@static if (is_linux() || is_bsd())
    libarb = String(split(libarb,".so")[1])
    libflint = String(split(libflint,".so")[1])
end
@static if is_apple()
    libarb = String(split(libarb,".dynlib")[1])
    libflint = String(split(libflint,".dynlib")[1])
end
@static if is_windows()
    libarb = String(split(libarb,".dll")[1])
    libflint = String(split(libflint,".dll")[1])
end

macro libarb(sym)
    (:($sym), libarb)
end

macro libflint(sym)
    (:($sym), libflint)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Preconditions must be fulfilled before an operation (update, delete, etc.) is carried out.

    IoK8sApimachineryPkgApisMetaV1Preconditions(;
        resourceVersion=nothing,
        uid=nothing,
    )

    - resourceVersion::String : Specifies the target ResourceVersion
    - uid::String : Specifies the target UID.
"""
mutable struct IoK8sApimachineryPkgApisMetaV1Preconditions <: SwaggerModel
    resourceVersion::Any # spec type: Union{ Nothing, String } # spec name: resourceVersion
    uid::Any # spec type: Union{ Nothing, String } # spec name: uid

    function IoK8sApimachineryPkgApisMetaV1Preconditions(;resourceVersion=nothing, uid=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1Preconditions, Symbol("resourceVersion"), resourceVersion)
        setfield!(o, Symbol("resourceVersion"), resourceVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1Preconditions, Symbol("uid"), uid)
        setfield!(o, Symbol("uid"), uid)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1Preconditions

const _property_map_IoK8sApimachineryPkgApisMetaV1Preconditions = Dict{Symbol,Symbol}(Symbol("resourceVersion")=>Symbol("resourceVersion"), Symbol("uid")=>Symbol("uid"))
const _property_types_IoK8sApimachineryPkgApisMetaV1Preconditions = Dict{Symbol,String}(Symbol("resourceVersion")=>"String", Symbol("uid")=>"String")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1Preconditions }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1Preconditions))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1Preconditions }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1Preconditions[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1Preconditions }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1Preconditions[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1Preconditions)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgApisMetaV1Preconditions }, name::Symbol, val)
end


text:
abstract type UnitaryOperator{T} <: Operator{T} end

inv(Q::UnitaryOperator) = Q'
transpose(Q::UnitaryOperator{T}) where {T<:Real} = Q'

ldiv_coefficients(Q::UnitaryOperator,v::AbstractVector;opts...) =
    mul_coefficients(Q',v;opts...)

struct HessenbergUnitary{uplo,T} <: UnitaryOperator{T}
    sign::Bool
    c::Vector{T}
    s::Vector{T}
    c∞::T
    s∞::T
    band::Int

    function HessenbergUnitary{uplo,T}(sgn::Bool,c::Vector{T},s::Vector{T},
                                       c∞::T,s∞::T,bnd::Int) where {uplo,T}
        @assert isapprox(s∞^2+c∞^2,1)
        @assert length(c)==length(s)+1

        @assert isapprox(abs(first(c)),1)

        @inbounds for k=2:length(c)
            @assert -100eps() ≤ c[k]^2+s[k-1]^2-1 ≤ 100eps()
        end

        new(sgn,c,s,c∞,s∞,bnd)
    end
end

function HessenbergUnitary(::Type{Val{uplo}},sign,c,s,c∞,s∞,band) where {uplo}
    @assert uplo=='L' || uplo=='U'
    HessenbergUnitary{uplo,promote_type(eltype(c),eltype(s),
                                           typeof(c∞),typeof(s∞))}(sign,c,s,c∞,s∞,band)
end

# we choose colstop based on s for simplicity
function hessenberg_colstop(s, s∞, j)
    n = length(s)
    tol = eps()

    # Compute the bandwidth of the matrix
    k=j
    cur = one(eltype(s))
    while abs(cur) > tol
        cur *= k ≤ n ? s[k] : s∞
        k += 1
    end
    k
end

function HessenbergUnitary(::Type{Val{uplo}},sign, c, s, c∞, s∞) where {uplo}
    @assert isapprox(s∞^2+c∞^2,1)
    @assert length(c)==length(s)+1

    @assert isapprox(abs(first(c)),1)

    @inbounds for k=2:length(c)
        @assert -100eps() ≤ c[k]^2+s[k-1]^2-1 ≤ 100eps()
    end

    band = mapreduce(j -> hessenberg_colstop(s, s∞, j)-j, max, 1:length(c))

    HessenbergUnitary(Val{uplo},sign,c,s,c∞,s∞,band)
end

adjoint(Q::HessenbergUnitary{'L',T}) where {T<:Real} =
    HessenbergUnitary(Val{'U'},Q.sign,Q.c,Q.s,Q.c∞,Q.s∞,Q.band)
adjoint(Q::HessenbergUnitary{'U',T}) where {T<:Real} =
    HessenbergUnitary(Val{'L'},Q.sign,Q.c,Q.s,Q.c∞,Q.s∞,Q.band)

bandwidths(Q::HessenbergUnitary{'L'}) = Q.band,1
bandwidths(Q::HessenbergUnitary{'U'}) = 1,Q.band

domainspace(::HessenbergUnitary) = ℓ⁰
rangespace(::HessenbergUnitary) = ℓ⁰

hc(c,c∞,k) = k ≤ length(c) ? c[k] : c∞
hs(s,s∞,k) = k ≤ length(s) ? s[k] : s∞

hc(Q::HessenbergUnitary,k) = hc(Q.c,Q.c∞,k)
hs(Q::HessenbergUnitary,k) = hs(Q.s,Q.s∞,k)

getindex(Q::HessenbergUnitary{'L'},k::Integer,j::Integer) =
    hessuni_getindex(Q.sign,Q.c,Q.s,Q.c∞,Q.s∞,j,k)

getindex(Q::HessenbergUnitary{'U'},k::Integer,j::Integer) =
    hessuni_getindex(Q.sign,Q.c,Q.s,Q.c∞,Q.s∞,k,j)

function hessuni_getindex(sgn::Bool,c::AbstractVector{T},s::AbstractVector{T},
                             c∞::T,s∞::T,
                             k::Integer,j::Integer) where T
    si = sgn ? 1 : -1

    if k>j+1
        zero(T)
    elseif k≥2 && j ==k-1
        -si*hs(s,s∞,k-1)
    else
        col0=hc(c,c∞,k)*hc(c,c∞,k+1)
        for p=k+1:j
            col0*=hs(s,s∞,p-1)*hc(c,c∞,p+1)/hc(c,c∞,p)
        end

        si*col0
    end
end

function mul_coefficients(Q::HessenbergUnitary{'U'},v::Vector;opts...)
    si=Q.sign ? 1 : -1

    ret = pad!(si*v,length(v)+1)
    # Compute each Givens rotation starting from the right

    for i = length(v):-1:1
        ret[i],ret[i+1] = hc(Q,i+1)*ret[i] + hs(Q,i)*ret[i+1],
                -hs(Q,i)*ret[i] + hc(Q,i+1)*ret[i+1]
    end
    ret[1]*=hc(Q,1)
    ret
end

function mul_coefficients(Q::HessenbergUnitary{'L'},v::Vector;tolerance=eps())
    N =  max(length(v),length(Q.s))+1
    si = Q.sign ? 1 : -1
    ret = pad!(si*v,N)

    # This part does the computation we are certain we have to do
    ret[1] *= hc(Q,1)
    for i = 1:N-1
        ret[i],ret[i+1] = hc(Q,i+1)*ret[i] -hs(Q,i)*ret[i+1],
                           hs(Q,i)*ret[i] + hc(Q,i+1)*ret[i+1]
    end

    # After this point, ret is monotonically decreasing to zero
    i = N
    while abs(ret[i]) > tolerance
        push!(ret,(Q.s∞)*ret[i])
        ret[i] *= Q.c∞
        i += 1
    end
    ret
end

-(Q::HessenbergUnitary{uplo}) where {uplo}=
    HessenbergUnitary{uplo,eltype(Q)}(!Q.sign,Q.c,Q.s,
                                          Q.c∞,
                                          Q.s∞,Q.band)

deflate(Q::HessenbergUnitary{uplo}) where {uplo}=HessenbergUnitary(Val{uplo},Q.sign,
                                                                  [(Q.sign ? 1 : -1)*sign(Q.c[1]);Q.c],
                                                                  [0;Q.s],Q.c∞,Q.s∞,Q.band)

deflate(Q::HessenbergUnitary,k::Integer) = k == 0 ? Q : deflate(deflate(Q),k-1)

struct BandedUnitary{uplo,T} <: UnitaryOperator{T}
    ops::Vector{HessenbergUnitary{uplo,T}}
end

adjoint(Q::BandedUnitary) = BandedUnitary(reverse!(map(adjoint,Q.ops)))

getindex(Q::BandedUnitary,k::Integer,j::Integer) = TimesOperator(Q.ops)[k,j]
bandwidths(Q::BandedUnitary) = bandwidths(TimesOperator(Q.ops))

domainspace(::BandedUnitary) = ℓ⁰
rangespace(::BandedUnitary) = ℓ⁰

function mul_coefficients(Q::BandedUnitary,v::Vector)
    ret=v
    for k=length(Q.ops):-1:1
        ret=mul_coefficients(Q.ops[k],ret)
    end
    ret
end


text:
using Azure
using Azure.StorageServices
using Test

function test_blobs(tenant_id, appid, password, subscription_id, resource_group_name, storage)
    creds = AzureCredentials(tenant_id, appid, password)
    ctx = AzureContext(creds)

    resource_group_name = "testresgroup"
    testcontainer = storage * "testcontainer"
    testblob1 = testcontainer * "/testblob1"
    testblob2 = testcontainer * "/testblob2"

    @test createContainer(ctx, subscription_id, resource_group_name, testcontainer)

    result = listContainers(ctx, subscription_id, resource_group_name, storage)
    @test haskey(result, "Containers")
    @test haskey(result["Containers"], "Container")
    containers = result["Containers"]["Container"]

    container = containers
    if isa(containers, Vector)
        # iterate through multiple containers
        for c in containers
            if c["Name"] == "testcontainer"
                container = c
            end
        end
    end

    @test haskey(container, "Name")
    @test container["Name"] == "testcontainer"

    # add a string as a blob
    @test putBlob(ctx, subscription_id, resource_group_name, testblob1, "BlockBlob"; block_blob_contents="hello world")

    # streaming a file as a blob
    mktemp() do path, wio
        println(wio, "hello world")
        close(wio)
        open(path) do rio
           @test putBlob(ctx, subscription_id, resource_group_name, testblob2, "BlockBlob"; block_blob_contents=rio, content_length=filesize(rio))
        end
    end

    # read blobs
    result = listBlobs(ctx, subscription_id, resource_group_name, testcontainer)
    @test haskey(result, "Blobs")
    @test haskey(result["Blobs"], "Blob")
    blobs = result["Blobs"]["Blob"]
    @test length(blobs) == 2
    blobnames = Set(b["Name"] for b in blobs)
    @test "testblob1" in blobnames
    @test "testblob2" in blobnames

    blobbytes1 = getBlob(ctx, subscription_id, resource_group_name, testblob1)
    @test String(blobbytes1) == "hello world"

    blobbytes2 = getBlob(ctx, subscription_id, resource_group_name, testblob2)
    blobstr2 = String(blobbytes2)
    @test strip(blobstr2) == "hello world"
    @test length(blobstr2) > length("hello world")

    # delete blobs
    @test deleteBlob(ctx, subscription_id, resource_group_name, testblob1)
    @test deleteBlob(ctx, subscription_id, resource_group_name, testblob2)

    @test deleteContainer(ctx, subscription_id, resource_group_name, testcontainer)
end


text:
using TEMPsolver_PDE
using Test

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


text:
using PackageCompiler
create_app(".", "./vchdfe", app_name="vchdfe", force=true, precompile_execution_file="deps/precompilation_execution_file.jl")