text:
module RedPitayaDAQServer

# package code goes here

using Sockets
import Sockets: send, connect

using Statistics
using LinearAlgebra

import Base: reset, iterate, length

export RedPitaya, receive, query, start, stop, disconnect, getLog

mutable struct RedPitaya
  host::String
  delim::String
  socket::TCPSocket
  dataSocket::TCPSocket
  decimation::Int64
  samplesPerPeriod::Int64
  periodsPerFrame::Int64
  isConnected::Bool
  isMaster::Bool
  destroyed::Bool
end

# Iterable Interface

length(rp::RedPitaya) = 1
iterate(rp::RedPitaya, state=1) = state > 1 ? nothing : (rp, state + 1)

"""
Send a command to the RedPitaya
"""
function send(rp::RedPitaya,cmd::String)
  @debug "send command: " cmd
  write(rp.socket,cmd*rp.delim)
end

"""
Receive a String from the RedPitaya
"""
function receive(rp::RedPitaya)
  readline(rp.socket)[1:end] #-2
end

"""
Perform a query with the RedPitaya. Return String
"""
function query(rp::RedPitaya,cmd::String)
  send(rp,cmd)
  return receive(rp)
end

"""
Perform a query with the RedPitaya. Parse result as type T
"""
function query(rp::RedPitaya,cmd::String,T::Type)
  a = query(rp,cmd)
  return parse(T,a)
end

