text:
# coding:utf-8

include("switch_depwarn.jl")
switch_depwarn!(false)

using DataFrames
using DataStructures
#using PyCall
#@pyimport polyglot.text as text

function input(prompt::AbstractString="")
  print(prompt)
  return chomp(readline())
end

# データ読み込み
d = open(deserialize, "./data/annotated_corpus.dic")
session_ids = filter(x -> !(in(x, [:unigram, :bigram])), collect(keys(d)))

# 記録用
if isfile("./data/simple_annotation.dic")
  log = open(deserialize, "./data/simple_annotation.dic")
else
  log = Dict()
end

# 一発話+segmentごとにデータを見ていく
for (i,sid) in enumerate(session_ids)
  # 登録済みならパス
  !haskey(log, sid) && (log[sid] = Dict())
  for j in collect(1:length(d[sid][:utterance]))
   ant, blob, spk, spkd = [], d[sid][:utterance][j],d[sid][:speaker][j], Dict(0=>"Librarian",1=>"Patron")
   dm, fc = d[sid][:domain_da_label][j],d[sid][:function_da_label][j]
   println("---------- Session $i: $j ----------")
   println(fc[2], " & ",dm[2])
   println(spkd[spk],": ",blob)
   while true
      println(ant)
      usrin = input(">>")
      usrin == "exit" && exit()
      if usrin == ""
        println(ant)
        chk = input("ok?")
        !in(chk, ["ok","y","yes",""]) && continue
        log[sid][i] = ant
        break
      elseif usrin == "back"
        pop!(ant)
      else
        push!(ant, usrin)
      end
    end
  end
  # 保存
  open(io -> serialize(io, log), "./data/simple_annotation.dic", "w")
end


text:
module PartialWaveFunctions

export f_logfact
include("factorials.jl")

export wignerd, wignerd_doublearg
export wignerD, wignerD_doublearg
export kronecker
include("wignerd.jl")

export CG, CG_doublearg
export clebschgordan, clebschgordan_doublearg
include("clebsch_gordan.jl")

end # module

# # # # # # # # # # # # # # # # # # # # # # # #
# created with PkgTemplates
# using PkgTemplates
# t = Template(;
#     user="mmikhasenko",
#     license="MIT",
#     authors="Misha Mikhasenko",
#     dir=joinpath(DEPOT_PATH[1], "dev"),
#     julia_version=v"1.2",
#     plugins=[
#         TravisCI(),
#         Codecov(),
#         AppVeyor(),
#         ],
#     )
# generate("PartialWaveFunctions", t)
# # # # # # # # # # # # # # # # # # # # # # # # # #

ext: jl


text:
#RSDeltaSigmaPort: Generate plots for noise transfer functions
#-------------------------------------------------------------------------------

#==Vector generators
===============================================================================#
"""`𝑓 = ds_freq(OSR=64, f0=0, quadrature=false)`

Frequency vector suitable for plotting the frequency response of an NTF.
"""
function ds_freq(;OSR::Int=64, f0::Float64=0, quadrature::Bool=false)
	local f_left, f_special
	if quadrature
		f_left = -0.5
		f_special = [f0 -f0]
	else
		f_left = 0
		f_special = f0
	end

	𝑓 = collect(range(f_left, stop=0.5, length=100))
	#Use finer spacing in the vicinity of the passband
	for fx = f_special
		f1 = max(f_left, fx-1/OSR)
		f2 = min(0.5,fx+2/OSR)
		deleteat!(𝑓, (𝑓 .<= f2) .& (𝑓 .>= f1))
		𝑓 = sort( vcat(𝑓, range(f1, stop=f2, length=100)) )
	end
	return 𝑓
end

#==NTF plots
===============================================================================#
function plotNTF(isbandpass)
	dBmin = -100 #Limit dB window to finite value
	𝑓min_log = 10^-2 #Limit min 𝑓 for log(𝑓-norm) plots
	nogrid = set(vmajor = false, vminor = false, hmajor=false, hminor=false)
	vgridonly = set(vmajor = true, vminor = true, hmajor=false, hminor=false)

	plot_PZ = plotPZ()
	plot_PZ.title = "NTF Poles and Zeros"

	plot_NTFlin = cons(:plot, legend=false,
		title = "NTF Magnitude Response",
		xyaxes=set(xscale=:lin, yscale=:lin, ymin=dBmin),
		labels=set(xaxis="Normalized frequency (1→f_s)", yaxis="dB"),
	)

	plot_NTFlog = cons(:plot, legend=false, grid = vgridonly,
#		title = "NTF Magnitude Response",
		xyaxes=set(xscale=:log, yscale=:lin, xmin=𝑓min_log, ymin=dBmin),
		labels=set(xaxis="Normalized frequency (1→f_B)", yaxis="dB"),
	)
	plot_NTFzoom = cons(:plot, legend=false, grid = vgridonly,
#		title = "NTF Magnitude Response",
		xyaxes=set(xscale=:lin, yscale=:lin, xmin=-0.6, ymax=0.6, ymin=dBmin),
		labels=set(xaxis="Normalized frequency offset", yaxis="dB"),
	)
	if !isbandpass; plot_NTFzoom = plot_NTFlog; end

	pcoll = push!(cons(:plot_collection, title="Sample Plot"),
		plot_PZ, plot_NTFlin, plot_NTFzoom
	)

	pcoll.bblist = [
		BoundingBox(0, 0.5, 0, 1), #Pole-zero
		BoundingBox(0.5, 1, 0, 0.5), #NTF linf
		BoundingBox(0.5, 1, 0.5, 1), #NTF zoom
	]

	return pcoll
end

function plotNTF!(pcoll, NTF::ZPKData, OSR::Int; color=:blue, f0::Real=0, STF=nothing)
	dfltline = cons(:a, line=set(style=:solid, color=color, width=2))
	gainline = cons(:a, line=set(style=:solid, color=:red, width=2))
	markerline = cons(:a, line=set(style=:dash, width=2.5))
	plot_PZ, plot_NTFlin, plot_NTFzoom = pcoll.plotlist
	isbandpass = f0!=0

	#Plot poles and zeros
	plotPZ!(plot_PZ, NTF, color=color)

	#Plot full-range NTF
	𝑓 = vcat(range(0, stop=0.75/OSR, length=100), range(0.75/OSR, stop=0.5, length=100))
	if isbandpass
		f1 = f0-1/(2*OSR); f2 = f0+1/(2*OSR)
		𝑓 = vcat(range(0, stop=f1, length=50), range(f1, stop=f2, length=100), range(f2, stop=0.5, length=50))
	end
	z = exp.(2j*π*𝑓)
	magNTF = dbv(evalTF(NTF, z))
	push!(plot_NTFlin,
		cons(:wfrm, DataF1(𝑓, magNTF), dfltline, label="|NTF|")
	)
	if STF != nothing
		plot_NTFlin.title = "NTF/STF Magnitude Response"
		magSTF = dbv(evalTF(STF,z))
		push!(plot_NTFlin,
			cons(:wfrm, DataF1(𝑓, magSTF), gainline, label="|STF|")
		)
	end

	#Plot "zoomed-in" NTF
	𝑓start = 0.01
	𝑓 = collect(range(𝑓start, stop=1.2, length=200)/(2*OSR))
	∫rng = [0, 0.5/OSR]
	if isbandpass
		f1 = f0-0.3/OSR; f2 = f0+0.3/OSR
		𝑓 = collect(range(f1, stop=f2, length=100))
		∫rng = f0 .+ [-.25, .25]/OSR
	end
	z = exp.(2j*π*𝑓)
	𝑓norm = (𝑓 .- f0)*(2*OSR)
	magNTF = dbv(evalTF(NTF, z))
	σ_NTF = dbv(rmsGain(NTF, ∫rng[1], ∫rng[2]))
	rmsgain_str = @sprintf("RMS Gain = %5.0fdB", σ_NTF)
	push!(plot_NTFzoom,
		cons(:wfrm, DataF1(𝑓norm, magNTF), dfltline, label="|NTF|"),
		cons(:atext, rmsgain_str, y=σ_NTF, offset=set(y=3), reloffset=set(x=0.5), align=:bc),
		cons(:hmarker, σ_NTF, markerline, strip=2),
	)

	return pcoll
end

plotNTF(NTF::ZPKData, args...; f0=0, kwargs...) =
	plotNTF!(plotNTF(f0!=0), NTF::ZPKData, args...; f0=f0, kwargs...)

#==documentNTF
===============================================================================#
"""`plot = documentNTF(NTF|ABCD|mod_struct; OSR=64, f0=0, quadrature=false, frespOnly=false, sizes=nothing)`

The first argument is either the NTF, ABCD matrix or a struct containing
NTF, OSR=64, f0=0, quadrature=0 and optionally stf.
If the first argument is a struct, then no other arguments should be supplied.
If the first argument is ABCD, the stf is also plotted.
"""
function documentNTF(arg1; OSR::Int=64, f0::Float64=0.0, quadrature::Bool=false, frespOnly::Bool=false, sizes=nothing)
	STF=nothing
	local NTF
	if isa(arg1, ZPKData)
		NTF = arg1
	else #if isa(arg1, Array) #Assume ABCD matrix
		ABCD = arg1
		NTF, STF = calculateTF(ABCD)
	end

	logplot = f0==0
	if quadrature
   	f_left = -0.5
   elseif logplot
   	f_left = 1e-3
   else
   	f_left = 0.0
	end

	pcoll = cons(:plot_collection, ncolumns=2, title="")
	if !frespOnly
		#Plot poles and zeros
		plot = plotPZ(NTF, color=:blue)
		push!(pcoll, plot)

		pcoll.bblist = [ #Format plot locations
			BoundingBox(0, 0.5, 0, 1), #Pole-zero plot
			BoundingBox(0.5, 1, 0, 1), #Frequency spectrum
		]
	end

	#Frequency response
	𝑓 = ds_freq(OSR=OSR, f0=f0, quadrature=quadrature)
	z = exp.(2π*j*𝑓)
	H = dbv.(evalTF(NTF, z))

	axes = logplot ? loglin : linlin
	plot = cons(:plot, axes, title="Frequency Response", legend=false,
		labels=set(xaxis="Normalized Frequency", yaxis=""),
		xyaxes=set(xmin=f_left, xmax=0.5, ymin=-100, ymax=15),
	)
	push!(pcoll, plot)

	H = waveform(𝑓, H)
	push!(plot,
		cons(:wfrm, H, line=set(style=:solid, color=:blue, width=2), label=""),
	)

	if !isnothing(STF)
#		plot.title = "NTF and STF"
		G = dbv.(evalTF(STF,z))
		G = waveform(𝑓, G)
		push!(plot,
			cons(:wfrm, G, line=set(style=:solid, color=:magenta, width=2), label=""),
		)
	end
	f1, f2 = ds_f1f2(OSR=OSR, f0=f0, iscomplex=quadrature)
	NG0 = dbv.(rmsGain(NTF, f1, f2))
	NG0w = waveform([max(f1,f_left), f2], NG0*[1,1])
	push!(plot,
		cons(:wfrm, NG0w, line=set(style=:dash, color=:black, width=2), label=""),
	)

	tstr = @sprintf("  %.0fdB", NG0)
	if f0==0 && logplot
		a = cons(:atext, tstr, x=sqrt(f_left*0.5/OSR), y=NG0+2, align=:bc)
	else
		a = cons(:atext, tstr, x=f2, y=NG0-1, align=:cl)
	end
	push!(plot, a)
	#set(h, 'FontSize',sizes.fs, 'FontWeight',sizes.fw);
	msg = @sprintf(" Inf-norm of H = %.2f\n 2-norm of H = %.2f", infnorm(NTF)[1], rmsGain(NTF,0,1));
	if f0<0.25
		a = cons(:atext, msg, x=0.48, y=0.0, align=:tr)
	else
		a = cons(:atext, msg, x=f_left, y=0.0, align=:tl)
	end
	push!(plot, a)
	#set(h, 'FontSize',sizes.fs, 'FontWeight',sizes.fw);
	if quadrature
		ING0 = dbv(rmsGain(ntf,-f1,-f2))
		ING0w = waveform(-1*[f1, f2], ING0*[1, 1])
		tstr = @sprintf("%.0fdB", ING0)
		push!(plot,
			cons(:wfrm, NG0w, line=set(style=:solid, color=:black, width=2), label=""),
			cons(:atext, tstr, x=-f0, y=ING0+1, align=:bc),
		)
	end
	#set(h, 'FontSize',sizes.fs, 'FontWeight',sizes.fw);

	return pcoll
end

"""`documentNTF(dsm, SYS=nothing; frespOnly=false, sizes=nothing)`
"""
function documentNTF(dsm::AbstractDSM, SYS=nothing; frespOnly::Bool=false, sizes=nothing)
	isnothing(SYS) && (SYS = synthesizeNTF(SYS))
	return documentNTF(SYS; OSR=dsm.OSR, f0=dsm.f0, quadrature=false, frespOnly=frespOnly, sizes=sizes)
end

#Last line


text:
function check_derivative(f, g, x::Number)
    auto_g = derivative(f)
    return maximum(abs.(g(x) - auto_g(x)))
end

function check_gradient(f, g, x::Vector{T}) where T <: Number
    auto_g = gradient(f)
    return maximum(abs.(g(x) - auto_g(x)))
end

function check_second_derivative(f, h, x::Number)
    auto_h = second_derivative(f)
    return maximum(abs.(h(x) - auto_h(x)))
end

function check_hessian(f, h, x::Vector{T}) where T <: Number
    auto_h = hessian(f)
    return maximum(abs.(h(x) - auto_h(x)))
end


text:
using OMEinsum
using Test
using LinearAlgebra
using CUDA

@testset "OMEinsum.jl" begin
    include("Core.jl")
    include("EinRule.jl")
    include("binaryrules.jl")
    include("utils.jl")
    include("einsum.jl")
    if CUDA.functional()
        include("cueinsum.jl")
    end
    include("autodiff.jl")
    include("einsequence.jl")
    include("interfaces.jl")

    include("contractionorder.jl")
end


text:
query= """
       {
         dog {
             name
               }
               }"""

@test Validatequery(Parse(query)) == "ok"

query="""
extend type Hello {
  world2: String
}
"""
try
  Validatequery(Parse(query))
catch e
  if e isa Diana.GraphQLError
    @test string(e.msg) =="{\"data\": null,\"errors\": [{\"message\": \"GraphQL cannot execute a request containing a TypeExtensionDefinition.\"}]}"
  end
end

query="""
type Hello {
  world2: String
}
"""
try
 Validatequery(Parse(query))
catch e
  if e isa Diana.GraphQLError
    @test string(e) =="Diana.GraphQLError(\"{\\\"data\\\": null,\\\"errors\\\": [{\\\"message\\\": \\\"GraphQL cannot execute a request containing a ObjectTypeDefinition.\\\"}]}\")"
  end
end

query="""
schema {
    query: RootQuery
  }
"""

try
  Validatequery(Parse(query))
catch e
  if e isa Diana.GraphQLError
    @test string(e) =="Diana.GraphQLError(\"{\\\"data\\\": null,\\\"errors\\\": [{\\\"message\\\": \\\"GraphQL cannot execute a request containing a SchemaDefinition.\\\"}]}\")"
  end
end

query="""
{
  dog {
    name
  }
}

query getName {
  dog {
    owner {
      name
    }
  }
}"""

try
  Validatequery(Parse(query))
catch e
  if e isa Diana.GraphQLError
    @test string(e) =="Diana.GraphQLError(\"{\\\"data\\\": null,\\\"errors\\\": [{\\\"message\\\": \\\"This anonymous operation must be the only defined operation.\\\"}]}\")"
  end
end

query="""
query getName {
  dog {
    name
  }
}

query getName {
  dog {
    owner {
      name
    }
  }
}"""

try
  Validatequery(Parse(query))
