text:
function make_model(InPath::AbstractString; OutPath="./autogeneratedmodel",
                    Host="bacteria", Model="Kinetics", Lang="julia")

  path_head = Base.@__DIR__
  # Error report: line number --> error message
  error_report = Dict{Int64, Array}()

  # include(joinpath(path_head, "preprocessor3.jl"))
  # include("preprocessor3.jl")
  println("\n---------loading sentences------------")
  # inputSentencesArray = Array of pairs of (line, line #)
  inputSentencesArray = getRidOfNewlineHashcommentBlankline(InPath)
  println("\n>>>>>tokenizing sentences" * ">"^10)
  # tokenizedInputSentencesArray = Array of pairs of (tokens, line #)
  tokenizedInputSentencesArray = sentenceTokenization2(inputSentencesArray)
  println("\n>>>>>normalizing and tagging" * ">"^10)
  reservedWordsPath = joinpath(path_head, "reservedWords.jl")
  # reservedWordsPath = "reservedWords.jl"
  # taggedSentencesArray = Array of pairs of (taggedTokens, line #)
  taggedSentencesArray = tokenClassification2(tokenizedInputSentencesArray, reservedWordsPath)
  # reshape for output observation
  printTagSen = [["$(y[1])/$(y[2])" for y in x] for (x, id) in taggedSentencesArray]
  # println(typeof(printTagSen))
  foreach(println, [join(ts, "  ") for ts in printTagSen])

  # include(joinpath(path_head, "biosymDecoder.jl"))
  # include("biosymDecoder.jl")
  println("\n------------information extraction------------")
  include(reservedWordsPath)
  # BioSymVerbInfo = Array of pairs of ((senType, line#), BioSym, paraSet)
  BioSymVerbInfo = extractBioSymVerbInformation(taggedSentencesArray,
                   reservedWords["SentenceType"], error_report)
  # println("print by foreach:")
  foreach(println, BioSymVerbInfo)
  println("\n>>>>>decoding bio symbols" * ">"^10)
  decodingBioSymGroups(BioSymVerbInfo, error_report)
  # BioSymVerbInfo = Array of pairs of ((senType, line#), BioSym, paraSet)
  println("\n>>>>>decoding bio symbols results" * ">"^10)
  printArrayOfTripleToken(BioSymVerbInfo, AbstractString)
  println("\n>>>>>type conversion dictionary" * ">"^10)
  typeConversionDict = setUpSymbolConvertionDict(BioSymVerbInfo, error_report)
  foreach(println, typeConversionDict)
  println("\n>>>>>replace each biosymbol string with generalBioSym
          composite" * ">"^10)
  # newVerbBiosymInfo = Array of ((senType, line#), BioSym, paraSet)
  newVerbBiosymInfo = replaceEachBioSymWithGeneralBioSym(BioSymVerbInfo,
                      typeConversionDict, error_report)
  printArrayOfTripleToken(newVerbBiosymInfo, generalBioSym)

  # include(joinpath(path_head, "semanticChecking.jl"))
  # include("semanticChecking.jl")
  println("\n------------right before IR generation-----------------")
  # preIR_DS = array of tuple ((verb, line #), reactant, product, catalyst)
  preIR_DS = semanticCheckingForEachTriplet(newVerbBiosymInfo,
             typeConversionDict, error_report)
  printArrayOfTripleToken(preIR_DS, generalBioSym)

  #######################
  # TO HERE, no more error report after this line
  #######################
  # print errors
  if length(error_report) != 0
    println("\n\n\n#################### Error report: Useful tips ########")
    for lineNum in sort(collect(keys(error_report)))
      println("in line $lineNum -->", )
      for err in error_report[lineNum]
        println("\t$(err)")
      end
    end
    exit()
  end

  println("\n------------model generation preparation------------------")
  # include(joinpath(path_head, "modelGeneration.jl"))
  # include("modelGeneration.jl")
  sys2userDict = Dict()
  for (key, val) in typeConversionDict
    sys2userDict[val] = key
  end
  # reform to match with NML_V1
  # preIR_DS = array of tuple ((verb, line #), [reactant, product, catalyst])
  (rnx, txtl, rnx_set, mRNA_set, m_protein_set) = preparing_rnxList_txtlDict(
    preIR_DS, sys2userDict, typeConversionDict)

  println(">>>>>rnx list" * ">"^10)
  foreach(x->println(x.rnxName), rnx)
  # Sorting, insert "BIOMASS" here
  (sorted_rnx_species_array, rnx_species2index_dict, rnx_index2species_dict,
   sorted_all_species_array, all_species2index_dict, all_index2species_dict, extra_species_num) =
    sorting_species_list(rnx_set, mRNA_set, m_protein_set, sys2userDict)
  println("------------")
  foreach(println, sorted_rnx_species_array)
  println("------------")
  foreach(println, sorted_all_species_array)

  # write program to disk
  output_file = Dict{String, String}()  # put all files in this dict
  # stoichiometric_matrix
  stoichiometric_matrix_buffer = build_stoichiometric_matrix_buffer(rnx, rnx_species2index_dict)
  output_file["stoichiometry.dat"] = stoichiometric_matrix_buffer

  println("\n------------Modeling Framework-----")
  # DISTRIBUTION file mapping
  DISTRIBUTION = Dict(
      "julia"   => ("distribution/dis.jl", "JuliaStrategy.jl"),
      "python2" => ("distribution/dis.py2", "Python2Strategy.jl"),
      "python"  => ("distribution/dis.py3", "Python3Strategy.jl"),
      "python3" => ("distribution/dis.py3", "Python3Strategy.jl"),
      "matlab"  => ("distribution/dis.m", "MATLABStrategy.jl")
  )
  target_lang = Lang
  model_frame = Model
  if target_lang == "julia"
    # include(joinpath(path_head, DISTRIBUTION["julia"][2]))
    fn_prefix = "jl_"
    file_suffix = ".jl"
    dis_file_folder = DISTRIBUTION["julia"][1]
  elseif ((target_lang == "python3") || (target_lang == "python"))
    # include(joinpath(path_head, DISTRIBUTION["python3"][2]))
    fn_prefix = "py3_"
    file_suffix = ".py"
    dis_file_folder = DISTRIBUTION["python3"][1]
  elseif target_lang == "python2"
    # include(joinpath(path_head, DISTRIBUTION["python2"][2]))
    fn_prefix = "py2_"
    file_suffix = ".py"
    dis_file_folder = DISTRIBUTION["python2"][1]
  elseif  target_lang == "matlab"
    # include(joinpath(path_head, DISTRIBUTION["matlab"][2]))
    fn_prefix = "m_"
    file_suffix = ".m"
    dis_file_folder = DISTRIBUTION["matlab"][1]
  else
    println("Generate model in julia as default")
    # include(joinpath(path_head, DISTRIBUTION["julia"][2]))
    fn_prefix = "jl_"
    target_lang = "julia"
    file_suffix = ".jl"
    dis_file_folder = DISTRIBUTION["julia"][1]
  end

  fnmap = Dict(
      "julia"   => (jl_build_kinetics_buffer, jl_build_data_dictionary_buffer,
                    jl_build_simulation_buffer, jl_build_solveODEBalances_buffer,
                    jl_generate_FBA_data_dictionary),
      "python2" => (py2_build_kinetics_buffer, py2_build_data_dictionary_buffer,
                    py2_build_simulation_buffer, py2_build_solveODEBalances_buffer,
                    py2_generate_FBA_data_dictionary),
      "python" => (py3_build_kinetics_buffer, py3_build_data_dictionary_buffer,
                    py3_build_simulation_buffer, py3_build_solveODEBalances_buffer,
                    py3_generate_FBA_data_dictionary),
      "python3" => (py3_build_kinetics_buffer, py3_build_data_dictionary_buffer,
                    py3_build_simulation_buffer, py3_build_solveODEBalances_buffer,
                    py3_generate_FBA_data_dictionary),
      "matlab"  => (m_build_kinetics_buffer, m_build_data_dictionary_buffer,
                    m_build_simulation_buffer, m_build_solveODEBalances_buffer,
                    m_generate_FBA_data_dictionary)
  )

  # Generate FBA model
  if model_frame in ["FBA", "FVA"]
    if target_lang == "matlab"
        data_dictionary_buffer, maximize_product_buffer = fnmap[target_lang][5](
               rnx, sorted_rnx_species_array, extra_species_num)
        output_file["DataDictionary" * file_suffix] = data_dictionary_buffer
        output_file["maximizeProductDictionary" * file_suffix] = maximize_product_buffer
    else
        data_dictionary_buffer = fnmap[target_lang][5](rnx,
                                   sorted_rnx_species_array, extra_species_num)
        output_file["DataDictionary" * file_suffix] = data_dictionary_buffer
    end
    # write & copy
    output_dir = OutPath
    if !(isdir(output_dir))
      mkpath(output_dir)
    end
    for (file_name, file) in output_file
      write(joinpath(output_dir, file_name), file)
    end
    for file in readdir(joinpath(path_head, dis_file_folder))
      if occursin(file_suffix, file)
          println(file)
          f = joinpath(path_head, dis_file_folder, file)
          cp(f, joinpath(output_dir, file); force=true)
      end
    end
    cp(InPath, joinpath(output_dir, basename(InPath)); force=true)
  # Generate Kinetic model
  else
    (kinetics_buffer, Monod_const, W_array, disassociation_const) =
      fnmap[target_lang][1](all_species2index_dict, rnx, txtl, sys2userDict)
    output_file["Kinetics" * file_suffix] = kinetics_buffer

    host_type = Host
    data_dictionary_buffer = fnmap[target_lang][2](host_type, sorted_all_species_array,
      all_species2index_dict,
      sorted_rnx_species_array, rnx, txtl, Monod_const, W_array, disassociation_const,
      collect(mRNA_set), collect(m_protein_set))
    output_file["DataDictionary" * file_suffix] = data_dictionary_buffer

    ODE_simulation_buffer = fnmap[target_lang][3](extra_species_num)
    output_file["Balances" * file_suffix] = ODE_simulation_buffer

    solveODEBalances_buffer = fnmap[target_lang][4](sorted_all_species_array,
      all_species2index_dict, collect(mRNA_set), collect(m_protein_set))
    output_file["SolveBalances" * file_suffix] = solveODEBalances_buffer

    output_dir = OutPath
    if !(isdir(output_dir))
      mkpath(output_dir)
    end
    if target_lang == "julia"
      cp(joinpath(path_head, "include.jl"), joinpath(output_dir, "include.jl"); force=true)
    end
    cp(InPath, joinpath(output_dir, basename(InPath)); force=true)
    for (file_name, file) in output_file
      write(joinpath(output_dir, file_name), file)
    end
  end

end


text:
export EqconstrainedOEDParam

import jInv.ForwardShare.ForwardProbType

"""
type EqconstrainedOEDParam <: ForwardProbType

Parameter for OED with equality constrained lower-level problem

	f(p) = argmin_f 0.5*| M(p)*f - d(p)|^2 + 0.5*|L*f|^2
	       subject to Ce*f = ce

STATUS: this is under development

Fields:

	Design  - design type (OEDTypeAParam)
	L       - Regularization Operator
	Ce      - equality constraint matrix
	ce      - equality constraints
	Xr      - field for reconstructions
	Ainv    - field for linear solvers
"""
type EqconstrainedOEDParam <: ForwardProbType
   Design::AbstractOEDParam
   L::AbstractArray{Float64}
   Ce::AbstractArray{Float64,2}
   ce::AbstractArray{Float64}
   Xr::AbstractArray{Float64,2}
   Ainv
end

function getKKT(H::SparseMatrixCSC{Float64}, Ce::AbstractArray)
	ne = size(Ce,1)
	if ne==0
		return H
	else
		return [H Ce';Ce spzeros(ne,ne)]
	end
end

function getKKT(H::Array{Float64}, Ce::AbstractArray)
	ne = size(Ce,1)
	if ne==0
		return H
	else
		return [H Ce';Ce 0]
	end
end

getNumberOfConstraints(pFor::EqconstrainedOEDParam) = size(pFor.Ce,1)

function Base.display(pFor::EqconstrainedOEDParam)
	(m,n,nrhs) = getProblemSize(pFor.Design)
	println("--- OED.EqconstrainedOEDParam <: ForwardProbType ---")
	println("Type of Design:                    $(typeof(pFor.Design))")
	println("Size of regularization operator:   $(size(pFor.L))")
	println("Number of constraints:             $(size(pFor.Ce,1))")
end

import jInv.ForwardShare.getData
function getData(p::Vector{Float64},pFor::EqconstrainedOEDParam,doClear=false)

	# get design
	M  = getM(p,pFor.Design,true)
	D  = getD(p,pFor.Design)
	Ce = pFor.Ce
	ce = pFor.ce
	L  = pFor.L
	n  = size(M,2)
	ne = size(Ce,1)

	# build and factorize KKT system
	H   = M'*M + L'*L
	KKT = getKKT(H,Ce)
    pFor.Ainv = lufact(KKT)

    xl  = pFor.Ainv\[M'*D;ce]
    Xr  = xl[1:n,:]

	pFor.Xr = Xr
    return Xr,pFor
end

ext: jl


text:
#shifts and memoizes the Gauss-Lobatto nodes for optical depth in a layer
@memoize function lobattonodes(n::Int64)::NTuple{2,Vector{Float64}}
    #gauss-lobatto nodes for integrating optical depth in pressure layer
    x, w = gausslobatto(n)
    #simply shift to [0,1]
    𝓍 = @. (x + 1)/2
    𝓌 = @. w/2
    return 𝓍, 𝓌
end

function lobattoevaluations(P, 𝒻T::Q, 𝒻μ::R, nlobatto::Int) where {Q,R}
    #setup new arrays
    np = length(P)
    U = typeof(𝒻T(P[1]))
    T = Array{U,2}(undef, nlobatto, np-1)
    μ = Array{U,2}(undef, nlobatto, np-1)
    #Gauss-Lobatto quadrature nodes, prescaled to [0,1]
    𝓍, _ = lobattonodes(nlobatto)
    for i ∈ 1:nlobatto
        @inbounds for j ∈ 1:np-1
            ΔP = P[j+1] - P[j]
            Pᵢ = P[j] + ΔP*𝓍[i]
            Tᵢ = 𝒻T(Pᵢ)
            μᵢ = 𝒻μ(Tᵢ,Pᵢ)
            T[i,j] = Tᵢ
            μ[i,j] = μᵢ
        end
    end
    return T, μ
end

function temperaturederivs(P, 𝒻T::Q, ϵ::Real=0.01) where {Q}
    n = length(P)
    ∂T = similar(P)
    δP = (P[2] - P[1])*ϵ
    ∂T[1] = (𝒻T(P[1] + δP) - 𝒻T(P[1]))/δP
    for i = 2:n-1
        δP = min(P[i+1] - P[i], P[i] - P[i-1])*ϵ
        ∂T[i] = (𝒻T(P[i] + δP) - 𝒻T(P[i] - δP))/(2*δP)
    end
    δP = (P[n] - P[n-1])*ϵ
    ∂T[n] = (𝒻T(P[n]) - 𝒻T(P[n] - δP))/δP
    return ∂T
end

function planckevaluations(P, ν, 𝒻T::Q) where {Q}
    U = typeof(planck(ν[1], 𝒻T(P[1])))
    B = Array{U,2}(undef, length(P), length(ν))
    #B′ = similar(B)
    @threads for i ∈ eachindex(P)
        Tᵢ = 𝒻T(P[i])
        @inbounds for j ∈ eachindex(ν)
            B[i,j] = planck(ν[j], Tᵢ)
            #B′[i,j] = dplanck(ν[j], Tᵢ)
        end
    end
    return B#, B′
end

function βevaluations(P, g, 𝒻T::Q, 𝒻μ::R, 𝒜::S) where {Q,R,S<:AbstractAbsorber}
    ν = 𝒜.ν
    C = 1e-4*𝐍𝐚/g
    β₀ = 𝜷(P[1], 𝒻T(P[1]), 𝒻μ(P[1], 𝒻T(P[1])), C, 1, 𝒜)
    β = Array{typeof(β₀),2}(undef, length(P), length(ν))
    @threads for j ∈ eachindex(ν)
        @inbounds for i ∈ eachindex(P)
            Pᵢ = P[i]
            Tᵢ = 𝒻T(Pᵢ)
            μᵢ = 𝒻μ(Tᵢ, Pᵢ)
            β[i,j] = 𝜷(Pᵢ, Tᵢ, μᵢ, C, j, 𝒜)
        end
    end
    return β
end

function 𝜷(P, T, μ, C, idx::Int, 𝒜::Q) where {Q<:AbstractAbsorber}
    #total cross-section [cm^2/molecule]
    σ = Σ(𝒜, idx, T, P)
    #absorption coefficient [m*s^2/kg], represents dτ/dP
    C*(σ/μ)
end

#the "linear in τ" approximation for planck emission of a discrete layer, see:
#  Clough, S. A., Iacono, M. J. & Moncet, J.-L. Line-by-line calculations of atmospheric fluxes and cooling rates: Application to water vapor. J. Geophys. Res. 97, 15761 (1992).
layerplanck(B₁, B₂, τ, t) = B₂*(1.0 - t) - (B₁ - B₂)*t + (1.0 - t)*(B₁ - B₂)/τ

layerplanck(B₁, B₂, τ) = layerplanck(B₁, B₂, τ, exp(-τ))

#-------------------------------------------------------------------------------
# optical depth paths

function 𝒹depth(P::AbstractVector{<:Real}, #ascending
                T::Matrix{<:Real}, #evaluated at Lobatto nodes
                μ::Matrix{<:Real}, #evaluated at Lobatto nodes
                𝒜::Q,
                idx::Int,
                C::Real,
                𝓂::Real, #1/cos(θ)
                nlobatto::Int) where {Q<:AbstractAbsorber}
    #setup
    np = length(P)
    𝓍, 𝓌 = lobattonodes(nlobatto)

    #first quadrature point is right at the surface
    @inbounds β₁ = 𝜷(P[1], T[1,1], μ[1,1], C, idx, 𝒜)
    #accumulating value
    τ = zero(β₁)
    #loop over pressure layers
    @inbounds for i ∈ 1:np-1
        #beginning and end of pressure layer
        P₁, P₂ = P[i], P[i+1]
        ΔP = P₂ - P₁
        #integrate optical depth over layer with Gauss-Lobotto quadrature
        τᵢ = zero(eltype(τ))
        #first β is reused from end of previous layer
        τᵢ += (ΔP*𝓌[1])*β₁
        for n ∈ 2:nlobatto-1
            #pressure coordinate of quadrature node, tricky because decreasing!
            Pₙ = P₁ + ΔP*𝓍[n]
            #absorption coefficient
            β = 𝜷(Pₙ, T[n,i], μ[n,i], C, idx, 𝒜)
            #contribute to optical depth integral
            τᵢ += (ΔP*𝓌[n])*β
        end
        #final point at end of pressure layer
        βₙ = 𝜷(P₂, T[end,i], μ[end,i], C, idx, 𝒜)
        τᵢ += (ΔP*𝓌[end])*βₙ
        #the end is now the beginning
        β₁ = βₙ
        #store τ
        τ += τᵢ*𝓂
    end
    return τ
end

function 𝒹depth!(τ::AbstractVector{Q}, #to be filled
                 P::AbstractVector{<:Real}, #ascending
                 T::Matrix{<:Real}, #evaluated at Lobatto nodes
                 μ::Matrix{<:Real}, #evaluated at Lobatto nodes
                 𝒜::R,
                 idx::Int,
                 C::Real,
                 nlobatto::Int)::Nothing where {Q<:Real,R<:AbstractAbsorber}
    #setup
    np = length(P)
    𝓍, 𝓌 = lobattonodes(nlobatto)
    τmin = 1e-6*one(Q)

    #first quadrature point
    β₁ = 𝜷(P[1], T[1,1], μ[1,1], C, idx, 𝒜)
    #loop over pressure layers
    @inbounds for i ∈ 1:np-1
        #beginning and end of pressure layer
        P₁, P₂ = P[i], P[i+1]
        ΔP = P₂ - P₁
        #integrate optical depth over layer with Gauss-Lobotto quadrature
        τᵢ = zero(Q)
        #first β is reused from end of previous layer
        τᵢ += (ΔP*𝓌[1])*β₁
        for n ∈ 2:nlobatto-1
            #pressure coordinate of quadrature node, tricky because decreasing!
            Pₙ = P₁ + ΔP*𝓍[n]
            #intra-layer absorption coefficient isn't pre-calculated
            βₙ = 𝜷(Pₙ, T[n,i], μ[n,i], C, idx, 𝒜)
            #contribute to optical depth integral
            τᵢ += (ΔP*𝓌[n])*βₙ
        end
        #final point at end of pressure layer
        βₙ = 𝜷(P₂, T[end,i], μ[end,i], C, idx, 𝒜)
        τᵢ += (ΔP*𝓌[end])*βₙ
        #the end is now the beginning for the next layer
        β₁ = βₙ
        #put a floor on τ and store it
        τ[i] = max(τᵢ, τmin)
    end
    nothing
end

#-------------------------------------------------------------------------------
# streams/fluxes up the atmosphere, no storage

function 𝒹streams(P::AbstractVector{<:Real},
                  T::Matrix{<:Real}, #evaluated at Lobatto nodes
                  μ::Matrix{<:Real}, #evaluated at Lobatto nodes
                  𝒜::Q,
                  idx::Int,
                  I₀::Real,
                  C::Real,
                  nstream::Int,
                  nlobatto::Int) where {Q<:AbstractAbsorber}
    #setup
    np = length(P)
    𝓂, 𝒲 = streamnodes(nstream)
    𝓍, 𝓌 = lobattonodes(nlobatto)
    ν = 𝒜.ν[idx]
    τmin = 1e-6*one(eltype(τ))

    #output flux [W/m^2/cm^-1]
    F = zero(I₀)
    #𝜷 coefficient at the surface
    βₛ = 𝜷(P[1], T[1,1], μ[1,1], C, idx, 𝒜)
    #iterate through streams of different angles

    @inbounds for k ∈ 1:nstream
        #initial planck emission from surface
        I = convert(typeof(F), I₀)
        #first quadrature point is right at the surface
        β₁ = βₛ
        #loop over pressure layers
        for i ∈ 1:np-1
            #beginning and end of pressure layer
            P₁, P₂ = P[i], P[i+1]
            ΔP = P₁ - P₂
            #integrate optical depth over layer with Gauss-Lobotto quadrature
            τ = zero(I)
            #first β is reused from end of previous layer
            τ += (ΔP*𝓌[1])*β₁*𝓂[k]
            for n ∈ 2:nlobatto-1
                #pressure coordinate of quadrature node, tricky because decreasing!
                Pₙ = P₁ - ΔP*𝓍[n]
                #absorption coefficient
                β = 𝜷(Pₙ, T[n,i], μ[n,i], C, idx, 𝒜)
                #contribute to optical depth integral
                τ += (ΔP*𝓌[n])*β*𝓂[k]
            end
            #final point at end of pressure layer
            βₙ = 𝜷(P₂, T[end,i], μ[end,i], C, idx, 𝒜)
            τ += (ΔP*𝓌[end])*βₙ*𝓂[k]
            τ = max(τ, τmin)
            #the end is now the beginning
            β₁ = βₙ
            #transmittance through layer
            t = exp(-τ)
            #CIM planck emission approximation
            B₁ = planck(ν, T[1,i])
            B₂ = planck(ν, T[end,i])
            Bₑ = layerplanck(B₁, B₂, τ, t)
            #update irradiance
            I = I*t + Bₑ
        end
        F += 𝒲[k]*I
    end
    return F
end

#-------------------------------------------------------------------------------
# core function for whole atmosphere upward and downward monochromatic fluxes

function 𝒹monoflux!(M⁺::AbstractVector{Q}, #downward monochromatic fluxes [W/m^2/cm^-1]
                    M⁻::AbstractVector{Q}, #upward monochromatic fluxes [W/m^2/cm^-1]
                    τ::AbstractVector{Q}, #optical depth of each layer [-] length(P) - 1
                    P::AbstractVector{<:Real}, #pressure coordinates of output
                    B::AbstractVector{<:Real}, #pre-evaluated Planck at P coordinates
                    ν::Real,
                    fS::R, #incoming stellar radiation fS(ν) [W/m^2]
                    fa::S, #surface albedo fa(ν)
                    θₛ::Real, #stellar radiation angle, corresponds to cos(θ) = 2/3
                    nstream::Int #number of streams to integrate in both directions
                    )::Nothing where {Q,R,S}
    #setup
    np = length(P)
    #number of layers
    L = np - 1
    #check sizes
    @assert np == length(M⁻) == length(M⁺)
    @assert L == length(τ)
    #gauss nodes for integrating streams
    𝓂, 𝒲 = streamnodes(nstream)
    #cosine of stellar zenith angle
    c = cos(θₛ)
    #avoid possible interpolator race conditions
    fS, fa = copyprofiles(fS, fa)
    #wipe any previous values
    @inbounds for i ∈ 1:np
        M⁺[i] = zero(Q)
        M⁻[i] = zero(Q)
    end

    #=============================
    #downward atmospheric emission
    =============================#
    @inbounds for k ∈ 1:nstream
        #initial irradiance
        I = zero(Q)
        𝓂ₖ = 𝓂[k]
        𝒲ₖ = 𝒲[k]
        for i ∈ 1:L
            τᵢ = τ[i]*𝓂ₖ #weight optical depth by 1/cos(θ)
            tᵢ = exp(-τᵢ)
            Bₑ = layerplanck(B[i], B[i+1], τᵢ, tᵢ)
            I = I*tᵢ + Bₑ
            M⁻[i+1] += 𝒲ₖ*I #quadrature weight applied to irradiance
        end
    end

    #==================================
    absorption of downward stellar flux
    ==================================#
    M⁻[1] += c*fS(ν) #initial flux scaled by angle
    Mₛ = M⁻[1]
    @inbounds for i ∈ 1:L
        Mₛ *= exp(-τ[i]/c) #apply layer transmittance
        M⁻[i+1] += Mₛ
    end

    #=================
    outgoing radiation
    =================#
    Iₛ⁺ = M⁻[end]*fa(ν)/π + B[end] #Lambertian reflection and planck emission
    M⁺[end] = Iₛ⁺*π #upward surface flux
    @inbounds for k ∈ 1:nstream
        I = Iₛ⁺
        𝓂ₖ = 𝓂[k]
        𝒲ₖ = 𝒲[k]
        for i ∈ L:-1:1
            τᵢ = τ[i]*𝓂ₖ #weight by 1/cos(θ)
            tᵢ = exp(-τᵢ)
            Bₑ = layerplanck(B[i+1], B[i], τᵢ, tᵢ)
            I = I*tᵢ + Bₑ
            M⁺[i] += 𝒲ₖ*I
        end
    end

    #stop at nothing!
    nothing
end


text:
function _RobustDelague_latex()
    return """
        ```math
        \\left\\{ r  \\; \\middle| \\begin{array}{ll}
        s.t.  \\quad (\\mathbb{E} [r] - \\hat{r}) ' \\Sigma^{-1} (\\mathbb{E} [r] - \\hat{r}) \\leq \\gamma_1 \\\\
        \\quad \\quad \\mathbb{E} [ (r - \\hat{r}) ' (r - \\hat{r}) ] \\leq \\gamma_2 \\Sigma \\\\
        \\end{array}
        \\right\\} \\\\
        ```
        """
end

"""
    RobustDelague <: AbstractMeanVariance

Delague's Ambiguity set:

$(_RobustDelague_latex())

Atributes:
- `predicted_mean::Array{Float64,1}` (latex notation ``\\hat{r}``): Predicted mean of returns.
- `predicted_covariance::Array{Float64,2}` (latex notation ``\\Sigma``): Predicted covariance of returns.
- `γ1::Float64`: Mean uncertainty parameter (has to be greater than 0).
- `γ2::Float64`: Covariance uncertainty parameter (has to be greater than 1).
- `utility_coeficients::Array{Float64,1}`: Piece-wise utility coeficients (default 1).
- `utility_intercepts::Array{Float64,1}`: Piece-wise utility intercepts (default 0).
"""
struct RobustDelague <: AbstractPortfolioFormulation # AbstractMeanVariance
    predicted_mean::Array{Float64,1}
    predicted_covariance::Array{Float64,2}
    γ1::Float64
    γ2::Float64
    utility_coeficients::Array{Float64,1}
    utility_intercepts::Array{Float64,1}
    number_of_utility_pieces::Int
    number_of_assets::Int
end

function RobustDelague(;
    predicted_mean::Array{Float64,1},
    predicted_covariance::Array{Float64,2},
    γ1::Float64,
    γ2::Float64,
    utility_coeficients::Array{Float64,1},
    utility_intercepts::Array{Float64,1},
)
    number_of_assets = size(predicted_mean, 1)
    number_of_utility_pieces = size(utility_coeficients, 1)
    if number_of_assets != size(predicted_covariance, 1)
        error("size of predicted mean ($number_of_assets) different to size of predicted covariance ($(size(predicted_covariance, 1)))")
    end
    if number_of_utility_pieces != size(utility_intercepts, 1)
        error("number of utility coeficients ($(number_of_utility_pieces)) different to number of utility intercepts ($(size(utility_intercepts,1)))")
    end

    @assert γ1 >= 0
    @assert γ2 >= 1
    @assert issymmetric(predicted_covariance)
    return RobustDelague(
        predicted_mean, predicted_covariance, γ1, γ2, utility_coeficients,
        utility_intercepts, number_of_utility_pieces, number_of_assets
    )
end

"""
    po_max_utility_return(formulation::RobustDelague)

Maximize expected utility of portfolio return under distribution uncertainty defined by
Delague's ambiguity set ([`RobustDelague`](@ref)).
"""
function po_max_utility_return(formulation::RobustDelague;
    current_wealth::Real = 1.0,
    model::JuMP.Model = base_model(formulation.number_of_assets; current_wealth=current_wealth),
    w = model[:w],
)
    # parameters
    r̄ = formulation.predicted_mean
    Σ = formulation.predicted_covariance
    numA = formulation.number_of_assets
    a = formulation.utility_coeficients
    b = formulation.utility_intercepts
    γ1 = formulation.γ1
    γ2 = formulation.γ2
    K = formulation.number_of_utility_pieces
    # dual variables
    @variable(model, P[i=1:numA, j=1:numA])
    @variable(model, p[i=1:numA])
    @variable(model, s)
    @variable(model, Q[i=1:numA, j=1:numA])
    @variable(model, q[i=1:numA])
    @variable(model, r)
    # constraints: from duality
    @constraint(model, p .== -q / 2 - Q * r̄)
    @SDconstraint(model, [[P p]; [p' s]] >= 0)
    for k in 1:K
        @SDconstraint(
            model, [[Q (q / 2 + a[k] * w / 2)]; [(q' / 2 + a[k] * w' / 2) (r + b[k])]] >= 0
        )
    end

    # objective
    @objective(
        model,
        Min, # Min because this is the dual problem
        γ2 * dot(Σ, Q) - sum(r̄'Q * r̄) + r + dot(Σ, P) - 2 * dot(r̄, p) + γ1 * s
    )

    return model
end


text:
@doc """
A simple bisection solver. Included for backwards compatability with previous scripts.
f is a function
a is one starting point
b is another starting point
f(a) and f(b) must have opposite signs.
n is optional and is the number of iterations
"""->function mybisect(f,a,b,n=50)
## Introduction
# Code from Introduction to Numerical Methods and Matlab Programming for
# Engineers( P.19)
# T. Young & M.J. Mohlenkamp
# Ohio University

# function [x e] = mybisect (f,a,b,n)
# Does n iterations of the bisection method for a function f
# Inputs : f -- a function
# a,b -- left and right edges of the interval
# n -- the number of bisections to do.
# Outputs : x -- the estimated solution of f(x) = 0
# e -- an upper bound on the error
# evaluate at the ends and make sure there is a sign change
## Code
c = f(a)
d = f(b)
if (c*d > 0.0)
    error("Function has same sign at both endpoints." )
end

for i in 1:n
    # find the middle and evaluate there
    x = (a + b )/2;
    y = f(x);

    if (y == 0.0) # solved the equation exactly
        e = 0;
        break # jumps out of the for loop
    end
    # decide which half to keep , so that the signs at the ends differ
    if (c*y < 0)
        b=x;
    else
        a=x;
    end
end
# set the best estimate for x and the error bound
x = (a + b )/2;
e = (b-a )/2;

return x
end


text:
using SweepOperator, LinearAlgebra, Test

# setup
n, p = 1000, 10
x = randn(n, p)
xtx = x'x

@testset "Sweep One By One" begin
    A = copy(xtx)
    B = copy(xtx)
    for j in 1:p
        sweep!(A, j)
        sweep!(A, j, true)
    end
    @test UpperTriangular(A) ≈ UpperTriangular(B)

    A = copy(xtx)
    B = copy(xtx)
    for j in 1:p
        sweep!(A, j)
    end
    for j in 1:p
        sweep!(A, j, true)
    end
    @test UpperTriangular(A) ≈ UpperTriangular(B)
end

@testset "Sweep All" begin
    A = copy(xtx)
    B = copy(xtx)
    sweep!(A, 1:p)
    sweep!(A, 1:p, true)
    @test A ≈ B
end

@testset "UpperTriangular" begin
    A = UpperTriangular(copy(xtx))
    B = UpperTriangular(copy(xtx))
    sweep!(A, 1:p)
    sweep!(A, 1:p, true)
    @test A ≈ B
end

@testset "Hermitian/Symmetric" begin
    A = Hermitian(copy(xtx))
    B = Symmetric(copy(xtx))
    sweep!(A, 1:p)
    sweep!(A, 1:p, true)
    sweep!(B, 1:p)
    sweep!(B, 1:p, true)
    @test A ≈ B ≈ xtx
end

@testset "Linear Regression" begin
    y = x * collect(1.:p) + randn(n)
    xy = [x y]
    xytxy = xy'xy
    sweep!(xytxy, 1:p)
    @test xytxy[1:p, end] ≈ x\y
end


text:
mutable struct MeasStack{G<:Number} # G = GreensEltype

  # ETPCs (equal-time pairing correlations)
  etpc_minus::Array{Float64, 2} # "P_(y,x)", i.e. d-wave
  etpc_plus::Array{Float64, 2} # "P+(y,x)", i.e. s-wave

  # ZFPCs (zero-frequency pairing correlations)
  zfpc_minus::Array{Float64, 2} # "P_(y,x)", i.e. d-wave
  zfpc_plus::Array{Float64, 2} # "P+(y,x)", i.e. s-wave

  # ETCDCs (equal-time charge density correlations)
  etcdc_minus::Array{ComplexF64, 2} # "C_(y,x)", i.e. d-wave
  etcdc_plus::Array{ComplexF64, 2} # "C+(y,x)", i.e. s-wave

  # ZFCDCs (zero-frequency charge density correlations)
  zfcdc_minus::Array{ComplexF64, 2} # "C_(y,x)", i.e. d-wave
  zfcdc_plus::Array{ComplexF64, 2} # "C+(y,x)", i.e. s-wave

  # ZFCCCs (zero-frequency current-current correlations)
  zfccc_xx::Array{ComplexF64, 2} # "Λxx(y,x)"
  zfccc_yy::Array{ComplexF64, 2} # "Λyy(y,x)"

  # TDGF
  Gt0::Vector{Matrix{G}}
  G0t::Vector{Matrix{G}}
  BT0Inv_u_stack::Vector{Matrix{G}}
  BT0Inv_d_stack::Vector{Vector{Float64}}
  BT0Inv_t_stack::Vector{Matrix{G}}
  BBetaT_u_stack::Vector{Matrix{G}}
  BBetaT_d_stack::Vector{Vector{Float64}}
  BBetaT_t_stack::Vector{Matrix{G}}
  BT0_u_stack::Vector{Matrix{G}}
  BT0_d_stack::Vector{Vector{Float64}}
  BT0_t_stack::Vector{Matrix{G}}
  BBetaTInv_u_stack::Vector{Matrix{G}}
  BBetaTInv_d_stack::Vector{Vector{Float64}}
  BBetaTInv_t_stack::Vector{Matrix{G}}

  MeasStack{G}() where G = new{G}()
end

mutable struct Stack{G<:Number} # G = GreensEltype
  n_elements::Int
  ranges::Array{UnitRange, 1}
  current_slice::Int # running internally over 0:p.slices+1, where 0 and p.slices+1 are artifcial to prepare next sweep direction.
  direction::Int

  eye_flv::Matrix{Float64} # UPGRADE: Replace by dynamically sized I everywhere
  eye_full::Matrix{Float64} # UPGRADE: Replace by dynamically sized I everywhere
  ones_vec::Vector{Float64}

  greens::Matrix{G}
  log_det::Float64 # contains logdet of greens_{p.slices+1} === greens_1
                            # after we calculated a fresh greens in propagate()

  C::Vector{G}
  S::Vector{G}
  R::Vector{G}
  eV::SparseMatrixCSC{G,Int64}

  meas::MeasStack{G}
  Bl::Matrix{G}

  # unsure about those -> always allocate them
  curr_U::Matrix{G}
  D::Vector{Float64}
  U::Matrix{G}
  T::Matrix{G}
  d::Vector{Float64}
  tmp::Matrix{G}
  tmp2::Matrix{G}

  # global update
  gb_u_stack::Array{G, 3}
  gb_d_stack::Matrix{Float64}
  gb_t_stack::Array{G, 3}

  gb_greens::Matrix{G}
  gb_log_det::Float64

  gb_hsfield::Array{Float64, 3}

  # etgreens stack
  u_stack::Array{G, 3}
  d_stack::Matrix{Float64}
  t_stack::Array{G, 3}

  # calc detratio
  delta_i::Matrix{G}
  M::Matrix{G}
  eVop1::Matrix{G}
  eVop2::Matrix{G}
  eVop1eVop2::Matrix{G}
  Mtmp::Matrix{G}
  Mtmp2::Matrix{G}

  # calc greens
  Ul::Matrix{G}
  Ur::Matrix{G}
  Dl::Vector{Float64}
  Dr::Vector{Float64}
  Tl::Matrix{G}
  Tr::Matrix{G}

  # update greens
  A::Matrix{G}
  B::Matrix{G}
  AB::Matrix{G}

  # propagate
  greens_temp::Matrix{G}

  Stack{G}() where G = new{G}()
end

@def stackshortcuts begin
  s = mc.s
  N = mc.l.sites
  flv = mc.p.flv
  safe_mult = mc.p.safe_mult
  G = geltype(mc)
end

function allocate_global_update!(mc)
  @stackshortcuts
  # Global update backup
  s.gb_u_stack = zero(s.u_stack)
  s.gb_d_stack = zero(s.d_stack)
  s.gb_t_stack = zero(s.t_stack)
  s.gb_greens = zero(s.greens)
  s.gb_log_det = 0.
  s.gb_hsfield = zero(mc.p.hsfield)
end

function allocate_etgreens_stack!(mc)
  @stackshortcuts
  s.u_stack = zeros(G, flv*N, flv*N, s.n_elements)
  s.d_stack = zeros(Float64, flv*N, s.n_elements)
  s.t_stack = zeros(G, flv*N, flv*N, s.n_elements)
end

function allocate_calc_detratio!(mc)
  @stackshortcuts
  ## calc_detratio
  s.M = zeros(G, flv, flv)
  s.Mtmp = s.eye_flv - s.greens[1:N:end,1:N:end]
  s.delta_i = zeros(G, size(s.eye_flv))
  s.Mtmp2 = zeros(G, size(s.eye_flv))
  s.eVop1eVop2 = zeros(G, size(s.eye_flv))
  s.eVop1 = zeros(G, flv, flv)
  s.eVop2 = zeros(G, flv, flv)
end

function allocate_calc_greens!(mc)
  @stackshortcuts
  s.Ul = Matrix{G}(I, flv*N, flv*N)
  s.Ur = Matrix{G}(I, flv*N, flv*N)
  s.Tl = Matrix{G}(I, flv*N, flv*N)
  s.Tr = Matrix{G}(I, flv*N, flv*N)
  s.Dl = ones(Float64, flv*N)
  s.Dr = ones(Float64, flv*N)
end

function allocate_update_greens!(mc)
  @stackshortcuts
  ## update_greens
  s.A = s.greens[:,1:N:end]
  s.B = s.greens[1:N:end,:]
  s.AB = s.A * s.B
end

function allocate_propagate!(mc)
  @stackshortcuts
  s.greens_temp = zeros(G, flv*N, flv*N)
end

function _initialize_stack(mc::AbstractDQMC)
  @stackshortcuts
  s.n_elements = convert(Int, mc.p.slices / safe_mult) + 1

  s.ranges = UnitRange[]
  for i in 1:s.n_elements - 1
    push!(s.ranges, 1 + (i - 1) * safe_mult:i * safe_mult)
  end

  s.eye_flv = Matrix{Float64}(I, flv, flv)
  s.eye_full = Matrix{Float64}(I, flv*N,flv*N)
  s.ones_vec = ones(flv*N)

  s.greens = zeros(G, flv*N, flv*N)

  # interaction matrix
  s.C = zeros(G, N)
  s.S = zeros(G, N)
  s.R = zeros(G, N)
  s.eV = spzeros(G, flv*N, flv*N)

  # init empty meas stack
  s.meas = MeasStack{G}()

  # slice matrix
  cbtype(mc) === CBFalse && (s.Bl = zeros(G, flv*N, flv*N))

  # unsure about those
  s.curr_U = zeros(G, flv*N, flv*N) # used in tdgf, add_slice_sequence_ and propagate
  s.D = zeros(Float64, flv*N) # used in calc_greens in mc (below) and in calc greens in fermion meas.
  s.d = zeros(Float64, flv*N) # same as above
  s.tmp = zeros(G, flv*N, flv*N)
  s.tmp2 = zeros(G, flv*N, flv*N)
  s.U = zeros(G, flv*N, flv*N)
  s.T = zeros(G, flv*N, flv*N)
end

function initialize_stack(mc::AbstractDQMC)
  @mytimeit mc.a.to "initialize_stack" begin
    _initialize_stack(mc)

    # allocate for dqmc
    allocate_etgreens_stack!(mc)
    allocate_global_update!(mc)
    allocate_calc_detratio!(mc)
    allocate_calc_greens!(mc)
    allocate_update_greens!(mc)
    allocate_propagate!(mc)

    # allocate for measurements during dqmc
    # allocate_etpc!(mc)

  end #timeit

  nothing
end

function build_stack(mc::AbstractDQMC)
  p = mc.p
  s = mc.s
  a = mc.a

  @mytimeit a.to "build_stack" begin

  s.u_stack[:, :, 1] = s.eye_full
  s.d_stack[:, 1] = s.ones_vec
  s.t_stack[:, :, 1] = s.eye_full

  @inbounds for i in 1:length(s.ranges)
    add_slice_sequence_left(mc, i)
  end

  s.current_slice = p.slices + 1
  s.direction = -1

  end #timeit

  nothing
end

"""
Updates stack[idx+1] based on stack[idx]
"""
function add_slice_sequence_left(mc::AbstractDQMC, idx::Int)
  s = mc.s
  a = mc.a

  copyto!(s.curr_U, s.u_stack[:, :, idx])

  # println("Adding slice seq left $idx = ", s.ranges[idx])
  for slice in s.ranges[idx]
    multiply_B_left!(mc, slice, s.curr_U)
  end

  @views rmul!(s.curr_U, Diagonal(s.d_stack[:, idx]))
  @mytimeit a.to "decompose_udt!" s.u_stack[:, :, idx + 1], T = @views decompose_udt!(s.curr_U, s.d_stack[:, idx + 1])
  @views mul!(s.t_stack[:, :, idx + 1],  T, s.t_stack[:, :, idx])
  nothing
end

"""
Updates stack[idx] based on stack[idx+1]
"""
function add_slice_sequence_right(mc::AbstractDQMC, idx::Int)
  s = mc.s
  a = mc.a

  copyto!(s.curr_U, s.u_stack[:, :, idx + 1])

  for slice in reverse(s.ranges[idx])
    multiply_daggered_B_left!(mc, slice, s.curr_U)
  end

  @views rmul!(s.curr_U, Diagonal(s.d_stack[:, idx + 1]))
  @mytimeit a.to "decompose_udt!" s.u_stack[:, :, idx], T = @views decompose_udt!(s.curr_U, s.d_stack[:, idx])
  @views mul!(s.t_stack[:, :, idx], T, s.t_stack[:, :, idx + 1])
  nothing
end

function wrap_greens!(mc::AbstractDQMC, gf::Matrix, curr_slice::Int=mc.s.current_slice, direction::Int=mc.s.direction)
  if direction == -1
    multiply_B_inv_left!(mc, curr_slice - 1, gf)
    multiply_B_right!(mc, curr_slice - 1, gf)
  else
    multiply_B_left!(mc, curr_slice, gf)
    multiply_B_inv_right!(mc, curr_slice, gf)
  end
  nothing
end

function wrap_greens(mc::AbstractDQMC, gf::Matrix,slice::Int=mc.s.current_slice,direction::Int=mc.s.direction)
  temp = copy(gf)
  wrap_greens!(mc, temp, slice, direction)
  return temp
end

"""
Calculates G(slice) using s.Ur,s.Dr,s.Tr=B(slice)' ... B(M)' and s.Ul,s.Dl,s.Tl=B(slice-1) ... B(1)
"""
function calculate_greens(mc::AbstractDQMC)
  # TODO: Use inv_one_plus_two_udts from linalg.jl here
  s = mc.s
  tmp = mc.s.tmp
  tmp2 = mc.s.tmp2
  a = mc.a
  @mytimeit a.to "calculate_greens" begin

  mul!(tmp, s.Tl, adjoint(s.Tr))
  rmul!(tmp, Diagonal(s.Dr))
  lmul!(Diagonal(s.Dl), tmp)
  @mytimeit a.to "decompose_udt!" s.U, s.T = decompose_udt!(tmp, s.D)

  mul!(tmp, s.Ul, s.U)
  s.U .= tmp
  mul!(tmp2, s.T, adjoint(s.Ur))
  myrdiv!(tmp, s.U', tmp2) # or tmp = s.U' / tmp2 or mul!(tmp, adjoint(s.U), inv(tmp2))
  tmp[diagind(tmp)] .+= s.D
  @mytimeit a.to "decompose_udt!" u, t = decompose_udt!(tmp, s.d)

  mul!(tmp, t, tmp2)
  s.T = inv(tmp)
  mul!(tmp, s.U, u)
  s.U = adjoint(tmp)
  s.d .= 1 ./ s.d

  copyto!(tmp2, s.U)
  lmul!(Diagonal(s.d), tmp2)
  mul!(s.greens, s.T, tmp2)
  end #timeit
  nothing
end

"""
Only reasonable immediately after calculate_greens() because it depends on s.U, s.d and s.T!
"""
function calculate_logdet(mc::AbstractDQMC)
  s = mc.s

  if mc.p.opdim == 1
    s.log_det = real(log(complex(det(s.U))) + sum(log.(s.d)) + log(complex(det(s.T))))
  else
    s.log_det = real(logdet(s.U) + sum(log.(s.d)) + logdet(s.T))
  end
end

################################################################################
# Propagation
################################################################################
function propagate(mc::AbstractDQMC)
  s = mc.s
  p = mc.p

  if s.direction == 1
    if mod(s.current_slice, p.safe_mult) == 0
      s.current_slice +=1 # slice we are going to
      # println("we are going to $(s.current_slice) with a fresh gf.")
      if s.current_slice == 1
        s.Ur[:, :], s.Dr[:], s.Tr[:, :] = s.u_stack[:, :, 1], s.d_stack[:, 1], s.t_stack[:, :, 1]
        s.u_stack[:, :, 1] = s.eye_full
        s.d_stack[:, 1] = s.ones_vec
        s.t_stack[:, :, 1] = s.eye_full
        s.Ul[:,:], s.Dl[:], s.Tl[:,:] = s.u_stack[:, :, 1], s.d_stack[:, 1], s.t_stack[:, :, 1]

        calculate_greens(mc) # greens_1 ( === greens_{m+1} )
        calculate_logdet(mc)

      elseif 1 < s.current_slice <= p.slices
        idx = Int((s.current_slice - 1)/p.safe_mult)

        s.Ur[:, :], s.Dr[:], s.Tr[:, :] = s.u_stack[:, :, idx+1], s.d_stack[:, idx+1], s.t_stack[:, :, idx+1]
        add_slice_sequence_left(mc, idx)
        s.Ul[:,:], s.Dl[:], s.Tl[:,:] = s.u_stack[:, :, idx+1], s.d_stack[:, idx+1], s.t_stack[:, :, idx+1]

        if p.all_checks
          copyto!(s.greens_temp, s.greens)
        end

        wrap_greens!(mc, s.greens_temp, s.current_slice - 1, 1)

        calculate_greens(mc) # greens_{slice we are propagating to}

        if p.all_checks
          diff = maximum(absdiff(s.greens_temp, s.greens))
          if diff > 1e-7
            @printf("->%d \t+1 Propagation instability\t %.4f\n", s.current_slice, diff)
          end
        end

      else # we are going to p.slices+1
        idx = s.n_elements - 1
        add_slice_sequence_left(mc, idx)
        s.direction = -1
        s.current_slice = p.slices+1 # redundant
        propagate(mc)
      end

    else
      # Wrapping
      wrap_greens!(mc, s.greens, s.current_slice, 1)

      s.current_slice += 1
    end

  else # s.direction == -1
    if mod(s.current_slice-1, p.safe_mult) == 0
      s.current_slice -= 1 # slice we are going to
      # println("we are going to $(s.current_slice) with a fresh+wrapped gf.")
      if s.current_slice == p.slices
        s.Ul[:, :], s.Dl[:], s.Tl[:, :] = s.u_stack[:, :, end], s.d_stack[:, end], s.t_stack[:, :, end]
        s.u_stack[:, :, end] = s.eye_full
        s.d_stack[:, end] = s.ones_vec
        s.t_stack[:, :, end] = s.eye_full
        s.Ur[:,:], s.Dr[:], s.Tr[:,:] = s.u_stack[:, :, end], s.d_stack[:, end], s.t_stack[:, :, end]

        calculate_greens(mc) # greens_{p.slices+1} === greens_1
        calculate_logdet(mc) # calculate logdet for potential global update

        # wrap to greens_{p.slices}
        wrap_greens!(mc, s.greens, s.current_slice + 1, -1)

      elseif 0 < s.current_slice < p.slices
        idx = Int(s.current_slice / p.safe_mult) + 1
        s.Ul[:, :], s.Dl[:], s.Tl[:, :] = s.u_stack[:, :, idx], s.d_stack[:, idx], s.t_stack[:, :, idx]
        add_slice_sequence_right(mc, idx)
        s.Ur[:,:], s.Dr[:], s.Tr[:,:] = s.u_stack[:, :, idx], s.d_stack[:, idx], s.t_stack[:, :, idx]

        if p.all_checks
          copyto!(s.greens_temp, s.greens)
        end

        calculate_greens(mc)

        if p.all_checks
          diff = maximum(absdiff(s.greens_temp, s.greens))
          if diff > 1e-7
            @printf("->%d \t-1 Propagation instability\t %.4f\n", s.current_slice, diff)
          end
        end

        wrap_greens!(mc, s.greens, s.current_slice + 1, -1)

      else # we are going to 0
        idx = 1
        add_slice_sequence_right(mc, idx)
        s.direction = 1
        s.current_slice = 0 # redundant
        propagate(mc)
      end

    else
      # Wrapping
      wrap_greens!(mc, s.greens, s.current_slice, -1)
      s.current_slice -= 1
    end
  end
  nothing
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
using ShowSet, Distributed, Factorions
@info "Setting up for parallel execution with $(Threads.nthreads()) threads"

addprocs(Threads.nthreads())
@everywhere using Factorions


text:
function SuperNode_PPLouvain(H::hypergraph,Ω::IntensityFunction;α,clusterpenalty=0,kmax=maximum(keys(H.E)),maxits::Int64=100,bigInt::Bool=true,verbose=true,scan_order="random",Z0 = collect(1:length(H.D)))
    randflag = !(scan_order == "lexical")
    cut_weights, vol_weights, e2n, n2e,w,d,elen = AON_Inputs(H,Ω.ω,α,kmax)
    Zset = SuperNode_PPLouvain(n2e,e2n,w,d,elen,cut_weights,vol_weights,kmax,randflag,maxits,verbose,Z0,clusterpenalty);
    Z = Zset[:,end];
    return Z
end

# function SuperNode_PPLouvain(H::hypergraph,Ω::IntensityFunction,kmax::Int64 = maximum(keys(H.E)),maxits::Int64=100,bigInt::Bool=true;α,verbose=true,scan_order="random", Z0 = collect(1:length(H.D)))
#    Z = SuperNode_PPLouvain(H,Ω;α=α,kmax=kmax,maxits=maxits,bigInt=bigInt,verbose=verbose,scan_order=scan_order,Z0 = Z0)
#    return Z
# end

function SuperNode_PPLouvain(node2edges::Vector{Vector{Int64}},
    edge2nodes::Vector{Vector{Int64}},w::Vector{Float64},
    d::Vector{Float64},elen::Vector{Int64},
    alp::Vector{Float64},bet::Vector{Float64},
    kmax::Int64,randflag::Bool=false,maxits::Int64=100,verbose::Bool=true,
    Zwarm::Vector{Int64}=Vector{Int64}(),clusterpenalty = 0)

    """
    Supernode version of Hypergraph Louvain for planted partition model
        (all-or-nothing cut penalties).

        * node2edges[i] = list of edge labels that node i is in
        * edge2nodes[e] = list of node labels that edge e contains
        * w[e] = weight of hyperedge e
        * d[i] = degree or weight of node i
        * elen[e] = number of nodes in hyperedge e
                     note that this may be > length(edge2nodes[e])
                     if some of the nodes are supernodes
        * alp = scaling values of for hyperedge cuts
        * bet = scaling values for volumes
        * kmax = maximum hyperedge size
        * Zwarm = warm start clustering (optional)
        * randflag = whether or not to permute node order
        * maxits = Maximum # of greedy passes over the node set
        * clusterpenalty = we include term clusterpenalty*log(k) in objective
            so if clusterpenalty > 0, there is incentive to form fewer clusters
    """

    n = length(d)
    # Step 1: greedy moves until no more improvement
    Z, improved = ANHL_Step(node2edges,edge2nodes,w,d,elen,alp,bet,kmax,
                            randflag,maxits,verbose,Zwarm,clusterpenalty)

    # Store all clusterings found
    if improved
        Zs = Z
        Z_old = copy(Z)
    else
        Zs = Z
    end

    # As long as something is still improving each time you call step 1, keep going
    while improved

        # Step 2: Collapse the clustering into supernodes
        @assert(minimum(Z_old) > 0)
        e2n = deepcopy(edge2nodes)
        Change_id_H!(e2n,Z_old)
        # Merge hyperedges that are the same and remove size 1 hyperedges
        e2n, wSuper, elenSuper =  Condense_H(e2n,w,elen)
        He2n = elist2incidence(e2n,maximum(Z_old))
        n2e = incidence2elist(He2n,true)
        dSuper = condense_d(d,Z_old)
        ########
        # Step 1: Go back to greedy local moves, this time on the reduced hypergraph
        Zsuper, improved = ANHL_Step(n2e,e2n,wSuper,dSuper,elenSuper,alp,bet,
                                kmax,randflag,maxits,verbose,Vector{Int64}(),clusterpenalty)
        N = length(Zsuper)    # N = number of supernodes = number clusters from last round

        @assert(minimum(Zsuper)>0)
        @assert(length(Zsuper)==length(dSuper))
        # println("$N, $improved")
        # Undo the procedure that merged nodes into super nodes, so that you get a clustering vector of length n
        # Do this only if the last call to Step 1 led to at least one greedy move.
        if improved
            # Extract what that new clustering is
            Z_new = zeros(Int64,n)

            # For each supernode, place all original node IDs that make it up
            # into a cluster of the supernodes label
            for i = 1:N

                # Get the cluster that supernode i is in
                SuperI_Cluster = Zsuper[i]

                # Get individual node ID that are in supernode i.
                SuperI_nodes = findall(x->x==i,Z_old)
                Z_new[SuperI_nodes] .= SuperI_Cluster
            end
            @assert(minimum(Z_new)>0)
            Zs = [Zs Z_new]
            Z_old = copy(Z_new)
        end
    end

    return Zs
end

function ANHL_Step(node2edges::Vector{Vector{Int64}},edge2nodes::Vector{Vector{Int64}},
    w::Vector{Float64},d::Vector{Float64},elen::Vector{Int64},
    alp::Vector{Float64},bet::Vector{Float64},kmax::Int64,
    randflag::Bool=false,maxits::Int64=100,verbose::Bool=true,Zwarm::Vector{Int64}=Vector{Int64}(),clusterpenalty=0)
    """
    Basic step Louvain algorithm: iterate through nodes and greedily move
    nodes to adjacent clusters. Does not form supernodes and does not recurse.

    NV: Should work also for degenerate hyperedges, though some bugs are possible

    Features:
    * has weight for size of hyperedge rather than computing the weights
    * allows you to warm start with a clustering
    * allows you to randomize node order
    """
    if verbose println("One step of all-or-nothing HyperLouvain") end
    m = length(edge2nodes)
    n = length(node2edges)

    # Randomize node visit order
    if randflag
        node_order = Random.randperm(n)
        # NOTE: may be faster to re-arrange
        #   all the data structures so that we march through
        #   Z in order.
    else
        node_order = 1:n
    end

    if verbose println("") end

    # Warm start clustering
    if length(Zwarm) > 0
        Z = renumber(Zwarm)
    else
        # Default is to put all nodes in their own cluster
        Z = collect(1:n)
    end

	K = maximum(Z)   			# number of clusters
    ClusVol = zeros(K)			# volume of clusters
	ClusSize = zeros(K)			# size of each cluster
	for i = 1:n
		ClusVol[Z[i]] += d[i]
		ClusSize[Z[i]] += 1
	end

    # Store cut penalties for each edge
    cutpenalty = zeros(m)
    for e = 1:m
        k = elen[e]      				# size of the edge
        cutpenalty[e] = alp[k]*w[e] 	# penalty for cutting it
    end

    # Store node neighbors of each node
    Neighbs = NeighborList(node2edges, edge2nodes)

    # Store all the edges for the node, minus that node id
    # This makes computing cut changes faster
    edgelists = Vector{Vector{Vector{Int64}}}()
    for i = 1:n
        iedges = Vector{Vector{Int64}}()    # each node has a list of hyperedges it's in, i.e. a list of lists of nodes.
        for e in node2edges[i]
            edge = edge2nodes[e]
            push!(iedges, setdiff(edge,i))  # don't save i itself
        end
        push!(edgelists,iedges)
    end

    improving = true
    iter = 0
    changemade = false

    # Initialize vols and cuts to track efficient local changes
    r = kmax
    toler = 1e-8

    mainstart = time()
    while improving && iter < maxits

        iter += 1
        if mod(iter,1) == 0
            if verbose println("Louvain Iteration $iter") end
        end
        improving = false

        # visit each node in turn
        for i = node_order

            # Cluster index for node i
            Ci_ind = Z[i]

            # Get the size of the cluster
			clus_size = ClusSize[Ci_ind]

            # Get the indices of i's neighbors--these define clusters we might move to
            Ni = Neighbs[i]

            # Get the neighboring clusters of i
            NC = unique(Z[Ni])

            # The default is to not move the node
            BestC = Z[i]
            BestImprove = 0

            # Set of edges that i is in
            Cv = node2edges[i]

            # Set of hyperedges that node i is in, not including i itself
            Cv_list = edgelists[i]

            # Volume of the set currently
            vS = ClusVol[Ci_ind]
            dv = d[i]

            for j = 1:length(NC)

				# Check how much it would improve to move i to to cluster j
                Cj_ind = NC[j]

                if Cj_ind == Ci_ind
                    continue
                end

                # Change in volume
                vJ = ClusVol[Cj_ind]
                Δvol = 0
                for k = 1:kmax
					# Better if this is smaller
                    Δvol += bet[k]*((vS-dv)^k + (vJ+dv)^k - vS^k - vJ^k)
                end

                # Change in cut
                Δcut = 0
                for eid = 1:length(Cv)
                    # change in cut for edge e when v moves from
                    # cluster Ci to cluster Cj
                    e = Cv[eid]
                    edge_noi = Cv_list[eid]

                    if elen[e] > 1
						we = cutpenalty[e]
                        Δcut += move_cut(i,Z,edge_noi,Ci_ind,Cj_ind,we)
                    end
                end

                # Change in objective due to change in cluster number
                Δclus = 0
                if clusterpenalty > 0 && clus_size == 1
                    Δclus = clusterpenalty*(log(K-1)-log(K))
                end
                change = Δcut + Δvol + Δclus # want this to be negative

                # Check if this is currently the best possible greedy move to make
                # The tolerance helps with making some behavior more stable:
                #   you only move if there's a numerically nonzero reason to move.
                if change < BestImprove - toler
                    BestImprove = change
                    BestC = Cj_ind
                    improving = true
                end
            end

            # Move i to the best new cluster, only if it strictly improves modularity
            if BestImprove < -toler

                # update clustering
                ci_old = Z[i]
                Z[i] = BestC

                # Remove i from its old cluster...
                ClusVol[ci_old] -= dv
				ClusSize[ci_old] -= 1

                # ...and add it to its new cluster
                ClusVol[BestC] += dv
				ClusSize[BestC] += 1
                changemade = true
                improving = true # we have a reason to keep iterating!

                if clus_size == 1
                    # we moved a singleton to another cluster
                    K -= 1
                end

            end
        end
    end
    mainloop = time()-mainstart
    if ~changemade
        improved = false
        if verbose println("No nodes moved clusters") end
    else
        improved = true
    end
    if verbose println("Main loop took $mainloop seconds") end
	Z = renumber(Z)
    return Z, improved

end

"""
Compute the change from moving node v from Ci to Cj
"""
function move_cut(v::Int64,Z::Vector{Int64},erest::Vector{Int64},Ci_ind::Int64,
    Cj_ind::Int64,we::Float64)

    # p = Z[erest]            # set of clusters in e, not counting e
    p1 = Z[erest[1]]
    na = notsame(erest,Z,p1)

    if na
        # if the edge is cut regardless of whether v is moved, return 0
        return 0
    else
        # p1 = p[1]   # this the the cluster all other nodes are in
        if Ci_ind == p1
            # then moving v cuts the edge
            return we
        elseif Cj_ind == p1
            # then moving v to Cj_ind "uncuts" the edge
            return -we
        else
            # moving doesn't matter
            return 0
        end
    end
end

"""
e is the set of nodes in the hyperedge. Z stores their clusters
"""
function notsame(e::Vector{Int64},Z::Vector{Int64},p1::Int64)
    for i = 2:length(e)
        if Z[e[i]] != p1
            return true
        end
    end
    return false
end

function Condense_H(edge2nodes::Vector{Vector{Int64}},w::Vector{Float64})
    """
    Condense Hypergraph.

        A hypergraph stored as an edgelist may contain some edges
        that are identical. This function finds them and condenses
        them instead into a weighted hypergraph edgelist
    """
    Hdict = Dict{Vector{Int64},Float64}()   # maps from set of nodes to a weight

    for e = 1:length(edge2nodes)
        edge = edge2nodes[e]
        sort!(edge)
        we = w[e]
        wt = get(Hdict,edge,0)
        Hdict[edge] = wt+we
    end

    newlist = Vector{Vector{Int64}}()
    for edge in keys(Hdict)
        push!(newlist,edge)
        push!(wnew,Hdict[edge])
    end

    return newlist, wnew
end

function Condense_H(edge2nodes::Vector{Vector{Int64}},w::Vector{Float64},label::Vector{Int64})
    """
    Condense Hypergraph.

        A hypergraph stored as an edgelist may contain some edges
        that are identical. This function finds them and condenses
        them instead into a weighted hypergraph edgelist.

        In this version, each hyperedge is additionally associated with a
        label, and we only merge nodes if they have the same node set
        AND the same label.

        The edge label is mainly meant to store original hyperedge size.
        Since some of the nodes may be supernodes, we need to keep track of
        the original hyperedge size.

        Also, get rid of size 1 hyperedges in the process
    """
    Hdict = Dict{Vector{Int64},Float64}()   # maps from set of nodes to a weight

    @assert(length(label)==length(edge2nodes))

    for e = 1:length(edge2nodes)
        edge = edge2nodes[e]
        lab = label[e]
        sort!(edge)
        ky = [lab; edge]    # the key is the label + edge
        we = w[e]
        wt = get(Hdict,ky,0)
        Hdict[ky] = wt+we
    end

    newlist = Vector{Vector{Int64}}()
    wnew = Vector{Float64}()
    labelsnew = Vector{Int64}() # indices of edges you keep
    for ky in keys(Hdict)
        edge = ky[2:end]
        lab = ky[1]
        if length(edge) > 1
            push!(newlist,edge)
            push!(wnew,Hdict[ky])
            push!(labelsnew,lab)
        end
    end

    return newlist, wnew, labelsnew
end

function Change_id_H!(edge2nodes::Vector{Vector{Int64}},Z::Vector{Int64},uniqueflag::Bool=true)
"""
Given a clustering vector Z, change node IDs in the hypergraph
so that node i becomes node Z[i]. This will potentially lead
to many degenerate hyperedges where the same node shows up multiple times.

This is useful as a super-node step in clustering heuristics like louvain.

If uniqueflag == true, don't store degenerate hyperedges. Instead, just save
one copy of the vector
"""

    for e = 1:length(edge2nodes)
        for i = 1:length(edge2nodes[e])
            edge2nodes[e][i] = Z[edge2nodes[e][i]]
        end
        if uniqueflag
            unique!(edge2nodes[e])
        end
    end

    return edge2nodes
end

function condense_d(d::Vector{Float64},Z::Vector{Int64})
    """
    Condense a degree vector to a weights vector by summing up volumes of
    clusters.
    """

    n = maximum(Z)
    dnew = zeros(n)
    for i = 1:n
        S = findall(x->x==i,Z)
        dnew[i] = sum(d[S])
    end
    return dnew
end

function AON_Inputs(H,ω,α,kmax)
    """
    Given a hyperedge H, ω function, and α which parameterizes ω,
    return all of the data structures you need to run the new
    all-or-nothing louvain algorithm
    """
    cut_weights = zeros(kmax)
    vol_weights = zeros(kmax)
    for k = 2:kmax
        cut_weights[k] = log(ω([1,k],α))-log(ω([0,k],α))
        vol_weights[k] = ω([1,k],α)-ω([0,k],α)
    end

    He2n, edge_weights = hypergraph2incidence(H)
    e2n = incidence2elist(He2n);
    n2e = incidence2elist(SparseArrays.sparse(He2n'));
    m = length(e2n)
    edge_len = zeros(Int64,m)
    for e = 1:m
        edge_len[e] = length(e2n[e])
    end
    deg = vec(sum(He2n,dims = 1))
    # @show H.D-deg
    # @assert(deg == H.D)

    return cut_weights, vol_weights, e2n, n2e, edge_weights,deg,edge_len

end

function Elist_to_Hypergraph(elist::Vector{Vector{Int64}}, maxsize::Int64=25)
    """
    Converts a binary hypergraph incidence matrix to type "hypergraph".
    Assumes edges are properly sorted already
    """
    E = Dict{Integer, Dict}()
    for edge in elist
        # sort!(edge)
        if length(edge) > maxsize; continue; end
        sz = length(edge)
        if !haskey(E, sz)
            E[sz] = Dict{}()
        end
        E[sz][edge] = 1
    end

    D = zeros(Int64, n)
    for (sz, edges) in E
        for (e, _) in edges
            D[e] .+= 1
        end
    end
    N = 1:n
    return hypergraph(N, E, D)
end

function Hmat_to_Hypergraph(H::SparseArrays.SparseMatrixCSC, maxsize::Int64=25)
    """
    Converts a binary hypergraph incidence matrix to type "hypergraph"
    """
    E = Dict{Integer, Dict}()
    elist = incidence2elist(H)
    for edge in elist
        sort!(edge)
        if length(edge) > maxsize; continue; end
        sz = length(edge)
        if !haskey(E, sz)
            E[sz] = Dict{}()
        end
        E[sz][edge] = 1
    end

    D = zeros(Int64, n)
    for (sz, edges) in E
        for (e, _) in edges
            D[e] .+= 1
        end
    end
    N = 1:n
    return hypergraph(N, E, D)
end

max_stars_repo_name: PhilChodrow/HypergraphModularity

max_issues_repo_name: PhilChodrow/HypergraphModularity

max_forks_repo_name: PhilChodrow/HypergraphModularity


text:
@testset "Auxiliary Functions Test" begin
    @testset "check constant columns" begin
		@test_throws Exception PLSRegressor.check_constant_cols([1.0 1;1 2;1 3])
		@test_throws Exception PLSRegressor.check_constant_cols([1.0;1;1][:,:])
		@test_throws Exception PLSRegressor.check_constant_cols([1.0 2 3])
		@test_throws Exception PLSRegressor.check_constant_cols([1.0; 1; 1][:,:])

		@test PLSRegressor.check_constant_cols([1.0 1;2 2;3 3])
		@test PLSRegressor.check_constant_cols([1.0;2;3][:,:])
	end

	@testset "checkparams" begin
	     #@test_logs PLSRegressor.check_params(2,1,"linear")
		 @test_throws Exception PLSRegressor.check_params(-1,2,"linear")
		 @test_throws Exception PLSRegressor.check_params(1,2,"x")

		 @test PLSRegressor.check_params(1,2,"linear")
	end

	@testset "checkdata" begin
		 @test_throws Exception PLSRegressor.check_data(zeros(0,0), 0)
		 @test_throws Exception PLSRegressor.check_data(zeros(1,1), 10)
		 @test PLSRegressor.check_data(zeros(1,1), 1)
	end
end


text:
# todo: wd, botd, topd, alt, ...

wd(w::WordType) = Ducer(length(w) + 1,
                        Set([1]), Set([length(w) + 1]),
                        [(k,k + 1) for k in 1:length(w)],
                        [(c,c) for c in w])

# wd(hi::WordType, lo::WordType) needs padding - - -

botd() = Ducer(0, Set{Int}[], Set{Int}[],
               Pair{Int,Int}[],
               Pair{eltype(ABET),eltype(ABET)}[])

# topd() = Ducer(...) can be done in different ways -
# what to do? two states, first hi,then lo?

topd() = Ducer(2, [1], [2],
               vcat([(1,1) for _ in ABET], [(1,2)],
                    [(2,2) for _ in ABET]),
               vcat([(c,NULL) for c in ABET], [(NULL,NULL)],
                    [(NULL,c) for c in ABET]))


text:
function supervise(strat::Strategy;
                   arg_values::Union{Vector, Nothing}=nothing,
                   px_trade::Symbol=:Open,
                   px_close::Symbol=:Settle,
                   limit::Int)::String
    if isnothing(arg_values)
        arg_values = convert(Vector, strat.indicator.paramset.arg_defaults)
    end
    if isempty(strat.backtest.trades)
        all_trades = generate_trades(strat, arg_values=arg_values)
    else
        all_trades = strat.backtest.trades
    end

    for asset in strat.universe.assets
        trades = all_trades[asset].values
        N = size(trades, 1)

        cnt = 1
        action = false
        signal = "hold"
        while(cnt < limit)
            for (i,rule) in enumerate(strat.rules)
                if trades[N-cnt-1,i] != 0
                    signal = string(rule.action)
                    action = true; break
                end
            end
            if action break end
            cnt+=1
        end
        return signal
    end
end

max_stars_repo_name: lithiferous/Strathreads

max_issues_repo_name: lithiferous/Strathreads

max_forks_repo_name: lithiferous/Strathreads


text:
#=
This is Julia implementation of MCTS by Tejan Karmali.
Reference: Python version of the implementation of MCTS by
(https://github.com/tensorflow/minigo)
=#

using DataStructures: DefaultDict
using Distributions: Dirichlet

# Exploration constant balancing priors vs. value net output
c_puct = 0.96
# How much to weight the priors vs. dirichlet noise when mixing
dirichlet_noise_weight = 0.25

struct MCTSRules
  max_game_length::Int
  # Concentrated-ness of the noise being injected into priors
  dirichlet_noise_alpha::Float32

  function MCTSRules(env)
    max_game_length = (env.N ^ 2 * 7) ÷ 5
    dirichlet_noise_alpha = 0.03 * env.max_action_space / env.action_space
    new(max_game_length, dirichlet_noise_alpha)
  end
end

mutable struct DummyNode
  #= A fake node of a MCTS search tree.

  This node is intended to be a placeholder for the root node, which would
  otherwise have no parent node. If all nodes have parents, code becomes
  simpler. =#
  parent::Void
  child_N::DefaultDict{Any, Float32, Float32}
  child_W::DefaultDict{Any, Float32, Float32}
  function DummyNode()
    new(nothing, DefaultDict{Any, Float32}(0.0f0), DefaultDict{Any, Float32}(0.0f0))
  end
end

mutable struct MCTSNode
  #= A node of a MCTS search tree.

  A node knows how to compute the action scores of all of its children,
  so that a decision can be made about which move to explore next. Upon
  selecting a move, the children dictionary is updated with a new node.

  position: A go.Position instance
  fmove: A move (coordinate) that led to this position, a flattened coord
          (raw number ∈ [1-N^2], with nothing a pass)
  parent: A parent MCTSNode.
  =#

  parent
  fmove
  position :: Position
  is_expanded::Bool
  losses_applied::Int
  child_N::Array{Float32, 1}
  child_W::Array{Float32, 1}
  original_prior::Array{Float32, 1}
  child_prior::Array{Float32, 1}
  children::Dict
  mcts_rules::MCTSRules

  function MCTSNode(position :: T, fmove = nothing, parent = nothing) where T <: Position
    if parent == nothing
      parent = DummyNode()
    end
    is_expanded = false
    losses_applied = 0  # number of virtual losses on this node
    child_N = zeros(Float32, position.env.action_space)
    child_W = zeros(Float32, position.env.action_space)
    # save a copy of the original prior before it gets mutated by d-noise.
    original_prior = zeros(Float32, position.env.action_space)
    child_prior = zeros(Float32, position.env.action_space)
    children = Dict()  # map of flattened moves to resulting MCTSNode
    mcts_rules = MCTSRules(position.env)
    new(parent, fmove, position, is_expanded, losses_applied,
    child_N, child_W, original_prior, child_prior, children, mcts_rules)
  end
end

legal_moves(x::MCTSNode) = all_legal_moves(x.position)

child_action_score(x::MCTSNode) = child_Q(x) * x.position.to_play .+
                                            child_U(x) - 1000 * (1 - legal_moves(x))

child_Q(x::MCTSNode) = x.child_W ./ (1 + x.child_N)

child_U(x::MCTSNode) = (c_puct * √(1 + N(x)) *
                x.child_prior ./ (1 + x.child_N))

Q(x::MCTSNode) = W(x) / (1 + N(x))

get_N(x::MCTSNode) = x.parent.child_N[x.fmove]
set_N!(x::MCTSNode, value) = x.parent.child_N[x.fmove] = value
N(x::MCTSNode) = get_N(x)

get_W(x::MCTSNode) = x.parent.child_W[x.fmove]
set_W!(x::MCTSNode, value) = x.parent.child_W[x.fmove] = value
W(x::MCTSNode) = get_W(x)

# Return value of position, from perspective of player to play
Q_perspective(x::MCTSNode) = Q(x) * x.position.to_play

function select_leaf(mcts_node::MCTSNode)
  current = mcts_node
  board_size = mcts_node.position.env.N
  pass_move = board_size * board_size + 1
  while true
    current_new_N = N(current) + 1
    set_N!(current, current_new_N)
    # if a node has never been evaluated, we have no basis to select a child.
    if !current.is_expanded
      break
    end
    # HACK: if last move was a pass, always investigate double-pass first
    # to avoid situations where we auto-lose by passing too early.
    if (length(current.position.recent) != 0
      && current.position.recent[end].move == nothing
      && current.child_N[pass_move] == 0)
      current = maybe_add_child!(current, pass_move)
      continue
    end
    cas = child_action_score(current)
    #max_score = maximum(cas)
    #possible_moves = find(x -> x == max_score, cas)
    #best_move = sample(possible_moves)
    best_move = findmax(cas)[2]
    current = maybe_add_child!(current, best_move)
  end
  return current
end

function maybe_add_child!(mcts_node::MCTSNode, fcoord)
  # Adds child node for fcoord if it doesn't already exist, and returns it
  if fcoord ∉ keys(mcts_node.children)
    new_position = play_move!(mcts_node.position, from_flat(fcoord, mcts_node.position))
    mcts_node.children[fcoord] = MCTSNode(new_position, fcoord, mcts_node)
  	end
  return mcts_node.children[fcoord]
end

function add_virtual_loss!(mcts_node::MCTSNode, up_to::MCTSNode)
  #= Propagate a virtual loss up to the root node.

    	Args:
        up_to: The node to propagate until. (Keep track of this! You'll
            need it to reverse the virtual loss later.)
  =#
  mcts_node.losses_applied += 1
  # This is a "win" for the current node; hence a loss for its parent node
  # who will be deciding whether to investigate this node again.
  loss = mcts_node.position.to_play
  set_W!(mcts_node, W(mcts_node) + loss)
  if (mcts_node.parent == nothing || mcts_node == up_to) return  end
  add_virtual_loss!(mcts_node.parent, up_to)
end

function revert_virtual_loss!(mcts_node::MCTSNode, up_to::MCTSNode)
  mcts_node.losses_applied -= 1
  revert = -mcts_node.position.to_play
	set_W!(mcts_node, W(mcts_node) + revert)
  if (mcts_node.parent == nothing || mcts_node == up_to) return end
  revert_virtual_loss!(mcts_node.parent, up_to)
end

function revert_visits!(mcts_node::MCTSNode, up_to::MCTSNode)
  #= Revert visit increments.

    Sometimes, repeated calls to select_leaf return the same node.
    This is rare and we're okay with the wasted computation to evaluate
    the position multiple times by the dual_net. But select_leaf has the
    side effect of incrementing visit counts. Since we want the value to
    only count once for the repeatedly selected node, we also have to
    revert the incremented visit counts.
  =#
  set_N!(mcts_node, N(mcts_node) - 1)
  if mcts_node.parent == nothing || mcts_node == up_to return end
  revert_visits!(mcts_node.parent, up_to)
end

function incorporate_results!(mcts_node::MCTSNode, move_probs, value, up_to)
  board_size = mcts_node.position.env.N
  @assert size(move_probs) == (board_size ^ 2 + 1,)
  # A finished game should not be going through this code path - should
  # directly call backup_value() on the result of the game.
	@assert !mcts_node.position.done
  if mcts_node.is_expanded
    revert_visits!(mcts_node, up_to)
    return
	end
  mcts_node.is_expanded = true
  mcts_node.original_prior .= mcts_node.child_prior .= move_probs
  # initialize child Q as current node's value, to prevent dynamics where
  # if B is winning, then B will only ever explore 1 move, because the Q
  # estimation will be so much larger than the 0 of the other moves.
  #
  # Conversely, if W is winning, then B will explore all go.N² + 1 moves before
  # continuing to explore the most favorable move. This is a waste of search.
  #
  # The value seeded here acts as a prior, and gets averaged into Q calculations
  mcts_node.child_W .= ones(Float32, mcts_node.position.env.action_space) * value
  backup_value!(mcts_node, value, up_to)
end

function backup_value!(mcts_node::MCTSNode, value, up_to::MCTSNode)
  #= Propagates a value estimation up to the root node.

  Args:
    value: the value to be propagated (1 = black wins, -1 = white wins)
    up_to: the node to propagate until.
	=#
  set_W!(mcts_node, W(mcts_node) + value)
  if mcts_node.parent == nothing || mcts_node == up_to return end
  backup_value!(mcts_node.parent, value, up_to)
end

#= true if the last two moves were Pass or if the position is at a move
  	greater than the max depth.
=#
is_done(mcts_node::MCTSNode) = mcts_node.position.done ||
mcts_node.position.n >= mcts_node.mcts_rules.max_game_length

function inject_noise!(mcts_node::MCTSNode)
  action_space = mcts_node.position.env.action_space
  dirch = rand(Dirichlet(mcts_node.mcts_rules.dirichlet_noise_alpha
                                                * ones(action_space)))
  mcts_node.child_prior .= mcts_node.child_prior * (1 - dirichlet_noise_weight) .+
                      dirch .* dirichlet_noise_weight
end

function children_as_π(mcts_node::MCTSNode, squash::Bool = false)
  #= Returns the child visit counts as a probability distribution, π
  If squash is true, exponentiate the probabilities by a temperature
  slightly larger than unity to encourage diversity in early play and
  hopefully to move away from 3-3s
  =#
  probs = mcts_node.child_N
  if squash
    probs = probs .^ 0.98
  end
  return probs ./ sum(probs)
end

#=
function most_visited_path_nodes(mcts_node::MCTSNode)
  node = mcts_node
  output = Array{MCTSNode, 1}()
  while node.children
    next_kid = findmax(node.child_N)[2]
    node = get(node.children, next_kid, nothing)
    if node == nothing break end
    push!(output, node)
  end
  return output
end

function most_visited_path(mcts_node::MCTSNode)
  node = mcts_node
  output = Array{String, 1}()
  while !isempty(node.children)
    next_kid = findmax(node.child_N)[2]
    node = get(node.children, next_kid, nothing)
    if node == nothing
      push!(output, "GAME END")
      break
    end
    push!(output, go.to_kgs(go.from_flat(node.fmove)) * " ($(N(node))) ==> ")
  end
  push!(output, @sprintf("Q: %.5f\n", Q(node)))
  return join(output)
end

function mvp_gg(mcts_node::MCTSNode)
  # Returns most visited path in go-gui VAR format e.g. 'b r3 w c17...
  node = mcts_node
  output = Array{String, 1}()
  while !isempty(node.children) && maximum(node.child_N) > 1
    next_kid = findmax(node.child_N)[2]
    node = node.children[next_kid]
    push!(output, go.to_kgs(go.from_flat(node.fmove)))
  end
  return join(output, " ")
end

function describe(mcts_node::MCTSNode)
  sort_order = collect(1:go.N * go.N + 1)
  sort!(sort_order, by = i -> (
      mcts_node.child_N[i], child_action_score(mcts_node)[i]), rev = true)
  soft_n = mcts_node.child_N / max(1, sum(mcts_node.child_N))
  prior = mcts_node.child_prior
  p_delta = soft_n - prior
  p_rel = zeros(p_delta)
  mask = prior .!= 0
  p_rel[mask] .= p_delta[mask] ./ prior[mask]

  # Dump out some statistics
  output = Array{String, 1}()
  push!(output, @sprintf("%.4f\n", Q(mcts_node)))
  push!(output, most_visited_path(mcts_node))
  push!(output,
      "move : action    Q     U     P   P-Dir    N  soft-N  p-delta  p-rel")
  for key in sort_order[1:15]
    if mcts_node.child_N[key] == 0 break end
    push!(output, @sprintf("\n%s   : % .3f % .3f %.3f %.3f %.3f %5d %.4f % .5f % .2f",
        go.to_kgs(go.from_flat(key)),
        child_action_score(mcts_node)[key],
        child_Q(mcts_node)[key],
        child_U(mcts_node)[key],
        mcts_node.child_prior[key],
        mcts_node.original_prior[key],
        trunc(Int, mcts_node.child_N[key]),
        soft_n[key],
        p_delta[key],
        p_rel[key]))
  end
  return join(output)
end
=#


text:
###############
# Composition #
###############

"""
    Composed(ts::A)

    ∘(b1::Bijector{N}, b2::Bijector{N})::Composed{<:Tuple}
    composel(ts::Bijector{N}...)::Composed{<:Tuple}
    composer(ts::Bijector{N}...)::Composed{<:Tuple}

where `A` refers to either
- `Tuple{Vararg{<:Bijector{N}}}`: a tuple of bijectors of dimensionality `N`
- `AbstractArray{<:Bijector{N}}`: an array of bijectors of dimensionality `N`

A `Bijector` representing composition of bijectors. `composel` and `composer` results in a
`Composed` for which application occurs from left-to-right and right-to-left, respectively.

Note that all the alternative ways of constructing a `Composed` returns a `Tuple` of bijectors.
This ensures type-stability of implementations of all relating methdos, e.g. `inv`.

If you want to use an `Array` as the container instead you can do

    Composed([b1, b2, ...])

In general this is not advised since you lose type-stability, but there might be cases
where this is desired, e.g. if you have a insanely large number of bijectors to compose.

# Examples
## Simple example
Let's consider a simple example of `Exp`:
```julia-repl
julia> using Bijectors: Exp

julia> b = Exp()
Exp{0}()

julia> b ∘ b
Composed{Tuple{Exp{0},Exp{0}},0}((Exp{0}(), Exp{0}()))

julia> (b ∘ b)(1.0) == exp(exp(1.0))    # evaluation
true

julia> inv(b ∘ b)(exp(exp(1.0))) == 1.0 # inversion
true

julia> logabsdetjac(b ∘ b, 1.0)         # determinant of jacobian
3.718281828459045
```

# Notes
## Order
It's important to note that `∘` does what is expected mathematically, which means that the
bijectors are applied to the input right-to-left, e.g. first applying `b2` and then `b1`:
```julia
(b1 ∘ b2)(x) == b1(b2(x))     # => true
```
But in the `Composed` struct itself, we store the bijectors left-to-right, so that
```julia
cb1 = b1 ∘ b2                  # => Composed.ts == (b2, b1)
cb2 = composel(b2, b1)         # => Composed.ts == (b2, b1)
cb1(x) == cb2(x) == b1(b2(x))  # => true
```

## Structure
`∘` will result in "flatten" the composition structure while `composel` and
`composer` preserve the compositional structure. This is most easily seen by an example:
```julia-repl
julia> b = Exp()
Exp{0}()

julia> cb1 = b ∘ b; cb2 = b ∘ b;

julia> (cb1 ∘ cb2).ts # <= different
(Exp{0}(), Exp{0}(), Exp{0}(), Exp{0}())

julia> (cb1 ∘ cb2).ts isa NTuple{4, Exp{0}}
true

julia> Bijectors.composer(cb1, cb2).ts
(Composed{Tuple{Exp{0},Exp{0}},0}((Exp{0}(), Exp{0}())), Composed{Tuple{Exp{0},Exp{0}},0}((Exp{0}(), Exp{0}())))

julia> Bijectors.composer(cb1, cb2).ts isa Tuple{Composed, Composed}
true
```

"""
struct Composed{A, N} <: Bijector{N}
    ts::A

    Composed(bs::C) where {N, C<:Tuple{Vararg{<:Bijector{N}}}} = new{C, N}(bs)
    Composed(bs::A) where {N, A<:AbstractArray{<:Bijector{N}}} = new{A, N}(bs)
end

isclosedform(b::Composed) = all(isclosedform, b.ts)

"""
    composel(ts::Bijector...)::Composed{<:Tuple}

Constructs `Composed` such that `ts` are applied left-to-right.
"""
composel(ts::Bijector{N}...) where {N} = Composed(ts)

"""
    composer(ts::Bijector...)::Composed{<:Tuple}

Constructs `Composed` such that `ts` are applied right-to-left.
"""
composer(ts::Bijector{N}...) where {N} = Composed(reverse(ts))

# The transformation of `Composed` applies functions left-to-right
# but in mathematics we usually go from right-to-left; this reversal ensures that
# when we use the mathematical composition ∘ we get the expected behavior.
# TODO: change behavior of `transform` of `Composed`?
@generated function ∘(b1::Bijector{N1}, b2::Bijector{N2}) where {N1, N2}
    if N1 == N2
        return :(composel(b2, b1))
    else
        return :(throw(DimensionMismatch("$(typeof(b1)) expects $(N1)-dim but $(typeof(b2)) expects $(N2)-dim")))
    end
end

# type-stable composition rules
∘(b1::Composed{<:Tuple}, b2::Bijector) = composel(b2, b1.ts...)
∘(b1::Bijector, b2::Composed{<:Tuple}) = composel(b2.ts..., b1)
∘(b1::Composed{<:Tuple}, b2::Composed{<:Tuple}) = composel(b2.ts..., b1.ts...)

# type-unstable composition rules
∘(b1::Composed{<:AbstractArray}, b2::Bijector) = Composed(pushfirst!(copy(b1.ts), b2))
∘(b1::Bijector, b2::Composed{<:AbstractArray}) = Composed(push!(copy(b2.ts), b1))
function ∘(b1::Composed{<:AbstractArray}, b2::Composed{<:AbstractArray})
    return Composed(append!(copy(b2.ts), copy(b1.ts)))
end

# if combining type-unstable and type-stable, return type-unstable
function ∘(b1::T1, b2::T2) where {T1<:Composed{<:Tuple}, T2<:Composed{<:AbstractArray}}
    error("Cannot compose compositions of different container-types; ($T1, $T2)")
end
function ∘(b1::T1, b2::T2) where {T1<:Composed{<:AbstractArray}, T2<:Composed{<:Tuple}}
    error("Cannot compose compositions of different container-types; ($T1, $T2)")
end

∘(::Identity{N}, ::Identity{N}) where {N} = Identity{N}()
∘(::Identity{N}, b::Bijector{N}) where {N} = b
∘(b::Bijector{N}, ::Identity{N}) where {N} = b

inv(ct::Composed) = Composed(reverse(map(inv, ct.ts)))

# # TODO: should arrays also be using recursive implementation instead?
function (cb::Composed{<:AbstractArray{<:Bijector}})(x)
    @assert length(cb.ts) > 0
    res = cb.ts[1](x)
    for b ∈ Base.Iterators.drop(cb.ts, 1)
        res = b(res)
    end

    return res
end

@generated function (cb::Composed{T})(x) where {T<:Tuple}
    @assert length(T.parameters) > 0
    expr = :(x)
    for i in 1:length(T.parameters)
        expr = :(cb.ts[$i]($expr))
    end
    return expr
end

function logabsdetjac(cb::Composed, x)
    y, logjac = forward(cb.ts[1], x)
    for i = 2:length(cb.ts)
        res = forward(cb.ts[i], y)
        y = res.rv
        logjac += res.logabsdetjac
    end

    return (rv = y, logabsdetjac = logjac)
end

@generated function logabsdetjac(cb::Composed{T}, x) where {T<:Tuple}
    N = length(T.parameters)

    expr = Expr(:block)
    push!(expr.args, :((y, logjac) = forward(cb.ts[1], x)))

    for i = 2:N - 1
        temp = gensym(:res)
        push!(expr.args, :($temp = forward(cb.ts[$i], y)))
        push!(expr.args, :(y = $temp.rv))
        push!(expr.args, :(logjac += $temp.logabsdetjac))
    end
    # don't need to evaluate the last bijector, only it's `logabsdetjac`
    push!(expr.args, :(logjac += logabsdetjac(cb.ts[$N], y)))

    push!(expr.args, :(return logjac))

    return expr
end

function forward(cb::Composed, x)
    rv, logjac = forward(cb.ts[1], x)

    for t in cb.ts[2:end]
        res = forward(t, rv)
        rv = res.rv
        logjac = res.logabsdetjac + logjac
    end
    return (rv=rv, logabsdetjac=logjac)
end

@generated function forward(cb::Composed{T}, x) where {T<:Tuple}
    expr = Expr(:block)
    push!(expr.args, :((y, logjac) = forward(cb.ts[1], x)))
    for i = 2:length(T.parameters)
        temp = gensym(:temp)
        push!(expr.args, :($temp = forward(cb.ts[$i], y)))
        push!(expr.args, :(y = $temp.rv))
        push!(expr.args, :(logjac += $temp.logabsdetjac))
    end
    push!(expr.args, :(return (rv = y, logabsdetjac = logjac)))

    return expr
end


text:
"""There are two subtypes of `AbstractDomain`:\n
* `FourierDomain{K}`: filtering is implemented as a product in the Fourier
domain. `K` is the dimension of the Fourier transform.\n
* `SpatialDomain{K}`: filtering is implemented as a convolution in the spatial
domain, or as a wavelet lifting scheme when it is possible. `K` is the number of
subscripts over which the convolution operates.
"""
abstract AbstractDomain

immutable FourierDomain{K} <: AbstractDomain
    dim::Val{K}
    function FourierDomain{K}(dim::Val{K})
        @assert (K > 0)
        @assert isa(K, Int)
        new(dim)
    end
end
FourierDomain(dim::Int) = FourierDomain{dim}(Val{dim}())

immutable SpatialDomain{K} <: AbstractDomain
    dim::Val{K}
    function SpatialDomain{K}(dim::Val{K})
        @assert (K > 0)
        @assert isa(K, Int)
        new(dim)
    end
end
SpatialDomain(dim::Int) = SpatialDomain{dim}(Val{dim}())

typealias LineDomains Union{FourierDomain{1},SpatialDomain{1}}
typealias PlaneDomains Union{FourierDomain{2},SpatialDomain{2}}


text:
using Plots

"""
  plot_training(acc, cost)

Plot the `Cost vs Epochs` and `Accuracy vs Epochs` graphs
for both the Training and Validation Set.

Arguments:
1. `acc`: A dictionary with keys `train` and `valid` and stores an
          array corresponding to the accuracy during training.
2. `cost`: A dictionary with keys `train` and `valid` and stores an
           array corresponding to the cost during training.
"""

function plot_training(acc, cost)
  gr()
  acc_train = acc["train"]
  acc_valid = acc["valid"]
  cost_train = acc["train"]
  cost_valid = acc["valid"]
  x = length(acc)

  p1 = plot(x, acc_train, xlabel = "Epochs", ylabel = "Accuracy", title = "Training Accuracy")

  p2 = plot(x, acc_valid, xlabel = "Epochs", ylabel = "Accuracy", title = "Validation Accuracy")

  p3 = plot(x, cost_train, xlabel = "Epochs", ylabel = "Cost", title = "Training Cost")

  p4 = plot(x, cost_valid, xlabel = "Epochs", ylabel = "Cost", title = "Validation Cost")

  plot(p1, p2, p3, p4, layout = (2, 2))
end

"""
  get_count(dataset)

Counts the number of abnormal images and the number of normal
images in the dataset.

Arguments:
1. `dataset`: For the input format see the output of `get_batched_images`.
              Any other input will not give proper results. Ideally pass
              the output from the `get_batched_images` function.
"""

function get_count(dataset)
  local categories = ["train", "valid"]
  total_abnormal_imgs = Dict("train" => 0, "valid" => 0)
  total_normal_imgs = Dict("train" => 0, "valid" => 0)
  for cat in categories
    for tup in dataset[cat]
      for i in tup[2]
        if i == 1
          total_abnormal_imgs[cat] += 1
        else
          total_normal_imgs[cat] += 1
        end
      end
    end
  end
  (total_abnormal_imgs, total_normal_imgs)
end

"""
  confusion_matrix(dataset)

Generates the confusion matrix based on the predicted values of the
model and the ground truth. The `x-axis` corresponds to the `ground truth`
and `y-axis` is the predicted value.

Arguments:
1. `dataset`: It must be an `array of tuples` or `array of arrays`
              where the first element of the `tuple` or `array` must
              be the `data` and the second element must be the `label`.
"""

function confusion_matrix(dataset)
  conf_mat = zeros(Int64, 2, 2)
  Flux.testmode!(model)
  for data in dataset
    pred = model(data[1] |> gpu)
    output = Int.(pred .>= threshold)
    for i in 1:length(output)
      conf_mat[output[i] + 1, Int(data[2][i] + 1)] += 1
    end
  end
  Flux.testmode!(model, false)
  conf_mat
end

"""
  Κ_score(conf_mat)

Return the kappa score from the confusion matrix

Arguments:
1. `conf_mat`: Confusion matrix. The orientation must be same as
               the output of `confusion_matrix` function.
"""

function Κ_score(conf_mat)
  tot_instances = sum(conf_mat)
  pₒ = (conf_mat[1, 1] + conf_mat[2, 2])/tot_instances
  pₑ = (sum(conf_mat[1, :]) * sum(conf_mat[:, 1]) + sum(conf_mat[:, 2]) * sum(conf_mat[2, :]))/(tot_instances ^ 2)
  (pₒ - pₑ)/(1 - pₑ)
end


text:
@plottest begin
    # simulation parameters
    T = 1000

    # neuron parameters
    vᵣ = 0
    τᵣ = 1.0
    vth = 1.0

    # create population
    # neuron 1, 2 excite neuron 3
    # neuron 3 inhibits neuron 1, 2
    weights = [ 0  0  1;
                0  0  1;
            -5 -5  0]
    pop = Population(weights; cell = () -> LIF(τᵣ, vᵣ),
                            synapse = Synapse.Alpha,
                            threshold = () -> Threshold.Ideal(vth))

    # create input currents
    low = ConstantRate(0.1)
    high = ConstantRate(0.99)
    switch(t; dt = 1) = (t < Int(T/2)) ? low(t; dt = dt) : high(t; dt = dt)
    n1synapse = QueuedSynapse(Synapse.Alpha())
    n2synapse = QueuedSynapse(Synapse.Alpha())
    excite!(n1synapse, filter(x -> x != 0, [low(t) for t = 1:T]))
    excite!(n2synapse, filter(x -> x != 0, [switch(t) for t = 1:T]))

    # simulate
    voltages = Dict([(i, Float64[]) for i in 1:3])
    cb = () -> begin
        for id in 1:size(pop)
            push!(voltages[id], getvoltage(pop[id]))
        end
    end
    outputs = simulate!(pop, T; cb = cb, inputs = [n1synapse, n2synapse, (t; dt) -> 0])

    rasterplot(outputs, label = ["Input 1", "Input 2", "Inhibitor"])
    title!("Raster Plot")
    xlabel!("Time (sec)")
end joinpath(datadir, "Winner-Take-All-Test.png") !isci


text:
module OpenTelemetry

using Reexport

@reexport using OpenTelemetryAPI
@reexport using OpenTelemetrySDK

end


text:
module RLAlgorithms
    include("envs/RLEnvironments.jl")
    include("helpers/RLHelpers.jl")
    using Reexport
    @reexport using .RLHelpers
    @reexport using .RLEnvironments
    using StatsBase: sample, Weights

    export episode!, episodes!

    function episodes!(f, memory, env, steps, maxt)
        j = 0
        while j < steps
            j += episode!(f, memory, env, min(maxt, steps-j))
        end
    end

    function episode!(f, memory::ReplayMemory{S, A, R},
            env::Environment{S, B, R}, maxt) where {S, A, B, R}
        reset!(env)
        done = false
        t = 0
        while t < maxt && !done
            s = copy(env.current)
            a = f(s)
            _, r, done = step!(env, a)
            append!(memory, s, a, r, done)
            t += 1
        end
        return t
    end

    function episode!(f, episode::PGEpisode{S, A, R},
            env::Environment{S, B, R}, maxt) where {S, A, B, R}
        reset!(env)
        done = false
        for _ in 1:maxt-1
            s = copy(env.current)
            p, v = f(s)
            a = sample(Weights(p))
            _, r, done = step!(env, a)
            append!(episode, s, a, r, v)
            done && break
        end
        if !done
            s = env.current
            p, v = f(s)
            append!(episode, s, argmax(p), v, v)
        end
        nothing
    end

    function episode!(f, memory::PGReplayMemory{S, A, R},
            env::Environment{S, B, R}, maxt) where {S, A, B, R}
        t = 0
        done = false
        reset!(env)
        while t < maxt && !done
            s = copy(env.current)
            a, loga, v = f(s)
            _, r, done = step!(env, a)
            append!(memory, s, a, loga, r, v)
            t += 1
        end
        lastval = done ? zero(R) : f(env.current)[3]
        finish!(memory, t, lastval)
        return t
    end
end


text:
Base.@deprecate ofstatictype(x, y) (z = oftype(unstatic(x), y); x isa Static ? static(z) : z)

macro staticnumbers(args::Expr, funs1::Expr, funs2::Expr, targets::Expr=:(()))
    @warn "The @staticnumbers macro has been renamed @generate_static_methods." # Not the right way to do this.
    args.head == :tuple || error("Expected a Tuple of numbers")
    funs1.head == :tuple || error("Expected a Tuple of 1-arg functions")
    funs2.head == :tuple || error("Expected a Tuple of 2-arg functions")
    targets.head == :tuple || error("Expected a Tuple of target numbers")
    return esc(Expr(:block,vcat(
       genstaticmethods1(funs1.args, args.args, targets.args),
       genstaticmethods2(funs2.args, args.args, targets.args)
       )...))
end


text:
using CSV
using DataFrames
using DelimitedFiles
using Printf
using Statistics

using SLiDE

"""
    make_uniform(df::DataFrame, cols::Array{Symbol,1})
This function is specific to the verify_datastream.jl file. It makes input DataFrames
uniform (same columns; should add support to standardize types)
"""
function make_uniform(df::DataFrame, cols::Array{Symbol,1})
    if size(propertynames(df)) == size(cols) && (length(setdiff(cols, propertynames(df))) > 0)
        df = edit_with(df, Rename.(propertynames(df), cols))[:,cols]
    end
    df = df[:,cols]

    # !!!! Assumptions about value column name and type: (1) values are stored in a column
    # named :value and (2) values are type FLOAT64. This is a common theme throughout SLiDE;
    # may need to address.
    df[!,:value] .= convert_type.(Float64, df[:,:value])

    df = edit_with(df, Drop(:value, 0.0, "=="));
    return unique(sort(df, propertynames(df)[1:end-1]));
end

# ******************************************************************************************

# PATHS TO DATA LOCATIONS:
# - path_side is general: the contents of this directory should have been created by
#     auto_standardize_data.
# - path_bluenote is specific to where Caroline stored WiNDC windc_datastream output.
#     This line must be changed to be consistent with user-specific file storage organizations.
path_slide = joinpath("data", "input")
path_bluenote = joinpath("data","windc_output","1b_stream_windc_base")

y = read_file(joinpath("dev", "verify_datastream.yml"));
lst = y["FileInput"];

df_attn = Dict()

# Iterate through list of datafiles to compare and shared column names.
for inp in lst[1:end-1]
    println("\n",uppercase(inp.f1));

    # For SLiDE and bluenote data sets, read and save the DataFrames.
    dfs = read_file(joinpath(path_slide, inp.f1));
    global dfs = make_uniform(dfs, inp.colnames);

    dfb = read_file(joinpath(path_bluenote, inp.f2));
    global dfb = make_uniform(dfb, inp.colnames);

    # If comparing SEDS data, compare unstandardized units.
    # if occursin("seds", inp.f1)
    #     x = [Rename.([:units_0, :value_0], [:units, :value]);
    #          Replace.([:source_code, :sector_code], "lower", "upper");
    #          Map("parse/units.csv", [:from], [:to], [:units], [:units])];
    #     global dfs = edit_with(dfs, x)
    #     global dfb = edit_with(dfb, x)
    # end

    global df_keys = compare_keys(copy.([dfs, dfb]), [:slide, :bluenote]);
    # size(df_keys,1) == 0 ? println("All keys are consistent.") : show(df_keys);
    println("");

    # Resolve MINOR discrepancy to compare values for heatrate calculations.
    x = Replace(:units, "btu per kWh generated", "btu per kilowatthour");
    global dfs = edit_with(dfs, x);
    global dfb = edit_with(dfb, x);

    # Compare values. Print a warning if there are discrepancies and
    # save the summary DataFrame in a list of DataFrames in need of attention.
    global df = compare_values(copy.([dfs, dfb]), [:slide, :bluenote]);
    size(df,1) > 0 && (df_attn[inp.f1[1:end-4]] = df)
end

# ******************************************************************************************
# cols = [:year, :state, :sgf_code, :units, :value];
# dfs = read_file(joinpath(path_slide, "sgf.csv"));
# dfs[!,:windc_code] .= ismissing.(dfs[:,:windc_code])
# dfs = unique(sort(dropmissing(dfs, :sgf_code), cols));
# display(first(dfs,4));

# dfb = read_file(joinpath(path_bluenote, "sgf_units.csv"));
# dfb = make_uniform(edit_with(dfb, Replace(:sgf_code, "othtax", "OTHTAX")), cols);
# display(first(dfb,4));

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
# for arithmetic etc. see AlgAssAbsOrd/Elem.jl

export trred

parent_type(::Type{AlgAssRelOrdElem{S, T, U}}) where {S, T, U} = AlgAssRelOrd{S, T, U}

parent_type(::AlgAssRelOrdElem{S, T, U}) where {S, T, U} = AlgAssRelOrd{S, T, U}

@doc Markdown.doc"""
    parent(x::AlgAssRelOrdElem) -> AlgAssRelOrd

> Returns the order containing $x$.
"""
@inline parent(x::AlgAssRelOrdElem) = x.parent

################################################################################
#
#  Parent check
#
################################################################################

function check_parent(x::AlgAssRelOrdElem{S, T, U}, y::AlgAssRelOrdElem{S, T, U}) where {S, T, U}
  return parent(x) === parent(y)
end

################################################################################
#
#  Parent object overloading
#
################################################################################

(O::AlgAssRelOrd{S, T, U})(a::AbsAlgAssElem{S}, check::Bool = true) where {S, T, U} = begin
  if check
    (x, y) = _check_elem_in_order(a, O)
    !x && error("Algebra element not in the order")
    return AlgAssRelOrdElem{S, T, U}(O, deepcopy(a), y)
  else
    return AlgAssRelOrdElem{S, T, U}(O, deepcopy(a))
  end
end

(O::AlgAssRelOrd{S, T, U})(a::AbsAlgAssElem{S}, arr::Vector{S}, check::Bool = false) where {S, T, U} = begin
  if check
    (x, y) = _check_elem_in_order(a, O)
    (!x || arr != y) && error("Algebra element not in the order")
    return AlgAssRelOrdElem{S, T, U}(O, deepcopy(a), y)
  else
    return AlgAssRelOrdElem{S, T, U}(O, deepcopy(a), deepcopy(arr))
  end
end

(O::AlgAssRelOrd{S, T, U})(arr::Vector{S}, check::Bool = true) where {S, T, U} = begin
  M = basis_matrix(O, copy = false)
  N = matrix(base_ring(algebra(O)), 1, degree(O), arr)
  NM = N*M
  x = elem_from_mat_row(algebra(O), NM, 1)
  return O(x, arr, check)
end

(O::AlgAssRelOrd{S, T, U})(a::AlgAssRelOrdElem{S, T, U}, check::Bool = true) where {S, T, U} = begin
  b = elem_in_algebra(a) # already a copy
  if check
    (x, y) = _check_elem_in_order(b, O)
    !x && error("Algebra element not in the order")
    return AlgAssRelOrdElem{S, T, U}(O, b, y)
  else
    return AlgAssRelOrdElem{S, T, U}(O, b)
  end
end

(O::AlgAssRelOrd)(a::T, check::Bool = true) where T = O(algebra(O)(a), check)

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(a::AlgAssRelOrdElem, dict::IdDict)
  b = parent(a)()
  b.elem_in_algebra = Base.deepcopy_internal(a.elem_in_algebra, dict)
  if a.has_coord
    b.has_coord = true
    b.coordinates = Base.deepcopy_internal(a.coordinates, dict)
  end
  return b
end

################################################################################
#
#  Special elements
#
################################################################################

(O::AlgAssRelOrd{S, T, U})() where {S, T, U} = AlgAssRelOrdElem{S, T, U}(O)

one(O::AlgAssRelOrd) = O(one(algebra(O)))

zero(O::AlgAssRelOrd) = O()

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, a::AlgAssRelOrdElem)
  print(io, elem_in_algebra(a, copy = false))
end

################################################################################
#
#  Representation matrices
#
################################################################################

@doc Markdown.doc"""
    representation_matrix(x::AlgAssRelOrdElem, action::Symbol = :left) -> MatElem

> Returns a matrix representing multiplication with $x$ with respect to the basis
> of `order(x)`.
> The multiplication is from the left if `action == :left` and from the right if
> `action == :right`.
"""
function representation_matrix(x::AlgAssRelOrdElem, action::Symbol = :left)

  O = parent(x)
  M = basis_matrix(O, copy = false)
  M1 = basis_mat_inv(O, copy = false)

  B = representation_matrix(elem_in_algebra(x, copy = false), action)
  B = M*B
  B = B*M1

  return B
end

################################################################################
#
#  Trace
#
################################################################################

@doc Markdown.doc"""
    tr(x::AlgAssRelOrdElem{S, T, U}) where { S, T, U } -> S

> Returns the trace of $x$.
"""
function tr(x::AlgAssRelOrdElem)
  return tr(elem_in_algebra(x, copy = false))
end

@doc Markdown.doc"""
    trred(x::AlgAssRelOrdElem{S, T, U}) where { S, T, U } -> S

> Returns the reduced trace of $x$.
"""
function trred(x::AlgAssRelOrdElem)
  return trred(elem_in_algebra(x, copy = false))
end

################################################################################
#
#  isone/iszero
#
################################################################################

iszero(a::AlgAssRelOrdElem) = iszero(elem_in_algebra(a, copy = false))

isone(a::AlgAssRelOrdElem) = isone(elem_in_algebra(a, copy = false))

isunit(a::AlgAssRelOrdElem) = isunit(norm(a))

################################################################################
#
#  Norm
#
################################################################################

function norm(a::AlgAssRelOrdElem)
  OK = base_ring(parent(a))
  return OK(norm(elem_in_algebra(a, copy = false)))
end

function normred(a::AlgAssRelOrdElem)
  OK = base_ring(parent(a))
  return OK(normred(elem_in_algebra(a, copy = false)))
end