# connect with timeout implementation
function Sockets.connect(host, port::Integer, timeout::Number)
  s = TCPSocket()
  t = Timer(_ -> close(s), timeout)
  try
    connect(s, host, port)
  catch e
    error("Could not connect to $(host) on port $(port)
            since the operations was timed out after $(timeout) seconds!")
  finally
    close(t)
  end
  return s
end

const _timeout = 1.0

function connect(rp::RedPitaya)
  if !rp.isConnected
    begin
    rp.socket = connect(rp.host, 5025, _timeout)
    send(rp, string("RP:Init"))
    connectADC(rp)
    rp.dataSocket = connect(rp.host, 5026, _timeout)
    rp.isConnected = true
    decimation(rp)
    end
  end
end

function disconnect(rp::RedPitaya)
  if rp.isConnected && !rp.destroyed
    @async begin
      close(rp.socket)
      close(rp.dataSocket)
      rp.isConnected = false
    end
  end
  return nothing
end

function getLog(rp::RedPitaya, log)
  command = "RP:LOG?"
  send(rp, command)
  chunk_size = 1024
  size = read(rp.dataSocket, Int64)
  recv = 0
  while (recv != size)
    buff = read(rp.dataSocket, min(chunk_size, size - recv))
    recv = recv + length(buff)
    write(log, buff)
  end
  close(log)
end

include("DAC.jl")
include("ADC.jl")
include("Cluster.jl")
include("SlowIO.jl")

function destroy(rp::RedPitaya)
  disconnect(rp)
  rp.destroyed = true
end

function RedPitaya(host, port=5025, isMaster=true)

  rp = RedPitaya(host,"\n", TCPSocket(), TCPSocket(), 1, 1, 1, false, isMaster, false)

  connect(rp)

  #rp.decimation = decimation(rp)
  #rp.samplesPerPeriod = samplesPerPeriod(rp)
  #rp.periodsPerFrame = periodsPerFrame(rp)

  finalizer(d -> destroy(d), rp)
  return rp
end

#=
{.pattern = "RP:DAC:CHannel#:COMPonent#:AMPlitude?", .callback = RP_DAC_GetAmplitude,},
 {.pattern = "RP:DAC:CHannel#:COMPonent#:AMPlitude", .callback = RP_DAC_SetAmplitude,},
 {.pattern = "RP:DAC:CHannel#:COMPonent#:FREQuency?", .callback = RP_DAC_GetFrequency,},
 {.pattern = "RP:DAC:CHannel#:COMPonent#:FREQuency", .callback = RP_DAC_SetFrequency,},
 {.pattern = "RP:DAC:CHannel#:COMPonent#:FACtor?", .callback = RP_DAC_GetModulusFactor,},
 {.pattern = "RP:DAC:CHannel#:COMPonent#:FACtor", .callback = RP_DAC_SetModulusFactor,},
 {.pattern = "RP:DAC:CHannel#:COMPonent#:PHAse?", .callback = RP_DAC_GetPhase,},
 {.pattern = "RP:DAC:CHannel#:COMPonent#:PHAse", .callback = RP_DAC_SetPhase,},
 {.pattern = "RP:DAC:MODe", .callback = RP_DAC_SetDACMode,},
 {.pattern = "RP:DAC:MODe?", .callback = RP_DAC_GetDACMode,},
 {.pattern = "RP:DAC:CHannel#:COMPonent#:MODulus", .callback = RP_DAC_ReconfigureDACModulus,},
 {.pattern = "RP:DAC:CHannel#:COMPonent#:MODulus?", .callback = RP_DAC_GetDACModulus,},
 {.pattern = "RP:ADC:DECimation", .callback = RP_ADC_SetDecimation,},
 {.pattern = "RP:ADC:DECimation?", .callback = RP_ADC_GetDecimation,},
 {.pattern = "RP:ADC:PERiod", .callback = RP_ADC_SetSamplesPerPeriod,},
 {.pattern = "RP:ADC:PERiod?", .callback = RP_ADC_GetSamplesPerPeriod,},
 {.pattern = "RP:ADC:FRAme", .callback = RP_ADC_SetPeriodsPerFrame,},
 {.pattern = "RP:ADC:FRAme?", .callback = RP_ADC_GetPeriodsPerFrame,},
 {.pattern = "RP:ADC:FRAmes:CURRent?", .callback = RP_ADC_GetCurrentFrame,},
 {.pattern = "RP:ADC:FRAmes:DATa", .callback = RP_ADC_GetFrames,},
 {.pattern = "RP:ADC:ACQCONNect", .callback = RP_ADC_StartAcquisitionConnection,},
 {.pattern = "RP:ADC:ACQSTATus", .callback = RP_ADC_SetAcquisitionStatus,},
 {.pattern = "RP:PDM:CHannel#:NextValue", .callback = RP_PDM_SetPDMNextValue,},
 {.pattern = "RP:PDM:CHannel#:NextValue?", .callback = RP_PDM_GetPDMNextValue,},
 {.pattern = "RP:PDM:CHannel#:CurrentValue?", .callback = RP_PDM_GetPDMCurrentValue,},
 {.pattern = "RP:XADC:CHannel#?", .callback = RP_XADC_GetXADCValueVolt,},
 {.pattern = "RP:WatchDogMode", .callback = RP_WatchdogMode,},
 {.pattern = "RP:RamWriterMode", .callback = RP_RAMWriterMode,},
 {.pattern = "RP:MasterTrigger", .callback = RP_MasterTrigger,},
 {.pattern = "RP:InstantResetMode", .callback = RP_InstantResetMode,},
 {.pattern = "RP:PeripheralAResetN?", .callback = RP_PeripheralAResetN,},
 {.pattern = "RP:FourierSynthAResetN?", .callback = RP_FourierSynthAResetN,},
 {.pattern = "RP:PDMAResetN?", .callback = RP_PDMAResetN,},
 {.pattern = "RP:WriteToRAMAResetN?", .callback = RP_WriteToRAMAResetN,},
 {.pattern = "RP:XADCAResetN?", .callback = RP_XADCAResetN,},
 {.pattern = "RP:TriggerStatus?", .callback = RP_TriggerStatus,},
 {.pattern = "RP:WatchdogStatus?", .callback = RP_WatchdogStatus,},
 {.pattern = "RP:InstantResetStatus?", .callback = RP_InstantResetStatus,},
 =#

end # module

max_stars_repo_name: Fypsilonn/RedPitayaDAQServer

max_issues_repo_name: Fypsilonn/RedPitayaDAQServer

max_forks_repo_name: Fypsilonn/RedPitayaDAQServer


text:
@testset "foreach" begin
    i = Indices([1,2,3,4,5])

    tmp = Ref(0)
    foreach(x -> tmp[] += x, i)
    @test tmp[] == 15

    tmp = Ref(0)
    foreach((x,y) -> tmp[] += x + y, i, i)
    @test tmp[] == 30

    tmp = Ref(0)
    foreach((x,y) -> tmp[] += x + y, i, copy(i))
    @test tmp[] == 30

    tmp = Ref(0)
    foreach((x,y,z) -> tmp[] += x + y + z, i, i, i)
    @test tmp[] == 45

    tmp = Ref(0)
    foreach((x,y,z) -> tmp[] += x + y + z, i, i, copy(i))
    @test tmp[] == 45

    d = i .+ 1

    tmp = Ref(0)
    foreach(x -> tmp[] += x, d)
    @test tmp[] == 20

    tmp = Ref(0)
    foreach((x,y) -> tmp[] += x + y, d, d)
    @test tmp[] == 40

    tmp = Ref(0)
    foreach((x,y) -> tmp[] += x + y, d, copy(d))
    @test tmp[] == 40

    tmp = Ref(0)
    foreach((x,y,z) -> tmp[] += x + y + z, d, d, d)
    @test tmp[] == 60

    tmp = Ref(0)
    foreach((x,y,z) -> tmp[] += x + y + z, d, d, copy(d))
    @test tmp[] == 60
end

ext: jl


text:
using EnvInfo: status, save, PkgStatusInfo
using Base.Test

@testset "status" begin
    @test status("EnvInfo") isa PkgStatusInfo
    @test all(psi isa PkgStatusInfo for psi in status(["EnvInfo", "JSON"]))
    @test all(psi isa PkgStatusInfo for psi in status())
end

@testset "save" begin
    @test save("/dev/null", ["EnvInfo", "JSON"]) === nothing
    @test save("/dev/null", ["EnvInfo", "JSON"],
               include_installed = true) === nothing
end


text:
#
###
# Gather-Scatter Operators - enforce continuity/ periodicity
###

# TODO write GatherScatterOp that calls NNlib.gather, scatter
struct GatherScatter{D} <: AbstractGatherScatterOperator{Bool,D}
    global_numbering
    implementation
end

"""
Q*Q'*u where Q: local -> global operator
"""
function DSS(u,l2g,g2l)

    Qu   = NNlib.scatter(+,u,l2g)
    QQtu = NNlib.gather(Qu,g2l)

    return v
end

function Qmatrix(n::Integer, periodic::Bool)
    Q = sparse(I,n, n-1)

    if periodic
        Q[end,1] = 1
    end

    Q
end

function GatherScatter(space::AbstractSpectralSpace{<:Number,D}) where{D}
    domain = get_domain(space)
    periodic = isperiodic(domain)
    npoints = get_numpoints(space)

    if !prod(periodic...)
        return IdentityOp{D}()
    end

    Qmats = Qmatrix.(npoints, periodic)

    Q = if D == 1
        MatrixOp(Qmats...)
    elseif D == 2
        TensorProductOp2D(Qmats...)
    elseif D == 3
        TensorProductOp3D(Qmats...)
    end

    QQt = Q * Q' # replace with call to NNlib gather-scatter
end

###
# Boundary Condition application
#
# Apply this boundary condition to that boundary tag
###

struct DirichletBC end

struct BoundaryCondition{T,D} <: AbstractBoundaryCondition{T,D}
    tag
    type # dirichlet, neumann
    dirichlet_func! # (ub, space) -> mul!(ub, I, false)
#   neumann_func!
    mask # implementation
end

function BoundaryCondition(tags, space::AbstractSpace<:Number,2;
                           dirichlet_func! =nothing, neumann_func! = nothing)

    mask = generate_mask(tags, space)

    BoundaryCondition()
end

"""
 bc = (:Dirichlet,:Neumann,:Dirichlet,:Dirichlet) at (rmin, rmax, smin, smax)

 :Dirichlet = Dirichlet = zeros ∂Ω data\n
 :Neumann   = Neumann   = keeps ∂Ω data

 A periodic mesh overwrites 'D' to 'N' in direction of periodicity.

 To achieve inhomogeneous Dirichlet condition, apply the formulation
 u = ub + uh, where uh is homogeneous part, and ub is an arbitrary
 smooth function on Ω. Then, solve for uh
"""
function generate_mask(tags, space::AbstractSpace{<:Number,2})
    (nr, ns,) = space.npoints

    periodic = isperiodic(space.domain)

    Ix = sparse(I,nr,nr)
    Iy = sparse(I,ns,ns)

    ix = collect(1:(nr))
    iy = collect(1:(ns))

    if(bc[1] == :Dirichlet) ix = ix[2:end]   end
    if(bc[2] == :Dirichlet) ix = ix[1:end-1] end
    if(bc[3] == :Dirichlet) iy = iy[2:end]   end
    if(bc[4] == :Dirichlet) iy = iy[1:end-1] end

    if(periodic[1]) ix = collect(1:(nr)); end
    if(periodic[2]) iy = collect(1:(ns)); end

    Rx = Ix[ix,:]
    Ry = Iy[iy,:]

    M = diag(Rx'*Rx) * diag(Ry'*Ry)'
    M = Array(M) .== true

    return M
end

function applyBC!(u::AbstractField{<:Number,D}, bc::BoundaryCondition{<:Number,D}) where{D}

    return u
end

function applyBC!()
end
#


text:
## Make an animation of a random walk
using Plots
using GraphRecipes
using SparseArrays
using DelimitedFiles
using LinearAlgebra
## Load the graph
ei,ej = copy.(eachcol(Int.(readdlm("undir-graph-14.edges"))))
xy = readdlm("undir-graph-14.xy")'
## A simple simulation
using Random
Random.seed!(1)
X = 1 # start on node 1
A = sparse(ei,ej,1,14,14) # make the sparse matrix
A = A+A' # make it symmetric (we only store half)
α = 0.85
for i=1:10
  global X
  if rand() <= α && length(findnz(A[:,X])[1]) >= 0
    Xn = rand(findnz(A[:,X])[1]) # random neighbor
    println("PageRank Walk moves from $X to $Xn")
  else
    Xn = rand(1:size(A,1)) # random node
    println("PageRank Walk jumps from $X to $Xn")
  end
  X = Xn
end

##
ei,ej = copy.(eachcol(Int.(readdlm("wiki-simple.edges"))))
xy = readdlm("wiki-simple.xy")'
pages = readlines("wiki-simple.nodes")
## A visual simulation of the random walk
gr(size=(300,300),dpi=300)
using Random
Random.seed!(1)
function prrw_on_graph(ei, ej, xy, alpha, start, nstep, nframes)
  n = size(xy,2)
  A = sparse(ei,ej,1,n,n)
  B = A.*A'
  U = A - B # we need A' so that the ei, ej plots

  p = graphplot(findnz(triu(B,1))[1:2]..., x = xy[1,:], y = xy[2,:],
    markercolor=:lightblue, markerstrokecolor=:white, curves=false,
    markersize=7, markerstrokewidth=1,
    linecolor=1, linealpha=0.8, linewidth=0.7, arrow=:both, shorten=0.85,
    axis_buffer=0.5, background=nothing)

  graphplot!(findnz(U)[1:2]..., x = xy[1,:], y = xy[2,:],
    markercolor=nothing, markerstrokecolor=:white, curves=false, curvature_scalar = -0.4,
    markersize=0, markerstrokewidth=0,
     linecolor=2, linealpha=0.8, linewidth=0.7, arrow = :tail, shorten=0.85,
    axis_buffer=0.5, background=:black)

  X = start
  # random walk step
  #p1 = scatter!([xy[1,X]],[xy[2,X]],markersize=4, color=1)
  p2 = scatter!([xy[1,X]],[xy[2,X]],markersize=10, color=:orange)

  xlims!(xlims().*1.2)
  ylims!((-0.5,7.5))
  for i=1:10
    annotate!(xy[1,i],xy[2,i],text("$i",9,:white))
  end

  Xn = X
  jump = false
  anim = @animate for i=0:(nframes*nstep)
    stepi = mod(i,nframes) # this is the frame number in the step
    if stepi == 0
      X = Xn
      # with prob alpha, follow a rand. neig
      if rand() <= alpha && length(ej[ei .== X]) > 0
        Xn = rand(ej[ei .== X])
        jump = false
      else # with prob. (1-alpha)
        Xn = rand(1:n) # pick a random node
        jump = true
        println("Jump")
      end
    end
    # interpolate between the two positions
    if jump
      pos = stepi/(nframes-2)
      if pos >= 1
        pos = 1
      else
        pos += 0.25*randn()
      end
      pos = min(max(pos,0),1) # cap in 0,1
    else
      pos = min(stepi/(nframes-2),1)
    end
    p[5] = [(1-pos)*xy[1,X] + pos*xy[1,Xn]],[(1-pos)*xy[2,X] + pos*xy[2,Xn]]
  end
  return anim
end
anim = prrw_on_graph(ei,ej,xy,0.85,1,100,20)
gif(anim, "pr-rw-on-wiki-simple.gif", fps=40)

## Make a simple convergence figure
gr()
theme(:dark)
n = max(maximum(ei),maximum(ej))
A = sparse(ei,ej,1,n,n)
d = vec(sum(A,dims=2 ))
P = A'*Diagonal(1.0./d)
v = ones(n)/n
α = 0.85
nsteps = 25
X = zeros(n,nsteps)
X[:,1] = v
for i=2:nsteps
  X[:,i] = α*(P*X[:,i-1]) + (1-α)*v
end
#plot(X',labels=pages,size=(300,300))
plot(X',size=(300,300),
  background_color_legend = nothing,foreground_color_legend = nothing,
  labels=pages,legend=:best)
xlabel!("Step")
ylabel!("Value")
savefig("pagerank-wiki-simple-convergence.pdf")

@show sum(X[:,end])


text:
function jacobian(m::HetDSGEGovDebt)
    reset_grids!(m)
    #truncate_distribution!(m)

    # Load in endogenous state and eq cond indices
    endo = augment_model_states(m.endogenous_states_original,#m.endogenous_states_unnormalized,
                         n_model_states_original(m))
    eq   = m.equilibrium_conditions

    # Load in parameters, steady-state parameters, and grids
    r::Float64     = m[:r].scaledvalue
    α::Float64     = m[:α].value
    H::Float64     = m[:H].value
    δ::Float64     = m[:δ].value
    μ_sp::Float64  = m[:μ_sp].value
#    ρ_sp::Float64  = m[:ρ_sp].value
#    σ_sp::Float64  = m[:σ_sp].value
    γ::Float64     = m[:γ].scaledvalue
    g::Float64     = m[:g].value
    η::Float64     = m[:η].value
    ρ_B::Float64    = m[:ρ_b].value
    ρ_G::Float64    = m[:ρ_g].value
    ρ_z::Float64    = m[:ρ_z].value
    ρ_μ::Float64    = m[:ρ_μ].value
    ρ_lamw::Float64 = m[:ρ_λ_w].value
    ρ_lamf::Float64 = m[:ρ_λ_f].value
    ρ_mon::Float64  = m[:ρ_rm].value
    spp::Float64   = m[:spp].value
    ϕh::Float64    = m[:ϕh].value
    ρ_R::Float64    = m[:ρR].value
    ψπ::Float64    = m[:ψπ].value
    ψy::Float64    = m[:ψy].value
    κ_p::Float64  = m[:κ_p].value
    κ_w::Float64  = m[:κ_w].value

    Tg::Float64 = m[:Tg].value
    δb::Float64 = m[:δb].value
    bg::Float64 = m[:bg].value

    R = 1 + r

    ell::Vector{Float64}  = m[:lstar].value
    c::Vector{Float64}    = m[:cstar].value
    μ::Vector{Float64}    = m[:μstar].value
    β::Float64            = m[:βstar].value

    T::Float64     = m[:Tstar].value
    ω::Float64     = m[:ωstar].value
    xstar::Float64 = m[:xstar].value
    ystar::Float64 = m[:ystar].value
    Rk::Float64    = m[:Rkstar].value
    kstar::Float64 = m[:kstar].value

    xgrid::Vector{Float64} = m.grids[:xgrid].points
    xwts::Vector{Float64}  = m.grids[:xgrid].weights
    sgrid::Vector{Float64} = m.grids[:sgrid].points
    swts::Vector{Float64}  = m.grids[:sgrid].weights
    fgrid::Matrix{Float64} = m.grids[:fgrid]
    xswts = kron(swts,xwts)

    zlo::Float64 = m[:zlo].value
    zhi::Float64 = m[:zhi].value

    nx::Int = get_setting(m, :nx)
    ns::Int = get_setting(m, :ns)
    nxns = nx*ns

    qp(z) = dmollifier_hetdsgegovdebt(z, zhi, zlo)
    qfunction_hetdsgegovdebt(x) = mollifier_hetdsgegovdebt(x, zhi, zlo) #/sumz

    unc = 1 ./ ell .<= repeat(xgrid,ns) .+ η

    dF1_dELL, dF1_dRZ, dF1_dELLP, dF1_dWHP, dF1_dTTP, ee =
        euler_equation_hetdsgegovdebt(nx, ns, qp, qfunction_hetdsgegovdebt, xgrid, sgrid, fgrid, unc, xswts,
                               R, γ, β, η, ell, T, ω, H)

    # KF Equation
    dF2_dWH, dF2_dRZ, dF2_dTT,dF2_dELL, bigΨ, dF2_dM =
        kolmogorov_fwd_hetdsgegovdebt(nx, ns, qfunction_hetdsgegovdebt, qp, xgrid, sgrid, fgrid, unc, xswts,
                               R, γ, ell, μ, η, T, ω, H, ee)

    # Market clearing, lambda function
    c = min.(1 ./ ell,repeat(xgrid,ns).+η)
    lam = (xswts.*μ)'*(1 ./ c) # average marginal utility which the union uses to set wages
    aggc = (xswts .* μ)'c
    ϕ = lam*ω/(H^ϕh) # now that we know lam in steady state, choose disutility to target hours H

    setup_indices!(m)
    normalize_model_state_indices!(m)

    nvars = get_setting(m, :nvars)
    # Make the Jacobian
    JJ = zeros(nvars, 2*nvars)

    # Euler equation
    JJ[eq[:eq_euler],endo[:l′_t]] = dF1_dELLP
    JJ[eq[:eq_euler],endo[:z′_t]]   = -dF1_dRZ
    JJ[eq[:eq_euler],endo[:w′_t]]   = dF1_dWHP
    JJ[eq[:eq_euler],endo[:L′_t]]  = dF1_dWHP
    JJ[eq[:eq_euler],endo[:t′_t]]  = dF1_dTTP
    JJ[eq[:eq_euler],endo[:b_t]]    = ell
    JJ[eq[:eq_euler],endo[:l_t]]  = dF1_dELL
    JJ[eq[:eq_euler],endo[:R_t]]   = dF1_dRZ

    # KF eqn
    JJ[eq[:eq_kolmogorov_fwd],endo[:kf′_t]]   = -Matrix{Float64}(I, nxns, nxns)
    JJ[eq[:eq_kolmogorov_fwd],endo[:kf_t]]   = dF2_dM
    JJ[eq[:eq_kolmogorov_fwd],endo[:l_t]] = dF2_dELL
    JJ[eq[:eq_kolmogorov_fwd],endo[:R_t]] = dF2_dRZ
    JJ[eq[:eq_kolmogorov_fwd],endo[:z_t]]    = -dF2_dM*dF2_dRZ
    JJ[eq[:eq_kolmogorov_fwd],endo[:w_t]]    = dF2_dM*dF2_dWH
    JJ[eq[:eq_kolmogorov_fwd],endo[:L_t]]   = dF2_dM*dF2_dWH
    JJ[eq[:eq_kolmogorov_fwd],endo[:t_t]]   = dF2_dM*dF2_dTT

    # aggregate consumption
    JJ[first(eq[:eq_agg_consumption]),first(endo[:C_t])]   = -aggc
    JJ[first(eq[:eq_agg_consumption]), endo[:l_t]] = -(μ .*unc.*xswts.*c)'
    JJ[first(eq[:eq_agg_consumption]), endo[:kf_t]] = (xswts.*c)' # note, now we linearize
    JJ[first(eq[:eq_agg_consumption]),first(endo[:z_t])]   = -(xswts.*c)'*dF2_dRZ
    JJ[first(eq[:eq_agg_consumption]),first(endo[:w_t])]   = (xswts.*c)'*dF2_dWH
    JJ[first(eq[:eq_agg_consumption]),first(endo[:L_t])]   = (xswts.*c)'*dF2_dWH
    JJ[first(eq[:eq_agg_consumption]),first(endo[:t_t])]   = (xswts.*c)'*dF2_dTT

    # lambda = average marginal utility
    JJ[first(eq[:eq_lambda]),first(endo[:margutil_t])] = lam
    JJ[first(eq[:eq_lambda]),endo[:kf_t]]   = -(xswts./c)' # note, now we linearize
    JJ[first(eq[:eq_lambda]),first(endo[:z_t])]   = (xswts./c)'*dF2_dRZ
    JJ[first(eq[:eq_lambda]),first(endo[:w_t])]   = -(xswts./c)'*dF2_dWH
    JJ[first(eq[:eq_lambda]),first(endo[:L_t])]   = -(xswts./c)'*dF2_dWH
    JJ[first(eq[:eq_lambda]),first(endo[:t_t])]   = -(xswts./c)'*dF2_dTT
    JJ[first(eq[:eq_lambda]),endo[:l_t]] = -(xswts.*unc.*μ./c)'

    # transfer
    JJ[first(eq[:eq_transfers]),first(endo[:t_t])]  = T
    JJ[first(eq[:eq_transfers]),first(endo[:capreturn_t])]  = -Rk*kstar
    JJ[first(eq[:eq_transfers]),first(endo[:k_t])]  = -Rk*kstar
    JJ[first(eq[:eq_transfers]),first(endo[:z_t])]  = Rk*kstar
    JJ[first(eq[:eq_transfers]),first(endo[:I_t])]   = xstar
    JJ[first(eq[:eq_transfers]),first(endo[:mc_t])]  = ystar
    #JJ[first(eq[:eq_transfers]),first(endo[:y_t])]   = (1-1/g)*ystar
    #JJ[first(eq[:eq_transfers]),first(endo[:g_t])]   = (ystar/g)
    JJ[first(eq[:eq_transfers]),first(endo[:tg_t])]   = Tg

    # investment
    JJ[first(eq[:eq_investment]),first(endo[:Q_t])]  = 1.
    JJ[first(eq[:eq_investment]),first(endo[:μ_t])] = 1.
    JJ[first(eq[:eq_investment]),first(endo[:I′_t])] = spp*(exp(3*γ))/R
    JJ[first(eq[:eq_investment]),first(endo[:z′_t])] = spp*(exp(3*γ))/R
    JJ[first(eq[:eq_investment]),first(endo[:I_t])]  = -spp*(exp(3*γ))/R - spp*exp(2*γ)
    JJ[first(eq[:eq_investment]),first(endo[:I_t1])] = spp*exp(2*γ)
    JJ[first(eq[:eq_investment]),first(endo[:z_t])]  = -spp*exp(2*γ)

    # tobin's q
    JJ[first(eq[:eq_tobin_q]),first(endo[:margutil_t])]  = 1.
    JJ[first(eq[:eq_tobin_q]),first(endo[:margutil′_t])] = -1.
    JJ[first(eq[:eq_tobin_q]),first(endo[:Q_t])]    = 1.
    JJ[first(eq[:eq_tobin_q]),first(endo[:z′_t])]   = 1.
    JJ[first(eq[:eq_tobin_q]),first(endo[:capreturn′_t])]  = -Rk/R
    JJ[first(eq[:eq_tobin_q]),first(endo[:Q′_t])]   = -(1-δ)/R

    # capital accumulation
    JJ[first(eq[:eq_capital_accumulation]),first(endo[:k′_t])] = 1.
    JJ[first(eq[:eq_capital_accumulation]),first(endo[:k_t])]  = -(1-δ)
    JJ[first(eq[:eq_capital_accumulation]),first(endo[:z_t])]   = (1-δ)
    JJ[first(eq[:eq_capital_accumulation]),first(endo[:μ_t])]  = -xstar/kstar
    JJ[first(eq[:eq_capital_accumulation]),first(endo[:I_t])]   = -xstar/kstar

    # wage phillips curve
    JJ[first(eq[:eq_wage_phillips]),first(endo[:π_w_t])]  = -1.
    JJ[first(eq[:eq_wage_phillips]),first(endo[:λ_w_t])] = 1. #(ϕ*H^ϕh)/Φw
    JJ[first(eq[:eq_wage_phillips]),first(endo[:L_t])]   = κ_w*ϕh #(ϕ*(H^ϕh)*(1+lamw)/lamw*Φw)*ϕh
    JJ[first(eq[:eq_wage_phillips]),first(endo[:margutil_t])]  = -κ_w #-(ϕ*(H^ϕh)*(1+lamw)/lamw*Φw)
    JJ[first(eq[:eq_wage_phillips]),first(endo[:w_t])]    = -κ_w #-(ϕ*(H^ϕh)*(1+lamw)/lamw*Φw)
    JJ[first(eq[:eq_wage_phillips]),first(endo[:π_w′_t])]  = β

    # price phillips curve
    JJ[first(eq[:eq_price_phillips]),first(endo[:π_t])]   = -1.
    JJ[first(eq[:eq_price_phillips]),first(endo[:mc_t])]   = κ_p #(1+lamf)/(lamf*Φp)
    JJ[first(eq[:eq_price_phillips]),first(endo[:λ_f_t])] = 1. #1/Φp
    JJ[first(eq[:eq_price_phillips]),first(endo[:π′_t])]  = 1/R

    # marginal cost
    JJ[first(eq[:eq_marginal_cost]),first(endo[:mc_t])] = 1.
    JJ[first(eq[:eq_marginal_cost]),first(endo[:w_t])]  = -(1-α)
    JJ[first(eq[:eq_marginal_cost]),first(endo[:capreturn_t])] = -α

    # gdp
    JJ[first(eq[:eq_gdp]),first(endo[:y_t])]  = 1.
    JJ[first(eq[:eq_gdp]),first(endo[:z_t])]  = α
    JJ[first(eq[:eq_gdp]),first(endo[:k_t])] = -α
    JJ[first(eq[:eq_gdp]),first(endo[:L_t])] = -(1-α)

    # optimal k/l ratio
    JJ[first(eq[:eq_optimal_kl]),first(endo[:capreturn_t])] = 1.
    JJ[first(eq[:eq_optimal_kl]),first(endo[:w_t])]  = -1.
    JJ[first(eq[:eq_optimal_kl]),first(endo[:L_t])] = -1.
    JJ[first(eq[:eq_optimal_kl]),first(endo[:k_t])] = 1.
    JJ[first(eq[:eq_optimal_kl]),first(endo[:z_t])]  = -1.

    # taylor rule
    JJ[first(eq[:eq_taylor]),first(endo[:i_t])]   = -1.
    JJ[first(eq[:eq_taylor]),first(endo[:i_t1])]  = ρ_R
    JJ[first(eq[:eq_taylor]),first(endo[:π_t])]   = (1-ρ_R)*ψπ
    JJ[first(eq[:eq_taylor]),first(endo[:y_t])]    = (1-ρ_R)*ψy
    JJ[first(eq[:eq_taylor]),first(endo[:y_t1])]  = -(1-ρ_R)*ψy
    JJ[first(eq[:eq_taylor]),first(endo[:z_t])]    = (1-ρ_R)*ψy
    JJ[first(eq[:eq_taylor]),first(endo[:rm_t])] = 1.

    # fisher eqn
    JJ[first(eq[:eq_fisher]),first(endo[:R_t])]  = 1.
    JJ[first(eq[:eq_fisher]),first(endo[:π′_t])] = 1.
    JJ[first(eq[:eq_fisher]),first(endo[:i_t])]  = -1.

    # wage inflation
    JJ[first(eq[:eq_nominal_wage_inflation]),first(endo[:π_w_t])] = 1.
    JJ[first(eq[:eq_nominal_wage_inflation]),first(endo[:π_t])]  = -1.
    JJ[first(eq[:eq_nominal_wage_inflation]),first(endo[:z_t])]   = -1.
    JJ[first(eq[:eq_nominal_wage_inflation]),first(endo[:w_t])]   = -1.
    JJ[first(eq[:eq_nominal_wage_inflation]),first(endo[:w_t1])]  = 1.

    # fiscal rule
    JJ[first(eq[:eq_fiscal_rule]),first(endo[:tg_t])]  = -Tg
    JJ[first(eq[:eq_fiscal_rule]),first(endo[:R_t])]   = δb*bg/R
    JJ[first(eq[:eq_fiscal_rule]),first(endo[:bg_t])]  = δb*bg*exp(-γ)
    JJ[first(eq[:eq_fiscal_rule]),first(endo[:z_t])]   = -δb*bg*exp(-γ)
    JJ[first(eq[:eq_fiscal_rule]),first(endo[:y_t])]   = δb*(1-(1/g))*ystar
    JJ[first(eq[:eq_fiscal_rule]),first(endo[:g_t])]   = δb*ystar/g

    # govt budget constraint
    JJ[first(eq[:eq_g_budget_constraint]),first(endo[:bg′_t])] = -(bg/R)
    JJ[first(eq[:eq_g_budget_constraint]),first(endo[:R_t])]   = bg/R
    JJ[first(eq[:eq_g_budget_constraint]),first(endo[:bg_t])]  = bg*exp(-γ)
    JJ[first(eq[:eq_g_budget_constraint]),first(endo[:z_t])]   = -bg*exp(-γ)
    JJ[first(eq[:eq_g_budget_constraint]),first(endo[:y_t])]   = (1-(1/g))*ystar
    JJ[first(eq[:eq_g_budget_constraint]),first(endo[:g_t])]   = ystar/g
    JJ[first(eq[:eq_g_budget_constraint]),first(endo[:tg_t])]  = -Tg

    # resource constraints
    JJ[first(eq[:eq_resource_constraint]),first(endo[:y_t])] = ystar/g
    JJ[first(eq[:eq_resource_constraint]),first(endo[:g_t])] = -ystar/g
    JJ[first(eq[:eq_resource_constraint]),first(endo[:I_t])] = -xstar
    JJ[first(eq[:eq_resource_constraint]),first(endo[:C_t])] = -aggc

    JJ[first(eq[:LI]),first(endo[:i′_t1])] = 1.
    JJ[first(eq[:LI]),first(endo[:i_t])]   = -1.

    JJ[first(eq[:LY]),first(endo[:y′_t1])] = 1.
    JJ[first(eq[:LY]),first(endo[:y_t])]   = -1.

    JJ[first(eq[:LW]),first(endo[:w′_t1])] = 1.
    JJ[first(eq[:LW]),first(endo[:w_t])]   = -1.

    JJ[first(eq[:LX]),first(endo[:I′_t1])] = 1.
    JJ[first(eq[:LX]),first(endo[:I_t])]   = -1.

    # discount factor shock
    JJ[first(eq[:eq_b]),first(endo[:b′_t])] = 1.
    JJ[first(eq[:eq_b]),first(endo[:b_t])]  = -ρ_B

    # g/y shock
    JJ[first(eq[:eq_g]),first(endo[:g′_t])] = 1.
    JJ[first(eq[:eq_g]),first(endo[:g_t])]  = -ρ_G

    # tfp growth shock
    JJ[first(eq[:eq_z]),first(endo[:z′_t])] = 1.
    JJ[first(eq[:eq_z]),first(endo[:z_t])]  = -ρ_z

    # investment shock
    JJ[first(eq[:eq_μ]),first(endo[:μ′_t])] = 1.
    JJ[first(eq[:eq_μ]),first(endo[:μ_t])]  = -ρ_μ

    # wage mkup shock
    JJ[first(eq[:eq_λ_w]),first(endo[:λ_w′_t])] = 1.
    JJ[first(eq[:eq_λ_w]),first(endo[:λ_w_t])]  = -ρ_lamw

    # price mkup shock
    JJ[first(eq[:eq_λ_f]),first(endo[:λ_f′_t])] = 1.
    JJ[first(eq[:eq_λ_f]),first(endo[:λ_f_t])]  = -ρ_lamf

    # monetary policy shock
    JJ[first(eq[:eq_rm]),first(endo[:rm′_t])] = 1.
    JJ[first(eq[:eq_rm]),first(endo[:rm_t])]  = -ρ_mon

    #consumption
    #=JJ[first(eq[:eq_consumption]), endo[:l_t]] = (μ .*unc.*xswts.*c)'
    JJ[first(eq[:eq_consumption]), endo[:kf_t]] = -(xswts.*c)' # note, now we linearize
    JJ[first(eq[:eq_consumption]),first(endo[:R_t])]   = -(xswts.*c)'*dF2_dRZ
    JJ[first(eq[:eq_consumption]),first(endo[:z_t])]   = (xswts.*c)'*dF2_dRZ
    JJ[first(eq[:eq_consumption]),first(endo[:w_t])]   = -(xswts.*c)'*dF2_dWH
    JJ[first(eq[:eq_consumption]),first(endo[:L_t])]   = -(xswts.*c)'*dF2_dWH
    JJ[first(eq[:eq_consumption]),first(endo[:t_t])]   = -(xswts.*c)'*dF2_dTT=#

    if !m.testing && get_setting(m, :normalize_distr_variables)
        JJ  = normalize(m, JJ)
    end
    return JJ
end

function euler_equation_hetdsgegovdebt(nx::Int, ns::Int,
                                qp::Function, qfunction::Function,
                                xgrid::Vector{Float64}, sgrid::Vector{Float64},
                                fgrid::Matrix{Float64},
                                unc::BitArray,
                                xswts::Vector{Float64},
                                R::Float64, γ::Float64, β::Float64,
                                η::Float64, ell::Vector{Float64}, T::Float64,
                                ω::Float64, H::Float64)
    nxns = nx*ns
    ee  = zeros(nxns,nxns) # ee[i,j] takes you from i to j
    ξ   = zeros(nxns,nxns)
    Ξ   = zeros(nxns,nxns)
    dF1_dELL = zeros(nxns,nxns)
    dF1_dRZ  = zeros(nxns)
    dF1_dELLP = zeros(nxns,nxns)
    dF1_dWHP = zeros(nxns)
    dF1_dTTP = zeros(nxns)
    for iss=1:ns
        for ia=1:nx
            i  = nx*(iss-1)+ia
            sumELL = 0.
            sumRZ  = 0.
            sumWH  = 0.
            sumTT  = 0.
            for isp=1:ns
                for iap=1:nx
                    ip = nx*(isp-1)+iap
                    ee[i,ip] = (xgrid[iap] - R*(exp(-γ))*max(xgrid[ia]-1/ell[i], -η) - T)/(ω*H*sgrid[isp])
                    ξ[i,ip] = ((β*R*xswts[i]*exp(-γ))/(ω*H*sgrid[isp])^2)*max(ell[ip],1/(xgrid[iap]+η))*qp(ee[i,ip])*fgrid[iss,isp]
                    Ξ[i,ip] = ((β*R*xswts[i]*exp(-γ))/(ω*H*sgrid[isp]))*max(ell[ip],1/(xgrid[iap]+η))*qfunction(ee[i,ip])*fgrid[iss,isp]
                    sumELL += ξ[i,ip]*R*(exp(-γ))*unc[i]/ell[i]
                    sumRZ  += ξ[i,ip]*R*(exp(-γ))*max(xgrid[ia] - 1/ell[i],-η)
                    dF1_dELLP[i,ip] = Ξ[i,ip]*unc[ip]
                    sumWH  += Ξ[i,ip] + ξ[i,ip]*ee[i,ip]*(ω*H*sgrid[isp])
                    sumTT  += ξ[i,ip]*T
                end
            end
            dF1_dELL[i,i] = -ell[i] - sumELL
            dF1_dRZ[i] = ell[i] - sumRZ
            dF1_dWHP[i] = -sumWH
            dF1_dTTP[i] = -sumTT
        end
    end
    return dF1_dELL, dF1_dRZ, dF1_dELLP, dF1_dWHP, dF1_dTTP, ee
end

function kolmogorov_fwd_hetdsgegovdebt(nx::Int, ns::Int,
                                qfunction::Function, qp::Function,
                                xgrid::Vector{Float64}, sgrid::Vector{Float64},
                                fgrid::Matrix{Float64}, unc::BitArray,
                                xswts::Vector{Float64},
                                R::Float64, γ::Float64,
                                ell::Vector{Float64}, μ::Vector{Float64},
                                η::Float64, T::Float64, ω::Float64, H::Float64, ee::Matrix{Float64})
    nxns = nx*ns
    bigΨ    = zeros(nxns,nxns)
    smallψ = zeros(nxns,nxns)
    dF2_dRZ = zeros(nxns)
    dF2_dM = zeros(nxns,nxns)
    dF2_dELL = zeros(nxns,nxns)
    dF2_dWH = zeros(nxns)
    dF2_dTT = zeros(nxns)

    for isp=1:ns
        for iap=1:nx
            ip  = nx*(isp-1)+iap
            sumWH = 0.
            sumRZ = 0.
            sumTT = 0.
            for iss=1:ns
                for ia=1:nx
                    i = nx*(iss-1)+ia
                    bigΨ[ip,i] = xswts[i]*μ[i]*qfunction(ee[i,ip])*fgrid[iss,isp]/(ω*H*sgrid[isp])
                    smallψ[ip,i] = xswts[i]*μ[i]*qp(ee[i,ip])*fgrid[iss,isp]/((ω*H*sgrid[isp])^2)
                    sumWH += bigΨ[ip,i] + smallψ[ip,i]*ee[i,ip]*(ω*H*sgrid[isp])
                    sumRZ += smallψ[ip,i]*(R*exp(-γ))*max(xgrid[ia] - 1/ell[i],-η)
                    dF2_dELL[ip,i] = smallψ[ip,i]*(R*exp(-γ))*(unc[i]/ell[i])
                    sumTT += smallψ[ip,i]*T
                    dF2_dM[ip,i] = xswts[i]*qfunction(ee[i,ip])*fgrid[iss,isp]/(ω*H*sgrid[isp]) # note, now we linearize
                end
            end
            dF2_dWH[ip] = -sumWH
            dF2_dRZ[ip] = -sumRZ
            dF2_dTT[ip] = -sumTT
        end
    end
    return dF2_dWH, dF2_dRZ, dF2_dTT, dF2_dELL, bigΨ, dF2_dM
end

function normalize(m::HetDSGEGovDebt, JJ::Matrix{Float64})

    Qx, Qy, Qleft, Qright  = compose_normalization_matrices(m)

    m <= Setting(:n_predetermined_variables, size(Qx, 1))
    m <= Setting(:Qx, Qx)
    m <= Setting(:Qy, Qy)
    #Qleft = get_setting(m, :Qleft) #m[:Qleft].value
    #Qright = get_setting(m, :Qright) #m[:Qright].value

	Jac1 = Qleft*JJ*Qright

    return Jac1
end

# use proxy distribuitions method - now `averaging' so that Qx'*Qx should not `add mass'
function avg_prox(n::Int, binsize::Int = 3, minn_noag::Int = 0)
    minn_noag_use = min(minn_noag,n)
    n_gps = Int(floor((n-minn_noag_use)/binsize)) # number of groups
    n_noag = n - binsize*n_gps
    Prox = cat(eye(n_noag), kron(eye(n_gps),ones(1,binsize)./sqrt(binsize)), dims = [1 2])
    n_new = n_noag + n_gps
    return (Prox, n_new, n_noag, n_gps)
end
# note: setting minn_noag = 0 seems to work the best

function make_S(n::Int, n_noag::Int = 0, binsize::Int = 1)
    P1 = ones(n,1)*sqrt(binsize)
    P1[1:n_noag] .= 1.
    Ptemp = eye(n)
    Ptemp = Ptemp[:,2:end]
    P2 = Ptemp
    P = [P1 P2]
    (QQQ,Rjunk)=qr(P)
    S         = QQQ[:,2:end]'
    return S
end

function compose_normalization_matrices(m::HetDSGEGovDebt)
    if get_setting(m, :poor_man_reduc)
        nx = get_setting(m, :nx1_state) #:nx)
        ns = get_setting(m, :ns)
        n = nx*ns #get_setting(m, :nx1) + get_setting(m, :nx2) #
        nscalars = get_setting(m, :nscalars)
        nyscalars = get_setting(m, :nyscalars)
        nxscalars = get_setting(m, :nxscalars)
        mindens = get_setting(m, :mindens)
        μ = m[:μstar].value

        nx1 = maximum(findall(μ[1:nx].>mindens))
        m <= Setting(:nx1_state, nx1)
        m <= Setting(:nx1_jump, nx1)

        setup_indices!(m)
        init_states_and_jumps!(m, get_setting(m, :states), get_setting(m, :jumps))
        normalize_model_state_indices!(m)
        endogenous_states_augmented = [:C_t1]
        for (i,k) in enumerate(endogenous_states_augmented); m.endogenous_states_augmented[k] = i + first(m.endogenous_states[get_setting(m, :jumps)[end]]) end #first(collect(values(m.endogenous_states))[end]) end

        m <= Setting(:n_model_states_augmented, get_setting(m, :n_model_states) +
                     length(m.endogenous_states_augmented))

        Reduc = eye(n)
        Reduc=Reduc[[1:nx1;nx+1:n],:]

        minn_noag = 0
        binsize = get_setting(m, :binsize) # this is the maximum binsize which seemed to leave the IRFs unchanged, you can experiment with this
        # setting binsize = 1 should return what we had before
        (ProxL, n_newL, n_noagL, n_gpsL) = avg_prox(nx1, binsize, minn_noag)
        (ProxH, n_newH, n_noagH, n_gpsH) = avg_prox(nx, binsize,minn_noag)
        m <= Setting(:nx1_state, n_newL)
        m <= Setting(:nx2_state, n_newH)
        if get_setting(m, :reduce_ell)
            m <= Setting(:nx1_jump, n_newL)
            m <= Setting(:nx2_jump, n_newH)
        end
        setup_indices!(m)
        init_states_and_jumps!(m, get_setting(m, :states), get_setting(m, :jumps))
        normalize_model_state_indices!(m)
        endogenous_states_augmented = [:C_t1]
        for (i,k) in enumerate(endogenous_states_augmented); m.endogenous_states_augmented[k] = i + first(m.endogenous_states[get_setting(m, :jumps)[end]]) end #first(collect(values(m.endogenous_states))[end]) end

        m <= Setting(:n_model_states_augmented, get_setting(m, :n_model_states) +
                     length(m.endogenous_states_augmented))

        Prox = cat(ProxL, ProxH, dims = [1 2])
        SL = make_S(n_newL, n_noagL, binsize)
        SH = make_S(n_newH, n_noagH, binsize)
        S = cat(SL, SH, dims = [1 2])

        if get_setting(m, :reduce_ell)
            Qleft     = cat(Prox*Reduc,S*Prox*Reduc,eye(nscalars), dims = [1 2])
            Qx        = cat(S*Prox*Reduc,eye(nxscalars), dims = [1 2])
            Qy        = cat(Prox*Reduc,eye(nyscalars), dims = [1 2])
        else
            Qleft     = cat(Reduc,S*Prox*Reduc,eye(nscalars), dims = [1 2])
            Qx        = cat(S*Prox*Reduc,eye(nxscalars), dims = [1 2])
            Qy        = cat(Reduc,eye(nyscalars), dims = [1 2])
        end

        Qright    = cat(Qx',Qy',Qx',Qy', dims = [1,2])

        return Qx, Qy, Qleft, Qright

    else
        nx = get_setting(m, :nx)
        ns = get_setting(m, :ns)
        nscalars = get_setting(m, :nscalars)
        nyscalars = get_setting(m, :nyscalars)
        nxscalars = get_setting(m, :nxscalars)

        # Create PPP matrix
        P1 = kron(Matrix{Float64}(I, ns,ns),ones(nx,1))
        Ptemp = Matrix{Float64}(I, nx, nx)
        Ptemp = Ptemp[:, 2:end]
        P2 = kron(Matrix{Float64}(I, ns, ns), Ptemp)
        P  = hcat(P1, P2)

        Q,R = qr(P)
        Q = Array(Q)
        S         = Q[:, ns+1:end]'

        nxns = nx*ns

        Qleft     = cat(Matrix{Float64}(I, nxns, nxns),S,Matrix{Float64}(I, nscalars, nscalars), dims = [1 2])
        Qx        = cat(S,Matrix{Float64}(I, nxscalars, nxscalars), dims = [1 2])
        Qy        = cat(Matrix{Float64}(I, nxns, nxns),Matrix{Float64}(I, nyscalars, nyscalars), dims = [1 2])
        Qright    = cat(Qx',Qy',Qx',Qy', dims = [1,2])

        return Qx, Qy, Qleft, Qright

    end
end

function truncate_distribution!(m::HetDSGEGovDebt)
    mindens = get_setting(m, :mindens)
    trunc_distr = get_setting(m, :trunc_distr)
    rescale_weights = get_setting(m, :rescale_weights)

    nx = get_setting(m, :nx)
    μ = m[:μstar].value
    ell = m[:lstar].value
    c = m[:cstar].value
    xgrid = m.grids[:xgrid].points
    #xlo = get_setting(m, :xlo)
    swts::Vector{Float64}  = m.grids[:sgrid].weights

    if trunc_distr
        oldnx = nx
        nx = maximum(findall(μ[1:nx]+μ[nx+1:2*nx] .> mindens)) # used to be 1e-8
        m[:μstar] = μ[[1:nx;oldnx+1:oldnx+nx]]
        m[:lstar] = ell[[1:nx;oldnx+1:oldnx+nx]]
        m[:cstar] = c[[1:nx;oldnx+1:oldnx+nx]]
        if rescale_weights
            xhi = xgrid[nx]
            xlo = xgrid[1]
            xscale = xhi-xlo
        end
        m <= Setting(:nx, nx)
        m <= Setting(:xhi, xhi)
        m <= Setting(:xscale, xscale)
        m.grids[:xgrid] = Grid(uniform_quadrature(xscale), xlo, xhi, nx, scale = xscale)
        m.grids[:weights_total] = kron(swts, m.grids[:xgrid].weights)
        nxns = nx*get_setting(m, :ns)
        m <= Setting(:n, nxns)
        m <= Setting(:nx, nx)

        setup_indices!(m)

        init_states_and_jumps!(m, get_setting(m, :states), get_setting(m, :jumps))

        normalize_model_state_indices!(m)

        endogenous_states_augmented = [:C_t1]
        for (i,k) in enumerate(endogenous_states_augmented); m.endogenous_states_augmented[k] = i + first(m.endogenous_states[get_setting(m, :jumps)[end]]) end #first(collect(values(m.endogenous_states))[end]) end

        m <= Setting(:n_model_states_augmented, get_setting(m, :n_model_states) +
                     length(m.endogenous_states_augmented))

    end
end


text:
# Methods for using Riesz analysis for building the demonstration warp map.

function higherordermonogenicanalysis(y::Array{T,D},
                                ω_set::Vector{T},
                                pass_band_factor::T,
                                scheme_cmd::String,
                                order::Int) where {T,D}

    N_bands = length(ω_set)

    #### Pre-filter.
    attenuation_factor_at_cut_off = 2 # factor of 2 attenuation at cut-off.
    reciprocal_cut_off_percentage = 1.0 # denominator is 0 to 1.
    σ_c = π/(reciprocal_cut_off_percentage*sqrt(2*log(attenuation_factor_at_cut_off)))

    LP, HP = getGaussianfilters(size(y),σ_c) # 3dB attenuation is approx. 2.355 std dev. for a Gaussian.

    Y = real.(ifft(fft(y).*LP)) # bandlimited version of y.
    residual = real.(ifft(fft(y).*HP))

    #### Split-band analysis.
    getfilterpairfunc = getGaussianfilters

    ϕY, ψY = runsplitbandanalysis(Y, ω_set, getfilterpairfunc)
    ηY = runbandpassanalysis(Y, ω_set, pass_band_factor, getfilterpairfunc)

    #### Riesz transform on the different filtered signals.
    H, ordering = gethigherorderRTfilters(Y,order)

    𝓡ϕY = collect( RieszAnalysisLimited(ϕY[s],H) for s = 1:N_bands)
    𝓡ψY = collect( RieszAnalysisLimited(ψY[s],H) for s = 1:N_bands)
    𝓡ηY = collect( RieszAnalysisLimited(ηY[s],H) for s = 1:N_bands)

    #### Generalized monogenic analysis.
    Hᵤ = Array{Array{Float64,D}}(N_bands)
    Aᵤ = Array{Array{Float64,D}}(N_bands)
    ϕᵤ = Array{Array{Float64,D}}(N_bands)

    # Metric tensor.
    metric_tensor_vec = collect( Combinatorics.multinomial(ordering[m]...) for m = 1:length(ordering) )

    # Perform analysis.
    for ss = 1:N_bands
        Hᵤ[ss] = directionalHilbert(metric_tensor_vec, 𝓡ηY[ss])
        Aᵤ[ss], ϕᵤ[ss] = monogenicanalysis(Hᵤ[ss], ηY[ss])
    end

    #### Compute ϕ.
    ϕ_set = collect( Hᵤ[s].*cos(ϕᵤ[s]) for s = 1:N_bands )

    return ϕ_set, 𝓡ϕY, 𝓡ψY, 𝓡ηY, H, ordering, Aᵤ, ϕᵤ, Hᵤ, Y, y, ϕY, ψY, ηY
end

function getRieszwarpmapsamples( y::Array{T,D},
                            ::Val{:simple},
                            ::Val{:uniform},
                            ω_set,
                            pass_band_factor) where {T,D}
    #
    N_bands = length(ω_set)

    Y = y

    #### Split-band analysis.
    ϕY, ψY = SignalTools.runsplitbandanalysis(Y, ω_set, SignalTools.getGaussianfilters)
    ηY = SignalTools.runbandpassanalysis(Y, ω_set, pass_band_factor, SignalTools.getGaussianfilters)

    # #### Riesz transform on the different filtered signals.
    # H, ordering = gethigherorderRTfilters(Y,order)
    #
    # 𝓡ϕY = collect( RieszAnalysisLimited(ϕY[s],H) for s = 1:N_bands)
    # 𝓡ψY = collect( RieszAnalysisLimited(ψY[s],H) for s = 1:N_bands)
    # 𝓡ηY = collect( RieszAnalysisLimited(ηY[s],H) for s = 1:N_bands)

    ϕ_set = ηY

    ϕ = reduce(+,ϕ_set)./N_bands

    return ϕ
    #return y #ηY[4]
end

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
using Roadways
import LibSpatialIndex
using NearestNeighbors
import PyPlot

filename = joinpath(Pkg.dir("Roadways"),"data","LA-UniversalStudios.hdf5")
open_map = load_map_from_HDF5(filename)
preprocess!(open_map)
open_map

seg1 = open_map.road_segments[3]
seg2 = open_map.road_segments[5]
r = 20
pt1 = rand(open_map.cached_curves[seg1.refline].pts)
pos1 = CartesianPos(
    x = pt1.x + r*(rand()),
    y = pt1.y + r*(rand()),
    θ = pt1.θ + r*(rand()),
)
pt2 = rand(open_map.cached_curves[seg2.refline].pts)
pos2 = CartesianPos(
    x = pt2.x + r*(rand()),
    y = pt2.y + r*(rand()),
    θ = pt2.θ + r*(rand()),
)

PyPlot.figure(figsize=[18,18])
ref_colors = ["blue","red","green","purple","cyan"]
for (i, segment_id) in enumerate(collect(keys(open_map.road_segments)))
    plot_segment(open_map, segment_id; ref_color=ref_colors[i],bound_color="white")
end
PyPlot.scatter(pos1.x,pos1.y,c="black")
PyPlot.scatter(pos2.x,pos2.y,c="black")
PyPlot.axis("equal")
PyPlot.legend()
PyPlot.show()

map_idx1 = project_to_map(open_map, pos1)[1]
map_idx2 = project_to_map(open_map, pos2)[1]
map_idx1, map_idx2

map_idx1.refline_index, map_idx2.refline_index

compute_longitudinal_offset_cached(open_map, map_idx1, map_idx2),
compute_lane_offset_cached(open_map, map_idx1, map_idx2)

seg_id, segment = rand(open_map.road_segments)
curve = open_map.curves[segment.refline]
cached_curve = open_map.cached_curves[segment.refline]
kdtree = open_map.spatial_index.kdtrees[segment.refline]
curve_pt = open_map.curve_points[curve.pts[Int(round(length(curve.pts)/2))]]

r = 10
pos = CartesianPos(
    x = curve_pt.x + r*(.5-rand()),
    y = curve_pt.y + r*(.5-rand()),
    θ = mod(curve_pt.θ + r*(.5-rand()),2π)
)

curve_proj = project_to_curve(cached_curve, pos);
proj_pt = curve_pt_from_projection(cached_curve, curve_proj)

lane_ids = get_lane_associations(open_map, seg_id, proj_pt)

PyPlot.figure(figsize=[12,12])
plot_segment(open_map, seg_id)
for lane_id in lane_ids
    curve = open_map.cached_curves[open_map.lane_segments[lane_id].centerline]
    PyPlot.plot([pt.x for pt in curve.pts],[pt.y for pt in curve.pts],c="lime")
end
# PyPlot.scatter(foot_pt.x, foot_pt.y,marker="o",c="blue",edgecolor="none",label="foot_pt")
PyPlot.scatter(proj_pt.x, proj_pt.y,marker="o",c="red",edgecolor="none",label="proj_pt")
PyPlot.legend()
PyPlot.axis("equal")

map_idxs = project_to_map(open_map, pos)

project_to_segment(open_map, seg_id, pos)

get_lane_associations(open_map, seg_id, proj_pt)

MapIndex()


text:
@testset "Constructors" begin
    using Distributions, BayesNN

    N = 1000  # number of observations
    D = 3     # number of covariates including intercept

    # generate data
    srand(1)
    X = rand(Normal(), N, D)
    y = rand(Bernoulli(0.5), N)

    @test_nowarn net = Nnet(y, X)
end


text:
export Polynomial

"""
    Polynomial{T<:Number, X}(coeffs::AbstractVector{T}, [var = :x])

Construct a polynomial from its coefficients `coeffs`, lowest order first, optionally in
terms of the given variable `var` which may be a character, symbol, or a string.

If ``p = a_n x^n + \\ldots + a_2 x^2 + a_1 x + a_0``, we construct this through
`Polynomial([a_0, a_1, ..., a_n])`.

The usual arithmetic operators are overloaded to work with polynomials as well as
with combinations of polynomials and scalars. However, operations involving two
polynomials of different variables causes an error except those involving a constant polynomial.

!!! note
    `Polynomial` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first
    index always corresponding to the constant term. In order to use the axis of `coeffs` as exponents,
    consider using a [`LaurentPolynomial`](@ref) or possibly a [`SparsePolynomial`](@ref).

# Examples
```jldoctest
julia> using Polynomials

julia> Polynomial([1, 0, 3, 4])
Polynomial(1 + 3*x^2 + 4*x^3)

julia> Polynomial([1, 2, 3], :s)
Polynomial(1 + 2*s + 3*s^2)

julia> one(Polynomial)
Polynomial(1.0)
```
"""
struct Polynomial{T <: Number, X} <: StandardBasisPolynomial{T, X}
    coeffs::Vector{T}
    function Polynomial{T, X}(coeffs::AbstractVector{S}) where {T <: Number, X, S}
        if Base.has_offset_axes(coeffs)
            @warn "ignoring the axis offset of the coefficient vector"
        end
        N = findlast(!iszero, coeffs)
        N == nothing && return new{T,X}(zeros(T,1))
        cs = T[coeffs[i] for i ∈ firstindex(coeffs):N]
        new{T,X}(cs)
    end
    # non-copying alternative assuming !iszero(coeffs[end])
    function Polynomial{T, X}(checked::Val{false}, coeffs::AbstractVector{T}) where {T <: Number, X}
        new{T, X}(coeffs)
    end
end

@register Polynomial

# scalar +,* faster  than standard-basis/common versions as it avoids a copy
function Base.:+(p::P, c::S) where {T, X, P <: Polynomial{T, X}, S<:Number}
    R = promote_type(T, S)
    Q = Polynomial{R,X}
    as = convert(Vector{R}, copy(coeffs(p)))
    as[1] += c
    iszero(as[end]) ? Q(as) : Q(Val(false), as)
end

function Base.:*(p::P, c::S) where {T, X, P <: Polynomial{T,X} , S <: Number}
    R = promote_type(T,S)
    Q = Polynomial{R, X}
    as = R[aᵢ * c for aᵢ ∈ coeffs(p)]
    iszero(as[end]) ? Q(as) : Q(Val(false), as)
end

# implement, as not copying speeds up multiplication by a factor of 2 or so
# over the default
function Base.:+(p1::P1, p2::P2) where {T,X, P1<:Polynomial{T,X},
                                        S,   P2<:Polynomial{S,X}}
    n1, n2 = length(p1), length(p2)
    R = promote_type(T,S)
    Q = Polynomial{R,X}
    if n1 == n2
        cs = ⊕(P1, p1.coeffs, p2.coeffs)
        return iszero(cs[end]) ? Q(cs) : Q(Val(false), cs)
    elseif n1 > n2
        cs = ⊕(P1, p1.coeffs, p2.coeffs)
    else
        cs = ⊕(P1, p2.coeffs, p1.coeffs)
    end

    Q(Val(false), cs)
end

function Base.:*(p::P, q::P) where {T,X, P<:Polynomial{T,X}}
    c = fastconv(p.coeffs, q.coeffs)
    return iszero(c[end]) ? P(c) : P(Val(false), c)
end


text:
"""
    savecheckpoint(dist, t)

Saves the distribution function to a binary file.

This file can be seen as a checkpoint or snapshot of the simulations internal state and be used to restart the simulation.

# Example
```jldoctest
julia> using JuSwalbe

julia> dist = JuSwalbe.DistributionD1Q3(f0 = ones(10), f1 = fill(0.1, 10), f2 = fill(-0.1, 10))
JuSwalbe.DistributionD1Q3{Array{Float64,1}}
  f0: Array{Float64}((10,)) [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
  f1: Array{Float64}((10,)) [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
  f2: Array{Float64}((10,)) [-0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1]

julia> savecheckpoint(dist)

julia> isfile("data/tmp/checkpoint_0.bson")
true

julia> load_dist = loadcheckpoint("data/tmp/checkpoint_0.bson")
JuSwalbe.DistributionD1Q3{Array{Float64,1}}
  f0: Array{Float64}((10,)) [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
  f1: Array{Float64}((10,)) [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
  f2: Array{Float64}((10,)) [-0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1]

julia> load_dist.f0 == dist.f0; load_dist.f1 == dist.f1; load_dist.f2 == dist.f2
true

julia> rm("data/tmp/checkpoint_0.bson")
```

See also: [`loadcheckpoint`](@ref)

# Additional information
The size of the checkpoint heavily depends on the system size.
In case of two spatial dimensions and a 4096^2 lattice a single checkpoint file for a Float64 simulation is about 1.2GB data.

"""
function savecheckpoint(dist::JuSwalbe.Distributionfunction; name::String="tmp", t::Int=0)
    # Generate a path for the checkpoint and a file name
    outpath = mkpath("data/$name")
    location = outpath * "/checkpoint_" * "$t" * ".bson"
    # BSON saves only dictonaries, so transform the distribution to a dictonary
    distdict = struct2dict(dist)
    # Remove all previous checkpoints
    rm.(glob("data/$name/checkpoint_*"))
    # Save the distribution and generate a checkpoint file
    bson(location, distdict)
end

"""
    savecheckpoints(dist, name, t)

Same as savecheckpoint but does not delete the previous checkpoint files.

See also [`savecheckpoint`](@ref), [`loadcheckpoint`](@ref)
"""
function savecheckpoints(dist::JuSwalbe.Distributionfunction; name::String="tmp", t::Int=0)
    # Generate a path for the checkpoint
    outpath = mkpath("data/$name")
    location = outpath * "/checkpoint_" * "$t" * ".bson"
    # BSON saves only dictonaries, so transform the distribution to a dictonary
    distdict = struct2dict(dist)

    # Save the distribution and generate a checkpoint file
    bson(location, distdict)
end

"""
    loadcheckpoint(file)

Loads a checkpoint file and generates a distribution function from it.

With the distribution function one can rerun a simulation from the point where the distribution function was saved.
This is helpful for the long runs, when it takes several hours to generate convergence.
Such it is adviced to generate a checkpoint roughly every three hours.

# Example
For an example see `savecheckpoint`.

See also: [`savecheckpoint`](@ref), [`savecheckpoints`](@ref)
"""
function loadcheckpoint(file)
    # Save the binary file to a dictonary
    distdict = BSON.load(file)
    # Check if it is one or two dimensional
    # and create a distribution
    if length(keys(distdict)) == 3
        dist = JuSwalbe.DistributionD1Q3(f0 = distdict[:f0], f1 = distdict[:f1], f2 = distdict[:f2])
    elseif length(keys(distdict)) == 9
        dist = JuSwalbe.DistributionD2Q9(f0 = distdict[:f0], f1 = distdict[:f1], f2 = distdict[:f2],
                                         f3 = distdict[:f3], f4 = distdict[:f4], f5 = distdict[:f5],
                                         f6 = distdict[:f6], f7 = distdict[:f7], f8 = distdict[:f8])
    end
    # Lastly return the distribution
    return dist
end

"""
    height2file(mom, name, time)

Saves the heightfield to a file in a binary format.

For data purposes and reproducibility one can save the heightfield to a file.
Such the file can be accessed with another software for e.g. post processing.

# Example
```jldoctest
julia> using JuSwalbe, BSON

julia> mom = JuSwalbe.Macroquant(height = ones(5,5), velocity=JuSwalbe.Twovector(x=zeros(5,5), y=zeros(5,5)), pressure=zeros(5,5), energy=zeros(5,5))
JuSwalbe.Macroquant{Array{Float64,2},JuSwalbe.Twovector{Array{Float64,2}}}
  height: Array{Float64}((5, 5)) [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]
  velocity: JuSwalbe.Twovector{Array{Float64,2}}
  pressure: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  energy: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

julia> height2file(mom)

julia> h = BSON.load("data/tmp/height_0.bson")
Dict{Symbol,Any} with 1 entry:
  :height => [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 …

julia> h[:height]
5×5 Array{Float64,2}:
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0

julia> mom.height == h[:height]
true

julia> rm("data/tmp/height_0.bson")

```

See also: [`savecheckpoint`](@ref), [`velocity2file`](@ref), [`velocityandheight2file`](@ref)
"""
function height2file(mom::Macroscopic_quantity; name::String="tmp", time::Int=0)
    outpath = mkpath("data/$name")
    location = outpath * "/height_" * "$time" * ".bson"
    # println("Saved checkpoint @ $location")
    bson(location, height = mom.height)
end

"""
    velocity2file(mom, name, time)

Saves the velocity field to a binary file called velocity_time.bson.

Independent of the dimensionality of the system, `D1Q3`, `D2Q9` saves the velocity to a file.
In case of `D2Q9` saves both velocity components to the file such `velocity.x` and `velocity.y`.

# Example
```jldoctest
julia> using JuSwalbe, BSON

julia> mom = simplemoment2d(5,5)
JuSwalbe.Macroquant{Array{Float64,2},JuSwalbe.Twovector{Array{Float64,2}}}
  height: Array{Float64}((5, 5)) [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]
  velocity: JuSwalbe.Twovector{Array{Float64,2}}
  pressure: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  energy: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

julia> velocity2file(mom)

julia> vel = BSON.load("data/tmp/velocity_0.bson")
Dict{Symbol,Any} with 2 entries:
  :velocity_x => [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; …
  :velocity_y => [-0.1 -0.1 … -0.1 -0.1; -0.1 -0.1 … -0.1 -0.1; … ; -0.1 -0.1 ……

julia> vel[:velocity_x]
5×5 Array{Float64,2}:
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1

julia> v = JuSwalbe.Twovector(x = vel[:velocity_x], y = vel[:velocity_y])
JuSwalbe.Twovector{Array{Float64,2}}
  x: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  y: Array{Float64}((5, 5)) [-0.1 -0.1 … -0.1 -0.1; -0.1 -0.1 … -0.1 -0.1; … ; -0.1 -0.1 … -0.1 -0.1; -0.1 -0.1 … -0.1 -0.1]

julia> rm("data/tmp/velocity_0.bson")
```

See also: [`height2file`](@ref), [`savecheckpoint`](@ref), [`velocityandheight2file`](@ref)
"""
function velocity2file(mom::Macroquant{Vector{T}, Vector{T}}; name::String="tmp", time::Int=0) where {T<:Number}
    outpath = mkpath("data/$name")
    location = outpath * "/velocity_" * "$time" * ".bson"

    bson(location, velocity = mom.velocity)

end
function velocity2file(mom::Macroquant{Matrix{T}, JuSwalbe.Twovector{Matrix{T}}}; name::String="tmp", time::Int=0) where {T<:Number}
    outpath = mkpath("data/$name")
    location = outpath * "/velocity_" * "$time" * ".bson"

    bson(location, velocity_x = mom.velocity.x, velocity_y = mom.velocity.y)

end

"""
    velocityandheight2file(mom, name, time)

Saves the velocity as well as the height field to a binary file called height_velocity_time.bson.

Independent of the dimensionality of the system, `D1Q3`, `D2Q9` saves the velocity and height to a file.
In case of `D2Q9` saves both velocity components to the file such `velocity.x` and `velocity.y`.

# Example
```jldoctest
julia> using JuSwalbe, BSON

julia> mom = simplemoment2d(5,5)
JuSwalbe.Macroquant{Array{Float64,2},JuSwalbe.Twovector{Array{Float64,2}}}
  height: Array{Float64}((5, 5)) [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]
  velocity: JuSwalbe.Twovector{Array{Float64,2}}
  pressure: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  energy: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

julia> velocityandheight2file(mom)

julia> h_vel = BSON.load("data/tmp/height_velocity_0.bson")
Dict{Symbol,Any} with 3 entries:
  :height     => [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; …
  :velocity_x => [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; …
  :velocity_y => [-0.1 -0.1 … -0.1 -0.1; -0.1 -0.1 … -0.1 -0.1; … ; -0.1 -0.1 ……

julia> h_vel[:height]
5×5 Array{Float64,2}:
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0

julia> h_vel[:velocity_x]
5×5 Array{Float64,2}:
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1

julia> rm("data/tmp/height_velocity_0.bson")

```

See also: [`height2file`](@ref), [`savecheckpoint`](@ref), [`velocityandheight2file`](@ref)
"""
function velocityandheight2file(mom::Macroquant{Vector{T}, Vector{T}}; name::String="tmp", time::Int=0) where {T<:Number}
    outpath = mkpath("data/$name")
    location = outpath * "/height_velocity_" * "$time" * ".bson"

    bson(location, height = mom.height, velocity = mom.velocity)

end
function velocityandheight2file(mom::Macroquant{Matrix{T}, JuSwalbe.Twovector{Matrix{T}}}; name::String="tmp", time::Int=0) where {T<:Number}
    outpath = mkpath("data/$name")
    location = outpath * "/height_velocity_" * "$time" * ".bson"

    bson(location, height = mom.height, velocity_x = mom.velocity.x, velocity_y = mom.velocity.y)

end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
export NEAT, HyperNEAT, GA_NEAT

"""
fitness::Function i::NEATIndiv -> computed_fitness(i)
"""

function core_NEAT(itype::Type, cfg::Dict, fitness::Function)
    evaluate!::Function = e::Evolution -> NEAT_evaluate!(e, fitness)
    selection::Function = i::Array{NEATIndiv} -> i[1]
    if cfg["selection_type"] == "tournament"
        selection =
            i::Array{NEATIndiv} -> NEAT_tournament(i, cfg["tournament_size"])
    elseif cfg["selection_type"] == "random_top"
        selection =
            i::Array{NEATIndiv} -> NEAT_random_top(i, cfg["survival_threshold"])
    else
        throw(ArgumentError("Wrong selection type: " + cfg["selection_type"]))
    end
    populate!::Function = e::Evolution -> NEAT_populate!(e, selection)
    Evolution(itype, cfg; evaluate = evaluate!, populate = populate!)
end

function NEAT(cfg::Dict, fitness::Function, fitness_args...)
    f::Function = i::NEATIndividual -> fitness(i, fitness_args...)
    core_NEAT(NEATIndividual, cfg, f)
end

function HyperNEAT(cfg::Dict, fitness::Function, fitness_args...)
    f::Function = i::HyperNEATIndividual -> fitness(i, fitness_args...)
    core_NEAT(HyperNEATIndividual, cfg, f)
end

## Genetic Algorithm with NEAT individuals

function GA_NEAT(itype::Type, cfg::Dict, fitness::Function; kwargs...)
    # Build before evaluation
    function evaluate(e::Evolution)
        build!.(e.population)
        Cambrian.fitness_evaluate!(e; fitness = fitness)
    end

    # Force NEAT mutation and crossover operators
    function neat_ga_populate!(e::Evolution)
        mut::Function = i::NEATIndiv -> NeuroEvolution.mutate(i, e.cfg)
        cross::Function = (p1, p2) -> NeuroEvolution.crossover(p1, p2, e.cfg)
        Cambrian.ga_populate!(e, mutation = mut, crossover = cross)
    end

    Evolution(itype, cfg; evaluate = evaluate, populate = neat_ga_populate!)
end


text: using Analytical, CSV, DataFrames


text:
function solve()
    n, x = [parse(Int, x) for x in split(readline())]
    s = readline()

    for c in s
        if c == 'o'
            x += 1
        elseif c == 'x' && x > 0
            x -= 1
        end
    end
    x
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct SimpleMarconatoMachine <: Machine
        R::Float64
        Xd::Float64
        Xq::Float64
        Xd_p::Float64
        Xq_p::Float64
        Xd_pp::Float64
        Xq_pp::Float64
        Td0_p::Float64
        Tq0_p::Float64
        Td0_pp::Float64
        Tq0_pp::Float64
        T_AA::Float64
        MVABase::Float64
        ext::Dict{String, Any}
        γd::Float64
        γq::Float64
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of 4-states synchronous machine: Simplified Marconato model
 The derivative of stator fluxes (ψd and ψq) is neglected and ωψd = ψd and
 ωψq = ψq is assumed (i.e. ω=1.0). This is standard when transmission network
 dynamics is neglected.

# Arguments
- `R::Float64`: Resistance after EMF in machine per unit, validation range: (0, nothing)
- `Xd::Float64`: Reactance after EMF in d-axis per unit, validation range: (0, nothing)
- `Xq::Float64`: Reactance after EMF in q-axis per unit, validation range: (0, nothing)
- `Xd_p::Float64`: Transient reactance after EMF in d-axis per unit, validation range: (0, nothing)
- `Xq_p::Float64`: Transient reactance after EMF in q-axis per unit, validation range: (0, nothing)
- `Xd_pp::Float64`: Sub-Transient reactance after EMF in d-axis per unit, validation range: (0, nothing)
- `Xq_pp::Float64`: Sub-Transient reactance after EMF in q-axis per unit, validation range: (0, nothing)
- `Td0_p::Float64`: Time constant of transient d-axis voltage, validation range: (0, nothing)
- `Tq0_p::Float64`: Time constant of transient q-axis voltage, validation range: (0, nothing)
- `Td0_pp::Float64`: Time constant of sub-transient d-axis voltage, validation range: (0, nothing)
- `Tq0_pp::Float64`: Time constant of sub-transient q-axis voltage, validation range: (0, nothing)
- `T_AA::Float64`: Time constant of d-axis additional leakage, validation range: (0, nothing)
- `MVABase::Float64`: Nominal Capacity in MVA, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `γd::Float64`
- `γq::Float64`
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct SimpleMarconatoMachine <: Machine
    "Resistance after EMF in machine per unit"
    R::Float64
    "Reactance after EMF in d-axis per unit"
    Xd::Float64
    "Reactance after EMF in q-axis per unit"
    Xq::Float64
    "Transient reactance after EMF in d-axis per unit"
    Xd_p::Float64
    "Transient reactance after EMF in q-axis per unit"
    Xq_p::Float64
    "Sub-Transient reactance after EMF in d-axis per unit"
    Xd_pp::Float64
    "Sub-Transient reactance after EMF in q-axis per unit"
    Xq_pp::Float64
    "Time constant of transient d-axis voltage"
    Td0_p::Float64
    "Time constant of transient q-axis voltage"
    Tq0_p::Float64
    "Time constant of sub-transient d-axis voltage"
    Td0_pp::Float64
    "Time constant of sub-transient q-axis voltage"
    Tq0_pp::Float64
    "Time constant of d-axis additional leakage"
    T_AA::Float64
    "Nominal Capacity in MVA"
    MVABase::Float64
    ext::Dict{String, Any}
    γd::Float64
    γq::Float64
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function SimpleMarconatoMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, T_AA, MVABase, ext=Dict{String, Any}(), )
    SimpleMarconatoMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, T_AA, MVABase, ext, ((Td0_pp*Xd_pp)/(Td0_p*Xd_p) )*(Xd-Xd_p), ((Tq0_pp*Xq_pp)/(Tq0_p*Xq_p) )*(Xq-Xq_p), [:eq_p, :ed_p, :eq_pp, :ed_pp], 4, InfrastructureSystemsInternal(), )