catch e
  if e isa Diana.GraphQLError
    @test string(e) =="Diana.GraphQLError(\"{\\\"data\\\": null,\\\"errors\\\": [{\\\"message\\\": \\\"There can only be one operation named 'getName'.\\\"}]}\")"
  end
end

query="""
query dogOperation {
  dog {
    name
  }
}

mutation dogOperation {
  mutateDog {
    id
  }
}"""

try
  Validatequery(Parse(query))
catch e
  if e isa Diana.GraphQLError
    @test string(e) =="Diana.GraphQLError(\"{\\\"data\\\": null,\\\"errors\\\": [{\\\"message\\\": \\\"There can only be one operation named 'dogOperation'.\\\"}]}\")"
  end
end

query="""
subscription sub {
  newMessage {
    body
    sender
  }
  __typename
}"""

try
  Validatequery(Parse(query))
catch e
  if e isa Diana.GraphQLError
    @test string(e) == "Diana.GraphQLError(\"{\\\"data\\\": null,\\\"errors\\\": [{\\\"message\\\": \\\"subscription must select only one top level field.\\\"}]}\")"
  end
end

query="""
subscription JAJA {
         Actor(
         filter: {
           mutation_in:[CREATED]
           }){
             updatedFields
             }
               File
               }"""

try
  Validatequery(Parse(query))
catch e
  if e isa Diana.GraphQLError
    @test string(e) =="Diana.GraphQLError(\"{\\\"data\\\": null,\\\"errors\\\": [{\\\"message\\\": \\\"subscription must select only one top level field.\\\"}]}\")"
  end
end

query= """
fragment multipleSubscriptions on Subscription {
  newMessage {
    body
    sender
  }
  disallowedSecondRootField
}"""

try
  Validatequery(Parse(query))
catch e
  if e isa Diana.GraphQLError
    @test string(e) =="Diana.GraphQLError(\"{\\\"data\\\": null,\\\"errors\\\": [{\\\"message\\\": \\\"subscription must select only one top level field\\\"}]}\")"
  end
end

query= """
 subscription{
         ...multipleSubscriptions
         ...multipleSubscriptions
         }"""

try
  Validatequery(Parse(query))
catch e
  if e isa Diana.GraphQLError
    @test string(e) =="Diana.GraphQLError(\"{\\\"data\\\": null,\\\"errors\\\": [{\\\"message\\\": \\\"subscription must select only one top level field.\\\"}]}\")"
  end
end

query = """
fragment nameFragment on Dog { # unused
  name
}

{
  dog {
    name
  }
}
                 """

try
  Validatequery(Parse(query))
catch e
  if e isa Diana.GraphQLError
    @test string(e) =="Diana.GraphQLError(\"{\\\"data\\\": null,\\\"errors\\\": [{\\\"message\\\": \\\"Fragment  nameFragment is not used.\\\"}]}\")"
  end
end

query = """
  {
         dog {
             ...fragmentOne
               ...fragmentotro
               }
               }
               fragment fragmentOne on Dog {
                 name
                 }
                 """

try
  Validatequery(Parse(query))
catch e
  if e isa Diana.GraphQLError
    @test string(e) =="Diana.GraphQLError(\"{\\\"data\\\": null,\\\"errors\\\": [{\\\"message\\\": \\\"Unknown fragment  fragmentotro.\\\"}]}\")"
  end
end

query = """
 fragment fragmentOne on Dog {
                 name
                 }

  {
         dog {
             ...fragmentOne
               ...fragmentotro
               }
               }

                 """

try
  Validatequery(Parse(query))
catch e
  if e isa Diana.GraphQLError
    @test string(e) == "Diana.GraphQLError(\"{\\\"data\\\": null,\\\"errors\\\": [{\\\"message\\\": \\\"Unknown fragment  fragmentotro.\\\"}]}\")"
  end
end

query = """
       {
         dog {
             ...dogFragment

               }
               }
               fragment dogFragment on Dog {
                 name
                   owner {
                       ...a
                         }
                         }

                         fragment dogFragment on Dog {
                           name
                             pets {
                                 ...dogFragment
                                   }
                                   }
"""

try
  Validatequery(Parse(query))
catch e
  if e isa Diana.GraphQLError
    @test string(e) =="Diana.GraphQLError(\"{\\\"data\\\": null,\\\"errors\\\": [{\\\"message\\\": \\\"There can only be one fragment named 'dogFragment'.\\\"}]}\")"
  end
end

query = """
       {
         dog {
             ...dogFragment

               }
               }
               fragment dogFragment on Dog {
                 name
                   owner {
                       ...a
                         }
                         }

                         fragment ownerFragment on Dog {
                           name
                             pets {
                                 ...dogFragment
                                   }
                                   }
           fragment a on Dog{
           ...b
           }
           fragment b on Dog{
           ...c
           }
           fragment c on Dog{
           ...ownerFragment
           }"""

try
  Validatequery(Parse(query))
catch e
  if e isa Diana.GraphQLError
    @test string(e) =="Diana.GraphQLError(\"{\\\"data\\\": null,\\\"errors\\\": [{\\\"message\\\": \\\"Cannot spread fragment dogFragment within itself via  ownerFragment a b c.\\\"}]}\")"
  end
end


text:
module Constructions

using GroupsCore
import AbstractAlgebra

include("direct_product.jl")
include("direct_power.jl")
include("wreath_product.jl")

end # of module Constructions


text:
# Autogenerated wrapper script for NetCDFF_jll for aarch64-linux-gnu-libgfortran3
export libnetcdff

using NetCDF_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("NetCDFF")
JLLWrappers.@declare_library_product(libnetcdff, "libnetcdff.so.7")
function __init__()
    JLLWrappers.@generate_init_header(NetCDF_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libnetcdff,
        "lib/libnetcdff.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Ex7_3c
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Example 7.3c - Nullspace-based synthesis

println("Example 7.3c")

A = diagm([ -2, -3, -1, 2, 3]);
Bu = [1 1 0 0 0]'; Bd = [0 0 2 0 0]';
Bf = [0 0 0 2 2;0 0 0 0 0]';
C = [-1 0 -1 1.5 0; 0 -1 0 0 2.5];
Du = [ 1 1]'; Dd = [1 0]'; Df = [1 0; 1 1];
p = 2; mu = 1; md = 1; mf = 2;       # enter dimensions

# setup the synthesis model
sysf = fdimodset(dss(A,[Bu Bd Bf],C,[Du Dd Df]),c =1:mu, d = mu .+ (1:md), f = (mu+md) .+ (1:mf));

# call of EFDSYN with the options for stability degree -3 and the synthesis
# of a scalar output filter
Qt, Rft, info = efdsyn(sysf, sdeg = -3, rdim = 1);

# compute Q and Rf; scale to meet example
scale = -1
Q = FDFilter(Qt.sys*scale,Qt.outputs,Qt.controls);
Rf = FDFilterIF(Rft.sys*scale,faults = Rft.faults);

# display results
println("Q = $(dss2rm(Q.sys))")
println("Rf = $(dss2rm(Rf.sys))")

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
@cache mutable struct DImplicitEulerCache{uType,uNoUnitsType,N} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  tmp::uType
  atmp::uNoUnitsType
  nlsolver::N
end

mutable struct DImplicitEulerConstantCache{N} <: OrdinaryDiffEqConstantCache
  nlsolver::N
end

function alg_cache(alg::DImplicitEuler,du,u,res_prototype,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  γ, c = 1, 1
  α = 1
  nlsolver = build_nlsolver(alg,u,uprev,p,t,dt,f,res_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,α,Val(false))

  DImplicitEulerConstantCache(nlsolver)
end

function alg_cache(alg::DImplicitEuler,du,u,res_prototype,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})

  γ, c = 1, 1
  α = 1
  nlsolver = build_nlsolver(alg,u,uprev,p,t,dt,f,res_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,α,Val(true))

  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)

  DImplicitEulerCache(u,uprev,uprev2,tmp,atmp,nlsolver)
end


text:
export mention,
    reply,
    plaintext,
    upload_file,
    set_game

"""
    mention(x::Union{DiscordChannel, Member, Role, User}) -> String

Get the mention string for an entity.
"""
mention(c::DiscordChannel) = "<#$(c.id)>"
mention(r::Role) = "<@&$(r.id)>"
mention(u::User) = "<@$(u.id)>"
function mention(m::Member)
    return ismissing(m.nick) || m.nick === nothing ? mention(m.user) : "<@!$(m.user.id)>"
end

"""
    reply(
        c::Client,
        m::Message,
        content::Union{AbstractString. AbstractDict, NamedTuple};
        at::Bool=false,
    )

Reply (send a message to the same [`DiscordChannel`](@ref)) to a [`Message`](@ref).
If `at` is set, then the message is prefixed with the sender's mention.
"""
function reply(c::Client, m::Message, content::AbstractString; at::Bool=false)
    content = at ? mention(m.author) * " " * content : content
    return create_message(c, m.channel_id; content=content)
end

function reply(c::Client, m::Message, embed::Union{AbstractDict, NamedTuple}; at::Bool=false)
    content = at ? mention(m.author) : ""
    return create_message(c, m.channel_id; content=content, embed=embed)
end

"""
    plaintext(m::Message) -> String
    plaintext(c::Client, m::Message) -> String

Get the [`Message`](@ref) contents with any [`User`](@ref) mentions replaced with their
plaintext. If a [`Client`](@ref) is provided, [`DiscordChannel`](@ref)s [`Role`](@ref) are
also replaced. However, only channels and roles stored in state are replaced; no API
requests are made.
"""
function plaintext(m::Message)
    content = m.content

    for u in coalesce(m.mentions, User[])
        name = "@$(u.username)"
        content = replace(content, "<@$(u.id)>" => name)
        content = replace(content, "<@!$(u.id)>" => name)
    end

    return content
end

function plaintext(c::Client, m::Message)
    content = m.content

    for u in coalesce(m.mentions, User[])
        member = get(c.state, Member; guild=m.guild_id, user=u.id)
        nick = if member !== nothing && member.nick isa String
            "@$(member.nick)"
        else
            "@$(u.username)"
        end
        content = replace(content, "<@$(u.id)>" => "@$(u.username)")
        content = replace(content, "<@!$(u.id)>" => "@$nick")
    end

    guild = get(c.state, Guild; guild=m.guild_id)
    if guild !== nothing
        for r in coalesce(m.mention_roles, Snowflake[])
            role = get(c.state, Role; guild=m.guild_id, role=r)
            if role !== nothing
                content = replace(content, "<@&$r>" => "@$(role.name)")
            end
        end

        for cap in unique(eachmatch(r"<#(\d+?)>", content))
            ch = get(c.state, DiscordChannel; channel=parse(Snowflake, first(cap.captures)))
            if ch !== nothing
                content = replace(content, cap.match => "#$(ch.name)")
            end
        end
    end

    return content
end

"""
    upload_file(c::Client, ch::DiscordChannel, path::AbstractString; kwargs...) -> Message

Send a [`Message`](@ref) with a file [`Attachment`](@ref). Any keywords are passed on to
[`create_message`](@ref).
"""
function upload_file(c::Client, ch::DiscordChannel, path::AbstractString; kwargs...)
    return create_message(c, ch.id; kwargs..., file=open(path))
end

"""
    set_game(
        c::Client,
        name::AbstractString,
        type::Union{ActivityType, Int}=AT_GAME,
        since::Union{Int, Nothing}=nothing,
        status::Union{PresenceStatus, AbstractString}=PS_ONLINE,
        afk::Bool=false,
        kwargs...,
    ) -> Bool

Shortcut for [`update_status`](@ref) to set the client's [`Activity`](@ref).
"""
function set_game(
    c::Client,
    game::AbstractString;
    type::Union{ActivityType, Int}=AT_GAME,
    since::Union{Int, Nothing}=nothing,
    status::Union{PresenceStatus, AbstractString}=PS_ONLINE,
    afk::Bool=false,
    kwargs...,
)
    activity = merge(Dict("name" => game, "type" => type), kwargs)

    return update_status(c, since, activity, status, afk)
end


text:
using Printf

function main_pre(n_x,n_y,s_x,e_x,s_y,e_y,s_z,e_z,n_z)
    xy_or_r = 1          # 1:x,y 2:r,t
    xnum = n_x           #       2:rnum
    ynum = n_y           #       2:tnum
    znum = n_z           #       2:tnum
    lenx = abs(e_x-s_x)           #       2:inr
    leny = abs(e_y-s_y)           #       2:outr
    lenz = abs(e_z-s_z)           #       2:outr
    st_ang  = -1/2*pi       # 2:angle
    en_ang  = -3/2*pi          # 2:angle

    outdir = "grid_morton"
    make_dir(outdir)
    result = "result_morton"
    make_dir(result)
    result = "post_result_morton"
    make_dir(result)
    nodes, xnum_max, ynum_max, znum_max = mk_gird(xnum,ynum,znum,lenx,leny,lenz,xy_or_r,s_x,s_y,s_z,st_ang,en_ang,outdir)
end

function mk_gird(xnum,ynum,znum,lenx,leny,lenz,xy_or_r,s_x,s_y,s_z,st_ang,en_ang,outdir)
    """
    nodes[i,j,k]
    i   : x,r方向の番号
    j   : y,theta方向の番号
    k=1 : 点のx座標
    k=2 : 点のy座標

    nodes[1,:]      : x方向境界
    nodes[:,1]      : y方向境界
    nodes[xnum+2,:] : x方向境界
    nodes[:,ynum+2] : y方向境界
    """
    xnum_max = xnum+1+2
    ynum_max = ynum+1+2
    znum_max = znum+1+2
    nodes = zeros(xnum_max, ynum_max, znum_max,3)
    if xy_or_r == 1
        for i in 2:xnum_max-1
            for j in 2:ynum_max-1
                for k in 2:znum_max-1
                    x = lenx/(xnum)*(i-2) + s_x
                    y = leny/(ynum)*(j-2) + s_y
                    z = lenz/(znum)*(k-2) + s_z
                    nodes[i,j,k,1] = x
                    nodes[i,j,k,2] = y
                    nodes[i,j,k,3] = z
                end
            end
        end

    end

    #=
    仮想セルの作成：現在あるセル境界線を延長して作成
    そのまま延長した際にセルが交差するのを防ぐため、仮想セルの大きさを100^-1倍にする.
    ベクトルで考えれば下記のような計算になる．たぶん
    また、[1,1]等の角にはダミー値が入っている.
    =#
    #=
    for j in 1:ynum_max
        for k in 1:2
            nodes[1,j,k]        =  1.01*nodes[2,j,k] - 0.01*nodes[3,j,k]
            nodes[xnum_max,j,k] =  1.01*nodes[xnum_max-1,j,k] - 0.01*nodes[xnum_max-2,j,k]
        end
    end
    for i in 1:xnum_max
        for k in 1:2
            nodes[i,1,k]        =  1.01*nodes[i,2,k] - 0.01*nodes[i,3,k]
            nodes[i,ynum_max,k] =  1.01*nodes[i,ynum_max-1,k] - 0.01*nodes[i,ynum_max-2,k]
        end
    end
    =#

    #=
    fff=outdir*"/nodes"
    open(fff,"w") do f
        write(f,"nodes: xnum, ynum , x, y\n")
        for i in 1:xnum_max
            for j in 1:ynum_max
                x = @sprintf("%8.8e", nodes[i,j,1])
                y = @sprintf("%8.8e", nodes[i,j,2])
                write(f,string(i)*" "*string(j)*" "*x*" "*y*"\n")
            end
        end
    end
    println("write "*fff)
    =#

    # nodes_forvtk
    nodes_num = zeros(Int,xnum_max, ynum_max, znum_max)
    fff=outdir*"/nodes_forvtk"
    open(fff,"w") do f
        write(f,"nodes: xnum, ynum , x, y\n")
        a=1
        for i in 2:xnum_max-1
            for j in 2:ynum_max-1
                for k in 2:znum_max-1
                    x = @sprintf("%8.8e", nodes[i,j,k,1])
                    y = @sprintf("%8.8e", nodes[i,j,k,2])
                    z = @sprintf("%8.8e", nodes[i,j,k,3])
                    write(f,string(a)*" "*x*" "*y*" "*z*"\n")
                    nodes_num[i,j,k] = a
                    a = a+1
                end
            end
        end
    end
    println("write "*fff)

    fff=outdir*"/nodesnum"
    open(fff,"w") do f
        write(f,"nodesnum: xnum_max, ynum_max\n")
        write(f,string(xnum_max)*" "*string(ynum_max)*" "*string(znum_max)*"\n")
    end

    ### element ###
    fff=outdir*"/element_forvtk"
    open(fff,"w") do f
        write(f,"elements:cell_xnum, l1,l2,l3,l4,u1,u2,u3,u4 \n")
        a=1
        for i in 2:xnum_max-2
            for j in 2:ynum_max-2
                for k in 2:znum_max-2
                    write(f,string(a))
                    for l in 1:2
                        d1 = @sprintf("%1.0f", nodes_num[i  ,j  ,k+(l-1)])
                        d2 = @sprintf("%1.0f", nodes_num[i  ,j+1,k+(l-1)])
                        d3 = @sprintf("%1.0f", nodes_num[i+1,j+1,k+(l-1)])
                        d4 = @sprintf("%1.0f", nodes_num[i+1,j  ,k+(l-1)])
                        write(f," "*d1*" "*d2*" "*d3*" "*d4)
                    end
                    write(f,"\n")
                    a = a+1
                end
            end
        end
    end
    println("write "*fff)

    return  nodes,xnum_max,ynum_max,znum_max
