text:
using CellListMap
using StaticArrays
using LinearAlgebra

export test
export compute_pairwise_velocity_moments
export compute_pairwise_velocity_los_pdf

global n_moments  = 9

function compute_pairwise_mean!(x,y,i,j,d2,hist,velocities, rbins,sides)
    d = x - y
    r = sqrt.(d2)
    ibin = searchsortedfirst(rbins, r) - 1
    dv = velocities[i] - velocities[j]
    v_r = LinearAlgebra.dot(dv,d)/r
    cos_theta = d[3]/r
    sin_theta = sqrt(d[1]*d[1]+d[2]*d[2])/r
    cos_phi = d[1]/sqrt(d[1]*d[1]+d[2]*d[2])
    sin_phi = d[2]/sqrt(d[1]*d[1]+d[2]*d[2])
    if(sqrt(d[1]*d[1]+d[2]*d[2]) < 1e-10)
        cos_phi = 1.0
        sin_phi = 0.0
    end
    v_t = dv[1] * cos_theta * cos_phi + dv[2] * cos_theta * sin_phi - dv[3] * sin_theta
    if ibin > 0
        hist[1][ibin] += 1
        hist[2][ibin] += v_r # v_r
        hist[3][ibin] += v_r * v_r # std_r
        hist[4][ibin] += v_r * v_r * v_r # skew_r
        hist[5][ibin] += v_r * v_r * v_r * v_r # kur_r
        hist[6][ibin] += v_t * v_t # std_t
        hist[7][ibin] += v_r * v_t * v_t # skew_rt
        hist[8][ibin] += v_t * v_t * v_t * v_t # kur_t
        hist[9][ibin] += v_r * v_r * v_t * v_t # kur_rt
    end
    return hist
end

function compute_pairwise_mean!(x,y,i,j,d2,hist,velocities_left, velocities_right, rbins,sides)
    d = x - y
    r = sqrt.(d2)
    ibin = searchsortedfirst(rbins, r) - 1
    dv = velocities_left[i] - velocities_right[j]
    v_r = LinearAlgebra.dot(dv,d)/r
    cos_theta = d[3]/r
    sin_theta = sqrt(d[1]*d[1]+d[2]*d[2])/r
    cos_phi = d[1]/sqrt(d[1]*d[1]+d[2]*d[2])
    sin_phi = d[2]/sqrt(d[1]*d[1]+d[2]*d[2])
    if(sqrt(d[1]*d[1]+d[2]*d[2]) < 1e-10)
        cos_phi = 1.0
        sin_phi = 0.0
    end
    v_t = dv[1] * cos_theta * cos_phi + dv[2] * cos_theta * sin_phi - dv[3] * sin_theta
    if ibin > 0
        hist[1][ibin] += 1
        hist[2][ibin] += v_r
        hist[3][ibin] += v_r * v_r
        hist[4][ibin] += v_r * v_r * v_r
        hist[5][ibin] += v_r * v_r * v_r * v_r
        hist[6][ibin] += v_t * v_t
        hist[7][ibin] += v_r * v_t * v_t
        hist[8][ibin] += v_t * v_t * v_t * v_t
        hist[9][ibin] += v_r * v_r * v_t * v_t
    end
    return hist
end

function convert_histogram_into_moments!(hist)
    n_pairs = hist[1]
    mask = (n_pairs .> 0)
    hist[2][mask] = hist[2][mask] ./ n_pairs[mask]
    hist[3][mask] = hist[3][mask] ./ n_pairs[mask] - hist[2][mask].^2
    hist[6][mask] = hist[6][mask] ./ n_pairs[mask]
    hist[4][mask] = (
                hist[4][mask] ./ n_pairs[mask]
                - 3. * hist[2][mask] .* hist[3][mask]
                - hist[2][mask].^3.
    )
    hist[7][mask] = (
            hist[7][mask] ./ n_pairs[mask]
            - hist[2][mask] .* hist[6][mask]
    )
    hist[5][mask] = (
            hist[5][mask] ./ n_pairs[mask]
            - 4. * hist[2][mask] .* hist[4][mask]
            - 6. * hist[2][mask].^2. .* hist[3][mask]
            - hist[2][mask] .^4.
    )
    hist[8][mask] = hist[8][mask] ./ n_pairs[mask]
    hist[9][mask] = (
        hist[9][mask] ./ n_pairs[mask]
        - 2. * hist[7][mask] .* hist[2][mask]
        - hist[2][mask].^2 .* hist[6][mask]
    )
    hist[3][:] = sqrt.(hist[3])
    hist[6][:] = sqrt.(hist[6])
    hist[4][:] = hist[4] ./ hist[3].^3
    hist[5][:] = hist[5] ./ hist[3].^4
    hist[8][:] = hist[8] ./ hist[6].^4
    hist[7][:] = hist[7] ./ hist[3].^2 ./ hist[6]
    hist[9][:] = hist[9] ./ hist[3].^2 ./ hist[6] .^2
end

function reduce_hist(hist,hist_threaded)
  hist = hist_threaded[1]
  for i in 2:length(hist_threaded)
    for moment in 1:n_moments
        hist[moment] .+= hist_threaded[i][moment]
    end
  end
  return hist
end

function reduce_hist_los(hist,hist_threaded)
  hist = hist_threaded[1]
  for i in 2:length(hist_threaded)
        hist .+= hist_threaded[i]
  end
  return hist
end