end

function SimpleMarconatoMachine(; R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, T_AA, MVABase, ext=Dict{String, Any}(), )
    SimpleMarconatoMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, T_AA, MVABase, ext, )
end

# Constructor for demo purposes; non-functional.
function SimpleMarconatoMachine(::Nothing)
    SimpleMarconatoMachine(;
        R=0,
        Xd=0,
        Xq=0,
        Xd_p=0,
        Xq_p=0,
        Xd_pp=0,
        Xq_pp=0,
        Td0_p=0,
        Tq0_p=0,
        Td0_pp=0,
        Tq0_pp=0,
        T_AA=0,
        MVABase=0,
        ext=Dict{String, Any}(),
    )
end

"""Get SimpleMarconatoMachine R."""
get_R(value::SimpleMarconatoMachine) = value.R
"""Get SimpleMarconatoMachine Xd."""
get_Xd(value::SimpleMarconatoMachine) = value.Xd
"""Get SimpleMarconatoMachine Xq."""
get_Xq(value::SimpleMarconatoMachine) = value.Xq
"""Get SimpleMarconatoMachine Xd_p."""
get_Xd_p(value::SimpleMarconatoMachine) = value.Xd_p
"""Get SimpleMarconatoMachine Xq_p."""
get_Xq_p(value::SimpleMarconatoMachine) = value.Xq_p
"""Get SimpleMarconatoMachine Xd_pp."""
get_Xd_pp(value::SimpleMarconatoMachine) = value.Xd_pp
"""Get SimpleMarconatoMachine Xq_pp."""
get_Xq_pp(value::SimpleMarconatoMachine) = value.Xq_pp
"""Get SimpleMarconatoMachine Td0_p."""
get_Td0_p(value::SimpleMarconatoMachine) = value.Td0_p
"""Get SimpleMarconatoMachine Tq0_p."""
get_Tq0_p(value::SimpleMarconatoMachine) = value.Tq0_p
"""Get SimpleMarconatoMachine Td0_pp."""
get_Td0_pp(value::SimpleMarconatoMachine) = value.Td0_pp
"""Get SimpleMarconatoMachine Tq0_pp."""
get_Tq0_pp(value::SimpleMarconatoMachine) = value.Tq0_pp
"""Get SimpleMarconatoMachine T_AA."""
get_T_AA(value::SimpleMarconatoMachine) = value.T_AA
"""Get SimpleMarconatoMachine MVABase."""
get_MVABase(value::SimpleMarconatoMachine) = value.MVABase
"""Get SimpleMarconatoMachine ext."""
get_ext(value::SimpleMarconatoMachine) = value.ext
"""Get SimpleMarconatoMachine γd."""
get_γd(value::SimpleMarconatoMachine) = value.γd
"""Get SimpleMarconatoMachine γq."""
get_γq(value::SimpleMarconatoMachine) = value.γq
"""Get SimpleMarconatoMachine states."""
get_states(value::SimpleMarconatoMachine) = value.states
"""Get SimpleMarconatoMachine n_states."""
get_n_states(value::SimpleMarconatoMachine) = value.n_states
"""Get SimpleMarconatoMachine internal."""
get_internal(value::SimpleMarconatoMachine) = value.internal