end

function make_dir(outdir)
    k = 0
    try rm(outdir,recursive=true)
    catch
        mkdir(outdir)
        k = 1
    end

    if k == 0
        mkdir(outdir)
    end
end


text:
using Pkg

packages = [
    "IJulia",
    "Plots",
    "PyPlot"]

for package=packages
    println("installing $package")
    Pkg.add(package)
end

Pkg.update()
Pkg.precompile()

max_stars_repo_name: MarshMue/math-tweets

max_issues_repo_name: MarshMue/math-tweets

max_forks_repo_name: MarshMue/math-tweets


text:
using UnPack, OrderedCollections, ForwardDiff, JLD2, NLsolve

mutable struct NKEZDisaster{T <: Real, S, N}
    β::T
    γ::T
    ψ::T
    ν::T
    ν̅::T
    χ::T
    δ::T
    α::T
    ϵ::T
    θ::T
    π_ss::T
    ϕ_r::T
    ϕ_π::T
    ϕ_y::T
    χ_y::T
    ρ_β::T
    ρ_l::T
    ρ_r::T
    σ_β::T
    σ_l::T
    σ_r::T
    μ_a::T
    σ_a::T
    κ_a::T
    disaster_occur_spec::Symbol
    disaster_intensity_spec::Symbol
    disaster_para::NamedTuple{S, NTuple{N, T}}
    N_approx::NamedTuple{(:q, :s₁, :s₂, :ω), NTuple{4, Int}}
    S::OrderedDict{Symbol, Int}
    J::OrderedDict{Symbol, Int}
    E::OrderedDict{Symbol, Int}
    SH::OrderedDict{Symbol, Int}
end

# Absent a better way, I assume (1) each specification of disaster risk
# has a unique name and (2) disaster_para has correctly named parameters,
# given the specification's name. To see the process implied by
# `disaster_occur_spec` and `disaster_intensity_spec`, see
# the functions `infer_ccgf` and `infer_X̅` at the end of this file.
function NKEZDisaster(disaster_occur_spec::Symbol = :PoissonNormalMixture,
                      disaster_intensity_spec::Symbol = :CoxIngersollRoss,
                      disaster_para::NamedTuple{S1, NTuple{N1, T}} =
                      (σ_k = .01, ρ_p = .08^(1. / 4.), p = .0355 / 4., σ_p = .0114 / 4. / (.02 / sqrt(4.)) / sqrt(.0355 / 4.));
                      β::T = .99, γ::T = 3.8, ψ::T = 1. / .75, ν::T = 1., ν̅ = 0.72,
                      χ::T = 4., δ::T = 0.025, α::T = 0.33, ϵ::T = 10., θ::T = 0.7,
                      π_ss::T = 0., ϕ_r::T = 0.5, ϕ_π::T = 1.3, ϕ_y::T = 0.25,
                      χ_y::T = 1.6, ρ_β::T = 0.1, ρ_l::T = 0.1,, ρ_r::T = 0.,
                      σ_β::T = sqrt((log(β) / 4.)^2 * (1. - ρ_β^2)),
                      σ_l::T = 0.01, σ_r::T = 0.01, μ_a::T = 0.0125,
                      σ_a::T = 0.01, κ_a::T = 1.,
                      N_approx::NamedTuple{(:q, :s₁, :s₂, :ω), NTuple{4, Int}} =
                      (q = 1, s₁ = 1, s₂ = 1, ω = 1)) where {T <: Real, S1, N1}

    @assert all(N_approx[k] > 0 for k in keys(N_approx)) "N_approx must be at least 1 for all variables."

    ## Create Indexing dictionaries.

    # Note that for the exogenous shock
    # state variables, instead of e.g. η_L and η_A, I use η_l and η_a
    # since the uppercase variable will not appear in the jumps/states.
    S_init  = [:k₋₁, :logΔ₋₁, :r₋₁, :output₋₁, :η_β, :η_l, :η_r, :a, :η_k] # State Variables
    J_init  = [:output, :c, :l, :v, :ce, :ω, :ℓ, :β̅, :w, :r, :π, :q, :x,
               :rk, :rq, :mc, :s₁, :s₂, :logΔ] # Jump variables
    E_init  = [:value_fnct, :certainty_equiv, :ez_fwd_diff,
               :eq_β̅, :wage, :labor_disutility, :euler, :cap_ret,
               :eq_mc, :kl_ratio, :eq_s₁, :eq_s₂,
               :tobin,, :eq_rq, :phillips_curve, :price_dispersion,
               :mp, :output_market_clear, :production] # Equations
    SH_init = [:ε_β, :ε_l, :ε_r, :ε_a, :ε_k, :ε_p] # Exogenous shocks

    # Add approximations for forward-difference equations
    for var in [:q, :s₁, :s₂, :ω]
        inds = (var == :q) ? (1:N_approx[var]) : (0:(N_approx[var] - 1))
        push!(J_init, [Symbol(:d, var, "$(i)") for i in inds]...)
        push!(J_init, [Symbol(:p, var, "$(i)") for i in 1:N_approx[var]]...)
        push!(E_init, [Symbol(:eq_d, var, "$(i)") for i in inds]...)
        push!(E_init, [Symbol(:eq_p, var, "$(i)") for i in 1:N_approx[var]]...)
    end

    # Specify random process(es) for whether a disaster occurs or not
    if disaster_occur_spec in [:PoissonNormalMixture, :Bernoulli]
        # Nothing need to be added
    end

    # Specify random process(es) for "intensity" (size or frequency) of a disaster.
    if disaster_intensity_spec in [:CoxIngersollRoss, :TwoStateMarkovChain,
                                   :TruncatedCoxIngersollRoss]
        push!(S_init, :p)
    elseif disaster_intensity_spec in [:LogAR1]
        push!(S_init, :logp)
    end

    S  = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(S_init))
    J  = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(J_init))
    E  = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(E_init))
    SH = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(SH_init))

    return NKEZDisaster(β, γ, ψ, ν,, ν̄,  χ, δ, α, ϵ, θ, π_ss, ϕ_r, ϕ_π, ϕ_y,
                        χ_y, ρ_β, ρ_l, ρ_r, σ_β, σ_l, σ_r, μ_a, σ_a, κ_a,
                        disaster_spec, disaster_para,
                        N_approx, S, J, E, SH)
end

function nk_ez_disaster(m::NKEZDisaster{T, SNK, NNK}) where {T <: Real, SNK, NNK}

    # Get parameters
    @unpack β, γ, ψ, ν, ν̅, χ, δ, α, ϵ, θ, π_ss, ϕ_r, ϕ_π, ϕ_y = m
    @unpack χ_y, ρ_β, ρ_l, ρ_r, σ_β, σ_l, σ_r, μ_a, σ_a, κ_a = m
    @unpack disaster_occur_spec, disaster_intensity_spec, disaster_para = m
    r_ss = infer_r_ss(m)
    X̅    = infer_X̅(m)
    𝔼η_k = infer_𝔼η_k(m)

    # Unpack indexing dictionaries
    @unpack N_approx, S, J, E, SH = m
    @unpack k₋₁, logΔ₋₁, r₋₁, output₋₁, η_β, η_l, η_r, a, η_k = S
    @unpack output, c, l, v, ce, ω, ℓ, β̅, w, r = J
    @unpack π, q, x, rk, rq, mc, s₁, s₂, logΔ  = J
    @unpack value_fnct, certainty_equiv, ez_fwd_diff = E
    @unpack eq_β̅, wage, labor_disutility, euler, cap_ret, eq_mc = E
    @unpack kl_ratio, eq_s₁, eq_s₂, tobin, eq_rq = E
    @unpack phillips_curve, price_dispersion, mp = E
    @unpack output_market_clear, production = E
    @unpack ε_β, ε_l, ε_r, ε_a, ε_k, ε_p = SH

    if disaster_intensity_spec in [:CoxIngersollRoss, :TwoStateMarkovChain,
                                   :TruncatedCoxIngersollRoss]
        p = m.S[:p]
        disaster_intensity_var = p
    elseif disaster_intensity_spec in [:LogAR1]
        logp = m.S[:logp]
        disaster_intensity_var = logp
    end

    Nz = length(S)
    Ny = length(J)
    Nε = length(SH)

    ## Define nonlinear equations

    # Some helper functions
    _Φ(Xin, Kin)  = X̅ ^ (1. / χ) / (1. - 1. / χ) * (Xin / Kin) ^ (1. - 1. / χ) - X̅ / (χ * (χ - 1.))
    _Φ′(Xin, Kin) = X̅ ^ (1. / χ) * (Xin / Kin) ^ (- 1. / χ)
    Φ(z, y)  = _Φ(exp(y[x]), exp(z[η_k] + z[k₋₁]))
    Φ′(z, y) = _Φ′(exp(y[x]), exp(z[η_k] + z[k₋₁]))
    m_ξ(z, y) = z[η_β] + log(β) - y[β̅] + γ * y[c] -
        (1. - γ) * y[ℓ] - (ψ - γ) * y[ce] - γ * μ_a
    μ_y_bgp(z, y) = μ_a + κ_a * 𝔼η_k # calculate growth rate of output along balanced growth path
    function m_fwd!(i, Γ₅, Γ₆)
        Γ₅[i, β̅] = 1.
        Γ₅[i, a] = -γ
        Γ₆[i, c] = -γ
        Γ₆[i, ℓ] = (1. - γ)
        Γ₆[i, v] = (ψ  - γ)
    end
    pstar(y) = log(ϵ / (ϵ - 1.)) + y[s₁] - y[s₂]
    μ_η_k    = infer_μ_disaster_occur(m)
    μ_disi   = infer_μ_disaster_intensity(m)

    function μ(F, z, y)
        # Expected value of η_k conditional on time t
        μ_η_k_v     = μ_η_k(z, y)

        F[k₋₁]      = log(1. - δ + Φ(z, y)) + z[η_k] + z[k₋₁]
        F[v₋₁]      = y[v]
        F[r₋₁]      = y[r]
        F[output₋₁] = y[output]
        F[η_β]      = ρ_β * z[η_β]
        F[η_l]      = ρ_l * z[η_l]
        F[η_r]      = ρ_r * z[η_r]
        F[a]        = κ_a * μ_η_k_v
        F[η_k]      = μ_η_k_v
        F[disaster_intensity_var] = μ_disi(z, y)
    end

    function ξ(F, z, y)
        F_type = eltype(F)

        ## Pre-evaluate (just once) some terms
        Φv     = Φ(z, y)
        Φ′v    = Φ′(z, y)
        pstarv = pstar(y)
        m_ξv   = m_ξ(z, y)

        ## Non-forward-difference equations
        F[value_fnct]          = 1. / (1. - ψ) * (y[β̅] + y[ω]) - y[v]
        F[certainty_equiv]     = 1. / (1. - ψ) * (y[β̅] - (z[η_β] + log(β)) + log(exp(y[ω]) - 1.)) - y[ce]
        F[wage]                = log(ψ) + z[η_l] + log(ν̅) + y[c] + ν * y[l] - (1. - ψ) / ψ * y[ℓ] - y[w]
        F[labor_disutility]    = ψ / (1. - ψ) * log(1. + (ψ - 1.) * exp(z[η_l]) * ν̅ *
                                                    exp((1. + ν) * y[l]) / (1. + ν)) - y[ℓ]
        F[euler]               = y[r] + m_ξv
        F[eq_mc]               = (1. - α) * (y[w] - log(1. - α)) + α * (y[rk] - log(α)) - y[mc]
        F[kl_ratio]            = log(α) - log(1. - α) + y[w] - y[rk] - (z[η_k] + z[k₋₁] - y[l])
        F[tobin]               = log(Φ′v) + y[q]
        F[eq_rq]               = log(1. - δ + Φv - Φ′v * exp(y[x] - (z[η_k] + z[k₋₁]))) - y[rq]
        F[phillips_curve]      = (1. - ϵ) * y[π] - log((1. - θ) * exp((1. - ϵ) * (pstarv + y[π])) + θ)
        F[price_dispersion]    = y[logΔ] - ϵ * y[π] - log((1. - θ) * exp(-ϵ * (pstarv + y[π])) + θ * exp(z[logΔ₋₁]))

        F[mp]                  = (1. - ϕ_r) * r_ss + ϕ_r * z[r₋₁] +
            (1. - ϕ_r) .* (ϕ_π * (y[π] - π_ss) + ϕ_y *
                           (y[output] - z[output₋₁] + (μ_a + z[a] - mp_μ_y_bgp(z, y)))) + z[η_r] - y[r]
        F[output_market_clear] = y[output] - log(exp(y[c]) + exp(y[x]))
        F[production]          = log(exp(α * z[k₋₁] + (1. - α) * y[l]) - χ_y) - y[logΔ] - y[output]
        F[eq_β̅]                = log(1. - exp(z[η_β])) - y[β̅]

        ## Forward-difference equations separately handled b/c recursions
        F[cap_ret]     = y[q]  - log(sum([exp(y[J[Symbol("dq$(i)")]]) for i in 1:N_approx[:q]]) +
                                     exp(y[J[Symbol("pq$(N_approx[:q])")]]))
        F[eq_s₁]       = y[s₁] - log(sum([exp(y[J[Symbol("ds₁$(i)")]]) for i in 0:(N_approx[:s₁] - 1)]) +
                                      exp(y[J[Symbol("ps₁$(N_approx[:s₁])")]]))
        F[eq_s₂]       = y[s₂] - log(sum([exp(y[J[Symbol("ds₂$(i)")]]) for i in 0:(N_approx[:s₂] - 1)]) +
                                      exp(y[J[Symbol("ps₂$(N_approx[:s₂])")]]))
        F[ez_fwd_diff] = y[ω]  - log(sum([exp(y[J[Symbol("dω$(i)")]]) for i in 0:(N_approx[:ω] - 1)]) +
                                     exp(y[J[Symbol("pω$(N_approx[:ω])")]]))

        # Set initial boundary conditions
        F[E[:eq_dq1]]  = -y[J[:dq1]] + m_ξv
        F[E[:eq_pq1]]  = -y[J[:pq1]] + m_ξv
        F[E[:eq_ds₁0]] = y[J[:ds₁0]] - y[mc] - y[output]
        F[E[:eq_ps₁1]] = μ_a + log(θ) - y[J[:ps₁1]] + m_ξv
        F[E[:eq_ds₂0]] = y[J[:ds₂0]] - y[output]
        F[E[:eq_ps₂1]] = μ_a + log(θ) - y[J[:ps₂1]] + m_ξv
        F[E[:eq_dω0]]  = y[J[:dω0]]
        F[E[:eq_pω1]]  = μ_a - y[c] - y[J[:pω1]] + m_ξv

        # Recursions for forward-difference equations
        for i in 2:N_approx[:q]
            F[E[Symbol("eq_dq$(i)")]]    = -y[J[Symbol("dq$(i)")]] + m_ξv
            F[E[Symbol("eq_pq$(i)")]]    = -y[J[Symbol("pq$(i)")]] + m_ξv
        end
        for i in 2:N_approx[:s₁]
            F[E[Symbol("eq_ds₁$(i-1)")]] = μ_a + log(θ) - y[J[Symbol("ds₁$(i-1)")]] + m_ξv
            F[E[Symbol("eq_ps₁$(i)")]]   = μ_a + log(θ) - y[J[Symbol("ps₁$(i)")]]   + m_ξv
        end
        for i in 2:N_approx[:s₂]
            F[E[Symbol("eq_ds₂$(i-1)")]] = μ_a + log(θ) - y[J[Symbol("ds₂$(i-1)")]] + m_ξv
            F[E[Symbol("eq_ps₂$(i)")]]   = μ_a + log(θ) - y[J[Symbol("ps₂$(i)")]]   + m_ξv
        end
        for i in 2:N_approx[:ω]
            F[E[Symbol("eq_dω$(i-1)")]] = μ_a - y[c] - y[J[Symbol("dω$(i-1)")]] + m_ξv
            F[E[Symbol("eq_pω$(i)")]]   = μ_a - y[c] - y[J[Symbol("pω$(i)")]]   + m_ξv
        end
    end

    # The cache is initialized as zeros so we only need to fill non-zero elements
    Λ = zeros(T, Nz, Ny)

    # The cache is initialized as zeros so we only need to fill non-zero elements
    function Σ(F, z)
        F_type = eltype(F)

        # AR(1) processes
        F[η_β, ε_β] = σ_β
        F[η_l, ε_l] = σ_l
        F[η_r, ε_r] = σ_r

        # Productivity process
        F[a, ε_a] = σ_a
        F[a, ε_k] = κ_a

        # Disaster risk
        F[η_k, ε_k] = one(F_type)
        F[disaster_intensity_var, ε_p] = Σ_disi(z)
    end

    ccgf = infer_ccgf(m)

    ## Forward-looking variables
    Γ₅ = zeros(T, Ny, Nz)
    Γ₆ = zeros(T, Ny, Ny)

    # Equations w/out SDF terms and are not forward-difference equations
    Γ₆[euler, π] = -one(T)

    # Equations with SDF terms but are not forward-difference equations
    m_fwd!(euler, Γ₅, Γ₆)

    # Forward difference equations: boundary conditions
    m_fwd!(E[:eq_dq1], Γ₅, Γ₆)
    Γ₅[E[:eq_dq1], η_k] = one(T)
    Γ₆[E[:eq_dq1], rk]  = one(T)

    m_fwd!(E[:eq_pq1], Γ₅, Γ₆)
    Γ₅[E[:eq_pq1], η_k] = one(T)
    Γ₆[E[:eq_pq1], q]   = one(T)
    Γ₆[E[:eq_pq1], rq]  = one(T)

    m_fwd!(E[:eq_ps₁1], Γ₅, Γ₆)
    Γ₅[E[:eq_ps₁1], a]  = one(T)
    Γ₆[E[:eq_ps₁1], π]  = convert(T, ϵ)
    Γ₆[E[:eq_ps₁1], s₁] = one(T)

    m_fwd!(E[:eq_ps₂1], Γ₅, Γ₆)
    Γ₅[E[:eq_ps₂1], a]  = one(T)
    Γ₆[E[:eq_ps₂1], π]  = convert(T, ϵ - 1.)
    Γ₆[E[:eq_ps₂1], s₂] = one(T)

    m_fwd!(E[:eq_pω1], Γ₅, Γ₆)
    Γ₆[E[:eq_pω1], c] = one(T)
    Γ₅[E[:eq_pω₁], a] = one(T)
    Γ₆[E[:eq_pω1], ω] = one(T)

    # Forward difference equations: recursions
    for i in 2:N_approx[:q]
        m_fwd!(E[Symbol("eq_dq$(i)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_dq$(i)")], η_k] = one(T)
        Γ₆[E[Symbol("eq_dq$(i)")], rq] = one(T)
        Γ₆[E[Symbol("eq_dq$(i)")], J[Symbol("dq$(i-1)")]] = one(T)

        m_fwd!(E[Symbol("eq_pq$(i)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_pq$(i)")], η_k] = one(T)
        Γ₆[E[Symbol("eq_pq$(i)")], rq] = one(T)
        Γ₆[E[Symbol("eq_pq$(i)")], J[Symbol("pq$(i-1)")]] = one(T)
    end

    for i in 2:N_approx[:s₁]
        m_fwd!(E[Symbol("eq_ds₁$(i-1)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_ds₁$(i-1)")], a] = one(T)
        Γ₆[E[Symbol("eq_ds₁$(i-1)")], π] = convert(T, ϵ)
        Γ₆[E[Symbol("eq_ds₁$(i-1)")], J[Symbol("ds₁$(i-2)")]] = one(T)

        m_fwd!(E[Symbol("eq_ps₁$(i)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_ps₁$(i)")], a] = one(T)
        Γ₆[E[Symbol("eq_ps₁$(i)")], π] = convert(T, ϵ)
        Γ₆[E[Symbol("eq_ps₁$(i)")], J[Symbol("ps₁$(i-1)")]] = one(T)
    end

    for i in 2:N_approx[:s₂]
        m_fwd!(E[Symbol("eq_ds₂$(i-1)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_ds₂$(i-1)")], a] = one(T)
        Γ₆[E[Symbol("eq_ds₂$(i-1)")], π] = convert(T, ϵ) - one(T)
        Γ₆[E[Symbol("eq_ds₂$(i-1)")], J[Symbol("ds₂$(i-2)")]] = one(T)

        m_fwd!(E[Symbol("eq_ps₂$(i)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_ps₂$(i)")], a] = one(T)
        Γ₆[E[Symbol("eq_ps₂$(i)")], π] = convert(T, ϵ) - one(T)
        Γ₆[E[Symbol("eq_ps₂$(i)")], J[Symbol("ps₂$(i-1)")]] = one(T)
    end

    for i in 2:N_approx[:ω]
        m_fwd!(E[Symbol("eq_dω$(i-1)")], Γ₅, Γ₆)
        Γ₆[E[Symbol("eq_dω$(i-1)")], c] = one(T)
        Γ₅[E[Symbol("eq_dω$(i-1)")], a] = one(T)
        Γ₆[E[Symbol("eq_dω$(i-1)")], J[Symbol("dω$(i-2)")]] = one(T)

        m_fwd!(E[Symbol("eq_pω$(i)")], Γ₅, Γ₆)
        Γ₆[E[Symbol("eq_pω$(i)")], c] = one(T)
        Γ₅[E[Symbol("eq_pω$(i)")], a] = one(T)
        Γ₆[E[Symbol("eq_pω$(i)")], J[Symbol("pω$(i-1)")]] = one(T)
    end

    ## Mapping from states to jump variables
    Ψ = zeros(T, Ny, Nz)

    ## Deterministic steady state as initial guess
    z, y = create_deterministic_ss_guess(m)

    return RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, ccgf, vec(z), vec(y), Ψ, Nε)