function get_pairwise_velocity_moments(
  positions, velocities,
  rbins,
  boxsize,
  cl, box
)
    hist = (
            zeros(Int,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
    )
    hist = map_pairwise!(
      (x,y,i,j,d2,hist) ->
      compute_pairwise_mean!(x,y,i,j,d2,hist,velocities, rbins, boxsize),
      hist, box, cl,
      reduce=reduce_hist,
      parallel=true,
      show_progress=false,
    )
    convert_histogram_into_moments!(hist)
    return hist
end

function get_pairwise_velocity_moments(
  positions_left, velocities_left,
  positions_right, velocities_right,
  rbins,
  boxsize,
  cl, box
)
    hist = (
            zeros(Int,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
            zeros(Float64,length(rbins)-1),
    )
    hist = map_pairwise!(
      (x,y,i,j,d2,hist) ->
      compute_pairwise_mean!(x,y,i,j,d2,hist,velocities_left, velocities_right,rbins, boxsize),
      hist, box, cl,
      reduce=reduce_hist
    )
    convert_histogram_into_moments!(hist)
    return hist
end

function compute_pairwise_velocity_moments(
    positions, velocities, rbins, boxsize
)
    println("Using n = ", Threads.nthreads() , " threads")
    Lbox = [boxsize,boxsize,boxsize]
    n = size(positions)[2]
    positions = reshape(reinterpret(SVector{3,Float64},positions),n)
    velocities = reshape(reinterpret(SVector{3,Float64},velocities),n)
    r_max = maximum(rbins)
    box = Box(Lbox, r_max, lcell=1)
    cl = CellList(positions,box)
    return get_pairwise_velocity_moments(
      positions,
      velocities,
      rbins,
      Lbox,
      cl,
      box,
    )

end

function compute_pairwise_velocity_moments(
    positions_left, velocities_left, positions_right, velocities_right, rbins, boxsize
)
    println("Using n = ", Threads.nthreads() , " threads")
    Lbox = [boxsize,boxsize,boxsize]
    positions_left = reshape(reinterpret(SVector{3,Float64},positions_left),size(positions_left)[2])
    velocities_left = reshape(reinterpret(SVector{3,Float64},velocities_left),size(velocities_left)[2])
    positions_right = reshape(reinterpret(SVector{3,Float64},positions_right),size(positions_right)[2])
    velocities_right = reshape(reinterpret(SVector{3,Float64},velocities_right),size(velocities_right)[2])
    r_max = maximum(rbins)
    box = Box(Lbox, r_max, lcell=1)
    cl = CellList(positions_left, positions_right, box)
    return get_pairwise_velocity_moments(
      positions_left,
      velocities_left,
      positions_right,
      velocities_right,
      rbins,
      Lbox,
      cl, box
    )

end

function compute_pairwise_los_pdf!(x,y,i,j,d2,hist,velocities,r_perp_bins, r_parallel_bins, vlos_bins,sides)
    d = x - y
    r_perp = sqrt(d[1]^2 + d[2]^2)
    r_parallel = d[3]
    vlos = (velocities[i][3] - velocities[j][3])*sign(r_parallel)
    if (r_perp > minimum(r_perp_bins)) & (r_perp < maximum(r_perp_bins)) & (r_parallel < maximum(r_parallel_bins)) & (r_parallel > minimum(r_parallel_bins)) & (vlos < maximum(vlos_bins)) & (vlos > minimum(vlos_bins))
        ibin_perp = searchsortedfirst(r_perp_bins, r_perp) - 1
        ibin_parallel = searchsortedfirst(r_parallel_bins, r_parallel) - 1
        ibin_vlos = searchsortedfirst(vlos_bins, vlos) - 1
        hist[ibin_perp, ibin_parallel, ibin_vlos] += 1
    end
    return hist
end

function compute_pairwise_los_pdf!(x,y,i,j,d2,hist,
        velocities_left,velocities_right,
        r_perp_bins, r_parallel_bins, vlos_bins,sides)
    d = x - y
    r_perp = sqrt(d[1]^2 + d[2]^2)
    r_parallel = d[3]
    vlos = (velocities_left[i][3] - velocities_right[j][3])*sign(r_parallel)
    if (r_perp < maximum(r_perp_bins)) & (r_parallel < maximum(r_parallel_bins)) & (r_parallel > minimum(r_parallel_bins)) & (vlos < maximum(vlos_bins)) & (vlos > minimum(vlos_bins))
        ibin_perp = searchsortedfirst(r_perp_bins, r_perp) - 1
        ibin_parallel = searchsortedfirst(r_parallel_bins, r_parallel) - 1
        ibin_vlos = searchsortedfirst(vlos_bins, vlos) - 1
        hist[ibin_perp, ibin_parallel, ibin_vlos] += 1
    end
    return hist
end

function get_pairwise_velocity_los_pdf(
  positions, velocities,
  r_perp_bins,
  r_parallel_bins,
  vlos_bins,
  boxsize,
  cl, box
)
    hist = zeros(Int,(
            length(r_perp_bins)-1,
            length(r_parallel_bins)-1,
            length(vlos_bins)-1
    )
    )
    return map_pairwise!(
      (x,y,i,j,d2,hist) ->
      compute_pairwise_los_pdf!(x,y,i,j,d2,hist,velocities,
        r_perp_bins, r_parallel_bins, vlos_bins, boxsize),
      hist, box, cl,
      reduce=reduce_hist_los,
      parallel=true,
      show_progress=false,
    )
end

function get_pairwise_velocity_los_pdf(
  positions_left, velocities_left,
  positions_right, velocities_right,
  r_perp_bins,
  r_parallel_bins,
  vlos_bins,
  boxsize,
  cl, box
)
    hist = zeros(Int,(
            length(r_perp_bins)-1,
            length(r_parallel_bins)-1,
            length(vlos_bins)-1
    )
    )
    return map_pairwise!(
      (x,y,i,j,d2,hist) ->
      compute_pairwise_los_pdf!(x,y,i,j,d2,hist,velocities_left, velocities_right,
        r_perp_bins, r_parallel_bins, vlos_bins, boxsize),
      hist, box, cl,
      reduce=reduce_hist_los,
      parallel=true,
      show_progress=false,
    )
end

function compute_pairwise_velocity_los_pdf(
    positions, velocities, r_perp_bins, r_parallel_bins, vlos_bins, boxsize
)
    println("Using n = ", Threads.nthreads() , " threads")
    Lbox = [boxsize,boxsize,boxsize]
    r_max = sqrt(maximum(r_perp_bins)^2 + maximum(r_parallel_bins)^2)
    positions = reshape(reinterpret(SVector{3,Float64},positions),size(positions)[2])
    velocities = reshape(reinterpret(SVector{3,Float64},velocities),size(velocities)[2])
    box = Box(Lbox, r_max, lcell=1)
    cl = CellList(positions,box)
    return get_pairwise_velocity_los_pdf(
      positions,
      velocities,
      r_perp_bins,
      r_parallel_bins,
      vlos_bins,
      Lbox,
      cl, box,
    )

end

function compute_pairwise_velocity_los_pdf(
    positions_left, velocities_left,
    positions_right, velocities_right,
    r_perp_bins, r_parallel_bins, vlos_bins, boxsize
)
    println("Using n = ", Threads.nthreads() , " threads")
    Lbox = [boxsize,boxsize,boxsize]
    r_max = sqrt(maximum(r_perp_bins)^2 + maximum(r_parallel_bins)^2)
    positions_left = reshape(reinterpret(SVector{3,Float64},positions_left),size(positions_left)[2])
    velocities_left = reshape(reinterpret(SVector{3,Float64},velocities_left),size(velocities_left)[2])
    positions_right = reshape(reinterpret(SVector{3,Float64},positions_right),size(positions_right)[2])
    velocities_right = reshape(reinterpret(SVector{3,Float64},velocities_right),size(velocities_right)[2])

    box = Box(Lbox, r_max, lcell=1)
    cl = CellList(positions_left,positions_right,box)
    return get_pairwise_velocity_los_pdf(
      positions_left,
      velocities_left,
      positions_right,
      velocities_right,
      r_perp_bins,
      r_parallel_bins,
      vlos_bins,
      Lbox,
      cl, box,
    )

end


text:
using Distributed
import Dagger: Context, Processor, OSProc, ThreadProc, get_parent, get_processors
import Dagger.Sch: ThunkOptions

@everywhere begin

struct UnknownStruct end

struct OptOutProc <: Dagger.Processor end

struct PathProc <: Dagger.Processor
    owner::Int
end
Dagger.get_parent(proc::PathProc) = OSProc(proc.owner)
Dagger.move(ctx, ::PathProc, ::OSProc, x::Float64) = x+1
Dagger.move(ctx, ::OSProc, ::PathProc, x::Float64) = x+2
Dagger.iscompatible(proc::PathProc, opts, f, args...) = true
Dagger.execute!(proc::PathProc, func, args...) = func(args...)

end

@testset "Processors" begin
    @testset "Parents/Children" begin
        tp = ThreadProc(1, 1)
        @test tp isa Processor
        op = get_parent(tp)
        @test op isa Processor
        @test op isa OSProc
        @test op.pid == 1
    end
    @testset "Function/argument compatability" begin
        unknown_func = () -> nothing
        tp = ThreadProc(1, 1)
        op = get_parent(tp)
        opts = ThunkOptions()
        us = UnknownStruct()
        for proc in (op, tp)
            @test Dagger.iscompatible_func(proc, opts, unknown_func)
            @test Dagger.iscompatible_arg(proc, opts, us)
            @test Dagger.iscompatible(proc, opts, unknown_func, us, 1, us, 2.0)
        end
    end
    @testset "Opt-in/Opt-out" begin
        @test Dagger.default_enabled(OSProc()) == true
        @test Dagger.default_enabled(ThreadProc(1,1)) == true
        @test Dagger.default_enabled(OptOutProc()) == false
    end
    @testset "Processor exhaustion" begin
        opts = ThunkOptions(proctypes=[OptOutProc])
        @test_throws CapturedException collect(delayed(sum; options=opts)([1,2,3]))
    end
    @testset "Roundtrip move()" begin
        ctx = Context()
        tp = ThreadProc(1, 1)
        op = get_parent(tp)
        value = rand()
        moved_value = Dagger.move(ctx, tp, op, Dagger.move(ctx, op, tp, value))
        @test value === moved_value
    end
    @testset "Generic path move()" begin
        @everywhere Dagger.add_callback!(proc->PathProc(myid()))
        ctx = Context()
        proc1 = first(filter(x->x isa PathProc, get_processors(OSProc(1))))
        proc2 = first(filter(x->x isa PathProc, get_processors(OSProc(2))))
        value = rand()
        moved_value = Dagger.move(ctx, proc1, proc2, value)
        @test moved_value == value+3
        @everywhere pop!(Dagger.PROCESSOR_CALLBACKS)
    end
    @testset "Add callback in same world" begin
        function addcb()
            @everywhere Dagger.add_callback!(cb)
            cb = eval(Dagger, :(proc->PathProc(myid())))
            opts = ThunkOptions(proctypes=[PathProc])
            collect(delayed(identity; options=opts)(1.0))
            @everywhere pop!(Dagger.PROCESSOR_CALLBACKS)
        end
    end
end

ext: jl


text:
__precompile__(false)

module Optics

import RefractiveIndex: Material

include("reflection.jl")
include("transmittance.jl")

export phaseangle, reflection, snellius, transmittance

end # module


text:
include("environment/blankenv.jl")
include("environment/iss.jl")
include("environment/iss_corner.jl")
include("environment/table.jl")
include("environment/forest.jl")


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/BoundingSphere.jl/blob/master/LICENSE

@testset "6 points 2d" begin

    pts = [
        [-2.740694122152086187327313382411375641822814941406250000000000000000000000000000, 8.010609127471919777008224627934396266937255859375000000000000000000000000000000],
        [-7.875078087682085836718215432483702898025512695312500000000000000000000000000000e-02, 8.253169461411594909350242232903838157653808593750000000000000000000000000000000],
        [-5.281301381389667426446976605802774429321289062500000000000000000000000000000000, 6.407720312422610753344542899867519736289978027343750000000000000000000000000000],
        [3.279069514617182434790265688207000494003295898437500000000000000000000000000000, -2.195389590186907824431727931369096040725708007812500000000000000000000000000000],
        [3.279069514617182434790265688207000494003295898437500000000000000000000000000000, -2.195389590186907824431727931369096040725708007812500000000000000000000000000000]]

    c,r = boundingsphere(pts, WelzlPivot())
    ball = MB.SqBall(c, r^2)
    @test MB.allinside(pts, ball, atol=10)
    @test_broken MB.allinside(pts, ball, rtol=1e-1)
    @test_broken MB.allinside(pts, ball, rtol=1e-6)

end


text:
## just takes realpart of operator
immutable ReOperator{O,T} <: Operator{T}
    op::O
end

ReOperator(op)=ReOperator{typeof(op),Float64}(op)
Base.convert{BT<:Operator}(::Type{BT},R::ReOperator)=ReOperator{typeof(R.op),eltype(BT)}(R.op)

for OP in (:rangespace,:domainspace,:bandinds)
    @eval $OP(R::ReOperator)=$OP(R.op)
end

getindex(RI::ReOperator,k::Integer,j::Integer) =
    real(RI.op[k,j])

choosedomainspace(R::ReOperator,sp::Space) = choosedomainspace(R.op,sp)
for OP in (:promotedomainspace,:promoterangespace)
    @eval begin
        $OP(R::ReOperator,sp::UnsetSpace) = ReOperator($OP(R.op,sp))
        $OP(R::ReOperator,sp::Space) = ReOperator($OP(R.op,sp))
    end
end

# TODO: can't do this because UnsetSpace might change type
#Base.real{T<:Real}(op::Operator{T})=op
Base.real(op::Operator)=ReOperator(op)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
Support for sparse arrays. Provides `AbstractSparseArray` and subtypes.
"""
module SparseArrays

using Base: ReshapedArray, promote_op, setindex_shape_check, to_shape, tail,
    has_offset_axes
using Base.Sort: Forward
using LinearAlgebra

import Base: +, -, *, \, /, &, |, xor, ==
import LinearAlgebra: mul!, ldiv!, rdiv!, cholesky, adjoint!, diag, eigen, dot,
    issymmetric, istril, istriu, lu, tr, transpose!, tril!, triu!,
    cond, diagm, factorize, ishermitian, norm, opnorm, lmul!, rmul!, tril, triu

import Base: @get!, acos, acosd, acot, acotd, acsch, asech, asin, asind, asinh,
    atan, atand, atanh, broadcast!, conj!, cos, cosc, cosd, cosh, cospi, cot,
    cotd, coth, count, csc, cscd, csch,
    exp10, exp2, findprev, findnext, floor, hash, argmin, inv,
    log10, log2, sec, secd, sech, show,
    sin, sinc, sind, sinh, sinpi, dropdims, sum, summary, tan,
    tand, tanh, trunc, abs, abs2,
    broadcast, ceil, complex, conj, convert, copy, copyto!, adjoint,
    exp, expm1, findall, findmax, findmin, float, getindex,
    vcat, hcat, hvcat, cat, imag, argmax, kron, length, log, log1p, max, min,
    maximum, minimum, one, promote_eltype, real, reshape, rot180,
    rotl90, rotr90, round, setindex!, similar, size, transpose,
    vec, permute!, map, map!, Array, diff, circshift!, circshift

using Random: GLOBAL_RNG, AbstractRNG, randsubseq, randsubseq!

export AbstractSparseArray, AbstractSparseMatrix, AbstractSparseVector,
    SparseMatrixCSC, SparseVector, blockdiag, droptol!, dropzeros!, dropzeros,
    issparse, nonzeros, nzrange, rowvals, sparse, sparsevec, spdiagm,
    sprand, sprandn, spzeros, nnz, permute, findnz

include("abstractsparse.jl")
include("sparsematrix.jl")
include("sparsevector.jl")
include("higherorderfns.jl")
include("linalg.jl")
include("deprecated.jl")

# temporarily moved here and commented out from from base/linalg/diagonal.jl, base/linalg/tridiag.jl
# and base/linalg/bidiag.jl due to their usage of spzeros
similar(B::Bidiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...)
similar(D::Diagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...)
similar(S::SymTridiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...)
similar(M::Tridiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...)

matprod(x, y) = x*y + x*y
const BiTriSym = Union{Bidiagonal,SymTridiagonal,Tridiagonal}
function *(A::BiTriSym, B::BiTriSym)
    TS = promote_op(matprod, eltype(A), eltype(B))
    mul!(similar(A, TS, size(A)...), A, B)
end

end

max_stars_repo_name: ashleysommer/julia

max_issues_repo_name: ashleysommer/julia

max_forks_repo_name: ashleysommer/julia


text:
# DESCRIPTION
# Horizontal stabilizer scheduling with airspeed and pitch rate.
#
# ------------------------------------------------------------------------------

function stabSched(Coll,Veq,ayg,q,StabSet)
    if (StabSet==0)
        aylim=max(min(ayg,0.13),-0.13)
        stab_ay=-0.0433*min((Veq-30.0),30.0)*aylim
        stab_q=0.16*180/pi*q
        VSchedule=-0.1119*min(max(Veq-80.0,0.0),67.0)
        PR50=min(max(Veq-30.0,0.0),50.0)
        C1=min(max(70.0-Coll,0.),20.0)
        C2=min(max(50.0-Coll,0.),50.0)
        CollSchedule=-1.0*PR50*(0.003286*(C1+C2)+0.67)
        IHT0=42.0+CollSchedule+VSchedule+stab_q+stab_ay
        # stabilator incidence [deg]
        StabInc=min(max(IHT0,-8.0),39.0)
    else
        StabInc=StabSet
    end
end

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text:
"""
    ww = rfweights(orofile, reffile, nf; weightsfn="", varname="", fsmooth=false)

Compute orographic weights from a fine-scale precipitation climatology file.

#Arguments
* `orofile`  : filename of input climatology
* `reffile`  : filename of reference file (for metadata, e.g. the file to downscale)
* `nf`       : refinement factor for spatial downscaling
* `weightsfn`: write weights to file weightsfn
* `varname`  : variable name in climatology
* `fsmooth`  : use smoothing instead of gp conservation

#Returns
* `ww`       : a weight matrix also saved to weightsfn

#Depends

This function uses external system calls using the "cdo" command (https://code.mpimet.mpg.de/projects/cdo/wiki/Cdo) which needs to be available on your system.
"""
function rfweights(orofile, reffile, nf; weightsfn="", varname="", fsmooth=false)

  # Create a reference gridrf.nc file (same grid as rainfarm output files)
  (pr,lon_mat,lat_mat)=read_netcdf2d(reffile, varname);
  # Create fine scale grid
  nss=size(pr)
  if (length(nss)>=3)
    pr=pr[:,:,1]
  end
  println(nss)
  ns=nss[1];
  (lon_f, lat_f)=lon_lat_fine(lon_mat, lat_mat,nf);

  rr=round.(Int,rand(1)*100000)[1]

  println("Output size: ",size(lon_f))
  if(varname=="")
    varname="pr"
    run(`cdo -s  setname,pr $reffile reffile_rr.nc`)
    reffile="reffile_rr.nc"
  end
  if(weightsfn=="")
    weightsfn="weights$rr.nc"
    fsavew=false
  else
    fsavew=true
  end

  # The rest is done in CDO
  println("Computing weights")
  write_netcdf2d("gridrf.nc",reshape(pr,ns,ns,1),lon_f,lat_f,varname,reffile)
  run(`cdo -s -b F32 timmean $orofile pr_orofile_$rr.nc`)
  run(`cdo -s -f nc copy gridrf.nc gridrf_2_$rr.nc`)
  run(`cdo -s -f nc remapbil,gridrf_2_$rr.nc pr_orofile_$rr.nc pr_remap_rr.nc`)
  if(fsmooth)
    (prr,lon,lat)=read_netcdf2d("pr_remap_rr.nc","")
    ww=prr./smoothconv(prr,ns);
    write_netcdf2d(weightsfn,ww,lon_f,lat_f,varname,reffile)
    run(`rm -f pr_remap_rr.nc pr_orofile_$rr.nc gridrf.nc gridrf_2_$rr.nc reffile_rr.nc`)
  else
    run(`cdo -s gridboxmean,$nf,$nf pr_remap_rr.nc pr_remap_gbm_$rr.nc`)
    run(`cdo -s remapnn,pr_remap_rr.nc pr_remap_gbm_$rr.nc pr_remap_nn_$rr.nc`)
    run(`cdo -s div pr_remap_rr.nc pr_remap_nn_$rr.nc $weightsfn`)
    run(`rm -f pr_remap_rr.nc pr_remap_gbm_$rr.nc pr_remap_nn_$rr.nc pr_orofile_$rr.nc gridrf.nc gridrf_2_$rr.nc reffile_rr.nc`)
  end
  (ww,lon_f,lat_f)=read_netcdf2d(weightsfn, "");
  if(!fsavew)
     run(`rm -f $weightsfn`)
  end
  return ww
end


text:
# draw histogram of sensor data 600 as time series

module DrawTimeSeries600
    using Plots, DataFrames, CSV
    pyplot()

    function main(is_test=false)
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_600.txt")
        df_600_mm = CSV.read(data_path, DataFrame,
                             header=["date", "time", "ir", "lidar"],
                             delim=' ')

        time_idx = Array(0:length(df_600_mm.lidar)-1)

        plot(time_idx, df_600_mm.lidar, lebel="time series")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/complex_dist/time_series_600mm/time_series_600_mm.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
immutable Color
    r::UInt8
    g::UInt8
    b::UInt8
    a::UInt8

    function Color(r::Integer, g::Integer, b::Integer)
        r=clamp(r, 0, 255); g=clamp(g, 0, 255); b=clamp(b, 0, 255);
        new(UInt8(r), UInt8(g), UInt8(b), UInt8(255))
    end

    function Color(r::Integer, g::Integer, b::Integer, a::Integer)
        r=clamp(r, 0, 255); g=clamp(g, 0, 255); b=clamp(b, 0, 255); a=clamp(a, 0, 255)
        new(UInt8(r), UInt8(g), UInt8(b), UInt8(a))
    end
end

const black = Color(0, 0, 0)
const white = Color(255, 255, 255)
const red = Color(255, 0, 0)
const green = Color(0, 255, 0)
const blue = Color(0, 0, 255)
const yellow = Color(255, 255, 0)
const magenta = Color(255, 0, 255)
const cyan = Color(0, 255, 255)
const transparent = Color(0, 0, 0, 0)


text:
using Test
using LinearAlgebra
using NiLang, NiLang.AD

@testset "i_norm2, dot" begin
    out = 0.0im
    vec = [1.0im, 2.0, 3.0]
    vec2 = [1.0, 2.0im, 5.0]
    @instr i_norm2(out, vec)
    @test out ≈ norm(vec)^2
    @test check_inv(i_norm2, (out, vec))

    out = 0.0im
    vec = [1.0im, 2.0, 3.0]
    vec2 = [1.0, 2.0im, 5.0]
    @instr i_dot(out, vec, vec2)
    @test out ≈ dot(vec, vec2)
    @test check_inv(i_dot, (out, vec, vec2))

    out = 0.0
    vec = [1.0, 2.0, 3.0]
    vec2 = [1.0, 2.0, 5.0]
    @test check_grad(i_norm2, (out, vec); verbose=true, iloss=1)

    out = 0.0
    @instr i_dot(out, vec, vec2)
    @test out ≈ dot(vec, vec2)
    @test check_inv(i_dot, (out, vec, vec2))

    @test check_grad(i_dot, (0.0, vec, vec2); verbose=true, iloss=1)

    m = randn(4,4)
    n = randn(4,4)
    out = 0.0
    @instr i_dot(out, m[:,2], n[:,4])
    @test out ≈ dot(m[:,2], n[:,4])
    @test check_inv(i_dot, (out, m[:,2], n[:,4]))

    @test check_grad(i_dot, (0.0, vec, vec2); verbose=true, iloss=1)
end

function naive_umm!(x, params)
    N = size(x, 1)
    k = 0
    for j=1:N
        for i=N-1:-1:j
            k += 1
            a, b = rot(x[i],x[i+1],params[k])
            x[i], x[i+1] = a, b
        end
    end
end

function inv_naive_umm!(x, params)
    N = size(x, 1)
    k = N*(N-1) ÷ 2
    for j=N:-1:1
        for i=j:N-1
            a, b = rot(x[i],x[i+1],-params[k])
            x[i], x[i+1] = a, b
            k -= 1
        end
    end
end
@testset "naive unitary" begin
    x = randn(200)
    params = randn(100*199).*2π

    x0 = copy(x)
    params0 = copy(params)
    naive_umm!(x, params)
    inv_naive_umm!(x, params)
    @test params ≈ params0
    @test x ≈ x0
end

@testset "unitary" begin
    x = randn(20)
    params = randn(10*19) * 2π

    x0 = copy(x)
    params0 = copy(params)
    Nx = length(x)
    @instr i_umm!(x, params)
    x1 = copy(x0)
    params1 = copy(params0)
    naive_umm!(x1, params1)
    @test params ≈ params1
    @test x ≈ x1
    @instr (~i_umm!)(x, params)
    @test params ≈ params0
    @test x ≈ x0
    @test check_inv(i_umm!, (x, params))
end


text:
using HTTP
using JSON3

HTTP.serve("127.0.0.1", 9174) do http
  obj = JSON3.read(HTTP.payload(http))
  println(obj)
  println("")
  println("token: $(obj.token)")

  resp = JSON3.write((content = obj.message.content, ))
  return HTTP.Response(resp)
end

max_stars_repo_name: Arkoniak/ZulipSnippetBot

max_issues_repo_name: Arkoniak/ZulipSnippetBot

max_forks_repo_name: Arkoniak/ZulipSnippetBot


text:
module HopfieldNets
    export HopfieldNet, DiscreteHopfieldNet, ContinuousHopfieldNet
    export update!, energy, settle!, train!, associate!, storkeytrain!

    include("generic.jl")
    include("discrete.jl")
    include("continuous.jl")
end


text:
using Sundials, Test, LinearAlgebra

# create mass matrix problems
function make_mm_probs(mm_A, ::Type{Val{iip}}) where iip
  # iip
  mm_b = vec(sum(mm_A; dims=2))
  function mm_f(du,u,p,t)
    mul!(du,mm_A,u)
    du .+= t * mm_b
    nothing
  end
  mm_g(du,u,p,t) = (@. du = u + t; nothing)

  # oop
  mm_f(u,p,t) = mm_A * (u .+ t)
  mm_g(u,p,t) = u .+ t

  mm_analytic(u0, p, t) = @. 2 * u0 * exp(t) - t - 1

  u0 = ones(3)
  tspan = (0.0, 1.0)

  prob = ODEProblem(ODEFunction{iip,true}(mm_f, analytic=mm_analytic, mass_matrix=mm_A), u0, tspan)
  prob2 = ODEProblem(ODEFunction{iip,true}(mm_g, analytic=mm_analytic), u0, tspan)

  prob, prob2
end

mm_A = Float64[-2 1 4; 4 -2 1; 2 1 3]
prob, prob2 = make_mm_probs(mm_A, Val{true})

sol = solve(prob,  ARKODE(),abstol=1e-8,reltol=1e-8)
sol2 = solve(prob2,ARKODE(),abstol=1e-8,reltol=1e-8)

@test norm(sol .- sol2) ≈ 0 atol=1e-7


text:
import GuillotineModels
const IR = GuillotineModels.InstanceReader

function kp_ub(e, d_, p_, a_, A)
	se = sortperm(e; rev = true)
	d, p, a = getindex.((d_, p_, a_), (se, se, se))
	rem_A = A
	ub = 0.0
	for (i, qt) in pairs(IndexLinear(), d)
		for _ in 1:qt
			if rem_A < a[i]
				ub += floor((rem_A  * p[i]) / a[i])
				@goto out
			end
			rem_A -= a[i]
			ub += p[i]
		end
	end
	@label out
	return ub
end

function extract_inst_data(inst_name)
	N, L, W, l, w, p, d = IR.read_from_file(inst_name)
	n_ = sum(d)
	a = l .* w
	A = L .* W
	csv_lb = sum((d .* a) ./ A)
	e = p ./ (l .* w)
	kp_ub_ = kp_ub(e, d, p, a, A)
	min_l, min_w, min_a = minimum.((l, w, a))
	return inst_name, N, n_, L, W, A, min_l, min_w, min_a, kp_ub_, csv_lb
end

function create_csv_for(inst_names)
	header = "instance;N;n_;L;W;A;min_l;min_w;min_a;kp_ub;csv_lb\n"
	lines = map(inst_names) do inst
		values = extract_inst_data(inst)
		join(map(string, values), ";")
	end
	body = join(lines, "\n")
	csv = header * body

	return csv
end

const DATASET_C = vcat(
	"A" .* string.(1:5),
	"CHL" .* string.([1, 2, 5, 6, 7]),
	"CU" .* ["1", "2"],
	"CW" .* string.(1:3),
	"Hchl" .* ["2", "9"],
	"Hchl" .* string.([3, 4, 6, 7, 8]) .* "s",
	"HH",
	"OF" .* ["1", "2"],
	"STS" .* ["2", "4"],
	"W",
	"2",
	"3"
)

println(create_csv_for("../instances/" .* DATASET_C))


text:
function naiveDFT(hk::Array{T}) where T
    N = length(hk);
    Wnk = Array{Complex{T}}(undef, N, N)

    prefactor = -2 * π * im / N

    for n = 0:(N - 1), k = 0:(N - 1)
        Wnk[k + 1, n + 1] = exp(prefactor * n * k)
    end

    return Wnk * hk
end

function naiveDFTIterative(hk::Array{T}) where T
    N = length(hk);
    Hk = zeros(Complex{T}, N)

    prefactor = -2 * π * im / N

    @inbounds @simd for k = 1:N
        s = zero(Complex{T})
        km1 = k - 1
        for j = 1:N
            s += exp(prefactor * (j - 1) * km1) * hk[j]
        end
        Hk[k] = s
    end

    return Hk
end

max_stars_repo_name: JulienKluge/ComputationalPhotonics

max_issues_repo_name: JulienKluge/ComputationalPhotonics

max_forks_repo_name: JulienKluge/ComputationalPhotonics


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

using Markdown
using InteractiveUtils

# ╔═╡ e00967d0-89b4-11eb-031d-c128221e86e1
begin
	using PlutoUI
	using NativeSVG
end

# ╔═╡ f2f3ba34-89b3-11eb-2577-83e1dcfdc1fc
md"""# Case Study: Data Structure Selection

At this point you have learned about Julia’s core data structures, and you have seen some of the algorithms that use them.

This chapter presents a case study with exercises that let you think about choosing data structures and practice using them.
"""

# ╔═╡ def6b6fc-89b4-11eb-2e66-b5cc3e72779d
md"""## Word Frequency Analysis

As usual, you should at least attempt the exercises before you read my solutions.

#### Exercise 13-1

Write a program that reads a file, breaks each line into words, strips whitespace and punctuation from the words, and converts them to lowercase.

!!! tip
    The function isletter tests whether a character is alphabetic.

#### Exercise 13-2

Go to [Project Gutenberg](https://www.gutenberg.org/) and download your favorite out-of-copyright book in plain text format.

Modify your program from the previous exercise to read the book you downloaded, skip over the header information at the beginning of the file, and process the rest of the words as before.

Then modify the program to count the total number of words in the book, and the number of times each word is used.

Print the number of different words used in the book. Compare different books by different authors, written in different eras. Which author uses the most extensive vocabulary?

#### Exercise 13-3

Modify the program from the previous exercise to print the 20 most frequently used words in the book.

#### Exercise 13-4

Modify the previous program to read a word list and then print all the words in the book that are not in the word list. How many of them are typos? How many of them are common words that should be in the word list, and how many of them are really obscure?
"""

# ╔═╡ 9016b5ea-89b5-11eb-062f-b59a784f61da
md"""## Random Numbers

Given the same inputs, most computer programs generate the same outputs every time, so they are said to be deterministic. Determinism is usually a good thing, since we expect the same calculation to yield the same result. For some applications, though, we want the computer to be unpredictable. Games are an obvious example, but there are more.

Making a program truly nondeterministic turns out to be difficult, but there are ways to make it at least seem nondeterministic. One of them is to use algorithms that gen‐ erate *pseudorandom* numbers. Pseudorandom numbers are not truly random because they are generated by a deterministic computation, but just by looking at the numbers it is all but impossible to distinguish them from random.

The function `rand` returns a random float between `0.0` and `1.0` (including `0.0` but not `1.0`). Each time you call rand, you get the next number in a long series. To see a sample, run this loop:

```julia
for i in 1:10
	x = rand()
	println(x)
end
```

The function `rand` can take an iterator or array as an argument and return a random element:

```julia
for i in 1:10
	x = rand(1:6)
	print(x, " ")
end
```

#### Exercise 13-5

Write a function named choosefromhist that takes a histogram as defined in “Dic‐ tionaries as Collections of Counters” and returns a random value from the histogram, chosen with probability in proportion to frequency. For example, for this histogram:

```julia
julia> t = ['a', 'a', 'b'];

julia> histogram(t)
Dict{Any,Any} with 2 entries:
  'a' => 2
  'b' => 1
```

your function should return `'a'` with probability `32` and `'b'` with probability `13`.
"""

# ╔═╡ 86793ed0-89b6-11eb-1759-cfa1ea349a6b
md"""## Word Histogram

You should attempt the previous exercises before you go on. You will also need the *emma.txt* file available from this book’s GitHub repository.

Here is a program that reads a file and builds a histogram of the words in the file:

```julia
function processfile(filename)
	hist = Dict()
	for line in eachline(filename)
		processline(line, hist)
	end
	hist
end

function processline(line, hist)
	line = replace(line, '-' => ' ')
	for word in split(line)
    	word = string(filter(isletter, [word...])...)
        word = lowercase(word)
        hist[word] = get(hist, word, 0) + 1
	end
end

hist = processfile("emma.txt");
```

This program reads *emma.txt*, which contains the text of *Emma* by Jane Austen.

`processfile` loops through the lines of the file, passing them one at a time to `processline`. The histogram `hist` is being used as an accumulator.

`processline` uses the function `replace` to replace hyphens with spaces before using `split` to break the line into an array of strings. It traverses the array of words and uses `filter`, `isletter`, and `lowercase` to remove punctuation and convert to lower case. (It is shorthand to say that strings are “converted”; remember that strings are immutable, so a function like `lowercase` returns new strings.)

Finally, `processline` updates the histogram by creating a new item or incrementing an existing one.

To count the total number of words in the file, we can add up the frequencies in the histogram:

```julia
totalwords(hist) = sum(values(hist))
```

The number of different words is just the number of items in the dictionary:

```julia
differentwords(hist) = length(hist)
```

Here is some code to print the results:

```julia
julia> println("Total number of words: ", totalwords(hist))
Total number of words: 162742
julia> println("Number of different words: ", differentwords(hist)) Number of different words: 7380
```
"""

# ╔═╡ 00bea97c-89b8-11eb-23a1-db0c7d6c7bd1
md"""## Most Common Words

To find the most common words, we can make an array of tuples, where each tuple contains a word and its frequency, and sort it. The following function takes a histogram and returns an array of word-frequency tuples:

```julia
function mostcommon(hist)
	t = []
	for (key, value) in hist
		push!(t, (value, key))
	end
	reverse(sort(t))
end
```

In each tuple, the frequency appears first, so the resulting array is sorted by frequency. Here is a loop that prints the 10 most common words:

```julia
t = mostcommon(hist)
println("The most common words are:")
for (freq, word) in t[1:10]
	println(word, "\t", freq)
end
```

I use a tab character (`'\t'`) as a “separator,” rather than a space, so the second column is lined up. Here are the results from *Emma*:

```
The most common words are:
to  5295
the 5266
and 4931
of  4339
i   3191
a   3155
it  2546
her 2483
was 2400
she 2364
```

!!! tip
    This code can be simplified using the `rev` keyword argument of the `sort` function. You can read about it in the documentation.
"""

# ╔═╡ 87a74886-89b8-11eb-29aa-a9c5075a1f95
md"""## Optional Parameters

We have seen built-in functions that take optional arguments. It is possible to write programmer-defined functions with optional arguments, too. For example, here is a function that prints the most common words in a histogram:

```julia
function printmostcommon(hist, num=10)
	t = mostcommon(hist)
	println("The most common words are: ")
	for (freq, word) in t[1:num]
		println(word, "\t", freq)
	end
end
```

The first parameter is required; the second is optional. The *default value* of `num` is `10`. If you only provide one argument:

```julia
printmostcommon(hist)
```

`num` gets the default value. If you provide two arguments:

```julia
printmostcommon(hist, 20)
```

`num` gets the value of the argument instead. In other words, the optional argument
overrides the default value.

If a function has both required and optional parameters, all the required parameters have to come first, followed by the optional ones.
"""

# ╔═╡ d5579180-89b8-11eb-27ec-d1b4e811c3cb
md"""## Dictionary Subtraction

Finding the words from a book that are not in the word list from *words.txt* is a problem you might recognize as set subtraction; that is, we want to find all the words from one set (the words in the book) that are not in the other (the words in the list).

`subtract` takes dictionaries `d1` and `d2` and returns a new dictionary that contains all the keys from `d1` that are not in `d2`. Since we don’t really care about the values, we set them all to `nothing`:

```julia
function subtract(d1, d2) res = Dict()
	for key in keys(d1)
		if key ∉ keys(d2)
			res[key] = nothing
		end
	end
	res
end
```

To find the words in the book you downloaded that are not in *words.txt*, you can use `processfile` to build a histogram for *words.txt*, and then `subtract`:

```julia
words = processfile("words.txt")
diff = subtract(hist, words)
println("Words in the book that aren't in the word list:")
for word in keys(diff)
	print(word, " ")
end
```

Here are some of the results from *Emma*:

```
Words in the book that aren't in the word list:
outree quicksighted outwardly adelaide rencontre jeffereys unreserved dixons
betweens ...
```

Some of these words are names and possessives. Others, like “rencontre,” are no longer in common use. But a few are common words that should really be in the list!

#### Exercise 13-6

Julia provides a data structure called `Set` that provides many common set operations. You can read about them in “Collections and Data Structures”, or read the documentation.

Write a program that uses set subtraction to find words in the book that are not in the word list.
"""

# ╔═╡ 4ee3d86e-89ba-11eb-1dd8-99e9a315b304
md"""## Random Words

To choose a random word from the histogram, the simplest algorithm is to build an array with multiple copies of each word, according to the observed frequency, and then choose from the array:

```julia
function randomword(h)
	t=[]
	for (word, freq) in h
		for i in 1:freq
			push!(t, word)
		end
	end
	rand(t)
end
```

This algorithm works, but it is not very efficient; each time you choose a random word it rebuilds the array, which is as big as the original book. An obvious improvement is to build the array once and then make multiple selections, but the array is still big.

An alternative is:

1. Use `keys` to get an array of the words in the book.
2. Build an array that contains the cumulative sum of the word frequencies (see “Exercise 10-2”). The last item in this array is the total number of words in the book, *n*.
3. Choose a random number from 1 to *n*. Use a bisection search (see “Exercise 10-10”) to find the index where the random number would be inserted in the cumulative sum.
4. Use the index to find the corresponding word in the word array.

#### Exercise 13-7

Write a program that uses this algorithm to choose a random word from the book.
"""

# ╔═╡ e697de38-89ba-11eb-0e7c-cb53a52eb709
md"""## Markov Analysis

If you choose words from the book at random, you can get a sense of the vocabulary, but you probably won’t get a sentence:

```
this the small regard harriet which knightley's it most things
```

A series of random words seldom makes sense because there is no relationship between successive words. For example, in a real sentence you would expect an article like “the” to be followed by an adjective or a noun, and probably not a verb or adverb.

One way to measure these kinds of relationships is Markov analysis, which character‐ izes, for a given sequence of words, the probability of the words that might come next. For example, the song “Eric, the Half a Bee” (by Monty Python) begins:

> Half a bee, philosophically,
>
> Must, ipso facto, half not be.
>
> But half the bee has got to be
>
> Vis a vis, its entity. D’you see?
>
> But can a bee be said to be
>
> Or not to be an entire bee
>
> When half the bee is not a bee
>
> Due to some ancient injury?

In this text, the phrase “half the” is always followed by the word “bee,” but the phrase “the bee” might be followed by either “has” or “is.”

The result of Markov analysis is a mapping from each prefix (like “half the” and “the bee”) to all possible suffixes (like “has” and “is”).

Given this mapping, you can generate random text by starting with any prefix and choosing at random from the possible suffixes. Next, you can combine the end of the prefix and the new suffix to form the next prefix, and repeat.

For example, if you start with the prefix “Half a,” then the next word has to be “bee,” because the prefix only appears once in the text. The next prefix is “a bee,” so the next suffix might be “philosophically,” “be,” or “due.”

In this example the length of the prefix is always 2, but you can do Markov analysis with any prefix length.

#### Exercise 13-8

Give Markov analysis a try.

1. Write a program to read text from a file and perform Markov analysis. The result should be a dictionary that maps from prefixes to a collection of possible suffixes. The collection might be an array, tuple, or dictionary; it is up to you to make an appropriate choice. You can test your program with prefix length 2, but you should write the program in a way that makes it easy to try other lengths.
2. Add a function to the previous program to generate random text based on the Markov analysis. Here is an example from *Emma* with prefix length 2:
   > “He was very clever, be it sweetness or be angry, ashamed or only amused, at such a stroke. She had never thought of Hannah till you were never meant for me?” “I cannot make speeches, Emma:” he soon cut it all himself.”
   For this example, I left the punctuation attached to the words. The result is almost syntactically correct, but not quite. Semantically, it almost makes sense, but not quite.
   What happens if you increase the prefix length? Does the random text make more sense?
3. Once your program is working, you might want to try a mash-up: if you combine text from two or more books, the random text you generate will blend the vocabulary and phrases from the sources in interesting ways.

*Credit*: This case study is based on an example from The Practice of Programming by Brian Kernighan and Rob Pike (Addison-Wesley).

!!! tip
    You should attempt this exercise before you go on.
"""

# ╔═╡ aa416bf8-89bb-11eb-00ce-1dadb76ecfc1
md"""## Data Structures

Using Markov analysis to generate random text is fun, but there is also a point to this exercise: data structure selection. In completing the previous exercise, you had to choose:

* How to represent the prefixes
* How to represent the collection of possible suffixes
* How to represent the mapping from each prefix to the collection of possible suffixes

The last one is easy: a dictionary is the obvious choice for a mapping from keys to corresponding values.

For the prefixes, the most obvious options are a string, array of strings, or tuple of strings.

For the suffixes, one option is an array; another is a histogram (dictionary).

How should you choose? The first step is to think about the operations you will need to implement for each data structure. For the prefixes, you need to be able to remove words from the beginning and add to the end. For example, if the current prefix is “Half a,” and the next word is “bee,” you need to be able to form the next prefix, “a bee.”

Your first choice might be an array, since it is easy to add and remove elements.

For the collection of suffixes, the operations you need to perform include adding a new suffix (or increasing the frequency of an existing one) and choosing a random suffix.

Adding a new suffix is equally easy for the array implementation or the histogram. Choosing a random element from an array is easy; choosing from a histogram is harder to do efficiently (see “Exercise 13-7”).

So far we have been talking mostly about ease of implementation, but there are other factors to consider in choosing data structures. One is runtime. Sometimes there is a theoretical reason to expect one data structure to be faster than another; for example, I mentioned that the in operator is faster for dictionaries than for arrays, at least when the number of elements is large.

But often you don’t know ahead of time which implementation will be faster. One option is to implement both of them and see which is better. This approach is called benchmarking. A practical alternative is to choose the data structure that is easiest to implement, and then see if it is fast enough for the intended application. If so, there is no need to go on. If not, there are tools, like the Profile module, that can identify the places in a program that take the most time.

The other factor to consider is storage space. For example, using a histogram for the collection of suffixes might take less space because you only have to store each word once, no matter how many times it appears in the text. In some cases, saving space can also make your program run faster, and in the extreme, your program might not run at all if you run out of memory. But for many applications, space is a secondary consideration after runtime.

One final thought: in this discussion, I have implied that you should use one data structure for both analysis and generation. But since these are separate phases, it would also be possible to use one structure for analysis and then convert to another structure for generation. This would be a net win if the time saved during generation exceeded the time spent in conversion.

!!! tip
    The Julia package `DataStructures` implements a variety of data structures that are tailored to specific problems, such as an ordered dictionary whose entries can be iterated over deterministically.
"""

# ╔═╡ 34c48de6-89bc-11eb-2f60-0de81d83b233
md"""## Debugging
When you are debugging a program, and especially if you are working on a hard bug, there are five things to try:

*Reading*:
Examine your code, read it back to yourself, and check that it says what you meant to say.

*Running*:
Experiment by making changes and running different versions. Often if you dis‐ play the right thing at the right place in the program, the problem becomes obvi‐ ous, but sometimes you have to build scaffolding.

*Ruminating*:
Take some time to think! What kind of error is it: syntax, runtime, or semantic? What information can you get from the error messages, or from the output of the program? What kind of error could cause the problem you’re seeing? What did you change last, before the problem appeared?

*Rubberducking*:
If you explain the problem to someone else, you sometimes find the answer before you finish asking the question. Often you don’t need the other person; you could just talk to a rubber duck. And that’s the origin of the well-known strategy called rubber duck debugging. I’m not making this up!

*Retreating*:
At some point, the best thing to do is back off, undoing recent changes, until you get back to a program that works and that you understand. Then you can start rebuilding.

Beginning programmers sometimes get stuck on one of these activities and forget the others. Each activity comes with its own failure mode.

For example, reading your code might help if the problem is a typographical error, but not if the problem is a conceptual misunderstanding. If you don’t understand what your program does, you can read it 100 times and never see the error, because the error is in your head.

Running experiments can help, especially if you run small, simple tests. But if you run experiments without thinking or reading your code, you might fall into a pattern I call “random walk programming,” which is the process of making random changes until the program does the right thing. Needless to say, random walk programming can take a long time.

You have to take time to think. As I’ve said already, debugging is like an experimental science. You should have at least one hypothesis about what the problem is. If there are two or more possibilities, try to think of a test that would eliminate one of them.

But even the best debugging techniques will fail if there are too many errors, or if the code you are trying to fix is too big and complicated. Sometimes the best option is to retreat, simplifying the program until you get to something that works and that you understand.

Beginning programmers are often reluctant to retreat because they can’t stand to delete a line of code (even if it’s wrong). If it makes you feel better, copy your program into another file before you start stripping it down. Then you can copy the pieces back one at a time.

Finding a hard bug requires reading, running, ruminating, and sometimes retreating. If you get stuck on one of these activities, try the others.
"""

# ╔═╡ 6a8e4a20-89bc-11eb-149c-ed4f0c5d3697
md"""## Glossary

*deterministic*:
Pertaining to a program that does the same thing each time it runs, given the same inputs.

*pseudorandom*:
Pertaining to a sequence of numbers that appears to be random, but is generated by a deterministic program.

*default value*:
The value given to an optional parameter if no argument is provided.

*override*:
To replace a default value with an argument.

*benchmarking*:
The process of choosing between data structures by implementing alternatives and testing them on a sample of the possible inputs.

*rubberducking*:
Debugging by explaining your problem to an inanimate object such as a rubber duck. Articulating the problem can help you solve it, even if the rubber duck doesn’t know Julia.
"""

# ╔═╡ 97b91802-89bc-11eb-1adc-856eeae921ab
md"""## Exercises Exercise 13-9

The “rank” of a word is its position in an array of words sorted by frequency: the most common word has rank 1, the second most common has rank 2, etc.

[Zipf’s law](https://en.wikipedia.org/wiki/Zipf%27s_law) describes a relationship between the ranks and frequencies of words in natural languages. Specifically, it predicts that the frequency, ``f`` , of the word with rank ``r`` is:

```math
f = cr^{−s}
```

where ``s`` and ``c`` are parameters that depend on the language and the text. If you take the logarithm of both sides of this equation, you get:

```math
\log f = \log c − s \log r
```

So if you plot ``\log f`` versus ``\log r``, you should get a straight line with slope −``s`` and intercept ``\log c``.

Write a program that reads a text from a file, counts word frequencies, and prints one line for each word, in descending order of frequency, with ``\log f`` and ``\log r``.

Install a plotting library:

```julia
julia> using Pkg

julia> pkg"add Plots"
```

Its usage is very easy:
```julia

julia> using Plots

julia> x = 1:10;

julia> y = x.^2;

julia> plot(x, y)
```

Use the Plots library to plot the results and check whether they form a straight line.
"""

# ╔═╡ Cell order:
# ╟─e00967d0-89b4-11eb-031d-c128221e86e1
# ╟─f2f3ba34-89b3-11eb-2577-83e1dcfdc1fc
# ╟─def6b6fc-89b4-11eb-2e66-b5cc3e72779d
# ╟─9016b5ea-89b5-11eb-062f-b59a784f61da
# ╟─86793ed0-89b6-11eb-1759-cfa1ea349a6b
# ╟─00bea97c-89b8-11eb-23a1-db0c7d6c7bd1
# ╟─87a74886-89b8-11eb-29aa-a9c5075a1f95
# ╟─d5579180-89b8-11eb-27ec-d1b4e811c3cb
# ╟─4ee3d86e-89ba-11eb-1dd8-99e9a315b304
# ╟─e697de38-89ba-11eb-0e7c-cb53a52eb709
# ╟─aa416bf8-89bb-11eb-00ce-1dadb76ecfc1
# ╟─34c48de6-89bc-11eb-2f60-0de81d83b233
# ╟─6a8e4a20-89bc-11eb-149c-ed4f0c5d3697
# ╟─97b91802-89bc-11eb-1adc-856eeae921ab


text:
struct Callbacks
    allocate_memory::Ptr{Cvoid}
    deallocate_memory::Ptr{Cvoid}
    get_os_pid::Ptr{Cvoid}
    enable_notify_shared_library::Ptr{Cvoid}
    disable_notify_shared_library::Ptr{Cvoid}
    get_symbol_address::Ptr{Cvoid}
    add_breakpoint::Ptr{Cvoid}
    remove_breakpoint::Ptr{Cvoid}
    set_breakpoint_state::Ptr{Cvoid}
    log_message::Ptr{Cvoid}
end

struct ArchitectureID
    handle::UInt64
end

struct ProcessID
    handle::UInt64
end

struct CodeObjectID
    handle::UInt64
end

struct AgentID
    handle::UInt64
end

struct QueueID
    handle::UInt64
end

struct DispatchID
    handle::UInt64
end

struct WaveID
    handle::UInt64
end

struct DisplacedSteppingID
    handle::UInt64
end

struct RegisterClassID
    handle::UInt64
end

struct RegisterID
    handle::UInt64
end

struct AddressClassID
    handle::UInt64
end

struct AddressSpaceID
    handle::UInt64
end

struct EventID
    handle::UInt64
end

struct SharedLibraryID
    handle::UInt64
end

struct BreakpointID
    handle::UInt64
end


text:
using Test, MendelPlots, TextParse, DataFrames, Compose, Cairo

const datadir = joinpath(dirname(@__FILE__), "data")
const datafile = datadir * "/gwasresults.txt"
data, colnames = csvread(datafile, ',', header_exists = true)
df = DataFrame([i for i in data], Symbol.(colnames))

@testset "Basic Manhattan Plots (pdf and png)" begin
    @time manhattan(df)
    @time manhattan(df, outfile = "manhattan.pdf")
    @time manhattan(df, outfile = "manhattan.svg")
    @test isfile("manhattan.svg")
    @test isfile("manhattan.pdf")
    @test isfile("manhattan.png")
    rm("manhattan.png")
    rm("manhattan.svg")
    rm("manhattan.pdf")
end

@testset "Advanced Manhattan Plot" begin
    @time manhattan(df[!, :pval], df[!, :chr]; titles = "Manhattan Plot", dpi = 300,
    ylabel = "expected", xlabel = "observed", outfile = "testmanhat.png",
    ymax = 20, linecolor = "green", signifline = 1e-7)
    @test isfile("testmanhat.png")
    rm("testmanhat.png")
end

@testset "Advanced Manhattan Plot 2" begin
    @time manhattan(df[!, :pval], df[!, :chr], df[!, :pos]; titles = "Manhattan Plot", dpi = 300,
    ylabel = "expected", xlabel = "observed", outfile = "testmanhat2.png",
    ymax = 20, linecolor = "green", signifline = 1e-7, fontsize = 5pt)
    @test isfile("testmanhat2.png")
    rm("testmanhat2.png")
end

@testset "Manhattan Annotation" begin
    @time manhattan(df;
    annotateinds = [4380
    5470
    6722
    7374
   12061
    4936
   14978
    7481],
    titles = "Manhattan Plot", dpi = 300,
    outfile = "testmanhat2.png",
    linecolor = "red", signifline = -log10(1e-5), fontsize = 15pt)
    @test isfile("testmanhat2.png")
    rm("testmanhat2.png")
end

@testset "Basic QQ Plots" begin
    @time qq(df; titles = "QQ Plot", outfile = "qqplot.svg")
    @time qq(df; titles = "QQ Plot", outfile = "qqplot.pdf")
    @test isfile("qqplot.svg")
    @test isfile("qqplot.pdf")
    rm("qqplot.pdf")
    rm("qqplot.svg")
end

@testset "Advanced QQ Plot" begin
    @time qq(df; titles = "QQ Plot", dpi = 300, outfile = "testqq.png",
    xmax = 10, xmin = 0, ymin = 1, ymax = 20)
    @test isfile("testqq.png")
    rm("testqq.png")
end

@testset "ArgumentErrors" begin
    @test_throws ArgumentError manhattan(df, pvalvar = "bad")
    @test_throws ArgumentError manhattan(df, chrvar = "chrs")
    @test_throws ArgumentError manhattan(df, outfile = "bad.jpeg")
    @test_throws ArgumentError qq(df, outfile = "bad.jpeg")
    @test_throws ArgumentError qq(df, pvalvar = "pvals", outfile = "bad.png")

end


text:
module Turing

# Code associated with running probabilistic programs as tasks
include("trace/trace.jl")

using Distributions
using Turing.Traces

# Turing essentials - modelling macros and inference algorithms
export @model, @assume, @observe, @predict, InferenceAlgorithm, IS, SMC, PG, sample

# Turing-safe data structures and associated functions
export TArray, tzeros

# Debugging helpers
export dprintln

# Inference code
include("distributions/distributions.jl")
include("core/util.jl")
include("core/compiler.jl")
include("core/intrinsic.jl")
include("core/conditional.jl")
include("core/container.jl")
include("core/io.jl")
include("samplers/sampler.jl")

## global data structures
const TURING = Dict{Symbol, Any}()
global sampler = nothing
global debug_level = 0

# debugging print function: The first argument controls the verbosity of message,
#  e.g. larger v leads to more verbose debugging messages.
dprintln(v, args...) = v > Turing.debug_level ? println(args...) : nothing

end


text:
using PyDSTool, PyCall

ics = Dict("x"=>1,"y"=>0.4)
pars = Dict("k"=>0.1,"m"=>0.5)
x_rhs = "y"
y_rhs = "-k*x/m"
vars = Dict("x"=>x_rhs,"y"=>y_rhs)
tspan = [0,30]
# keys()...

dsargs = build_ode("SHM",ics,pars,vars,tspan)
d = solve_ode(dsargs)

#using Plots
#plot(d[:t],[d[:x],d[:y]])


text:
mutable struct DescriptorPoolCreateInfo
    mSizes::Vector{vk.VkDescriptorPoolSize}
    mFlags::vk.VkDescriptorPoolCreateFlags
    mMaxSets::UInt32

    mHandleRef::Ref{vk.VkDescriptorPoolCreateInfo}

    function DescriptorPoolCreateInfo()
        this = new()
        this.mSizes = Vector{vk.VkDescriptorPoolSize}()
        this.mFlags = 0
        this.mMaxSets = 0
        return this
    end
end

function handleRef(this::DescriptorPoolCreateInfo)
    this.mHandleRef = Ref(vk.VkDescriptorPoolCreateInfo(
        vk.VK_STRUCTURE_TYPE_DESCRIPTOR_POOL_CREATE_INFO, #sType::VkStructureType
        C_NULL, #pNext::Ptr{Cvoid}
        this.mFlags, #flags::VkDescriptorPoolCreateFlags
        this.mMaxSets, #maxSets::UInt32
        length(this.mSizes), #poolSizeCount::UInt32
        pointer(this.mSizes) #pPoolSizes::Ptr{VkDescriptorPoolSize}
    ))
end

function allowFreeing(this::DescriptorPoolCreateInfo, val::Bool = true)
    if val
        this.mFlags |= vk.VK_DESCRIPTOR_POOL_CREATE_FREE_DESCRIPTOR_SET_BIT
    else
        this.mFlags &= (~Int(vk.VK_DESCRIPTOR_POOL_CREATE_FREE_DESCRIPTOR_SET_BIT))
    end
end

function setMaxSets(this::DescriptorPoolCreateInfo, count::UInt32)
    this.mMaxSets = count
end

function addSize(this::DescriptorPoolCreateInfo, type::vk.VkDescriptorType, count::UInt32)
    push!(this.mSizes, vk.VkDescriptorPoolSize(type, count))
end


text:
"""
```
forecast(m, system, z0; enforce_zlb = false, shocks = Matrix{S}(undef, 0,0))

forecast(m, altpolicy, z0, shocks)

forecast(system, z0, shocks; enforce_zlb = false)
```

The first method produces a forecast, given a state space system,
initial state, and shocks, using information about the desired forecast
contained in `m`. It enforces the ZLB by using monetary policy shocks.

The second method is similar but differs in two ways. First,
it produces forecasts specifically when an alternative policy
is used. Second, it enforces the ZLB by treating it as a temporary
alternative policy.

The third method is an internal function used by the first two methods.

### Inputs

- `system::System{S}`: state-space system matrices
- `z0::Vector{S}`: state vector in the final historical period
- `shocks::Matrix{S}`: `nshocks` x `nperiods` matrix of shocks to use when
  forecasting. Note that in the first method, `nperiods` doesn't necessarily
  have to equal `forecast_horizons(m)`; it will be truncated or padded with
  zeros appropriately

**Method 1 and 2 only:**

- `m::AbstractDSGEModel`

**Method 2 only:**

- `altpolicy::Symbol`: Which alternative policy is being used
- `obs::Matrix{S <: Real}`: matrix of forecasted observables

where `S<:AbstractFloat`.

### Keyword Arguments

- `cond_type::Symbol`: one of `:none`, `:semi`, or `:full`, used to determine
  how many periods to forecast ahead. If `cond_type in [:semi, :full]`, the
  forecast horizon is reduced by the number of periods of conditional
  data. Defaults to `:none`.
- `enforce_zlb::Bool`: whether to enforce the zero lower bound. Defaults to
  `false`.
- `shocks::Matrix{S}`: matrix of size `nshocks` x `shock_horizon` of shock
  innovations under which to forecast. If `shock_horizon > horizon`, the extra
  periods of shocks will be ignored; if `shock_horizon < horizon`, zeros will be
  filled in for the shocks hitting the remaining forecasted periods.
- `draw_shocks::Bool`: if `isempty(shocks)`, indicates whether to draw shocks
  according to:

  1. If `forecast_tdist_shocks(m)`, draw `horizons` many shocks from a
     `Distributions.TDist(forecast_tdist_df_val(m))`
  2. Otherwise, draw `horizons` many shocks from a
     `DegenerateMvNormal(zeros(nshocks), sqrt(system[:QQ]))`

  or to set `shocks` to a `nshocks` x `horizon` matrix of zeros. Defaults to
  `false`. If `shocks` is provided as a keyword argument, this flag has no
  effect.

**Method 2 only:**

- `temporary_altpolicy_max_iter::Int`: maximum number of iterations for which the
    function attempts to enforce the ZLB as a temporary alternative policy. Defaults to 10.

- `set_zlb_regime_vals::Function`: user-provided function that adds additional regimes to
    regime-switching parameters if not enough regimes exist to impose the ZLB
    as a temporary alternative policy. Defaults to `identity`, and nothing will happen
    if this is the case.

- `tol::{<: Real}`: Tolerance for the smallest permissible value for the nominal interest rate.
    Defaults to -1e-14.

### Outputs

- `states::Matrix{S}`: matrix of size `nstates` x `horizon` of forecasted states
- `obs::Matrix{S}`: matrix of size `nobs` x `horizon` of forecasted observables
- `pseudo::Matrix{S}`: matrix of size `npseudo` x `horizon` of forecasted
  pseudo-observables
- `shocks::Matrix{S}`: matrix of size `nshocks` x `horizon` of shock innovations
"""
function forecast(m::AbstractDSGEModel, system::Union{RegimeSwitchingSystem{S}, System{S}},
                  z0::Vector{S}; cond_type::Symbol = :none, enforce_zlb::Bool = false,
                  shocks::AbstractMatrix{S} = Matrix{S}(undef, 0, 0), draw_shocks::Bool = false) where {S<:AbstractFloat}

    # Numbers of things
    nshocks = n_shocks_exogenous(m)
    horizon = forecast_horizons(m; cond_type = cond_type)

    if isempty(shocks)
        # Populate shocks matrix
        if draw_shocks
            # If Regime Switching System, need to get proper QQ in each forecast period.
            if isa(system, RegimeSwitchingSystem)
                shocks = zeros(nshocks, horizon)
                μ = zeros(nshocks)
                # Forecast regime indices (starting with 1 in 1st fcast period, going up to horizon)
                regime_inds = get_fcast_regime_inds(m, horizon, cond_type)
                # The regime in which forecast starts (indexing from presample_start as reg=1)
                sys_ind = cond_type == :none ? get_setting(m, :reg_forecast_start) :
                    max(get_setting(m, :reg_forecast_start), get_setting(m, :reg_post_conditional_end))
                for ts in regime_inds
                    σ = sqrt.(system[sys_ind, :QQ])
                    dist = if forecast_tdist_shocks(m)
                        # Use t-distributed shocks
                        ν = forecast_tdist_df_val(m)
                        DegenerateDiagMvTDist(μ, σ, ν)
                    else
                        # Use normally distributed shocks
                        DegenerateMvNormal(μ, σ)
                    end
                    shocks[:, ts] = rand(dist, length(ts))
                    sys_ind += 1
                end
            else
                μ = zeros(S, nshocks)
                σ = sqrt.(system[:QQ])
                dist = if forecast_tdist_shocks(m)
                    # Use t-distributed shocks
                    ν = forecast_tdist_df_val(m)
                    DegenerateDiagMvTDist(μ, σ, ν)
                else
                    # Use normally distributed shocks
                    DegenerateMvNormal(μ, σ)
                end
                shocks = rand(dist, horizon)
            end

            # Forecast without anticipated shocks
            if n_mon_anticipated_shocks(m) > 0
                ind_ant1 = m.exogenous_shocks[Symbol(get_setting(m, :monetary_policy_shock), :l1)]
                ind_antn = m.exogenous_shocks[Symbol(get_setting(m, :monetary_policy_shock), "l$(n_mon_anticipated_shocks(m))")]
                ant_shock_inds = ind_ant1:ind_antn
                shocks[ant_shock_inds, :] .= 0
            end
        else
            shocks = zeros(S, nshocks, horizon)
        end
    else
        # Adjust size of shocks matrix, padding with zeros or cutting off
        # periods of shocks if necessary
        shock_horizon = size(shocks, 2)
        if shock_horizon <= horizon
            shocks0 = zeros(nshocks, horizon - shock_horizon)
            shocks = hcat(shocks, shocks0)
        else
            shocks = shocks[:, 1:horizon]
        end
    end

    if haskey(get_settings(m), :n_periods_no_shocks)
        if get_setting(m, :n_periods_no_shocks) > 0
            shocks[:, 1:get_setting(m, :n_periods_no_shocks)] = zeros(size(shocks, 1), get_setting(m, :n_periods_no_shocks))
        end
    end

    # Populate shocks matrix under alternative policy, if
    # user has specified a function to do so
    alt_policy = alternative_policy(m)

    # TODO: check that these handling of the altpolicy cases work properly
    if (alt_policy.solve != identity &&
        alt_policy.forecast_init != identity)
        shocks, z0 = alt_policy.forecast_init(m, shocks, z0, cond_type = cond_type)
    end

#=
    # This code block has been moved to forecast_one_draw in drivers.jl
    # b/c we want to actually alter the smoothed states so that pgap or ygap
    # starts from the user-specified gap value and that this start is actually recorded.
    # TODO: check that these handling of the altpolicy cases work properly
    # Added separately to cover case where you're not using altpolicy but are using gensys2
    if haskey(m.settings, :pgap_type) && haskey(get_settings(m), :pgap_value)
        if get_setting(m, :pgap_type) == :ngdp
            _, z0 = ngdp_forecast_init(m, shocks, z0, cond_type = cond_type)
        end
    end

    if haskey(m.settings, :ygap_type) && haskey(get_settings(m), :ygap_value)
        if get_setting(m, :ygap_type) == :smooth_ait_gdp
            _, z0 = smooth_ait_gdp_forecast_init(m, shocks, z0, cond_type = cond_type)
        end
    end
=#

    # Get variables necessary to enforce the zero lower bound in the forecast
    ind_r = m.observables[get_setting(m, :nominal_rate_observable)]
    ind_r_sh = m.exogenous_shocks[get_setting(m, :monetary_policy_shock)]
    zlb_value = forecast_zlb_value(m)

    if isa(system, RegimeSwitchingSystem)
        forecast(m, system, z0, shocks; cond_type = cond_type, enforce_zlb = enforce_zlb,
                 ind_r = ind_r, ind_r_sh = ind_r_sh, zlb_value = zlb_value)
    else
        forecast(system, z0, shocks; enforce_zlb = enforce_zlb,
                 ind_r = ind_r, ind_r_sh = ind_r_sh, zlb_value = zlb_value)
    end
end

function forecast(system::System{S}, z0::Vector{S},
    shocks::Matrix{S}; enforce_zlb::Bool = false, ind_r::Int = -1,
    ind_r_sh::Int = -1, zlb_value::S = 0.13/4) where {S<:AbstractFloat}

    # Unpack system
    T, R, C = system[:TTT], system[:RRR], system[:CCC]
    Q, Z, D = system[:QQ], system[:ZZ], system[:DD]
    Z_pseudo, D_pseudo = system[:ZZ_pseudo], system[:DD_pseudo]

    # Setup
    nshocks = size(R, 2)
    nstates = size(T, 2)
    nobs    = size(Z, 1)
    npseudo = size(Z_pseudo, 1)
    horizon = size(shocks, 2)

    # Define our iteration function
    function iterate(z_t1, ϵ_t)
        z_t = C + T*z_t1 + R*ϵ_t

        # Change monetary policy shock to account for 0.13 interest rate bound
        if enforce_zlb
            interest_rate_forecast = getindex(D + Z*z_t, ind_r)
            if interest_rate_forecast < zlb_value
                # Solve for interest rate shock causing interest rate forecast to be exactly ZLB
                ϵ_t[ind_r_sh] = 0.
                z_t = C + T*z_t1 + R*ϵ_t
                ϵ_t[ind_r_sh] = getindex((zlb_value - D[ind_r] - Z[ind_r, :]'*z_t) / (Z[ind_r, :]' * R[:, ind_r_sh]), 1)

                # Forecast again with new shocks
                z_t = C + T*z_t1 + R*ϵ_t

                # Confirm procedure worked
                interest_rate_forecast = getindex(D + Z*z_t, ind_r)
                @assert interest_rate_forecast >= zlb_value - 0.01 "interest_rate_forecast = $interest_rate_forecast must be >= zlb_value - 0.01 = $(zlb_value - 0.01)"
            end
        end
        return z_t, ϵ_t
    end

    # Iterate state space forward
    states = zeros(S, nstates, horizon)
    states[:, 1], shocks[:, 1] = iterate(z0, shocks[:, 1])
    for t in 2:horizon
        states[:, t], shocks[:, t] = iterate(states[:, t-1], shocks[:, t])
    end

    # Apply measurement and pseudo-measurement equations
    obs    = D .+ Z*states
    pseudo = D_pseudo .+ Z_pseudo * states

    # Return forecasts
    return states, obs, pseudo, shocks
end

function forecast(m::AbstractDSGEModel, system::RegimeSwitchingSystem{S}, z0::Vector{S},
    shocks::Matrix{S}; cond_type::Symbol = :none, enforce_zlb::Bool = false, ind_r::Int = -1,
    ind_r_sh::Int = -1, zlb_value::S = 0.13/4) where {S<:AbstractFloat}

    # Determine how many regimes occur in the forecast, depending
    # on whether we need to subtract conditional forecast regimes or not
    n_fcast_reg = get_setting(m, :n_fcast_regimes)
    if cond_type != :none
        n_fcast_reg -= get_setting(m, :reg_post_conditional_end) - # Remove number of regimes switched since starting
            get_setting(m, :reg_forecast_start)                    # the conditional forecast
        if n_fcast_reg == 0 # Then no new regimes after conditional forecast ends
            n_fcast_reg = 1 # So we use the last conditional forecast regime
        end
    end

    Ts = Vector{Matrix{Float64}}(undef, n_fcast_reg)
    Rs = Vector{Matrix{Float64}}(undef, n_fcast_reg)
    Cs = Vector{Vector{Float64}}(undef, n_fcast_reg)
    Qs = Vector{Matrix{Float64}}(undef, n_fcast_reg)
    Zs = Vector{Matrix{Float64}}(undef, n_fcast_reg)
    Ds = Vector{Vector{Float64}}(undef, n_fcast_reg)
    Z_pseudos = Vector{Matrix{Float64}}(undef, n_fcast_reg)
    D_pseudos = Vector{Vector{Float64}}(undef, n_fcast_reg)

    # Determine in which regime the forecast starts, after accounting for
    # conditional forecast regimes, if applicable
    reg_fcast_cond_start = (cond_type == :none) ? get_setting(m, :reg_forecast_start) :
        max(get_setting(m, :reg_forecast_start), get_setting(m, :reg_post_conditional_end))

    # n_regimes is incorrectly being computed, as is n_fcast_reg
    # Unpack system
    for (ss_ind, sys_ind) in enumerate(reg_fcast_cond_start:get_setting(m, :n_regimes))
        Ts[ss_ind], Rs[ss_ind], Cs[ss_ind] = system[sys_ind, :TTT], system[sys_ind, :RRR], system[sys_ind, :CCC]
        Qs[ss_ind], Zs[ss_ind], Ds[ss_ind] = system[sys_ind, :QQ], system[sys_ind, :ZZ], system[sys_ind, :DD]
        Z_pseudos[ss_ind], D_pseudos[ss_ind] = system[sys_ind, :ZZ_pseudo], system[sys_ind, :DD_pseudo]
    end

    # Setup
    nshocks = size(Rs[1], 2)
    nstates = size(Ts[1], 2)
    nobs    = size(Zs[1], 1)
    npseudo = size(Z_pseudos[1], 1)
    horizon = size(shocks, 2)

    # Define our iteration function
    function iterate(z_t1, ϵ_t, T, R, C, Q, Z, D)
        z_t = C + T*z_t1 + R*ϵ_t

        # Change monetary policy shock to account for 0.13 interest rate bound
        if enforce_zlb
            interest_rate_forecast = getindex(D + Z*z_t, ind_r)
            if interest_rate_forecast < zlb_value
                # Solve for interest rate shock causing interest rate forecast to be exactly ZLB
                ϵ_t[ind_r_sh] = 0. # get forecast when MP shock
                z_t = C + T*z_t1 + R*ϵ_t # is zeroed out
                z_t_old = C + T*z_t1 + R*ϵ_t
                ϵ_t[ind_r_sh] = getindex((zlb_value - D[ind_r] - Z[ind_r, :]'*z_t) / (Z[ind_r, :]' * R[:, ind_r_sh]), 1)

                # Forecast again with new shocks
                z_t = C + T*z_t1 + R*ϵ_t

                # Confirm procedure worked
                interest_rate_forecast = getindex(D + Z*z_t, ind_r)

                # Subtract a small number to deal with numerical imprecision
                @assert interest_rate_forecast >= zlb_value - 0.01 "interest_rate_forecast = $interest_rate_forecast must be >= zlb_value - 0.01 = $(zlb_value - 0.01)."
            end
        end
        return z_t, ϵ_t
    end

    # These are the fcast regime indices
    regime_inds = get_fcast_regime_inds(m, horizon, cond_type)

    # Iterate state space forward
    states = zeros(nstates, horizon)
    obs = zeros(nobs, horizon)
    pseudo = zeros(npseudo, horizon)

    states[:, 1], shocks[:, 1] = iterate(z0, shocks[:, 1], Ts[1], Rs[1], Cs[1], Qs[1], Zs[1], Ds[1])
    obs[:, 1] = Ds[1] .+ Zs[1] * states[:, 1]
    pseudo[:, 1] = D_pseudos[1] .+ Z_pseudos[1] * states[:, 1]

    # If there's multiple regimes in forecast period, go through each set of indices. Otherwise, just take the first set
    if length(regime_inds) > 1
        for i in 1:length(regime_inds)
            ts = regime_inds[i]
            if length(ts) > 0 # Check in case the regime is empty
                if ts[1] == 1 # We already did this period in the initialization
                    ts = 2:ts[end]
                end
            end
            for t in ts
                states[:, t], shocks[:, t] = iterate(states[:, t - 1], shocks[:, t], Ts[i], Rs[i], Cs[i], Qs[i], Zs[i], Ds[i])
                obs[:, t] = Ds[i] .+ Zs[i] * states[:, t]
                pseudo[:, t] = D_pseudos[i] .+ Z_pseudos[i] * states[:, t]
            end
        end
    else
        for t in 2:horizon
            states[:, t], shocks[:, t] = iterate(states[:, t - 1], shocks[:, t], Ts[1], Rs[1], Cs[1], Qs[1], Zs[1], Ds[1])
            obs[:, t] = Ds[1] .+ Zs[1] * states[:, t]
            pseudo[:, t] = D_pseudos[1] .+ Z_pseudos[1] * states[:, t]
        end
    end

    # Return forecasts
    return states, obs, pseudo, shocks
end

# Automatic enforcing of the ZLB as a temporary alternative policy
function forecast(m::AbstractDSGEModel, altpolicy::Symbol, z0::Vector{S}, states::AbstractMatrix{S},
                  obs::AbstractMatrix{S}, pseudo::AbstractMatrix{S}, shocks::AbstractMatrix{S};
                  cond_type::Symbol = :none, uncertain_altpolicy::Bool = false,
                  temporary_altpolicy_max_iter::Int = 10, set_zlb_regime_vals::Function = identity,
                  tol::S = -1e-14) where {S <: Real}

    # Grab "original" settings" so they can be restored later
    is_regime_switch = haskey(get_settings(m), :regime_switching) ? get_setting(m, :regime_switching) : false
    is_replace_eqcond = haskey(get_settings(m), :replace_eqcond) ? get_setting(m, :replace_eqcond) : false
    is_gensys2 = haskey(get_settings(m), :gensys2) ? get_setting(m, :gensys2) : false

    # Grab some information about the forecast
    n_hist_regimes = haskey(DSGE.get_settings(m), :n_hist_regimes) ? get_setting(m, :n_hist_regimes) : 1
    has_reg_dates = haskey(get_settings(m), :regime_dates) # calculate dates of first & last regimes we need to add for the alt policy
    if has_reg_dates && length(get_setting(m, :regime_dates)) > 1
        start_altpol_date = get_setting(m, :regime_dates)[2]
        end_altpol_date   = get_setting(m, :regime_dates)[get_setting(m, :reg_post_conditional_end)] +
            Dates.Month(3) * (forecast_horizons(m; cond_type = cond_type) + 1)
    else
        start_altpol_date = date_forecast_start(m)
        end_altpol_date   = iterate_quarters(start_altpol_date, forecast_horizons(m; cond_type = cond_type) + 1)
    end

    first_zlb_regime = findfirst(obs[get_observables(m)[:obs_nominalrate], :] .<
                                 get_setting(m, :forecast_zlb_value))
    if !isnothing(first_zlb_regime) # Then there are ZLB regimes to enforce
        first_zlb_regime += n_hist_regimes
        if cond_type != :none
            if is_regime_switch
                first_zlb_regime += get_setting(m, :reg_post_conditional_end) - get_setting(m, :reg_forecast_start)
            else
                first_zlb_regime += subtract_quarters(get_setting(m, :date_conditional_end),
                                                      get_setting(m, :date_forecast_start)) + 1
            end
        end

        # Setup for loop enforcing zero rate
        states = similar(states)
        pseudo = similar(pseudo)
        obs = copy(obs) # This way, we don't overwrite the underlying obs matrix

        # Now iteratively impose ZLB periods until there are no negative rates in the forecast horizon
        orig_regimes      = haskey(get_settings(m), :n_regimes) ? get_setting(m, :n_regimes) : 1
        orig_regime_dates = haskey(get_settings(m), :regime_dates) ? get_setting(m, :regime_dates) : Dict{Int, Date}()
        for iter in 0:(size(obs, 2) - 1)
            # Calculate the number of ZLB regimes. For now, we add in a separate regime for every
            # period b/n the first and last ZLB regime in the forecast horizon. It is typically the case
            # that this is necessary anyway but not always, especially depending on the drawn shocks
            n_total_regimes = first_zlb_regime + iter + 1 # plus 1 for lift off

            # Set up parameters if there are switching parameter values.
            #
            # User needs to provide a function which takes in the model object `m`
            # and the total number of regimes (after adding the required temporary regimes),
            # and sets up regime-switching parameters for these new additional regimes.
            if set_zlb_regime_vals != identity
                set_zlb_regime_vals(m, n_total_regimes)
            end

            # Set up regime dates
            altpol_regime_dates = Dict{Int, Date}(1 => date_presample_start(m))
            for (regind, date) in zip(2:n_total_regimes, # We want to have enough regimes to fill up altpol_regime_dates properly
                                      start_altpol_date:Dates.Month(3):end_altpol_date)
                altpol_regime_dates[regind] = date
            end
            m <= Setting(:regime_dates, altpol_regime_dates)
            m <= Setting(:regime_switching, true)
            setup_regime_switching_inds!(m; cond_type = cond_type)

            # Set up replace_eqcond entries
            m <= Setting(:replace_eqcond, true)
            m <= Setting(:gensys2, true)
            replace_eqcond = Dict{Int, Function}()     # Which rule to replace with in which periods
            for regind in first_zlb_regime:(n_total_regimes - 1)
                replace_eqcond[regind] = zero_rate_replace_eq_entries # Temp ZLB rule in this regimes
            end

            replace_eqcond[n_total_regimes] = if altpolicy == :ait
                ait_replace_eq_entries # Add permanent AIT regime
            elseif altpolicy == :ngdp
                ngdp_replace_eq_entries # Add permanent NGDP regime
            elseif altpolicy == :smooth_ait_gdp
                smooth_ait_gdp_replace_eq_entries # Add permanent smooth AIT-GDP regime
            elseif altpolicy == :smooth_ait_gdp_alt
                smooth_ait_gdp_alt_replace_eq_entries # Add permanent smooth AIT-GDP (alternative specification) regime
            else # Default to historical policy
                eqcond
            end
            m <= Setting(:replace_eqcond_func_dict, replace_eqcond)

            # Recompute to account for new regimes
            system = compute_system(m; apply_altpolicy = true)

            # Forecast!
            states[:, :], obs[:, :], pseudo[:, :] = forecast(m, system, z0; cond_type = cond_type, shocks = shocks)

            # Delete extra regimes added to implement the temporary alternative policy, or else updating the parameters
            # in forecast_one will not work.
            if set_zlb_regime_vals != identity
                for p in m.parameters
                    if haskey(p.regimes, :value)
                        if length(p.regimes[:value]) > orig_regimes
                            for i in (orig_regimes + 1):length(p.regimes[:value])
                                delete!(p.regimes[:value], i)
                            end
                        end
                    end
                end
            end

            if all(obs[get_observables(m)[:obs_nominalrate], :] .> tol)
                # Restore the original number of regimes and regime dates
                m <= Setting(:n_regimes, orig_regimes)
                if isempty(orig_regime_dates)
                    delete!(get_settings(m), :regime_dates)
                else
                    m <= Setting(:regime_dates, orig_regime_dates)
                    setup_regime_switching_inds!(m; cond_type = cond_type)
                    n_reg_p1 = get_setting(m, :n_regimes) + 1
                end
                break
            end
        end
    end

    if !all(obs[get_observables(m)[:obs_nominalrate], :] .> tol)
        @warn "Unable to enforce the ZLB. Throwing NaN for forecasts"
        states .= fill(NaN, size(states))
        obs    .= fill(NaN, size(obs))
        pseudo .= fill(NaN, size(pseudo))
    end

    # Restore original settings
    m <= Setting(:regime_switching, is_regime_switch)
    m <= Setting(:replace_eqcond,   is_replace_eqcond)
    m <= Setting(:gensys2,          is_gensys2)

    return states, obs, pseudo
end

function get_fcast_regime_inds(m::AbstractDSGEModel, horizon::Int, cond_type::Symbol)

    fcast_post_cond_date = (cond_type == :none) ? date_forecast_start(m) : # If cond forecast, then we want to start forecasting
        max(date_forecast_start(m), iterate_quarters(date_conditional_end(m), 1)) # from the period after the conditional end period
    # last_date = iterate_quarters(fcast_post_cond_date, -1) # minus 1 to get one date before the forecast start
    last_date = fcast_post_cond_date
    last_ind = 1

    # Construct vector of time periods for each regime in the forecast periods after the conditional forecast
    regime_inds = Vector{UnitRange{Int}}(undef, 0)
    start_reg = (cond_type == :none) ? get_setting(m, :reg_forecast_start) :
        max(get_setting(m, :reg_forecast_start), get_setting(m, :reg_post_conditional_end))
    for i in start_reg:(get_setting(m, :n_regimes) - 1) # Last regime handled separately
        qtr_diff = subtract_quarters(get_setting(m, :regime_dates)[i + 1], last_date)
        push!(regime_inds, last_ind:(last_ind + qtr_diff - 1))
        last_ind += qtr_diff
        last_date = get_setting(m, :regime_dates)[i + 1]
    end
    regime_inds = push!(regime_inds, last_ind:horizon) # Covers case where final hist regime is also first (and only) forecast regime.

    return regime_inds
end


text:
"""
"""
function delete_nothings!(x::AbstractVector)::Nothing
    filter!(e->e≠nothing, x)
    return nothing
end

is_nothing(x::Nothing)::Bool = true
is_nothing(x::Any)::Bool = false

"""
"""
is_nothing


text:
hsep(io::IO, c='-') = println(io, repeat(c, 40))
export hsep
hsep(args...; kwargs...) = hsep(stdout::IO, args...; kwargs...)


text:
import Images

const writedir = joinpath(tempdir(), "Images")

if !isdir(writedir)
    mkdir(writedir)
end

a = rand(2,2)
fn = joinpath(writedir, "2by2.png")
Images.imwrite(a, fn)
sleep(0.2)
b = Images.imread(fn)
@assert convert(Array, b) == scale(scaleinfo(a), a)


text:
module PointwiseKDEs

export PointwiseKDE

using Distributions
using LinearAlgebra
using StatsFuns
using Random

struct PointwiseKDE{T <: Real} <: ContinuousMultivariateDistribution
    pts::Matrix{T}
    CF::Cholesky{T, Matrix{T}}
end

"""
    PointwiseKDE(points::Matrix{T <: Real})

Construct a Gaussian KDE from the given matrix of points.

`points` should have size `(ndim, npts)`.  Bandwidth selection is Scott's rule.

The result is a `MultivariateDistribution` object suitable for use with `rand`
or `logpdf` or ....
"""
function PointwiseKDE(pts::Matrix{T}) where T <: Real
    nd, np = size(pts)
    C = cov(pts, dims=2) ./ np^(2/(4+nd))
    CF = cholesky(C)

    PointwiseKDE(pts, CF)
end

function Distributions.length(d::PointwiseKDE)
    size(d.pts, 1)
end

function Distributions.sampler(d::PointwiseKDE)
    d
end

function Distributions.eltype(d::PointwiseKDE{T}) where T <: Number
    T
end

function Distributions._rand!(rng::AbstractRNG, d::PointwiseKDE{T}, x::AbstractArray{T}) where T <: Real
    nd, np = size(d.pts)
    i = rand(rng, 1:np)

    pt = vec(d.pts[:,i])
    cpt = d.CF.L*randn(rng, nd)

    for i in eachindex(x)
        x[i] = pt[i] + cpt[i]
    end
    x
end

function Distributions._rand!(rng::AbstractRNG, d::PointwiseKDE{T}, x::AbstractMatrix{T}) where T <: Real
    nd, np = size(x)

    for j in 1:np
        x[:,j] = Distributions._rand!(rng, d, x[:,j])
    end

    x
end

function Distributions._logpdf(d::PointwiseKDE{T}, x::AbstractArray{T}) where T <: Real
    nd, np = size(d.pts)

    logp::T = zero(T) - Inf
    for i in 1:np
        r = x .- vec(d.pts[:,i])
        ru = d.CF.L \ r

        logp = logaddexp(logp, -0.5*dot(ru, ru))
    end

    logp - sum(0.5*log(2*pi) .+ log.(diag(d.CF.L))) - log(np)
end

end # module


text:
# included in runtests.jl
using Base.Test
using HTTP
using HttpServer
using WebSockets
import WebSockets:  generate_websocket_key,
    websocket_handshake,
    upgrade
import HttpCommon:Request
import HttpServer:  is_websocket_handshake

function templaterequests()
    chromeheaders = Dict{String, String}( "Connection"=>"Upgrade",
                                            "Upgrade"=>"websocket")
    firefoxheaders = Dict{String, String}("Connection"=>"keep-alive, Upgrade",
                                            "Upgrade"=>"websocket")
    chromerequest = HttpCommon.Request("GET", "", chromeheaders, "")
    firefoxrequest= Request("GET", "", firefoxheaders, "")
    chromerequest_HTTP = HTTP.Messages.Request("GET", "/", collect(chromeheaders))
    firefoxrequest_HTTP = HTTP.Messages.Request("GET", "/", collect(firefoxheaders))
    return [chromerequest, firefoxrequest, chromerequest_HTTP, firefoxrequest_HTTP]
end

sethd(r::Request, pa::Pair) = push!(r.headers, pa)
sethd(r::HTTP.Messages.Request, pa::Pair) = HTTP.Messages.setheader(r, HTTP.Header(pa))

takefirstline(buf::IOBuffer) = strip(split(buf |> take! |> String, "\r\n")[1])
takefirstline(buf::BufferStream) = strip(split(buf |> read |> String, "\r\n")[1])
function handshakeresponse(request::Request)
    cli = HttpServer.Client(2, IOBuffer())
    websocket_handshake(request, cli)
    strip(takefirstline(cli.sock))
end
function handshakeresponse(request::HTTP.Messages.Request)
    buf = BufferStream()
    c = HTTP.ConnectionPool.Connection(buf)
    t = HTTP.Transaction(c)
    s = HTTP.Streams.Stream(request, t)
    upgrade(dummywshandler, s)
    close(buf)
    takefirstline(buf)
end

"""
The dummy websocket don't use TCP. Close won't work, but we can manipulate the contents
using otherwise the same functions as TCP sockets.
"""
dummyws(server::Bool)  = WebSocket(BufferStream(), server)

function dummywshandler(req, dws::WebSockets.WebSocket{BufferStream})
    close(dws.socket)
    close(dws)
end

# Test generate_websocket_key"
@test generate_websocket_key("dGhlIHNhbXBsZSBub25jZQ==") == "s3pPLMBiTxaQ9kYGzzhZRbK+xOo="

# Test reject / switch format"
io = IOBuffer()
const REJECT = "HTTP/1.1 400 Bad Request"
Base.write(io, Response(400))
@test takefirstline(io) == REJECT
Base.write(io, HTTP.Messages.Response(400))
@test takefirstline(io) == REJECT

const SWITCH = "HTTP/1.1 101 Switching Protocols"
Base.write(io, Response(101))
@test takefirstline(io) == SWITCH
Base.write(io, HTTP.Messages.Response(101))
@test takefirstline(io) == SWITCH

chromerequest, firefoxrequest, chromerequest_HTTP, firefoxrequest_HTTP = Tuple(templaterequests())
wshandler = WebSocketHandler((x,y)->nothing)
for request in [chromerequest, firefoxrequest]
    @test is_websocket_handshake(wshandler, request) == true
end

# "Test simple handshakes that are unacceptable"

for r in templaterequests()
    @test handshakeresponse(r) == REJECT
    sethd(r, "Sec-WebSocket-Version"    => "11")
    @test handshakeresponse(r) == REJECT
    sethd(r, "Sec-WebSocket-Version"    => "13")
    @test handshakeresponse(r) == REJECT
    sethd(r,   "Sec-WebSocket-Key"      => "mumbojumbobo")
    @test handshakeresponse(r) == REJECT
    sethd(r,  "Sec-WebSocket-Key"       => "17 bytes key is not accepted")
    @test handshakeresponse(r) == REJECT
    sethd(r,  "Sec-WebSocket-Key"       => "16 bytes key this is surely")
    sethd(r,  "Sec-WebSocket-Protocol"  => "unsupported")
    @test handshakeresponse(r) == REJECT
end

#  Test simple handshakes, acceptable
for r in templaterequests()
    sethd(r, "Sec-WebSocket-Version"    => "13")
    sethd(r,  "Sec-WebSocket-Key"       => "16 bytes key this is surely")
    @test handshakeresponse(r) == SWITCH
end

#  Test unacceptable subprotocol handshake subprotocol
for r in templaterequests()
    sethd(r, "Sec-WebSocket-Version"    => "13")
    sethd(r, "Sec-WebSocket-Key"        => "16 bytes key this is surely")
    sethd(r, "Sec-WebSocket-Protocol"       => "my.server/json-zmq")
    @test handshakeresponse(r) == REJECT
end

#  add simple subprotocol to acceptable list
@test true == WebSockets.addsubproto("xml")

# add subprotocol with difficult name
@test true == WebSockets.addsubproto("my.server/json-zmq")

# "Test handshake subprotocol now acceptable"
for r in templaterequests()
    sethd(r, "Sec-WebSocket-Version"    => "13")
    sethd(r,  "Sec-WebSocket-Key"        => "16 bytes key this is surely")
    sethd(r, "Sec-WebSocket-Protocol"       => "xml")
    @test handshakeresponse(r) == SWITCH
    sethd(r, "Sec-WebSocket-Protocol"       => "my.server/json-zmq")
    @test handshakeresponse(r) == SWITCH
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LIBMINC_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LIBMINC")
JLLWrappers.@generate_main_file("LIBMINC", UUID("e365f29a-713b-59c4-b64a-7ac3c1ec6045"))
end  # module LIBMINC_jll