text:
## Provides a very basic SMC sampler using random walk Metropolis proposals
## the first only allows one-dimensional states, the second allows any
## state of fixed-dimension

## This could clearly be improved to allow proposal covariance matrices to be
## supplied; currently only multipliers for an identity covariance are
## permitted. Perhaps more importantly, I would like to implement a
## mechanism for adaptively choosing the temperatures, e.g. by using the
## conditional ESS of:
## Zhou, Y., Johansen, A.M. and Aston, J.A., 2016. Toward Automatic Model
## Comparison: An Adaptive Sequential Monte Carlo Approach. Journal of
## Computational and Graphical Statistics, 25(3), pp. 701--726.
## A functional, generic SMC sampler could then be included in
## SequentialMonteCarlo.jl

module SMCSampler

using SequentialMonteCarlo
using RNGPool
using StaticArrays
using LinearAlgebra
using Random

## one-dimensional state
mutable struct SMCSP
  x::Float64
  lpibar0::Float64
  lpibar1::Float64
  SMCSP() = new()
end

## for one-dimensional states
function makeRWSMCSampler(mu, lpibar0, lpibar1, betas::Vector{Float64},
  taus::Vector{Float64}, iterations::Vector{Int64})
  n = length(betas)
  @assert n == length(taus) == length(iterations)

  @inline function lG(p::Int64, particle::SMCSP, ::Nothing)
    if p == n
      return 0.0
    end
    @inbounds return (betas[p + 1] - betas[p]) *
      (particle.lpibar1 - particle.lpibar0);
  end
  @inline function M!(newParticle::SMCSP, rng::RNG, p::Int64,
    particle::SMCSP, ::Nothing)
    if p == 1
      newParticle.x = mu(rng)
      newParticle.lpibar0 = lpibar0(newParticle.x)
      newParticle.lpibar1 = lpibar1(newParticle.x)
    else
      y::Float64 = particle.x
      lv0::Float64 = particle.lpibar0
      lv1::Float64 = particle.lpibar1
      for i = 1:iterations[p]
        @inbounds z::Float64 = y + taus[p] * randn(rng)
        z_lv0::Float64 = lpibar0(z)
        z_lv1::Float64 = lpibar1(z)
        @inbounds vy::Float64 = (1.0 - betas[p]) * lv0 + betas[p] * lv1
        @inbounds vz::Float64 = (1.0 - betas[p]) * z_lv0 + betas[p] * z_lv1
        lu::Float64 = -randexp(rng)
        if lu < vz - vy
          y = z
          lv0 = z_lv0
          lv1 = z_lv1
        end
      end
      newParticle.x = y
      newParticle.lpibar0 = lv0
      newParticle.lpibar1 = lv1
    end
  end
  return SMCModel(M!, lG, n, SMCSP, Nothing)