end

function create_deterministic_ss_guess(m::NKEZDisaster{T, SNK, NNK}) where {T <: Real, SNK, NNK}

    ## Set up

    # Get parameters
    @unpack β, γ, ψ, ν, ν̅, χ, δ, α, ϵ, θ, π_ss, ϕ_r, ϕ_π, ϕ_y = m
    @unpack χ_y, ρ_β, ρ_l, ρ_r, σ_β, σ_l, σ_r, μ_a, σ_a, κ_a = m
    @unpack disaster_occur_spec, disaster_intensity_spec, disaster_para = m
    r_ss = infer_r_ss(m)
    X̅    = infer_X̅(m)
    𝔼η_k = infer_𝔼η_k(m)

    # Unpack indexing dictionaries
    @unpack N_approx, S, J, E, SH = m
    @unpack k₋₁, logΔ₋₁, r₋₁, output₋₁, η_β, η_l, η_r, a, η_k = S
    @unpack output, c, l, v, ce, ω, ℓ, β̅, w, r = J
    @unpack π, q, x, rk, rq, mc, s₁, s₂, logΔ  = J

    ## Create guesses for deterministic steady state
    z = Vector{T}(undef, Nz)
    y = Vector{T}(undef, Ny)

    # AR(1) start at 0
    η_β0 = 0.
    η_l0 = 0.
    η_r0 = 0.

    # Disaster shock assumed to occur deterministically
    # and equals the unconditional expected value
    η_k0 = 𝔼η_k
    A0   = exp(κ_a * η_k0)

    # Variables known outright
    Ω0  = 1. / (1. - (β * A0 * exp(μ_a)) ^ (1. - ψ))
    V0  = ((1. - β) * Ω0) ^ (1. / (1. - ψ))
    𝒞ℰ0 = ((1. - β) / β * (Ω0 - 1.)) ^ (1. / (1. - ψ))
    M0  = β * (β * Ω0 / (Ω0 - 1.)) ^ ((ψ - γ) / (1. - ψ)) * (A0 * exp(μ_a)) ^ (-γ)
    R0  = exp(r_ss)
    Q0  = 1.
    Rq0 = 1 / η_k0 - X̅
    Rk0 = 1. / (M * exp(η_k0)) - Rq0
    expβ̅ = 1. - exp(η_β0) * β

    # Guesses
    L0 = .5548
    Δ0 = 1. # true if π_ss = 0, otherwise this is only a reasonable guess
    ℒ0 = (1. + (ψ - 1.) * exp(η_l0) * ν̅ * L0^(1. + ν) / (1. + ν))^(ψ / (1. - ψ))

    # Implied values given guesses
    C0_fnct = Cin -> Cin[1] + X̅ * (α / (1. - α) * ψ * ν̅ * C0 * L0^ν / ℒ0 / RK0 * L0) -
        ((α / (1. - α) * * ψ * ν̅ * C0 * L0^ν / ℒ0 / RK0) ^ α * L0 - χ_y) / Δ0
    C0_guess = NaN
    for theguess in .5:.5:10.
        try
            C0_fnct([theguess])
            C0_guess = theguess
        catch e
        end
    end
    C0 = nlsolve(C0_fnct, [C0_guess]).zero[1]
    W0 = ψ * exp(η_l0) * ν̅ * C0 * L0^ν / ℒ0^((1. - ψ) / ψ)
    MC0 = (1. / (1. - α)) ^ (1. - α) * (1. / α) ^ α * W0 ^ (1. - α) * RK0 ^ α
    K0  = (α / (1. - α) * W0 / RK0 * L0) / η_k0
    X0  = X̅ * η_k0 * K0
    Y0  = ((η_k0 * K0) ^ α * L0 ^ (1. - α) - χ_y) / Δ0
    Π0  = exp(π_ss)
    S₁0 = MC0 * Y0 / (1. - exp(μ_a) * θ * M0 * A0 * Π0 ^ ϵ)
    S₂0 = Y0 / (1. - exp(μ_a) * θ * M0 * A0 * Π0 ^ (ϵ - 1.))
    z .= [convert(T, x) for x in log.([K0, Δ0, R0, Y0, exp.([η_β0, η_l0, η_r0, log(A0), η_k0])...])]
    y[1:19] = [convert(T, x) for x in log.([Y0, C0, L0, V0, 𝒞ℰ0, Ω0, ℒ0, expβ̅, W0, R0, Π0, Q0, X0, Rk0, Rq0,
                                            MC0, S₁0, S₂0, Δ0])]

    y[J[:dq1]] = convert(T, log(M0 * Rk0))
    y[J[:pq1]] = convert(T, log(Rq0 * M0 * Q0))
    y[J[:ds₁0]] = convert(T, log(MC0 * Y0))
    y[J[:ps₁1]] = convert(T, log(exp(μ_a) * θ * M0 * A0 * Π0^ϵ * S₁0))
    y[J[:ds₂0]] = convert(T, log(Y0))
    y[J[:ps₂1]] = convert(T, log(exp(μ_a) * θ * M0 * A0 * Π0^(ϵ - 1.) * S₂0))

    # NEED TO ADD GUESSES FOR omega

    for i in 2:N_approx
        y[J[Symbol("dq$(i)")]] = convert(T, log(M0) + η_k0 + log(Rq0) + y[J[Symbol("dq$(i-1)")]])
        y[J[Symbol("pq$(i)")]] = convert(T, log(M0) + η_k0 + log(Rq0) + y[J[Symbol("pq$(i-1)")]])
        y[J[Symbol("ds₁$(i-1)")]] = convert(T, μ_a + log(θ) + log(M0) + log(A0) + ϵ * π_ss + y[J[Symbol("ds₁$(i-2)")]])
        y[J[Symbol("ps₁$(i)")]] = convert(T, μ_a + log(θ) + log(M0) + log(A0) + ϵ * π_ss + y[J[Symbol("ps₁$(i-1)")]])
        y[J[Symbol("ds₂$(i-1)")]] = convert(T, μ_a + log(θ) + log(M0) + log(A0) + (ϵ - 1.) * π_ss + y[J[Symbol("ds₂$(i-2)")]])
        y[J[Symbol("ps₂$(i)")]] = convert(T, μ_a + log(θ) + log(M0) + log(A0) + (ϵ - 1.) * π_ss + y[J[Symbol("ps₂$(i-1)")]])
    end
end

# Infer the value of η_k in the stochastic steady state
function infer_𝔼η_k(m::NKEZDisaster)
    d = m.disaster_para
    𝔼η_k = if m.disaster_occur_spec == :PoissonNormalMixture
        # η_{k, t} ∼ N(-jₜ, jₜ σ_k^2), jₜ ∼ Poisson(pₜ₋₁)
        # ⇒ 𝔼[ η_{k, t} ] = 𝔼[ 𝔼[η_{k, t} ∣ jₜ] ] = 𝔼[ -jₜ ]
        # = -𝔼[ 𝔼[ jₜ ∣ pₜ₋₁] ] = -𝔼[ pₜ₋₁ ].
        if m.disaster_intensity_spec == :CoxIngersollRoss
            # pₜ₋₁ ∼ discretized CIR process w/unconditional mean p
            # -𝔼[ pₜ₋₁ ] = -p
            -d[:p]
        elseif m.disaster_intensity_spec == :TwoStateMarkovChain
            # pₜ₋₁ ∼ Markov Chain with states p_ and p̅;
            # (respective) persistence probabilities ρ_ and ρ̅
            # -𝔼[ pₜ₋₁ ] = -(ergodic mean)
            -((1. - d[:ρ̅ₚ]) * d[:p_] + (1. - d[:ρ_ₚ] * d[:p̅])) / (2. - (d[:ρ_ₚ] + d[:ρ̅ₚ]))
        end
    elseif m.disaster_occur_spec == :Bernoulli
        # η_{k, t} ∼ Bernoulli(pₜ₋₁) taking values η_ w/probability pₜ₋₁ and zero otherwise.
        # ⇒ 𝔼[ η_{k, t} ] = 𝔼[ 𝔼[ η_{k, t} ∣ pₜ₋₁ ] ] = η_ 𝔼[ pₜ₋₁]
        if m.disaster_intensity_spec == :CoxIngersollRoss
            # pₜ₋₁ ∼ discretized CIR process w/unconditional mean p
            # η_ 𝔼[ pₜ₋₁ ] = η_ p
            d[:η_] * d[:p]
        elseif m.disaster_intensity_spec == :TwoStateMarkovChain
            # pₜ₋₁ ∼ Markov Chain with states p_ and p̅;
            # (respective) persistence probabilities ρ_ and ρ̅
            # η_ 𝔼[ pₜ₋₁ ] = η_ (ergodic mean)
            d[:η_] * ((1. - d[:ρ̅ₚ]) * d[:p_] + (1. - d[:ρ_ₚ] * d[:p̅])) / (2. - (d[:ρ_ₚ] + d[:ρ̅ₚ]))
        end
    end

    return 𝔼η_k
end

# Infer steady state investment rate given the disaster shock specification
function infer_X̅(m::NKEZDisaster)
    return m.χ / (m.χ + 1.) * (1. / exp(infer_𝔼η_k(m)) + m.δ - 1.)
end

# Figure out the steady state interest rate
# given EZ preferences and the disaster shock specification
function infer_r_ss(m::NKEZDisaster)

    # a = σ_a ε_a + κ_a * η_k
    # ⇒ in stochastic steady state, a = κ_a 𝔼[ η_k ]
    𝔼η_k = infer_𝔼η_k(m)
    A    = exp(m.κ_a * 𝔼η_k)

    # Stochastic steady state is the expected state,
    # conditional on shocks always equaling zero.
    Ω̃    = 1. / (1. - (m.β * A * exp(m.μ_a))^(1. - m.ψ))
    M    = m.β * (m.β * Ω̃ / (Ω̃ - 1.))^((ψ - γ) / (1. - ψ))
    return m.π_ss - log(M)
end

# Figure out the ccgf given the disaster shock specification
function infer_ccgf(m::NKEZDisaster)
    function ccgf(F, α, z)
        # F .= .5 * RiskAdjustedLinearizations.diag(α * α') # slower but this is the underlying math
        F .= vec(.5 * sum(α.^2, dims = 2)) # faster implementation
    end
end

# Infer state transition equations for
# the disaster shock occurrence η_k
function infer_μ_disaster_occur(m::NKEZDisaster)
    d = m.disaster_para

    # Define expected disaster shock proportion to pₜ b/c conditionally linear in pₜ
    𝔼ₜη_k_div_pₜ = if m.disaster_occur_spec == :PoissonNormalMixture
        # η_{k, t} ∼ N(-jₜ, jₜ σ_k^2), jₜ ∼ Poisson(pₜ₋₁)
        # ⇒ 𝔼ₜ[ η_{k, t + 1} ] = 𝔼ₜ[ 𝔼[η_{k, t + 1} ∣ jₜ₊₁] ] = 𝔼ₜ[ -jₜ₊₁ ]
        # = -𝔼ₜ[ 𝔼[ jₜ₊₁ ∣ pₜ] ] = -𝔼ₜ[ pₜ ] = -pₜ.
        -1.
    elseif m.disaster_occur_spec == :Bernoulli
        # η_{k, t} ∼ Bernoulli(pₜ₋₁) taking values η_ w/probability pₜ₋₁ and zero otherwise.
        # ⇒ 𝔼ₜ[ η_{k, t + 1} ] = 𝔼ₜ[ 𝔼[ η_{k, t + 1} ∣ pₜ ] ] = η_ 𝔼ₜ[ pₜ] = η_ pₜ
        d[:η_]
    end

    𝔼ₜη_k = if m.disaster_intensity_spec in [:CoxIngersollRoss, :TwoStateMarkovChain, :TruncatedCoxIngersollRoss]
        state_i = m.S[:p]
        _𝔼ₜη_k_linear(z, y) = 𝔼ₜη_k_div_pₜ * z[state_i]
    elseif m.disaster_intensity_spec in [:LogAR1]
        state_i = m.S[:logp]
        _𝔼ₜη_k_loglinear(z, y) = 𝔼ₜη_k_div_pₜ * exp(z[state_i])
    end

    return 𝔼ₜη_k
end

# Infer state transition equations for
# the disaster shock intensity p
function infer_μ_disaster_intensity(m::NKEZDisaster)
    d = m.disaster_para
    mdisi = m.disaster_intensity_spec

    μ_p = if mdisi == :CoxIngersollRoss
        state_i = m.S[:p]
        @inline _μ_p_cir(z, y) = (1. - d[:ρ_p]) * d[:p] + d[:ρ_p] * z[state_i]
    elseif mdisi == :TwoStateMarkovChain
        state_i = m.S[:p]
        @inline function _μ_p_2mc(z, y)
            if z[state_i] == d[:p̅]
                d[:ρ̅_p] * d[:p̅] + (1. - d[:ρ̅_p]) * d[:p_]
            else
                d[:ρ_p] * d[:p_] + (1. - d[:ρ_p]) * d[:p̅]
            end
        end
    elseif mdisi == :TruncatedCoxIngersollRoss
        error("TruncatedCoxIngersollRoss not implemented yet")
    elseif mdisi == :LogAR1
        state_i = m.S[:logp]
        @inline _μ_p_logar1(z, y) = (1 - d[:ρ_p]) * log(d[:p]) + d[:ρ_p] * z[state_i]
    end

    return μ_p
end

function infer_Σ_disaster_intensity(m::NKEZDisaster)
    d = m.disaster_para
    mdisi = m.disaster_intensity_spec

    Σ_p = if mdisi in [:CoxIngersollRoss, :TruncatedCoxIngersollRoss]
        state_i = m.S[:p]
        @inline _Σ_p_cir(z) = sqrt(z[state_i]) * d[:σ_p]
    elseif mdisi == :TwoStateMarkovChain
        @inline _Σ_p_2mc(z) = one(eltype(z))
    elseif mdisi == :LogAR1
        state_i = m.S[:logp]
        @inline _Σ_p_logar1(z) = d[:σ_p]
    end

    return Σ_p
end

# Infer the desired CCGF function
function infer_ccgf(m::NKEZDisaster)
    d = m.disaster_para
    S = m.S
    SH = m.SH
    not_dis_keys = setdiff(collect(keys(SH)), [:ε_k, :ε_p])
    not_dis_inds = [SH[i] for i in not_dis_keys]
    ccgf = if m.disaster_occur_spec == :PoissonNormalMixture
        # apply Poisson mgf to C_2(A) = -A + σ_j^2 A^2 / 2
        # Poisson mgf w/intensity pₜ₋₁ is exp((exp(s) - 1) pₜ₋₁)
        # and then subtract s * E_t[\eta_{k, t + 1}]

        if m.disaster_intensity_spec == :CoxIngersollRoss
            function _ccgf_poissonnormalmixture_cir(F, A, z)
                F  .= sum(A[:, vcat(not_dis_inds, SH[:ε_p])].^2, dims = 2) .* .5 # Gaussian parts
                A_k = @view A[:, SH[:ε_k]]
                F .+= ((exp.(-A_k + A_k.^2 .* (d[:σ_k] ^ 2 / 2.)) .- 1.) + A_k) .* z[S[:p]] # ε_k
            end
        elseif m.disaster_intensity_spec == :TwoStateMarkovChain
            function _ccgf_poissonnormalmixture_2smc(F, A, z)
                F .= sum((@view A[:, not_dis_inds]).^2, dims = 2) .* .5 # Gaussian parts
                A_k = @view A[:, SH[:ε_k]]
                F .+= ((exp.(-A_k + A_k.^2 .* (d[:σ_k] ^ 2 / 2.)) .- 1.) + A_k) .* z[S[:p]] # ε_k

                # ε_p
                A_p = @view A[:, SH[:ε_p]]
                if z[S[:p]] == d[:p̅]
                    F .+= log.((1. - d[:ρ̅_p]) .* exp.(A_p * d[:p_]) + d[:ρ̅_p] .* exp.(A_p * d[:p̅])) -
                        A_p .* (d[:ρ̅_p] * d[:p̅] + (1. - d[:ρ̅_p]) * d[:p_])
                else
                    F .+= log.((1. - d[:ρ_p]) .* exp.(A_p .* d[:p̅]) + d[:ρ_p] .* exp.(A_p .* d[:p_])) -
                        A_p .* (d[:ρ_p] * d[:p_] + (1. - d[:ρ_p]) * d[:p̅])
                end
            end
        elseif m.disaster_intensity_spec == :LogAR1
            function _ccgf_poissonnormalmixture_logar1(F, A, z)
                F .= sum(A[:, vcat(not_dis_inds, SH[:ε_p])].^2, dims = 2) .* .5 # Gaussian parts
                A_k = @view A[:, SH[:ε_k]]
                F .+= ((exp.(-A_k + A_k.^2 .* (d[:σ_k] ^ 2 / 2.)) .- 1.) + A_k) .* exp(z[S[:logp]]) # ε_k
            end
        end
    elseif m.disaster_occur_spec == :Bernoulli
        if m.disaster_intensity_spec == :CoxIngersollRoss
            function _ccgf_bernoulli_cir(F, A, z)
                F .= sum(A[:, vcat(not_dis_inds, SH[:ε_p])].^2, dims = 2) .* .5 # Gaussian parts
                A_k = @view A[:, SH[:ε_k]]
                F .+= log((1. - z[S[:p]] + z[S[:p]]) .* exp.(A_k * d[:η_])) .-
                    A_k .* (d[:η_] * z[S[:p]]) # ε_k
            end
        elseif m.disaster_intensity_spec == :TwoStateMarkovChain
            function _ccgf_bernoulli_2smc(F, A, z)
                F .= sum((@view A[:, not_dis_inds]).^2, dims = 2) .* .5 # Gaussian parts
                A_k = @view A[:, SH[:ε_k]]
                F .+= log.((1. - z[S[:p]] + z[S[:p]]) .* exp.(A_k * d[:η_])) .-
                    A_k .* (d[:η_] * z[S[:p]]) # ε_k

                # ε_p
                A_p = @view A[:, SH[:ε_p]]
                if z[S[:p]] == d[:p̅]
                    F .+= log.((1. - d[:ρ̅_p]) .* exp.(A_p * d[:p_]) + d[:ρ̅_p] .* exp.(A_p * d[:p̅])) -
                        A_p .* (d[:ρ̅_p] * d[:p̅] + (1. - d[:ρ̅_p]) * d[:p_])
                else
                    F .+= log.((1. - d[:ρ_p]) .* exp.(A_p * d[:p̅]) + d[:ρ_p] .* exp.(A_p * d[:p_])) -
                        A_p .* (d[:ρ_p] * d[:p_] + (1. - d[:ρ_p]) * d[:p̅])
                end
            end
        elseif m.disaster_intensity_spec == :LogAR1
            function _ccgf_bernoulli_logar1(F, A, z)
                F .= sum(A[:, vcat(not_dis_inds, SH[:ε_p])].^2, dims = 2) .* .5 # Gaussian parts
                A_k = @view A[:, SH[:ε_k]]
                F .+= log((1. - exp(z[S[:logp]]) + exp(z[S[:logp]])) .* exp.(A_k * d[:η_])) .-
                    A_k .* (d[:η_] * exp(z[S[:logp]])) # ε_k
            end
        end
    end

    if isnothing(ccgf)
        error("Either the specification of the disaster shock's occurrence ($(m.disaster_occur_spec)) or intensity " *
              "$(m.disaster_intensity_spec) is not recognized.")
    else
        return ccgf
    end
end


text:
module Transforms
export ColTrans, datasource, mergecols, onehot, onecold, replacevals, standardize, fit!, transform, invert, collecttrans,
    deserialize, tonumerical, copytransform, inids, outids, intype, outtype, getdatasource, invtransform

using DataFrames
using DataStructures: counter, OrderedDict, OrderedSet
using JSON
using MLDataUtils
using Statistics

import ..MLGoodies: toposort, dview, ncols

# TODO different meta and IO options?
customshow(io, e) = show(io, e)
customshow(io, e::AbstractString) = print(io, "\"$e\"")
customshow(io, p::Pair) = begin customshow(io, first(p)); print(io, " => "); customshow(io, last(p)) end

function customshow(io, v::AbstractVector, brackets = true)
    brackets && print(io, "[")
    for (i, e) in enumerate(v)
        i > 1 && print(io, ", ")
        customshow(io, e)
    end
    brackets && print(io, "]")
end

customshow(io, v::AbstractDict, brackets = true) = customshow(io, collect(v), brackets)

_smap_process(smap) = smap isa AbstractDict ? smap : Dict(:default => smap)

# ----- ColTrans -----
# --------------------
"Abstract Column transform."
abstract type ColTrans end

# ----- Connections -------
# -------------------------

"""
Describes input connection by its source transform `trn` and possibly its `type`.
"""
mutable struct Inp
    trn::ColTrans
    type::Type
    Inp(trn, type = Nothing) = new(trn, type)
end

"""
Describes output connection assigning it a `type`.
"""
mutable struct Out
    type::Type
    Out(type = Nothing) = new(type)
end

"""
Maps id to a transformation.
"""
const IdColTrans = Pair{Symbol,<:ColTrans}

"""
Maps input id to the input connection (`Inp`).
"""
const IdInp = Pair{Symbol,Inp}

IdInp(i::IdInp) = i
IdInp(i::IdColTrans) = first(i) => Inp(last(i))

"""
Maps output id to the output connection (`Out`).
"""
const IdOut = Pair{Symbol,Out}

IdOut(o::IdOut) = o
IdOut(o::IdColTrans) = first(o) => Out()

"""
Manage transform node input and output connections.
"""
mutable struct Connections
    ins::OrderedDict{Symbol,Inp}
    outs::OrderedDict{Symbol,Out}
    maxins::Int
    maxouts::Int
    function Connections(ins::OrderedDict{Symbol,Inp}, outs::OrderedDict{Symbol,Out}; maxins = typemax(Int), maxouts = typemax(Int))
        length(ins) > maxins && throw(ArgumentError("too many inputs $(length(ins)) > $(maxins)!"))
        length(outs) > maxouts && throw(ArgumentError("too many outputs $(length(outs)) > $(maxouts)!"))
        new(ins, outs, maxins, maxouts)
    end
end

function Connections(ins = IdInp[], outs = IdOut[]; kwargs...)
    length(ins) == length(unique(map(first, ins))) || error("duplicate input ids!")
    length(outs) == length(unique(map(first, outs))) || error("duplicate output ids!")
    Connections(OrderedDict(IdInp[IdInp(i) for i in ins]), OrderedDict(IdOut[IdOut(o) for o in outs]); kwargs...)
end

function Base.push!(c::Connections, id::Symbol, inp::Inp)
    length(c.ins) >= c.maxins && throw(ArgumentError("more than $(c.maxins) inputs not allowed!"))
    id ∈ keys(c.ins) && throw(ArgumentError("input with id $(id) already exists!"))
    push!(c.ins, id => inp)
end

function Base.push!(c::Connections, id::Symbol, out::Out)
    length(c.outs) >= c.maxouts && throw(ArgumentError("more than $(c.maxouts) outputs not allowed!"))
    id ∈ keys(c.outs) && throw(ArgumentError("output with id $(id) already exists!"))
    push!(c.outs, id => out)
end

inids(c::Connections) = collect(keys(c.ins))
intrns(c::Connections) = [v.trn for v in values(c.ins)]
outids(c::Connections) = collect(keys(c.outs))

function Base.show(io::IO, c::Connections)
    inrule = (r)->"$(first(r))(#$(last(r).trn.id))"
    if length(c.ins) > 0
        print(io, join(inrule.(collect(c.ins)), ", "))
    else
        print(io, "empty")
    end
end

function Base.getindex(c::Connections, e::Symbol)
    for (col, tran) in c.in
        col == e && return tran
    end
    throw(KeyError(e))
end

# ----- ColTrans -----
# --------------------
coltranscounter = Threads.Atomic{Int}(1)
_newid() = Threads.atomic_add!(coltranscounter, 1)

inids(c::ColTrans) = inids(c.con)
intrns(c::ColTrans) = intrns(c.con)
intrn(c::ColTrans, id::Symbol) = c.con.ins[id].trn
intype(c::ColTrans, id::Symbol) = c.con.ins[id].type
outids(c::ColTrans) = outids(c.con)
outtype(c::ColTrans, id::Symbol) = c.con.outs[id].type