end

function defaultSMCSampler1D()
  function makenormlogpdf(mu::Float64, var::Float64)
    c1::Float64 = - 0.5 / var
    c2::Float64 = - 0.5 * log(2 * π * var)
    @inline function lpdf(x::Float64)
      v::Float64 = x - mu
      return c1*v*v + c2
    end
    return lpdf
  end

  betas::Vector{Float64} = [0, 0.0005, 0.001, 0.0025, 0.005, 0.01, 0.025, 0.05,
    0.1, 0.25, 0.5, 1]
  taus::Vector{Float64} = [0.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0,
    1.0, 1.0]
  iterations::Vector{Int64} = [0, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10]

  x0::Float64 = 0.0
  v0::Float64 = 100.0
  sqrtv0::Float64 = sqrt(v0)

  lpibar0 = makenormlogpdf(x0, v0)
  @inline function mu(rng::RNG)
    return x0 + sqrtv0 * randn(rng)
  end

  lw1::Float64 = log(0.3)
  lw2::Float64 = log(0.7)
  lpdf1 = makenormlogpdf(-10.0, 0.01)
  lpdf2 = makenormlogpdf(10.0, 0.04)

  @inline function lpibar1(x::Float64)
    v1::Float64 = lw1 + lpdf1(x)
    v2::Float64 = lw2 + lpdf2(x)
    m::Float64 = max(v1, v2)
    return m + log(exp(v1 - m) + exp(v2 - m));
  end

  return makeRWSMCSampler(mu, lpibar0, lpibar1, betas, taus, iterations),
    lpibar1
end

## fixed-dimensional state
mutable struct SMCSamplerParticle{d}
  x::SVector{d, Float64}
  lpibar0::Float64
  lpibar1::Float64
  SMCSamplerParticle{d}() where d = new()
end

## scratch space for computations
struct SMCSScratch{d}
  tmp::MVector{d, Float64}
end
SMCSScratch{d}() where d = SMCSScratch{d}(MVector{d, Float64}(undef))

## for fixed-dimensional states
function makeRWSMCSampler(d::Int64, mu, lpibar0, lpibar1, betas::Vector{Float64},
  taus::Vector{Float64}, iterations::Vector{Int64})
  n = length(betas)
  @assert n == length(taus) == length(iterations)

  @inline function lG(p::Int64, particle::SMCSamplerParticle{d1},
    scratch::SMCSScratch{d1}) where d1
    if p == n
      return 0.0
    end
    return (betas[p + 1] - betas[p]) * (particle.lpibar1 - particle.lpibar0);
  end
  @inline function M!(newParticle::SMCSamplerParticle{d1}, rng::RNG, p::Int64,
    particle::SMCSamplerParticle{d1}, scratch::SMCSScratch{d1}) where d1
    if p == 1
      newParticle.x = mu(rng, scratch)
      newParticle.lpibar0 = lpibar0(newParticle.x)
      newParticle.lpibar1 = lpibar1(newParticle.x)
    else
      y::SVector{d1, Float64} = particle.x
      lv0::Float64 = particle.lpibar0
      lv1::Float64 = particle.lpibar1
      for i = 1:iterations[p]
        randn!(rng, scratch.tmp)
        z::SVector{d1, Float64} = y + taus[p] * scratch.tmp
        z_lv0::Float64 = lpibar0(z)
        z_lv1::Float64 = lpibar1(z)
        vy::Float64 = (1.0 - betas[p]) * lv0 + betas[p] * lv1
        vz::Float64 = (1.0 - betas[p]) * z_lv0 + betas[p] * z_lv1
        lu::Float64 = -randexp(rng)
        if lu < vz - vy
          y = z
          lv0 = z_lv0
          lv1 = z_lv1
        end
      end
      newParticle.x = y
      newParticle.lpibar0 = lv0
      newParticle.lpibar1 = lv1
    end
  end
  return SMCModel(M!, lG, n, SMCSamplerParticle{d}, SMCSScratch{d})