idstr(t::ColTrans) = "#$(t.id)"
instancestr(t::ColTrans) = string(typeof(t).name.name) * "($(idstr(t)))"

"""
    isfit(t::ColTrans)

Find if the transform `t` is fit.
"""
isfit(t::ColTrans) = true

"""
    isfitrec(t::ColTrans)

Find recursively if the transform `t` is fit.
"""
isfitrec(t::ColTrans) = all((isfit(t), isfitrec.(intrns(t))...))

"""
     checkfit(t::ColTrans)

Throw error if the transform `t` is not fit.
"""
checkfit(t::ColTrans) = isfit(t) || error("$(t) not initialized, use fit!")

"""
     checkfitrec(t::ColTrans)

Throw error if the transform `t` or any of its predecessors is not fit.
"""
checkfitrec(t::ColTrans) = isfitrec(t) || error("$(t) not initialized, use fit!")

"""
     check(t::ColTrans)

Check transform `t` validity and throw error if needed.

Currently checking for more than one `DataSource`. More to come...
"""
function check(t::ColTrans)
    ts = typeof.(collecttrans(t))
    cnts = counter(ts)
    if cnts[DataSource] > 1
        error("More than one DataSource!")
    end
end

function _fit_helper!(c::Connections, t::ColTrans, smap)
    smap = _smap_process(smap)
    #TODO optimize! Some calls are now performed multiple times!
    data = OrderedDict() # maps predecessor transform to it output
    cache = Dict{ColTrans,Any}()
    for (id, input) in c.ins # fit all inputs recursively storing results to the cache
        if input.trn ∉ keys(cache)
            cache[input.trn] = fit!(input.trn, smap)
        end
        data[id] = cache[input.trn][id]
    end
    ret = fit!(t, data, smap) # compute output of this transform
    for (i, d) in zip(values(c.ins), values(data)) # fix input types based on output of predecessors
        if i.type == Nothing
            i.type = eltype(d)
        end
    end
    for (o, r) in zip(values(c.outs), values(ret)) # set ouput types
        if o.type == Nothing
            o.type = eltype(r)
        end
    end
    ret
end

"""
    fit!(t::ColTrans, smap)

Fit all transforms rooted at `t` using the source map `smap`.
"""
fit!(t::ColTrans, smap) = _fit_helper!(t.con, t, smap)

function _checktransform(t::ColTrans, smap)
    checkfit(t)
    length(smap) > 0 || error("no source!")
    ns = [nobs(s) for s in values(smap)]
    minimum(ns) != maximum(ns) && throw(ArgumentError("source sizes not same; $ns"))
end

"""
    transform(t::ColTrans, smap, ::Type{<:AbstractDataFrame})

Apply transform `t` to source map `smap` to create a new DataFrame.

# Example
```jldoctest
julia> using DataFrames;

julia> df = DataFrame(name = ["car", "bus", "cat"], category = [1, 1, 2], weight = [1000, 10000, 5.0])
3×3 DataFrames.DataFrame
│ Row │ name   │ category │ weight  │
│     │ String │ Int64    │ Float64 │
├─────┼────────┼──────────┼─────────┤
│ 1   │ car    │ 1        │ 1000.0  │
│ 2   │ bus    │ 1        │ 10000.0 │
│ 3   │ cat    │ 2        │ 5.0     │

julia> transform(tonumerical(df), df, DataFrame)
3×5 DataFrames.DataFrame
│ Row │ name_bus │ name_car │ name_cat │ category │ weight    │
│     │ Float64  │ Float64  │ Float64  │ Float64  │ Float64   │
├─────┼──────────┼──────────┼──────────┼──────────┼───────────┤
│ 1   │ -0.57735 │ 1.1547   │ -0.57735 │ -0.57735 │ -0.484631 │
│ 2   │ 1.1547   │ -0.57735 │ -0.57735 │ -0.57735 │ 1.14998   │
│ 3   │ -0.57735 │ -0.57735 │ 1.1547   │ 1.1547   │ -0.665346 │
```
"""
function transform(t::ColTrans, smap, ::Type{<:AbstractDataFrame})
    smap = _smap_process(smap)
    _checktransform(t, smap)
    os = outids(t)
    data = [t[o](smap) for o in os]
    DataFrame(data, os)
end

"""
    transform(t::ColTrans, smap, type::Type{<:AbstractMatrix})

Apply transform `t` to source map `smap` to create a new Matrix.

# Example
```jldoctest
julia> using DataFrames;

julia> df = DataFrame(name = ["car", "bus", "cat"], category = [1, 1, 2], weight = [1000, 10000, 5.0])
3×3 DataFrames.DataFrame
│ Row │ name   │ category │ weight  │
│     │ String │ Int64    │ Float64 │
├─────┼────────┼──────────┼─────────┤
│ 1   │ car    │ 1        │ 1000.0  │
│ 2   │ bus    │ 1        │ 10000.0 │
│ 3   │ cat    │ 2        │ 5.0     │

julia> transform(tonumerical(df), df, Matrix{Float32})
5×3 Array{Float32,2}:
 -0.57735    1.1547   -0.57735
  1.1547    -0.57735  -0.57735
 -0.57735   -0.57735   1.1547
 -0.57735   -0.57735   1.1547
 -0.484631   1.14998  -0.665346
```
"""
function transform(t::ColTrans, smap, type::Type{<:AbstractMatrix})
    smap = _smap_process(smap)
    _checktransform(t, smap)
    n = nobs(first(values(smap)))
    os = outids(t)
    M = type(undef, length(os), n)
    for (i, o) in enumerate(os)
        M[i,:] .= t[o](smap)
    end
    M
end

# ----- DataSource -----
# ----------------------
struct DataSource <: ColTrans
    id::Int
    con::Connections # no inputs here
    cols::OrderedDict{Symbol,Function} # column id to function which reads the data based on source map
    sids::OrderedDict{Symbol,Symbol} # column id to sid
    function DataSource(con::Connections, cols::OrderedDict{Symbol,Function}, sids::OrderedDict{Symbol,Symbol})
        con.maxins != 0 && throw(ArgumentError("maxins != 0!"))
        new(_newid(), con, cols, sids)
    end
end

Base.getindex(t::DataSource, col::Symbol) = t.cols[col]

function Base.show(io::IO, t::DataSource)
    print(io, instancestr(t) * "[" * join(["$col($(t.sids[col]))" for col in keys(t.cols)], ", ") * "]")
end

"""
    datasource(s::AbstractDataFrame; sid::Symbol=:default)

Create a `DataSource` with source id `sid` based on a DataFrame `s`.
"""
datasource(s::AbstractDataFrame; sid::Symbol=:default) = datasource(s, Symbol.(names(s)); sid=sid)

"""
    datasource(s, names::Vector{Symbol}; sid::Symbol=:default)

Create a `DataSource` with source id `sid` based on a source object `s`.

The column names are given by `names`.  Method `ncols(s)` must be defined for the object.
"""
function datasource(s, names::Vector{Symbol}; sid::Symbol=:default)
    ncols(s) != length(names) && error("source $(typeof(s)) size $(size(s)) does not match the number of col names ($(length(names)))")
    colsid = OrderedDict{Symbol,Symbol}(col => sid for col in names)
    datasource(colsid)
end

"""
    datasource(colsid::OrderedDict)

Create a `DataSource` based on ordered dictionary of column id => source id pairs.

Data container must support `dview` method and access by column indices.
"""
function datasource(colsid::OrderedDict)
    con = Connections(IdInp[], [col => Out() for col in keys(colsid)]; maxins = 0)
    colsid = OrderedDict{Symbol,Symbol}(Symbol(col) => Symbol(sid) for (col, sid) in colsid)
    colfunc = OrderedDict{Symbol,Function}(col => ((i)->dview(i[sid], :, colidx)) for (colidx, (col, sid)) in enumerate(colsid))
    DataSource(con, colfunc, colsid)
end

function fit!(t::DataSource, data, smap)
    # nothing to fit here just give column id => column data dict
    OrderedDict(col => t[col](smap) for col in outids(t))
end

"""
    getdatasource(t::ColTrans)

Finds the only DataSource in the transformation DAG rooted in `t`.
"""
function getdatasource(t::ColTrans)
    filter(trn->typeof(trn) == DataSource, collecttrans(t)) |> first
end

# ----- MergeCols -----
# ---------------------
struct MergeCols <: ColTrans
    id::Int
    con::Connections
    function MergeCols(con::Connections)
        new(_newid(), con)
    end
end

"""
     mergecols(ins...)

Merges all input id to input connection mappings in `ins`.

No id duplicates are allowed.
"""
function mergecols(ins...)
    outs = IdOut[first(i) => Out() for i in ins]
    MergeCols(Connections(ins, outs))
end

"""
     mergecols(cs::ColTrans...)

Merges outputs of all transforms `cs`.

All transforms must be fitted or error is thrown. Also no output id duplicates are allowed.
"""
function mergecols(cs::ColTrans...)
    ins = IdInp[]
    cols = Set{Symbol}()
    for c in cs
        checkfitrec(c)
        for outid in outids(c)
            outid ∈ cols && error("duplicate symbols: $(outid)")
            push!(ins, outid => Inp(c))
            push!(cols, outid)
        end
    end
    mergecols(ins...)
end

function Base.show(io::IO, t::MergeCols)
    print(io, instancestr(t) * "[" * repr(t.con) * "]")
end

"""
    Base.getindex(t::MergeCols, col)

Give a column evaluation function for `col`.

The column evaluation function (CEF) expects a source map parameter (`smap`) to give actual results.
"""
function Base.getindex(t::MergeCols, col)
    fsrc = intrn(t, col)[col]
    (i)->func(t, col, OrderedDict(col => fsrc(i)))
end

"""
    func(t::MergeCols, col, data)

Get column evaluation function for `col` based on input `data`.
"""
func(t::MergeCols, col, data) = data[col]

"""
    fit!(t::MergeCols, data, smap)

Nothing to fit for `mergecols`. Just return a dict of column evaluation functions.
"""
function fit!(t::MergeCols, data, smap)
    OrderedDict(col => func(t, col, data) for col in outids(t))
end

# ----- One Hot -----
# -------------------
struct OneHot <: ColTrans
    id::Int
    con::Connections
    tbl::OrderedDict{Symbol,Any} # output column id to value
    function OneHot(con::Connections, tbl::OrderedDict{Symbol,Any})
        con.maxins > 1 && throw(ArgumentError("OneHot: only single input allowed!"))
        new(_newid(), con, tbl)
    end
end

"""
    onehot(in, tblpairs::Pair...)

Define one hot encoding transformation.

Input is given as: `in = icol => itrans`, where `icol` is the output column of the input
transform `itrans`. Transform table `tblpairs` can be either given or automatically generated later using by `fit!()`.
If generated, the output columns are denoted `icol_v1`, `icol_v2`, ..., where `v`i denotes possible values of `icol`
(sorted when fitting). Default output type is `Bool`.
"""
function onehot(in, tblpairs::Pair...)
    OneHot(Connections([in], []; maxins = 1), OrderedDict{Symbol,Any}(tblpairs...))
end

isfit(t::OneHot) = length(t.tbl) > 0

function Base.show(io::IO, t::OneHot)
    print(io, instancestr(t) * "[" * repr(t.con) * "|")
    if length(t.tbl) == 0
        print(io, "not initialized!]")
    else
        customshow(io, t.tbl, false)
    end
    print(io, "]")
end

function Base.getindex(t::OneHot, col)
    checkfit(t)
    icol = first(inids(t)) # only possible input id
    fsrc = intrn(t, icol)[icol] # get CEF for this column
    (i)->func(t, col, fsrc(i)) # compute output, if generated the the col is one of icol_VAL1, icol_VAL2, ...
end

function func(t::OneHot, col, data)
    col ∉ keys(t.tbl) && ArgumentError("unknown output column $col")
    T = outtype(t, col)
    T = T == Nothing ? Bool : T
    T.(data .== t.tbl[col])
end

function fit!(t::OneHot, data, smap)
    icol = first(inids(t))
    if !isfit(t)
        unq = sort(unique(data[icol]))
        for val in unq
            newcol = Symbol(string(icol) * "_" * string(val))
            t.tbl[newcol] = val
            push!(t.con, newcol, Out())
        end
    end
    OrderedDict(col => func(t, col, data[icol]) for col in outids(t))
end

# ----- One Cold ----
# -------------------
struct OneCold <: ColTrans
    id::Int
    con::Connections
    tbl::OrderedDict{Symbol,Any} # input column id to output value
    function OneCold(con::Connections, tbl::OrderedDict{Symbol,Any})
        con.maxouts > 1 && throw(ArgumentError("OneCold: only single output allowed!"))
        new(_newid(), con, tbl)
    end
end

"""
    onecold(ins, outid::Symbol, tbl::OrderedDict{Symbol,Any})

Define one cold transform.

# Examples
The following transform finds maximum in `itrans` columns `:I1` and `:I2` and outputs
either `"a"` (for maximum in `:I1`) or `2` (for maximum in `:I2`) to column `:O`.

```julia
onecold([:I1 => itrans, :I2 => itrans], :O, OrderedDict(:I1 => "a", :I2 => 2))
```
"""
function onecold(ins, outid::Symbol, tbl::OrderedDict{Symbol,Any})
    OneCold(Connections(ins, [outid => Out()]; maxouts = 1), tbl)
end

function Base.show(io::IO, t::OneCold)
    print(io, instancestr(t) * "[" * repr(t.con) * "|$(outids(t)[1]): ")
    customshow(io, t.tbl)
    print(io, "]")
end

function Base.getindex(t::OneCold, col)
    col != first(outids(t)) && error("OneCold: unknown output column $(col) != $(first(outids(t)))")
    fsrcs = [intrn(t, icol)[icol] for icol in inids(t)]
    (i)->func(t, [fsrc(i) for fsrc in fsrcs])
end

function func(t::OneCold, data)
    vals = [t.tbl[icol] for icol in inids(t)]
    [vals[argmax(obs)] for obs in zip(data...)]
end

function fit!(t::OneCold, data, smap)
    OrderedDict(col => func(t, data) for col in outids(t))
end

# ----- ReplaceVals -----
# -----------------------
struct ReplaceVals <: ColTrans
    id::Int
    con::Connections
    tbl::OrderedDict
    function ReplaceVals(con::Connections, tbl::OrderedDict)
        con.maxins > 1 && throw(ArgumentError("ReplaceVals: only single input allowed!"))
        con.maxouts > 1 && throw(ArgumentError("ReplaceVals: only single output allowed!"))
        new(_newid(), con, tbl)
    end
end

"""
    replacevals(in, tblpairs::Pair...)

Replace values of a column.

Input `in` is given as: `in = icol => itrans`, where `icol` is the output column of the input
transform `itrans`. Rewriting rules are given as `Pair`s. Nonlisted values are kept untouched.

# Examples

The following transforms `itrans` column `:C` values `X` and `Z` to the respective lowercase letters.
```julia
replacevals(:C => itrans, "X" => "x", "Z" => "z")
```
"""
function replacevals(in, tblpairs::Pair...)
    ReplaceVals(Connections([in], [first(in) => Out()]; maxins = 1, maxouts = 1), OrderedDict(tblpairs...))
end

function Base.show(io::IO, t::ReplaceVals)
    print(io, instancestr(t) * "[" * repr(t.con) * "|")
    customshow(io, t.tbl, false)
    print(io, "]")
end

function func(t::ReplaceVals, data)
    map(e->get(t.tbl, e, e), data)
end

function Base.getindex(t::ReplaceVals, col)
    fsrc = intrn(t, col)[col] # same name of icol and col
    (i)->func(t, fsrc(i))
end

function fit!(t::ReplaceVals, data, smap)
    OrderedDict(col => func(t, data[col]) for col in outids(t))
end

# ----- Standardize ----
# ----------------------
mutable struct Standardize <: ColTrans
    id::Int
    con::Connections
    μ::AbstractFloat
    σ::AbstractFloat
    function Standardize(con::Connections, μ::AbstractFloat, σ::AbstractFloat)
        con.maxins > 1 && throw(ArgumentError("Standardize: only single input allowed!"))
        con.maxouts > 1 && throw(ArgumentError("Standardize: only single output allowed!"))
        isnan(σ)  || σ > 0 || throw(ArgumentError("Standardize: σ ≦ 0!: $(σ)"))
        new(_newid(), con, μ, σ)
    end
end

"""
    standardize(in, μ::AbstractFloat = NaN, σ::AbstractFloat = NaN)

Standardize a column.

Default output type is `Float64`.

The following standardizes `itrans` column `:B`.
```julia
standardize(:B => itrans)
```
"""
function standardize(in, μ::AbstractFloat = NaN, σ::AbstractFloat = NaN)
    Standardize(Connections([in], [first(in) => Out()]; maxins = 1, maxouts = 1), μ, σ)
end

isfit(t::Standardize) = !(isnan(t.μ) || isnan(t.σ))

Base.show(io::IO, t::Standardize) = print(io, instancestr(t) * "[" * repr(t.con) * "::{$(intype(t, inids(t)[1]))}|μ = $(t.μ) , σ = $(t.σ)::{$(outtype(t, outids(t)[1]))}]")

function func(t::Standardize, col, data)
    T = outtype(t, col)
    T = T == Nothing ? Float64 : T
    try
        return T.((data .- t.μ) ./ t.σ)
    catch e
        e isa InexactError && return T.(round.((data .- t.μ) ./ t.σ))
        rethrow(e)
    end
end

function Base.getindex(t::Standardize, col)
    checkfit(t)
    fsrc = intrn(t, col)[col] # col is same as icol
    (i)->func(t, col, fsrc(i))
end

function fit!(t::Standardize, data, smap)
    icol = first(inids(t))
    din = data[icol]
    if !isfit(t)
        t.μ = mean(din)
        t.σ = std(din)
        if t.σ <= 0
            t.σ = 1
        end
    end
    OrderedDict(col => func(t, col, data[icol]) for col in outids(t))
end

# ----- INVERSE -----
# -------------------
_inverse(o::DataSource, outs::Vector{IdOut}) = MergeCols(Connections(IdInp[], outs))

_inverse(o::OneHot, outs::Vector{IdOut}) = OneCold(Connections(IdInp[], outs; maxouts = 1), copy(o.tbl))

_inverse(o::OneCold, outs::Vector{IdOut}) = OneHot(Connections(IdInp[], outs; maxins = 1), copy(o.tbl))

function _inverse(o::ReplaceVals, outs::Vector{IdOut})
    length(values(o.tbl)) > length(unique(values(o.tbl))) &&
        throw(ArgumentError("cannot invert ReplaceVals $(idstr(o)), rules: $(o.tbl)"))
    revtbl = OrderedDict(v => k for (k, v) in o.tbl)
    ReplaceVals(Connections(IdInp[], outs; maxins = 1, maxouts = 1), revtbl)
end

_inverse(o::Standardize, outs::Vector{IdOut}) = Standardize(Connections(IdInp[], outs; maxins = 1, maxouts = 1), -o.μ / o.σ, 1 / o.σ)

function _invhelper(o::ColTrans, oi::OrderedDict{ColTrans,ColTrans})
    if o ∉ keys(oi)
        outs = IdOut[id => Out(intype(o, id)) for id in inids(o)]
        oi[o] = _inverse(o, outs)
    end
    for pid in inids(o)
        ptrn = intrn(o, pid)
        _invhelper(ptrn, oi)
        icon = oi[ptrn].con
        otype = ptrn.con.outs[pid].type
        pid ∈ keys(icon.ins) || push!(icon, pid, Inp(oi[o], otype))
    end
end

"""
    invert(odst::ColTrans, isrc::ColTrans)

Invert original transformation `odst`, the source for the inverted transform will be `isrc`.
"""
function invert(odst::ColTrans, isrc::ColTrans)
    if !(odst isa MergeCols)
        odst = mergecols(odst)
    end
    check(odst)
    checkfit(odst)
    osrc = getdatasource(odst)

    isrc.con.outs = OrderedDict(id => Out(intype(odst, id)) for id in outids(odst))
    oi = OrderedDict{ColTrans,ColTrans}(odst => isrc)
    _invhelper(odst, oi)
    idst = oi[osrc]
    idst.con.outs = OrderedDict(id => Out(intype(idst, id)) for id in inids(idst))
    idst
end

"""
    invert(odst::ColTrans; sid=:default)

Invert original transformation `odst`, feed it with a new `datasource` with given source id `sid`
based on `odst outputs`.
"""
function invert(odst::ColTrans; sid=:default)
    invert(odst, datasource(OrderedDict{Symbol,Symbol}(col => sid for col in outids(odst))))
end

"""
    invtransform(t::ColTrans, smap, type::Type)

Invert original transformation `odst` and immediatelly run `transform`.

Shortcut for `transform(invert(t), smap, type)`.
"""
invtransform(t::ColTrans, smap, type::Type) = transform(invert(t), smap, type)

# ----- COLLECTING/SORTING -----
# ------------------------------
function collecttrans(root::ColTrans)
    set = OrderedSet{ColTrans}()
    function helper(t::ColTrans)
        t ∈ set && return
        push!(set, t)
        for i in intrns(t)
            helper(i)
        end
    end
    helper(root)
    set
end

function toposort(root::ColTrans)
    nodes = collect(collecttrans(root))
    toposort(nodes, v->OrderedSet(intrns(v)))
end

# ----- SERIALIZATION -----
# -------------------------

function JSON.lower(t::Connections)
    ins = OrderedDict(id => OrderedDict("from" => string(v.trn.id), "type" => v.type) for (id, v) in t.ins)
    outs = OrderedDict(id => v.type for (id, v) in t.outs)
    OrderedDict("ins" => ins, "outs" => outs, "maxins" => t.maxins, "maxouts" => t.maxouts)
end
JSON.lower(t::DataSource) = OrderedDict("type" => "DataSource", "con" => t.con, "colsid" => t.sids)
JSON.lower(t::OneHot) = OrderedDict("type" => "OneHot", "con" => t.con, "tbl" => t.tbl, "tbltype" => typeof(t.tbl))
JSON.lower(t::OneCold) = OrderedDict("type" => "OneCold", "con" => t.con, "tbl" => t.tbl, "tbltype" => typeof(t.tbl))
JSON.lower(t::ReplaceVals) = OrderedDict("type" => "ReplaceVals", "con" => t.con, "tbl" => t.tbl, "tbltype" => typeof(t.tbl))
JSON.lower(t::Standardize) = OrderedDict("type" => "Standardize", "con" => t.con, "μ" => t.μ, "σ" => t.σ)

function JSON.lower(ddst::MergeCols)
    ts = toposort(ddst)
    tdict = OrderedDict()
    for t in ts
        if t == ddst
            tdict[t.id] = OrderedDict("type" => "MergeCols", "con" => t.con)
        else
            tdict[t.id] = JSON.lower(t)
        end
    end
    tdict
end

Base.parse(::Type{String}, s::String) = s
Base.parse(::Type{Symbol}, s::String) = Symbol(s)

function _converttable(tbl, tbltype)
    T = eval(Meta.parse(tbltype))
    T(parse(T.parameters[1], k) => v for (k, v) in tbl)
end

deserializetype(t::Type{DataSource}, p::OrderedDict, newnodes) = datasource(OrderedDict(Symbol(k) => v for (k, v) in p["colsid"]))
deserializetype(t::Type{OneHot}, p::OrderedDict, newnodes) = OneHot(deserializecons(p["con"], newnodes),
    _converttable(p["tbl"], p["tbltype"]))
deserializetype(t::Type{OneCold}, p::OrderedDict, newnodes) = OneCold(deserializecons(p["con"], newnodes),
    OrderedDict{Symbol,Any}(Symbol(k) => v for (k, v) in p["tbl"]))
deserializetype(t::Type{ReplaceVals}, p::OrderedDict, newnodes) = ReplaceVals(deserializecons(p["con"], newnodes),
    _converttable(p["tbl"], p["tbltype"]))
deserializetype(t::Type{Standardize}, p::OrderedDict, newnodes) = Standardize(deserializecons(p["con"], newnodes),
    p["μ"], p["σ"])
deserializetype(t::Type{MergeCols}, p::OrderedDict, newnodes) = MergeCols(deserializecons(p["con"], newnodes))

function deserializecons(d, newnodes)
    ins = OrderedDict{Symbol,Inp}(Symbol(id) => Inp(newnodes[info["from"]],
        eval(Meta.parse(info["type"]))) for (id, info) in d["ins"])
    outs = OrderedDict{Symbol,Out}(Symbol(id) => Out(eval(Meta.parse(type))) for (id, type) in d["outs"])
    #TODO fix maxins/maxouts in corresponding ColTrans
    Connections(ins, outs; maxins = d["maxins"], maxouts = d["maxouts"])
end

function deserializehelper(dict, newnodes)
    (dict isa OrderedDict && "type" ∈ keys(dict)) || return dict
    t = eval(Meta.parse(dict["type"]))
    deserializetype(t, dict, newnodes)
end

function deserialize(dict)
    function predf(v)
        "in" ∉ keys(last(v)) && return []
        string.(unique(values(last(v)["in"]))) # JSON keys as strings
    end

    orgnodes = toposort(dict, predf, first)
    newnodes = OrderedDict{String,ColTrans}()

    dst = nothing
    for n in orgnodes
        newnodes[n] = deserializehelper(dict[n], newnodes)
        if newnodes[n] isa MergeCols
            dst !== nothing && error("multiple MergeCols!")
            dst = newnodes[n]
        end
    end
    dst === nothing && error("no MergeCols in serialized data!")
    dst
end

# ----- BUILDERS -----
# --------------------
#TODO make these work for more than DataFrames!

"""
    nominal(df::AbstractDataFrame; maxunique::Union{Int,Nothing} = nothing)

Select nominal attributes from DataFrame.

Returns DataFrame column ids. These are ids corresponding to other than Real column types. Even Real types can be included
if less then or equal to `maxunique` values are present in data.
"""
function nominal(df::AbstractDataFrame; maxunique::Union{Int,Nothing} = nothing)
    ret = OrderedSet{Symbol}()
    for (name, type) in zip(names(df), eltype.(eachcol(df)))
        if type <: Real
            if maxunique != nothing && length(unique(df[name])) <= maxunique
                push!(ret, name)
            end
        else
            push!(ret, name)
        end
    end
    ret
end

"""
    tonumerical(df::AbstractDataFrame; nominal = nothing, standardize = true, sid = :default)

Automaticaly create a DataFrame transform to preprocess data as expected by many ML models such as neural networks.

All columns are transformed to numerical types. One may select columns (the `nominal` parameter) which should be
transformed using one hot encoding or `nominal()` is used to determine them automatically. All outputs are standardized
by default - this behavior can be altered using the `standardize` parameter.
"""
function tonumerical(df::AbstractDataFrame; nominal = nothing, standardize = true, sid = :default)
    if nominal === nothing
        nominal = Transforms.nominal(df)
    end
    nominal = OrderedSet(nominal)
    sid = Symbol(sid)
    smap = Dict(sid => df)

    src = datasource(df; sid=sid)
    transforms = Pair{Symbol,ColTrans}[]
    for id in outids(src)
        if id ∈ nominal
            oh = onehot(id => src)
            fit!(oh, smap)
            for outid in outids(oh)
                push!(transforms, outid => oh)
            end
        else
            push!(transforms, id => src)
        end
    end
    if standardize
        stransforms = Pair{Symbol,ColTrans}[]
        for (id, trn) in transforms
            st = Transforms.standardize(id => trn)
            fit!(st, smap)
            push!(stransforms, id => st)
        end
        transforms = stransforms
    end
    me = mergecols(transforms...)
end

"""
    copytransform(df::AbstractDataFrame, sid = :default)

Create identity transform for a DataFrame.
"""
function copytransform(df::AbstractDataFrame, sid = :default)
    sid = Symbol(sid)
    src = datasource(df; sid=sid)
    mergecols(src)
end

end

max_stars_repo_name: dhonza/MLGoodies


text:
function +(x::SMCg{N,T},y::SMCg{N,T}) where {N,T<:AbstractFloat}
	return SMCg{N,T}(x.cc+y.cc, x.cv+y.cv, x.cc_grad+y.cc_grad, x.cv_grad+y.cv_grad,
						 (x.Intv+y.Intv),(x.cnst && y.cnst),x.IntvBox,x.xref)
end
function -(x::SMCg{N,T},c::SMCg{N,T}) where {N,T<:AbstractFloat}
	return x + (-c)
end
function -(x::SMCg{N,T}) where {N,T<:AbstractFloat}
	return SMCg{N,T}(-x.cv, -x.cc, -x.cv_grad, -x.cc_grad, -x.Intv,
							 x.cnst, x.IntvBox, x.xref)
end
function /(x::SMCg{N,T},y::SMCg{N,T}) where {N,T<:AbstractFloat}
	pos_orth::Bool = (x.Intv.lo >= 0) && (y.Intv.lo >= 0)
	if (x==y)
		#println("ran me to 1")
		return one(T)
	elseif (MC_param.multivar_refine) && (~ (MC_param.mu >= 1)) && (pos_orth)
		#println("ran me to 2")
		Intv::Interval{T} = x.Intv/y.Intv
		cv1::T,pos1::Int64 = mid3(x.cv,x.cc,x.Intv.lo)
		cv1t::T = (cv1+sqrt(x.intv.lo*x.intv.hi))
				  /(sqrt(x.intv.lo)+sqrt(x.intv.hi))
		cv2t::T,pos2::Int64 = mid3(y.cv,y.cc,y.Intv.hi)
		cv::T = sqr(cv1t)/cv2t
		cv_grad::SVector{N,T} = 2.0*(cv1t/cv2t)/(sqrt(x.Intv.lo)+sqrt(x.Intv.hi))*
								  mid_grad(x.cc_grad, x.cv_grad, pos1) - ((cv1t/cv2t)^2)*
									mid_grad(y.cc_grad, y.cv_grad, pos2)
		cc1::T,pos1 = mid3(x.cv,x.cc,x.Intv.hi)
		cc2::T,pos2 = mid3(y.cv,y.cc,y.Intv.lo)
		gcc1::T = y.Intv.hi*cc1 - x.Intv.lo*cc2 + x.Intv.lo*y.Intv.lo
		gcc2::T = y.Intv.lo*cc1 - x.Intv.hi*cc2 + x.Intv.hi*y.Intv.hi
		if gcc1 <= gcc2
			cc::T = gcc1/(y.Intv.lo*y.Intv.hi)
			cc_grad::SVector{N,T} = one(T)/y.Intv.lo*mid_grad(x.cc_grad, x.cv_grad, pos1)
							- x.Intv.lo/(y.Intv.lo*y.Intv.hi)*mid_grad(y.cc_grad, y.cv_grad, pos2)
		else
			cc = gcc2/(y.Intv.lo*y.Intv.hi)
			cc_grad = one(T)/y.Intv.hi*mid_grad(x.cc_grad, x.cv_grad, pos1)
							- x.Intv.hi/(y.Intv.lo*y.Intv.hi)*mid_grad(y.cc_grad, y.cv_grad, pos2)
		end
		cnst::Bool = y.cnst ? x.cnst : (x.cnst ? y.cnst : x.cnst || y.cnst)
		return SMCg{T}(cc,cv,cc_grad,cv_grad,Intv,cnst,x.IntvBox,x.xref)
	else
		#println("ran me to 3")
		return x*inv(y)
	end
end