end

function defaultSMCSampler()
  function makelogMVN(μ::SVector{2, Float64}, Σ::SMatrix{2, 2, Float64})
    invΣ = inv(Σ)
    lognc = - 0.5 * 2 * log(2 * π) - 0.5 * logdet(Σ)
    function lpdf(x::SVector{d, Float64}) where d
      v = x - μ
      return lognc - 0.5*dot(v, invΣ * v)
    end
    return lpdf
  end

  betas = [0, 0.0005, 0.001, 0.0025, 0.005, 0.01, 0.025, 0.05,
    0.1, 0.25, 0.5, 1]
  taus = [0.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0, 1.0]
  iterations = [0, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10]

  μ0 = SVector{2, Float64}(0.0, 0.0)
  Σ0 = SMatrix{2, 2, Float64}(100.0, 0.0, 0.0, 100.0)
  A0 = SMatrix{2, 2, Float64}(cholesky(Σ0).L)

  lpibar0 = makelogMVN(μ0, Σ0)
  @inline function mu(rng::RNG, scratch::SMCSScratch{2})
    randn!(rng, scratch.tmp)
    return μ0 + A0 * scratch.tmp
  end

  lw1 = log(0.3)
  lw2 = log(0.7)
  μ11 = SVector{2, Float64}(-5.0, -5.0)
  Σ11 = SMatrix{2, 2, Float64}(0.01, 0.005, 0.005, 0.01)
  μ12 = SVector{2, Float64}(5.0, 5.0)
  Σ12 = SMatrix{2, 2, Float64}(0.04, -0.01, -0.01, 0.04)
  lpdf1 = makelogMVN(μ11, Σ11)
  lpdf2 = makelogMVN(μ12, Σ12)

  @inline function lpibar1(x::SVector{2, Float64})
    v1::Float64 = lw1 + lpdf1(x)
    v2::Float64 = lw2 + lpdf2(x)
    m::Float64 = max(v1, v2)
    return m + log(exp(v1 - m) + exp(v2 - m));
  end

  return makeRWSMCSampler(2, mu, lpibar0, lpibar1, betas, taus, iterations),
    lpibar1
end

end