for numtype in union(int_list,float_list)
eval(quote
function +(x::SMCg{N,T},y::$numtype) where {N,T<:AbstractFloat}
	return SMCg{N,T}(x.cc+y, x.cv+y, x.cc_grad, x.cv_grad, (x.Intv+y),
	            x.cnst, x.IntvBox,x.xref)
end
function +(x::$numtype,y::SMCg{N,T}) where {N,T<:AbstractFloat}
	return SMCg{N,T}(x+y.cc, x+y.cv, y.cc_grad, y.cv_grad, (x+y.Intv),
	            y.cnst, y.IntvBox,y.xref)
end
function -(x::SMCg{N,T},c::$numtype) where {N,T<:AbstractFloat}
	return x + (-c)
end
function -(c::$numtype,x::SMCg{N,T}) where {N,T<:AbstractFloat}
	return c + (-x)
end
function *(x::SMCg{N,T},c::$numtype) where {N,T<:AbstractFloat}
	if (c>=zero(c))
		temp1::T = convert(T,c*x.cc)
		temp2::T = convert(T,c*x.cv)
		return SMCg{N,T}(temp1,temp2,c*x.cc_grad,c*x.cv_grad,c*x.Intv,x.cnst,x.IntvBox,x.xref)
	elseif (c<zero(c))
		temp1 = convert(T,c*x.cv)
		temp2 = convert(T,c*x.cc)
		return SMCg{N,T}(temp1,temp2,c*x.cv_grad,c*x.cc_grad,c*x.Intv,x.cnst,x.IntvBox,x.xref)
	end
end
function *(c::$numtype,x::SMCg{N,T}) where {N,T<:AbstractFloat}
	return x*c
end
@inline sqr(x::$numtype) = x*x
function /(x::SMCg{N,T},y::$numtype) where {N,T<:AbstractFloat}
	x*inv(y)
end
function /(x::$numtype,y::SMCg{N,T}) where {N,T<:AbstractFloat}
	x*inv(y)
end
end)
end
# convex relaxation of square (Khan2016)
function cv_sqr(xMC2::T,xL::T,xU::T) where {T<:AbstractFloat}
  if (zero(xMC2)<=xL||xU<=zero(xMC2))
		  return xMC2^2,2*xMC2
  elseif ((zero(xMC2)<=xMC2)&&((xL<zero(xL))&&(zero(xU)<xU)))
      return (xMC2^3)/xU,(3*xMC2^2)/xU
  elseif (zero(xMC2)>xMC2)&&((xL<zero(xL))&&(zero(xU)<xU))
      return (xMC2^3)/xL,(3*xMC2^2)/xL
  end
end
# concave relaxation of square (Khan2016)
function cc_sqr(xMC1::T,xL::T,xU::T) where {T<:AbstractFloat}
  return xL^2+(xL+xU)*(xMC1-xL),(xL+xU)
end

function sqr_cc_NS(x::T,lo::T,hi::T) where {T<:AbstractFloat}
  return line_seg(x,lo,lo^2,hi,hi^2),dline_seg(x,lo,lo^2,hi,hi^2,2*x)
end
function sqr_cv_NS(x::T,lo::T,hi::T) where {T<:AbstractFloat}
  return x^2,2*x
end
function sqr(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  	eps_min,blank = mid3(x.Intv.lo,x.Intv.hi,zero(x.Intv.lo))
  	eps_max = ifelse((abs(x.Intv.lo)>=abs(x.Intv.hi)),x.Intv.lo,x.Intv.hi)
	midcc,cc_id = mid3(x.cc,x.cv,eps_max)
	midcv,cv_id = mid3(x.cc,x.cv,eps_min)
	if (MC_param.mu >= 1)
  	cc,dcc = cc_sqr(midcc,x.Intv.lo,x.Intv.hi)
  	cv,dcv = cv_sqr(midcv,x.Intv.lo,x.Intv.hi)
		gcc1,gdcc1 = cc_sqr(x.cv,x.Intv.lo,x.Intv.hi)
		gcv1,gdcv1 = cv_sqr(x.cv,x.Intv.lo,x.Intv.hi)
		gcc2,gdcc2 = cc_sqr(x.cc,x.Intv.lo,x.Intv.hi)
		gcv2,gdcv2 = cv_sqr(x.cc,x.Intv.lo,x.Intv.hi)
		cv_grad = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
		cc_grad = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
	else
		cc,dcc = sqr_cc_NS(midcc,x.Intv.lo,x.Intv.hi)
		cv,dcv = sqr_cv_NS(midcv,x.Intv.lo,x.Intv.hi)
		cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
		cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
	end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, x.Intv^2, x.cnst, x.IntvBox,x.xref)
end

########### Defines functions required for linear algebra packages

one(x::SMCg{N,T}) where {N,T<:AbstractFloat} = SMCg{N,T}(one(Float64),one(Float64),zero.(x.cc_grad),zero.(x.cv_grad),Interval(one(Float64)),x.cnst,x.IntvBox,x.xref)

zero(x::SMCg{N,T}) where {N,T<:AbstractFloat} = SMCg{N,T}(zero(x.cc),zero(x.cv),zeros(SVector{N}),zeros(SVector{N}),Interval(zero(x.Intv.lo)),x.cnst,x.IntvBox,x.xref)

real(x::SMCg{N,T}) where {N,T<:AbstractFloat} = x

dist(x1::SMCg{N,T}, x2::SMCg{N,T}) where {N,T<:AbstractFloat} = max(abs(x1.cc-x2.cc), abs(x1.cv-x2.cv))

eps(x::SMCg{N,T}) where {N,T<:AbstractFloat} = max(eps(x.cc), eps(x.cv))

#=
###### Defines boolean operators

"""==(x::SMC,y::SMC) defines == for SMC type
"""
function ==(x::SMCg{N,T},y::SMCg{N,T})
	x.Intv == y.Intv && x.cv == y.cv && x.cc == y.cc
end
"""!=(x::SMC,y::SMC) defines != for SMC type
"""
function !=(x::SMCg{N,T},y::SMCg{N,T})
	x.Intv != y.Intv || x.cv != y.cv || x.cc != y.cc
end
"""<=(x::SMC,y::SMC) defines <= for SMC type
"""
function <=(x::SMCg{N,T},y::SMCg{N,T})
	x.Intv <= y.Intv && x.cv <= y.cv && x.cc <= y.cc
end
"""=>(x::SMC,y::SMC) defines => for SMC type
"""
function >=(x::SMCg{N,T},y::SMCg{N,T})
	x.Intv >= y.Intv && x.cv >= y.cv && x.cc >= y.cc
end
""">(x::SMC,y::SMC) defines > for SMC type
"""
function >(x::SMCg{N,T},y::SMCg{N,T})
	x.Intv > y.Intv && x.cv > y.cv && x.cc > y.cc
end
"""<(x::SMC,y::SMC) defines < for SMC type
"""
function <(x::SMCg{N,T},y::SMCg{N,T})
	x.Intv < y.Intv && x.cv < y.cv && x.cc < y.cc
end
=#


text:
# this is a script to train a basic VAE on the protobuffer pdf data
using GenerativeModels
using GMExtensions
using ConditionalDists
using PdfFraud
using JLD2
using FileIO
using Flux
using CuArrays
using Dates
using ValueHistories
using ProgressMeter
using ArgParse
using MLDataPattern
using DrWatson
quickactivate(dirname(dirname(pathof(PdfFraud)))) # this is so that DrWatson saves all the info

# this is for running the script from REPL - default values
if length(ARGS) == 0
	push!(ARGS, "2")
	push!(ARGS, "2")
	push!(ARGS, "4")
	push!(ARGS, "--test")
end

# parse arguments
s = ArgParseSettings()
@add_arg_table s begin
    "ldimsize"
    	required = true
    	arg_type = Int
    	help = "size of latent dimension"
   "channels"
    	required = true
    	arg_type = Int
    	help = "a list of channel numbers"
    	nargs = '+'
    "--scaling"
    	default = [2]
    	help = "a scalar or a vector of scaling factors"
    	arg_type = Int
    	nargs = '+'
    "--kernelsize"
    	default = [5]
    	help = "a scalar or a vector of kernelsizes"
    	arg_type = Int
    	nargs = '+'
    "--eta"
    	default = 0.001f0
    	arg_type = Float32
    	help = "learning rate"
    "--nepochs"
    	default = 10
    	arg_type = Int
    	help = "number of training epochs"
    "--batchsize"
    	default = 128
    	arg_type = Int
    	help = "batch size"
    "--beta"
    	default = 1.0f0
    	arg_type = Float32
    	help = "elbo β"
    "--savepath"
    	default = ""
    	arg_type = String
    	help = "where should the model be saved"
    "--save-frequency"
    	default = 1
    	arg_type = Int
    	help = "how often (in a number of epochs) should a model be saved"
    "--test"
    	action = :store_true
    	help = "test run"
end
parsed_args = parse_args(ARGS, s)

# extract settings from the parsed arguments
vectorize(x) = (length(x) == 1) ? repeat(x, nconv) : x
ldim = parsed_args["ldimsize"]
channels = parsed_args["channels"]
nconv = length(channels)
scalings = vectorize(parsed_args["scaling"])
kernelsizes = vectorize(parsed_args["kernelsize"])
η = parsed_args["eta"]
nepochs = parsed_args["nepochs"]
batchsize = parsed_args["batchsize"]
β = parsed_args["beta"]
savepath = parsed_args["savepath"]
save_freq = parsed_args["save-frequency"]
test = parsed_args["test"]
if test
	batchsize = 2
	savepath = "test"
end

# prepare the savepath
savepath = (length(savepath) > 0 && savepath[1] == "/") ? savepath :
	joinpath(PdfFraud.experiment_path, savepath, "models")
mkpath(savepath)

# get data
X, labels, sha256, page_nums = PdfFraud.load_pdf_data()
X = X |> gpu;

# construct the convolutional VAE model
(h, w, c) = size(X)[1:3]
densedims = []
vecdim = reduce(*, (h,w,c))
T = Float32

encoder = GMExtensions.conv_encoder((h, w, c), ldim*2, kernelsizes, channels, scalings;
	densedims = densedims)
enc_dist = CMeanVarGaussian{DiagVar}(encoder)

decoder = GMExtensions.conv_decoder((h, w, c), ldim, reverse(kernelsizes),
	reverse(channels), reverse(scalings); densedims = densedims, vec_output = true,
	vec_output_dim = vecdim + 1)
dec_dist = CMeanVarGaussian{ScalarVar}(decoder)

model = VAE(ldim, enc_dist, dec_dist) |> gpu

# create stuff for training
loss(x) = PdfFraud.elbo(model, x, β=β)
opt = ADAM(η)
if !test
	data = [(gpu(Array(x)),) for x in RandomBatches(X, batchsize, ceil(Int, size(X,4)/batchsize))];
	validation_data = gpu(Array(collect(RandomBatches(X, batchsize, 1))[1]));
else
	data = [(X[:,:,:,1:batchsize],) for _ in 1:100];
	validation_data = data[1][1];
end
@info "Training with $(length(data)) batches of $(batchsize) samples."
h = MVHistory()
tstart = "$(now())"
iepoch = 0

# train
function train_epoch(model, history, experiment_args)
	# setup callback
	prog = Progress(length(data))
	global iepoch += 1
	filename = joinpath(savepath, tstart*"_epoch-$(iepoch).bson")
	cb = PdfFraud.history_progress_cb(prog, model, validation_data, history, loss)

	# train
	Flux.train!(loss, params(model), data, opt, cb = cb)

	# save model
	if iepoch%save_freq == 0
		PdfFraud.save_model(filename, model, history, experiment_args)
	end
end
Flux.@epochs nepochs train_epoch(model, h, parsed_args)


text:
#=
Example use of XGBoost module

=#

#Modules
using XGBoost
using Distributions

#Includes
include "readlibsvm.jl"

# Get the data
#xTrain, yTrain = readlibsvm("../data/agaricus.txt.train", (6513, 126))
#xTest, yTest = readlibsvm("../data/agaricus.txt.test", (1611, 126))

# Number of rounds
#num_rounds = 2

#=WAYS OF TRAINING=#
# Dense matrix training of xgboost
#dense = xgboost(xTrain, num_round, label=yTrain, eta=1, max_depth=2, objective="binary:logistic")

# Sparse matrix training
#param = ["eta"=>1, "max_depth"=>2, "objective"=>"binary:logistic")
#sparse = xgboost(sptrain(xTrain), num_round, label=yTrain, param=param)

# Dmatrix
#dmat = DMatrix(xTrain, label=yTrain)
#model = xgboost(dmat, num_rounds, param=param)

# Alternatively, you can specify a libSVM formatted txt file as the data input

#=PREDICTION=#
#predict = predict(dense, xTest)

# Open the data file
filename = "/home/potterzot/data/uci_ml/agaricus-lepiota.data"
f = open(filename, "r")

k = 1
nfold = 5

if rand(0:nfold)==k
  # write test
else
  # write train

max_stars_repo_name: potterzot/jstats

max_issues_repo_name: potterzot/jstats

max_forks_repo_name: potterzot/jstats


text:
function inner_dual(thePortfolio::SparsePortfolioData, indices::Array{Int64})

    n=thePortfolio.f
    m=thePortfolio.m
    f=size(indices,1)

    innerQP = Model(solver=MosekSolver(MSK_DPAR_INTPNT_QO_TOL_PFEAS=1e-6, MSK_DPAR_INTPNT_QO_TOL_DFEAS=1e-6, MSK_IPAR_LOG=0, MSK_IPAR_MAX_NUM_WARNINGS=0)) #
    #

    @variable(innerQP, α[1:n])
    @variable(innerQP, λ)
    @variable(innerQP, w[1:f])
    @variable(innerQP, ρ[1:f]>=0.0)
    @variable(innerQP, βl[1:m]>=0.0)
    @variable(innerQP, βu[1:m]>=0.0)

    @constraint(innerQP, w.>=(thePortfolio.X[:, indices]'*α+thePortfolio.A[:, indices]'*(βl-βu)+λ*ones(f)+ρ-thePortfolio.d[indices]))
    @objective(innerQP, Max, (-0.5*α'*α-sum((thePortfolio.γ[indices]/2.0).*w.*w)+thePortfolio.Y'*α+λ+βl'*thePortfolio.l-βu'*thePortfolio.u+ρ'*thePortfolio.min_investment[indices])')

    status=solve(innerQP; suppress_warnings=true)
    return Dual(getvalue(α), getvalue(λ), getvalue(βl), getvalue(βu), getvalue(ρ), getvalue(w), getobjectivevalue(innerQP), status)
end

function inner_dual2(thePortfolio::SparsePortfolioData, indices::Array{Int64}, s::Array{Float64})

    n=thePortfolio.f
    m=thePortfolio.m
    f=size(indices,1)

    innerQP = Model(solver=MosekSolver(MSK_DPAR_INTPNT_QO_TOL_PFEAS=1e-6, MSK_DPAR_INTPNT_QO_TOL_DFEAS=1e-6, MSK_IPAR_LOG=0, MSK_IPAR_MAX_NUM_WARNINGS=0)) #
    #

    @variable(innerQP, α[1:n])
    @variable(innerQP, λ)
    @variable(innerQP, w[1:f])
    @variable(innerQP, ρ[1:f]>=0.0)
    @variable(innerQP, βl[1:m]>=0.0)
    @variable(innerQP, βu[1:m]>=0.0)

    @constraint(innerQP, w.>=(thePortfolio.X[:, indices]'*α+thePortfolio.A[:, indices]'*(βl-βu)+λ*ones(f)+ρ-thePortfolio.d[indices]))
    @objective(innerQP, Max, (-0.5*α'*α-sum((thePortfolio.γ[indices]/2.0).*s[indices].*w.*w)+thePortfolio.Y'*α+λ+βl'*thePortfolio.l-βu'*thePortfolio.u+ρ'*(s[indices].*thePortfolio.min_investment[indices]))')

    status=solve(innerQP; suppress_warnings=true)
    return Dual(getvalue(α), getvalue(λ), getvalue(βl), getvalue(βu), getvalue(ρ), getvalue(w), getobjectivevalue(innerQP), status)
end