text:
"""
`module  JAC.BascisGenerate`
    ... a submodel of JAC that contains methods that support tasks related to spectroscopic computation.
"""
module BascisGenerate

    using Printf, ..AngularMomentum, ..Basics, ..Continuum, ..Defaults, ..Einstein, ..ManyElectron,
                  ..PhotoEmission, ..Radial

    export generate

    """
    `Basics.generate("condensed multiplet: by single weight", multiplet::Multiplet)`
        ... to condense/reduce the number of CSF in the basis of the given multiplet due to a single 'weight';
            a multiplet::Multiplet is returned.  **Not yet implemented !**
    """
    function Basics.generate(sa::String, multiplet::Multiplet)
        !(sa == "condensed multiplet: by single weight")   &&   error("Unsupported keystring = $sa")
        error("Not yet implemented !")
    end

    """
    `Basics.generate("configuration list: NR, from basis", basis::Basis)`
        ... to (re-) generate the list of NR configurations from the given basis; a confList::Array{Configuration,1} is returned.
    """
    function Basics.generate(sa::String, basis::Basis)
        !(sa == "configuration list: NR, from basis")   &&   error("Unsupported keystring = $sa")
        confList    = Configuration[]
        NoElectrons = sum( basis.csfs[1].occupation )
        for  csf in basis.csfs
            shellList = Dict{Shell,Int64}[]
            for  i = 1:length(csf.occupation)
                n = csf.subshell[i].n;    l = Basics.subshell_l( csf.subshell[i] );    sh = Shell(n, l)
                if     haskey(shellList, sh)    shellList(sh) = shellList(sh) + csf.occupation[i]
                else   shellList = merge( shellList, Dict(sh => csf.occupation[i]) )
                end
            end
            confNew = Configuration( shellList, NoElectrons )
            push!( confList, confNew )
        end
        confListNew = ConfigurationsExcludeDoubles(confList)

        return( confListNew )
    end

    """
    `Basics.generate("configuration list: NR, single-configuration", refConf::Configuration, NoExcitations::Int64, fromShells::Array{Shell,1},`
                     toShells::Array{Shell,1})
        ... to generate a non-relativistic configuration list, including the given reference configuration (refConf) and with
            all configurations that differ by NoExcitations from the fromShells into the toShells; an Array{Configuration,1}
            is returned.
    """
    function Basics.generate(sa::String, refConf::Configuration, NoExcitations::Int64, fromShells::Array{Shell,1}, toShells::Array{Shell,1})
        !(sa == "configuration list: NR, single-configuration")   &&   error("Unsupported keystring = $sa")
        confList = [refConf]
        # First prepare a proper reference configuration that also includes all fromShells and toShells with zero occupation
        shellDict = deepcopy(refConf.shells)
        allShells = deepcopy(fromShells)
        for  sh in toShells
            if     sh in allShells
            else   append!(allShells, [sh] )
            end
        end

        for  sh in allShells
        if  !haskey(shellDict, sh )    shellDict = merge( shellDict, Dict( sh => 0 ) )    end
        end

        rconf = Configuration(shellDict, refConf.NoElectrons)
        println("aa", rconf)

        # Add single excitations
        if      NoExcitations == 0
        elseif  NoExcitations == 1
            for  fromsh in fromShells
                for  tosh in toShells
                    confnew = deepcopy(rconf)
                    confnew.shells[fromsh] = confnew.shells[fromsh] - 1
                    confnew.shells[tosh]   = confnew.shells[tosh]   + 1
                    println("bb", confnew)

                    # Check that this is valid configuration and append, if appropriate
                    add = true
                    for (k,v) in  confnew.shells
                        if  v <= -1  ||  Parity(rconf) != Parity(confnew)    add = false    end
                    end
                    if  add
                        for (k,v) in confnew.shells
                        if  v == 0   delete!(confnew.shells, k)   end
                        end
                        push!(confList, confnew);   println("fromsh = $fromsh  tosh = $tosh")
                    end
                end
            end

        # Add double excitations
        elseif  NoExcitations == 2
            for  fromsha in fromShells
            for  fromshb in fromShells
                for  tosha in toShells
                for  toshb in toShells
                    confnew = deepcopy(rconf)
                    confnew.shells[fromsha] = confnew.shells[fromsha] - 1
                    confnew.shells[fromshb] = confnew.shells[fromshb] - 1
                    confnew.shells[tosha]   = confnew.shells[tosha]   + 1
                    confnew.shells[toshb]   = confnew.shells[toshb]   + 1
                    println("bb", confnew)

                    # Check that this is valid configuration and append, if appropriate
                    add = true
                    for (k,v) in  confnew.shells
                        if  v <= -1  ||  Parity(rconf) != Parity(confnew)    add = false    end
                    end
                    if  add
                        for (k,v) in confnew.shells
                        if  v == 0   delete!(confnew.shells, k)   end
                        end
                        push!(confList, confnew)
                    end
                end
                end
            end
            end

        else
            error("stop a")
        end

        # Exclude all configurations that appear twice (and more) in the list
        confList = Basics.ConfigurationsExcludeDoubles(confList)

        return( confList )
    end

    """
    `Basics.generate("configuration list: relativistic", conf::Configuration)`
        ... to split/decompose a non-relativistic configuration into an list of relativistic ConfigurationR[]. The proper
            occupuation of the relativistic subshells is taken into account.
    """
    function Basics.generate(sa::String, conf::Configuration)
        subshellList = Subshell[]
        confList     = ConfigurationR[]

        !(sa == "configuration list: relativistic")  &&   error("Unsupported keystring = $sa")

        initialize = true;    NoElectrons = 0
        for (k, occ)  in conf.shells
            NoElectrons     = NoElectrons + occ
            wa              = Basics.shellSplitOccupation(k, occ)
            subshellListNew = Dict{Subshell,Int64}[]

            if  initialize
                subshellListNew = wa;    initialize = false
            else
                for  s in 1:length(subshellList)
                    for  a in 1:length(wa)
                        wb = Base.merge( subshellList[s], wa[a] )
                        push!(subshellListNew, wb)
                    end
                end
            end
            subshellList = deepcopy(subshellListNew)
        end

        for subsh in subshellList
        wa = ConfigurationR(subsh, NoElectrons)
        push!(confList, wa)
        end

        return( confList )
    end

    """
    `Basics.generate("CSF list: from single ConfigurationR", conf::ConfigurationR, subshellList::Array{Subshell,1})`
        ... to construct from a given (relativistic) configuration all possible CSF with regard to the subshell order as specified
            by subshellList; a list::Array{CsfR,1} is returned.
    """
    function Basics.generate(sa::String, conf::ConfigurationR, subshellList::Array{Subshell,1})
        parity  = Basics.determineParity(conf)
        csfList = CsfR[];   useStandardSubshells = true;    first = true;    previousCsfs = CsfR[]
        # subhshellList = Subshell[];
        for  subsh in subshellList
            if   subsh in keys(conf.subshells)    occ = conf.subshells[subsh]    else    occ = 0    end
            if   first
                stateList   = ManyElectron.provideSubshellStates(subsh, occ)
                currentCsfs = CsfR[]
                for  state in stateList
                    push!( currentCsfs, CsfR( true, AngularJ64(state.Jsub2//2), parity, [state.occ], [state.nu],
                                                [AngularJ64(state.Jsub2//2)], [AngularJ64(state.Jsub2//2)], Subshell[]) )
                end
                previousCsfs = copy(currentCsfs)
                first        = false
            else
                # Now support also all couplings of the subshell states with the CSFs that were built-up so far
                stateList   = ManyElectron.provideSubshellStates(subsh, occ)
                currentCsfs = CsfR[]
                ##x for  i = 1:length(previousCsfs)    println("generate-aa: ", previousCsfs[i])    end
                for  csf in  previousCsfs
                    for  state in stateList
                        occupation = deepcopy(csf.occupation);    seniority = deepcopy(csf.seniority);
                        subshellJ  = deepcopy(csf.subshellJ);     subshells = deepcopy(csf.subshells)
                        push!(occupation, state.occ);   push!(seniority, state.nu);   push!(subshellJ, AngularJ64(state.Jsub2//2) )
                        push!(subshells, subsh)
                        newXList = oplus( csf.subshellX[end], AngularJ64(state.Jsub2//2) )
                        for  newX in newXList
                            subshellX = deepcopy(csf.subshellX);   push!(subshellX, newX)
                            push!( currentCsfs, CsfR( true, subshellX[end], parity, occupation, seniority, subshellJ, subshellX, Subshell[]) )
                        end
                    end
                end
                previousCsfs = copy(currentCsfs)
            end
        end

        return( previousCsfs )
    end

    """
    `Basics.generate("shells: ordered list for NR configurations", confs::Array{Configuration,1})`
        ... to generate for confs, i.e. all the given (non-relativistic) configurations, a common and ordered shell list;
            a list::Array{Shell,1} is returned.
    """
    function Basics.generate(sa::String, confs::Array{Configuration,1})
        shells = Shell[]

        !(sa == "shells: ordered list for NR configurations")  &&   error("Unsupported keystring = $sa")

        wa = Defaults.getDefaults("ordered shell list: non-relativistic", 7)
        for  a in wa
            for  cf in 1:length(confs)
                ks = keys(confs[cf].shells)
                if  a in ks   push!(shells, a);    break    end
            end
        end

        return( shells )
    end

    """
    `Basics.generate("subshells: ordered list for relativistic configurations", confs::Array{ConfigurationR,1})`
        ... to generate for confs, i.e. all the given (relativistic) configurations, common and ordered subshell list;
            a list::Array{Subshell,1} is returned.
    """
    function Basics.generate(sa::String, confs::Array{ConfigurationR,1})
        subshells = Subshell[]

        !(sa == "subshells: ordered list for relativistic configurations")  &&   error("Unsupported keystring = $sa")

        wa = Defaults.getDefaults("ordered subshell list: relativistic", 7)
        for  a in wa
            for  cf in 1:length(confs)
                ks = keys(confs[cf].subshells)
                ## Do include 'empty' subshells into the subshell list if they are specified by the given configurations
                ## if  a in ks   &&   confs[cf].subshells[a]  !=  0     push!(subshells, a);    break    end
                if  a in ks     push!(subshells, a);    break    end
            end
        end
        ##x println("***subshells = $subshells")
        return( subshells )
    end

    """
    `Basics.generate("subshells: ordered list for two bases", basisA::Basis,  basisB::Basis)`
        ... to generate common and ordered subshell list for the two basis A and B; a list::Array{Subshell,1} is returned.
    """
    function Basics.generate(sa::String, basisA::Basis,  basisB::Basis)
        function areEqual(nn::Int64, sha::Array{Subshell,1}, shb::Array{Subshell,1})
            # Determines whether the first nn subshells are equal in sha and shb (true) or not (false)
            for  i = 1:nx
                if    sha[i] != shb[i]    return( false )   end
            end
            return( true )
        end

        subshells = Subshell[]

        !(sa == "subshells: ordered list for two bases")  &&   error("Unsupported keystring = $sa")

        nx = min(length(basisA.subshells), length(basisB.subshells))
        if  areEqual(nx, basisA.subshells, basisB.subshells)
            # If subshell order is equal, any subshell order is accepted for those subshells that occur in both basis
            for  i = 1:nx
                if    basisA.subshells[i] != basisB.subshells[i]   error("Inconsistent subshells of two bases.")
                else  push!( subshells, basisA.subshells[i])
                end
            end
            #
            if       length(basisA.subshells) > nx
                for  i = nx+1:length(basisA.subshells)    push!(subshells, basisA.subshells[i])    end
            elseif   length(basisB.subshells) > nx
                for  i = nx+1:length(basisB.subshells)    push!(subshells, basisB.subshells[i])    end
            end
        else
            println("basisA.subshells = $(basisA.subshells)")
            println("basisB.subshells = $(basisB.subshells)")
            # If subshell order is NOT equal, all subshell in basisA and basisB must follow standard order
            standardList = Defaults.getDefaults("ordered subshell list: relativistic", 7)
            na = 0;  nb = 0
            for  sh in standardList
                if  sh in basisA.subshells  ||   sh in basisB.subshells   push!( subshells, sh)    end
            end
            nn = 0;
            for sh in basisA.subshells
                if    !(sh in subshells) push!( subshells, sh)
                else  wb = findall(x->x==sh, subshells);    if wb[1] <= nn   error("stop a")   else   nn = wb[1]  end
                end
            end
            nn = 0;
            for sh in basisB.subshells
                if    !(sh in subshells) push!( subshells, sh)
                else  wb = findall(x->x==sh, subshells);    if wb[1] <= nn   error("stop a")   else   nn = wb[1]  end
                end
            end
            println("*** extended subshells from two basis = $subshells")
        end
        ##x println("subshells from two basis = $subshells")

        return( subshells )
    end

    """
    `Basics.generate("single-electron spectrum: STO", N::Int64, potential::Radial.Potential, grid::Radial.Grid; N_0::Int64=30, alpha_0::Float64=1.0,
                     beta_0::Float64=1.1)`
        ... to generate a complete one-electron spectrum with N positive and N negative states, and by using even-tempered Slater-type
            orbitals (STO) with parameters ``\alpha_i = \alpha_0 \beta_0^i``; a spectrum::SingleElecSpectrum is returned where just
            N0 positive and N_0 negative are kept for later use.  **Not yet implemented !**

    `Basics.generate("single-electron spectrum: STO, positive", N::Int64, potential::Radial.Potential, grid::Radial.Grid; N_0::Int64=30, alpha_0::Float64=1.0,
                     beta_0::Float64=1.1)`
        ... to generate the same but to return only the N_0 positive states.  **Not yet implemented !**
    """
    function Basics.generate(sa::String, N::Int64, potential::Radial.Potential, grid::Radial.Grid; N_0::Int64=30, alpha_0::Float64=1.0, beta_0::Float64=1.1)
        !(sa == "single-electron spectrum: STO"  &&  sa == "single-electron spectrum: STO, positive")   &&   error("Unsupported keystring = $sa")
        error("Not yet implemented !")

        #= * define N normalized states in stoplus::Array{Vector{Float64},1} and stominus::Array{Vector{Float64},1}
            on the given grid and with given parameters
        * calculate one-electron (NxN) Hamiltonian matrix with potential h_pq and overlap S_pq, p,q = 1..N
        * solve generalized eigenvalue problem h c_vec = epsilon S c_vec
        * set spectrum::SingleElecSpectrum with just N0 functions
        * check normalization and orthogonality of the functions   =#

        return( nothing )
    end

    """
    `Basics.generateConfigurationsForExcitationScheme(confs::Array{Configuration,1}, exScheme::Basics.NoExcitationScheme,
                                                    nMax::Int64, lValues::Array{Int64,1})`
        ... generates a list of non-relativistic configurations for the given (reference) confs and the excitation scheme.
            All orbitals in standard order are considered up to maxShell = (n_max, l_max).
    """
    function Basics.generateConfigurationsForExcitationScheme(confs::Array{Configuration,1}, exScheme::Basics.NoExcitationScheme,
                                                            nMax::Int64, lValues::Array{Int64,1})
        @warn(  "No excitations are included if exScheme::NoExcitationScheme. ")
        println("No excitations are included if exScheme::NoExcitationScheme. ")
        #
        newConfList = confs
        return( newConfList )
    end

    """
    `Basics.generateConfigurationsForExcitationScheme(confs::Array{Configuration,1}, exScheme::Basics.DeExciteSingleElectron,
                                                    nMax::Int64, lValues::Array{Int64,1})`
        ... generates a list of non-relativistic configurations for the given (reference) confs and the excitation scheme.
            All orbitals in standard order are considered up to maxShell = (n_max, l_max).
    """
    function Basics.generateConfigurationsForExcitationScheme(confs::Array{Configuration,1}, exScheme::Basics.DeExciteSingleElectron,
                                                            nMax::Int64, lValues::Array{Int64,1})
        confList = Configuration[];    NoElectrons = confs[1].NoElectrons
        shellList   = Basics.generateShellList(confs, nMax, lValues)

        # Create all configurations with single excitations from the given list of configurations
        for  conf in confs
            # Take one electron fromShell and add toShell
            for  (fromShell,occ)  in  conf.shells
                for  toShell  in  shellList
                    newShells = deepcopy( conf.shells )
                    ##x print("newShells = $newShells   fromShell = $fromShell    toShell = $toShell")
                    if      fromShell == toShell
                    elseif  haskey(conf.shells, toShell )   fromOcc = occ;   toOcc = conf.shells[toShell]
                            if  fromOcc - 1 < 0                        println("");    continue    end
                            if  toOcc   + 1 > 2*(2*toShell.l + 1)      println("");    continue    end
                            newShells[fromShell] = newShells[fromShell] - 1
                            newShells[toShell]   = newShells[toShell] + 1
                    else    fromOcc = occ
                            if  fromOcc - 1 < 0                        println("");    continue    end
                            newShells[fromShell] = newShells[fromShell] - 1
                            newShells = Base.merge( newShells, Dict( toShell => 1))
                    end
                    # Add a new configuration
                    push!( confList, Configuration( newShells, NoElectrons))
                    println("  newconf = $(Configuration( newShells, NoElectrons))")
                end
            end
        end

        nbefore     = length(confList)
        newConfList = Configuration[]
        for  conf  in  confList
            if  conf in newConfList    continue;    else    push!( newConfList,  conf)    end
        end
        nafter      = length(newConfList)
        println("Basics.generateConfigurationsForExcitationScheme():: Number of configs is $nbefore (before) and $nafter (after).")

        #
        ##x push!( newConfList, confs[1])
        ##x push!( newConfList, Configuration("[Ne] 3s 3p^6 4s") )
        ##x push!( newConfList, Configuration("[Ne] 3s 3p^6 5s") )
        return( newConfList )
    end

    """
    `Basics.generateLevelWithExtraElectron(newOrbital, symt::LevelSymmetry, level::Level)`
        ... generates a (new) level with one extra electron in subshell newOrbital.subshell and with overall symmetry symt. The function
            assumes that all CSF in the basis of level have the same symmetry as the level itself, that the new subshell is not yet part
            of the basis and that CSF with the total symmetry can be constructed; it terminates with an error if one of this assumptions
            is violated. A newLevel::Level with the same number of CSF and the same representation (eigenvector) as the given level is
            returned. From a physics viewpoint, newLevel refers to a antisymmetrized product state of [level x orbital(sh, energy)] J^P
            (symt) and with an total energy = level.energy + newOrbital.energy.
    """
    function Basics.generateLevelWithExtraElectron(newOrbital, symt::LevelSymmetry, level::Level)
        basis = level.basis;    newSubshells = copy(basis.subshells);    newCsfs = CsfR[];   J = level.J;   parity = level.parity

        push!(newSubshells, newOrbital.subshell)
        Defaults.setDefaults("relativistic subshell list", newSubshells; printout=false)
        newOrbitals = Base.merge( copy(basis.orbitals), Dict( newOrbital.subshell => newOrbital ))

        for  i = 1:length(basis.csfs)
            stateList   = ManyElectron.provideSubshellStates(newOrbital.subshell, 1)
            if  basis.csfs[i].J != J   ||   basis.csfs[i].parity != parity    error("Improper symmetry of CSF.")             end
            if  length(stateList) != 1                                        error("Improper number of subshell states.")   end
            if  AngularMomentum.triangularDelta(J, Basics.subshell_j(newOrbital.subshell), symt.J) != 1
                                                                            error("Improper coupling of (new) subshell & total J.")   end
            substate   = stateList[1]
            occupation = copy(basis.csfs[i].occupation);    push!(occupation, substate.occ)
            seniority  = copy(basis.csfs[i].seniority);     push!(seniority,  substate.nu)
            subshellJ  = copy(basis.csfs[i].subshellJ);     push!(subshellJ,  AngularJ64(substate.Jsub2//2) )
            subshellX  = copy(basis.csfs[i].subshellX);     push!(subshellX,  symt.J)
            push!(newCsfs, CsfR(true, symt.J, symt.parity, occupation, seniority, subshellJ, subshellX, Subshell[] ) )
        end

        newBasis = Basis(true, basis.NoElectrons+1, newSubshells, newCsfs, copy(basis.coreSubshells), newOrbitals)
        newLevel = Level(symt.J, AngularM64(symt.J), symt.parity, level.index, level.energy+newOrbital.energy, level.relativeOcc,
                        true, newBasis, level.mc)
        return( newLevel )
    end

    """
    `Basics.generateLevelWithExtraSubshell(sh::Subshell, level::Level)`
        ... generates a (new) level with one extra subshell sh but with the same overall symmetry as before. The function assumes that the
            new subshell is not yet part of the basis; it terminates with an error if one of this assumptions is violated.
            A newLevel::Level with the same symmetry, number of CSF, energy and the same representation (eigenvector) as the given
            level is returned.
    """
    function Basics.generateLevelWithExtraSubshell(sh::Subshell, level::Level)
        basis = level.basis;    newSubshells = copy(basis.subshells);    newCsfs = CsfR[];   J = level.J;   parity = level.parity
        newOrbitals = copy(basis.orbitals)

        push!(newSubshells, sh)
        Defaults.setDefaults("relativistic subshell list", newSubshells; printout=false)

        for  i = 1:length(basis.csfs)
            stateList   = ManyElectron.provideSubshellStates(sh, 0)
            if  length(stateList) != 1                                        error("Improper number of subshell states.")   end
            substate   = stateList[1]
            occupation = copy(basis.csfs[i].occupation);    push!(occupation, substate.occ)
            seniority  = copy(basis.csfs[i].seniority);     push!(seniority,  substate.nu)
            subshellJ  = copy(basis.csfs[i].subshellJ);     push!(subshellJ,  AngularJ64(substate.Jsub2//2) )
            subshellX  = copy(basis.csfs[i].subshellX);     push!(subshellX,  subshellX[end])
            push!(newCsfs, CsfR(true, basis.csfs[i].J, basis.csfs[i].parity, occupation, seniority, subshellJ, subshellX, Subshell[] ) )
        end

        newBasis = Basis(true, basis.NoElectrons, newSubshells, newCsfs, copy(basis.coreSubshells), newOrbitals)
        newLevel = Level(level.J, level.M, level.parity, level.index, level.energy, level.relativeOcc, true, newBasis, level.mc)
        return( newLevel )
    end

    """
    `Basics.generateLevelWithSymmetryReducedBasis(level::Level)`
        ... generates a level with a new basis and representation that only includes the CSF with the same symmetry as the level
            itself. Both, the underlying CSF basis and the eigenvectors are adopted accordingly. A (new) level::levelR is returned.
    """
    function Basics.generateLevelWithSymmetryReducedBasis(level::Level)
        basis = level.basis;    newCsfs = CsfR[];    newMc = Float64[]

        for  i = 1:length(basis.csfs)
            if  basis.csfs[i].J == level.J  &&   basis.csfs[i].parity == level.parity
                push!(newCsfs, basis.csfs[i]);    push!(newMc, level.mc[i])
            end
        end
        newBasis = Basis(true, basis.NoElectrons, basis.subshells, newCsfs, basis.coreSubshells, basis.orbitals)
        newLevel = Level(level.J, level.M, level.parity, level.index, level.energy, level.relativeOcc, true, newBasis, newMc)
        return( newLevel )
    end

    """
    `Basics.generateShellList(confs::Array{Configuration,1}, nMax::Int64, lValues::Array{Int64,1})`
        ... generates a list of (non-relativistic) shells in standard order that contains all shells from the given list of
            configurations as well as the shells up to the principal quantum number nMax and the list of orbital quantum
            numbers lValues. A shellList::Array{Shell,1} is returned.
    """
    function Basics.generateShellList(confs::Array{Configuration,1}, nMax::Int64, lValues::Array{Int64,1})
        shellList = Shell[]
        # Add all shells from the given configurations
        for  conf in confs
            for (k,v) in conf.shells    push!( shellList, k)    end
        end
        # Add all shells for the given quantum numbers
        for  n = 1:nMax
            for  ll in  lValues
                if  n >= ll + 1     push!( shellList, Shell(n,ll))      end
            end
        end
        ##x println("shellList = $shellList ")

        # Now bring the shells in standard order
        newShellList = Shell[]
        for  n = 1:nMax
            for  l = 0:30
                if  n >= l + 1  && Shell(n,l) in shellList    push!( newShellList, Shell(n,l))     end
            end
        end
        println("newShellList = $newShellList ")

        return( newShellList )
    end

end # module


text:
using MCMCBenchmarks,CSV

ProjDir = @__DIR__
cd(ProjDir)

# Load results
folder = "2019_10_14T07_44_00/"
results = CSV.read(ProjDir * "/results/" * folder * "results.csv")

# Make plots
include("primary_plots.jl")