text:
using Test
using DFTK

include("testcases.jl")

# TODO Once we have converged SCF densities in a file it would be better to instead / also
#      test the energies of these densities and compare them directly to the reference
#      energies obtained in the data files

@testset "Evaluate energies of guess density" begin
    Ecut = 15
    n_bands = 8
    fft_size = [27, 27, 27]

    Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.psp))
    model = model_DFT(silicon.lattice, [Si => silicon.positions], [:lda_x, :lda_c_vwn])
    basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)

    ρ0 = guess_density(basis, [Si => silicon.positions])
    E, H = energy_hamiltonian(basis, nothing, nothing; ρ=ρ0)

    @test E["Hartree"] ≈  0.3527293727197568  atol=5e-8
    @test E["Xc"]      ≈ -2.3033165870558165  atol=5e-8

    # Run one diagonalization and compute energies
    res = diagonalize_all_kblocks(lobpcg_hyper, H, n_bands, tol=1e-9)
    occupation = [[2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0, 0.0]
                  for i in 1:length(basis.kpoints)]
    ρnew = compute_density(H.basis, res.X, occupation)
    E, H = energy_hamiltonian(basis, res.X, occupation; ρ=ρnew)

    @test E["Kinetic"]        ≈  3.291847293270256   atol=5e-8
    @test E["AtomicLocal"]    ≈ -2.367978663117999   atol=5e-8
    @test E["AtomicNonlocal"] ≈  1.6527493682542034  atol=5e-8
    @test E["Hartree"]        ≈  0.6477025793366571  atol=5e-8
    @test E["Xc"]             ≈ -2.4375329416161162  atol=5e-8
    @test E["Ewald"]          ≈ -8.397893578467201   atol=5e-8
    @test E["PspCorrection"]  ≈ -0.294622067031369   atol=5e-8

    # Now we have a reasonable set of ψ, we make up a crazy model, and check the energies
    model = model_DFT(silicon.lattice,
                      [Si => silicon.positions],
                      [:gga_x_pbe, :gga_c_pbe],
                      extra_terms=[ExternalFromReal(X -> cos(1.2*(X[1]+X[3]))),
                                   ExternalFromFourier(X -> cos(1.3*(X[1]+X[3]))),
                                   PowerNonlinearity(1.2, 2.4),
                                   Magnetic(X -> [1, cos(1.4*X[2]), exp(X[3])])]
                      )
    basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)
    ρ = from_real(basis, ρnew.real)

    E, H = energy_hamiltonian(basis, res.X, occupation; ρ=ρ)

    @test E["Kinetic"]             ≈  3.291847293270256   atol=5e-8
    @test E["AtomicLocal"]         ≈ -2.367978663117999   atol=5e-8
    @test E["AtomicNonlocal"]      ≈  1.6527493682542034  atol=5e-8
    @test E["Hartree"]             ≈  0.6477025793366571  atol=5e-8
    @test E["Xc"]                  ≈ -2.456212919662419   atol=5e-8
    @test E["Ewald"]               ≈ -8.397893578467201   atol=5e-8
    @test E["PspCorrection"]       ≈ -0.294622067031369   atol=5e-8
    @test E["ExternalFromReal"]    ≈  0.139216686139006   atol=5e-8
    @test E["ExternalFromFourier"] ≈  0.057896835498415   atol=5e-8
    @test E["PowerNonlinearity"]   ≈  0.142649748399169   atol=5e-8
    @test E["Magnetic"]            ≈ -0.447218096875610   atol=5e-7
end


text:
@testset "MetNetUtils.jl" begin

    include("iders_tests.jl")
    include("getters_tests.jl")
    include("set_mets_tests.jl")
    include("set_rxn_tests.jl")
    include("invert_bkwds_tests.jl")
    include("del_rxn_tests.jl")
    include("well_scaled_model_tests.jl")

end

ext: jl


text:
##################################
# Operador de cruce  BLX-alpha
##############################
using Random

"""
    BLX(C1, C2, 𝛂)
Operador de cruce BLX, dados dos cromosomas
devuele dos descencientes.
Si C1 == C2 entoneces el cruce es la identidad
"""
function BLX(C1, C2, alpha=0.3)
    l = length(C1)
    H = Array{Float64, 2}(undef, 2, l)
    for j in 1:l
        c_max = maximum([C1[j], C2[j]])
        c_min = minimum([C1[j], C2[j]])
        I = c_max - c_min
        # auxiliales
        longitud_intervalo = I + 2*I* alpha
        for i in 1:2
            H[i,j] = min(1.0,max(0,rand()*longitud_intervalo + (c_min-I*alpha)))
        end
    end
    return [H[1,:], H[2,:]]
end


text:
# ------------------ #
# --- Structures --- #
# ------------------ #

"""
    iio_context_info

Empty struct to manage pointers to the `iio_context_info` C struct.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/structiio__context__info.html)
"""
mutable struct iio_context_info
end

"""
    iio_scan_block

Empty struct to manage pointers to the `iio_scan_block` C struct.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Scan.html#gad7fd2ea05bf5a8cebaff26b60edb8a13)
"""
mutable struct iio_scan_block
end

"""
    iio_scan_context

Empty struct to manage pointers to the `iio_scan_context` C struct.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Scan.html#gaa333dd2e410a2769cf5685019185d99c)
"""
mutable struct iio_scan_context
end

"""
    iio_context

Empty struct to manage pointers to the `iio_context` C struct.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/structiio__context.html)
"""
mutable struct iio_context
end

"""
    iio_device

Empty struct to manage pointers to the `iio_device` C struct.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/structiio__device.html)
"""
mutable struct iio_device
end

"""
    iio_channel

Empty struct to manage pointers to the `iio_channel` C struct.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/structiio__channel.html)
"""
mutable struct iio_channel
end

"""
    iio_buffer

Empty struct to manage pointers to the `iio_buffer` C struct.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/structiio__buffer.html)
"""
mutable struct iio_buffer
end

"""
    iio_data_format

Empty struct to manage pointers to the `iio_data_format` C struct.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/structiio__data__format.html)
"""
mutable struct iio_data_format
    length::Cuint;
    bits::Cuint;
    shift::Cuint;
    is_signed::Cuchar;
    is_fully_defined::Cuchar;
    is_be::Cuchar;
    with_scale::Cuchar;
    scale::Cdouble;
    repeat::Cuint;
end

# -------------- #
# --- Macros --- #
# -------------- #
"""
    @assert_null_pointer expression

    Throws an error if the expression returns a null pointer.
    If the expression return a `Ptr{type}` the error message is "type* null pointer".
"""
macro assert_null_pointer(ex)
    quote
        local val = $(esc(ex));
        typeassert(val, Ptr) != Ptr{Any}(0) || NO_ASSERT || throw(AssertionError(string(typeof(val))[5:end-1] * "* null pointer"));
        val
    end
end

# other C values
macro assert_Cstring(ex)
    :($(esc(ex)) != C_NULL || NO_ASSERT || error("Invalid Cstring: NULL pointer"));
end

# ------------- #
# --- Enums --- #
# ------------- #

"""
    iio_modifier

Copy of the equivalent C enum.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/iio_8h.html#a944ad22f426e09cdbb493081a05472e5)
"""
@enum iio_modifier begin
    IIO_NO_MOD
	IIO_MOD_X
	IIO_MOD_Y
	IIO_MOD_Z
	IIO_MOD_X_AND_Y
	IIO_MOD_X_AND_Z
	IIO_MOD_Y_AND_Z
	IIO_MOD_X_AND_Y_AND_Z
	IIO_MOD_X_OR_Y
	IIO_MOD_X_OR_Z
	IIO_MOD_Y_OR_Z
	IIO_MOD_X_OR_Y_OR_Z
	IIO_MOD_LIGHT_BOTH
	IIO_MOD_LIGHT_IR
	IIO_MOD_ROOT_SUM_SQUARED_X_Y
	IIO_MOD_SUM_SQUARED_X_Y_Z
	IIO_MOD_LIGHT_CLEAR
	IIO_MOD_LIGHT_RED
	IIO_MOD_LIGHT_GREEN
	IIO_MOD_LIGHT_BLUE
	IIO_MOD_QUATERNION
	IIO_MOD_TEMP_AMBIENT
	IIO_MOD_TEMP_OBJECT
	IIO_MOD_NORTH_MAGN
	IIO_MOD_NORTH_TRUE
	IIO_MOD_NORTH_MAGN_TILT_COMP
	IIO_MOD_NORTH_TRUE_TILT_COMP
	IIO_MOD_RUNNING
	IIO_MOD_JOGGING
	IIO_MOD_WALKING
	IIO_MOD_STILL
	IIO_MOD_ROOT_SUM_SQUARED_X_Y_Z
	IIO_MOD_I
	IIO_MOD_Q
	IIO_MOD_CO2
	IIO_MOD_VOC
	IIO_MOD_LIGHT_UV
	IIO_MOD_LIGHT_DUV
	IIO_MOD_PM1
	IIO_MOD_PM2P5
	IIO_MOD_PM4
	IIO_MOD_PM10
	IIO_MOD_ETHANOL
	IIO_MOD_H2
end

"""
    iio_chan_type

Copy of the equivalent C enum.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/iio_8h.html#a29714c3a5add6b599e29be0485ca548b)
"""
@enum iio_chan_type begin
    IIO_VOLTAGE
	IIO_CURRENT
	IIO_POWER
	IIO_ACCEL
	IIO_ANGL_VEL
	IIO_MAGN
	IIO_LIGHT
	IIO_INTENSITY
	IIO_PROXIMITY
	IIO_TEMP
	IIO_INCLI
	IIO_ROT
	IIO_ANGL
	IIO_TIMESTAMP
	IIO_CAPACITANCE
	IIO_ALTVOLTAGE
	IIO_CCT
	IIO_PRESSURE
	IIO_HUMIDITYRELATIVE
	IIO_ACTIVITY
	IIO_STEPS
	IIO_ENERGY
	IIO_DISTANCE
	IIO_VELOCITY
	IIO_CONCENTRATION
	IIO_RESISTANCE
	IIO_PH
	IIO_UVINDEX
	IIO_ELECTRICALCONDUCTIVITY
	IIO_COUNT
	IIO_INDEX
	IIO_GRAVITY
	IIO_POSITIONRELATIVE
	IIO_PHASE
	IIO_MASSCONCENTRATION
	IIO_CHAN_TYPE_UNKNOWN = C_INT_MAX
end


text:
using Revise
import LPDM, D3Trees, LightDarkLPDM, execute1d
# execute1d.execute(solver=:lpdm)
steps = 1
# execute1d.execute(s0=3/2*π, solver_mode = :lpdm, action_mode = :standard, obs_mode = :discrete, steps=steps, vis=Int64[], output=1)
# Juno.@enter execute1d.execute(s0=3/2*π, solver_mode = :lpdm, action_mode = :standard, obs_mode = :continuous, steps=steps, vis=Int64[1], output=3)
Juno.@run execute1d.execute(s0=3/2*π, solver_mode = :lpdm, action_mode = :standard, obs_mode = :discrete, steps=steps, vis=Int64[1], output=3)

# execute1d.batch_execute(n=50, debug=0, reward_func=:quadratic)
# execute1d.batch_execute(n=50, debug=0, reward_func=:quadratic)
# execute1d.batch_execute(n=30, debug=0)

max_stars_repo_name: ebalaban/LightDarkLPDM

max_issues_repo_name: tomerarnon/LightDarkLPDM

max_forks_repo_name: tomerarnon/LightDarkLPDM


text:
using StructuresKit

shape_name = "W14X90"

shape_info = StructuresKit.CrossSection.AISC(shape_name)

#A Ix Iy J Cw xc yc xs ys
section_properties = [(shape_info.A, shape_info.Ix, shape_info.Iy, shape_info.J, shape_info.Cw, 0.0, shape_info.d/2, 0.0, shape_info.d/2)]

#E  ν
material_properties = [(29000,0.30)]

#kx ky kϕ hx hy
springs = [0.0, 0.0, 0.0, 0.0, 0.0]

#member information
#L dL SectionProperties MaterialProperties Springs
member_definitions = [(13.0*12, 1.0,1,1,1)]

num_nodes = Int(member_definitions[1][1]/member_definitions[1][2]) + 1

#P
loads = 10000.0 * ones(num_nodes)

#end boundary conditions
#type=1 u''=v''=ϕ''=0 (simply supported), type=2 u'=v'=ϕ'=0  (fixed), type=3 u''=v''=ϕ''=u'''=v'''=ϕ'''=0 (free end, e.g., a cantilever)
end_boundary_conditions = [1 1]

#supports
#location where u=v=ϕ=0
supports = [0.0, 13.0*12]

#imperfection
L = member_definitions[1][1]
dL = member_definitions[1][2]
z = 0: dL :L
Δo = L / 1000 * sin.(π*z/L)

uo = Δo  #in the weak axis
vo = zeros(Float64, num_nodes)
ϕo = zeros(Float64, num_nodes)

imperfections = [uo, vo, ϕo]

#Calculate weak axis flexural buckling load
E = material_properties[1][1]
Iy = shape_info.Iy
Pcry = π^2*E*Iy/L^2   #weak axis buckling load

#calculate yield load
Py = 50*shape_info.A

#define load range for a VR test
P = 0:0.99*Pcry/10:0.99*Pcry

#inelasticity
inelasticity_flag = 0

# mutable struct Properties

#     A::Array{Float64,1}
#     Ix::Array{Float64,1}
#     Iy::Array{Float64,1}
#     J::Array{Float64,1}
#     Cw::Array{Float64,1}
#     xc::Array{Float64,1}
#     yc::Array{Float64,1}
#     xs::Array{Float64,1}
#     ys::Array{Float64,1}
#     xo::Array{Float64,1}
#     yo::Array{Float64,1}
#     Io::Array{Float64,1}
#     E::Array{Float64,1}
#     ν::Array{Float64,1}
#     G::Array{Float64,1}
#     kx::Float64
#     ky::Float64
#     kϕ::Float64
#     hx::Float64
#     hy::Float64
#     uo::Array{Float64,1}
#     vo::Array{Float64,1}
#     ϕo::Array{Float64,1}
#     P::Array{Float64,1}
#     dz::Array{Float64,1}
#     z::Array{Float64,1}
#     dm::Array{Float64,1}
#     Azz::Array{Float64,2}
#     Azzzz::Array{Float64,2}
#     supports::Array{Float64,1}

# end

# dz, z, dm = Mesh.define_line_element(member_definitions)

#    NumberOfNodes=length(dz)+1

#    #define property vectors
#    A =  Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 1)
#    Ix = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 2)
#    Iy = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 3)
#    J = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 4)
#    Cw = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 5)
#    xc = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 6)
#    yc = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 7)
#    xs = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 8)
#    ys = Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 9)

#    xo = -(xc .- xs)
#    yo = yc .- ys

#    Io = Ix .+ Iy .+ A .* (xo.^2 + yo.^2)

#    E = Mesh.create_line_element_property_array(member_definitions, dm, dz, material_properties, 4, 1)
#    ν = Mesh.create_line_element_property_array(member_definitions, dm, dz, material_properties, 4, 2)
#    G = E./(2 .*(1 .+ ν))

#    kx = springs[1]
#    ky = springs[2]
#    kϕ = springs[3]

#    hx = springs[4]
#    hy = springs[5]

#    #define imperfections
#    uo = imperfections[1]
#    vo = imperfections[2]
#    ϕo = imperfections[3]

#    #define load
#    P =  loads

#    #calculate derivative operators
#    Azzzz,Azz = Column.calculate_derivative_operators(dz)

#    #apply left and right end boundary condition stencils to derivative operators
#    NthDerivative = 4
#    Azzzz = Column.apply_end_boundary_conditions(Azzzz, end_boundary_conditions, NthDerivative, dz)

#    NthDerivative = 2
#    Azz = Column.apply_end_boundary_conditions(Azz, end_boundary_conditions, NthDerivative, dz)

#    column_properties = Properties(A, Ix, Iy, J, Cw, xc, yc, xs, ys, xo, yo, Io, E, ν, G, kx, ky, kϕ, hx, hy, uo, vo, ϕo, P, dz, z, dm, Azz, Azzzz, supports)

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

# dz, z, dm = StructuresKit.Mesh.define_line_element(member_definitions)
# A =  StructuresKit.Mesh.create_line_element_property_array(member_definitions, dm, dz, section_properties, 3, 1)

def_coordinates = []
connectivity = []

for i = 1:length(P)

    loads = P[i] * ones(num_nodes)

    #define column information
    column_properties = Column.initialize(member_definitions, section_properties, material_properties, loads, springs, end_boundary_conditions, supports, imperfections)

    #solve for column deformation
    u, v, ϕ = Column.solve(column_properties, inelasticity_flag)

    #add initial imperfection to solution
    u_total = uo .+ u
    v_total = vo .+ v
    ϕ_total = ϕo .+ ϕ

    #create surface mesh of undeformed column

    #define discretization along member
    zcoords = 0:dL:L
    n=(4, 4, 4, 4, 4)

    xcoords, ycoords = StructuresKit.CrossSection.wshape_nodes(shape_info, n)

    #undeformed shape
    coordinates, connectivity = StructuresKit.Mesh.surface(xcoords, ycoords, zcoords)

    #deformed shape

    σ_total_normal = Column.normal_stresses(xcoords, ycoords, zcoords, u, v, column_properties)
    ϵ_total_normal = Column.normal_stains(column_properties, σ_total_normal)
    w = Column.warping_displacements(zcoords, ϵ_total_normal)

    scale_u = 1.0
    scale_v = 1.0
    scale_ϕ = 1.0
    scale_w = 1.0

    x_deformed, y_deformed, z_deformed = Column.deformed_shape(xcoords, ycoords, zcoords, u, v, ϕ, w, scale_u, scale_v, scale_ϕ, scale_w)

    filename = string(@__DIR__, "/", "W14x90", "_SS_", string(i-1), ".ply")

    def_coordinates = [x_deformed y_deformed z_deformed]

    # Visualize.create_ply_file(def_coordinates, connectivity, filename)

end

using Makie
scene = poly(def_coordinates, connectivity, color=:lightgray, shading=:true, show_axis=:true, overdraw=:false, strokecolor = (:black, 0.6), strokewidth = 1)


text:
"""
    DatasourceMetadata(dataSource::DataSource; context=nothing)

Returns a Dict{String, String} with "timestamp" and "maxIngestedEventTime" as
keys.

context should be a Dict if provided.
"""
mutable struct DatasourceMetadata <: Query
    queryType::String
    dataSource::DataSource
    context
    function DatasourceMetadata(dataSource; context=nothing)
        nothing_or_type(context, Dict)
        new("dataSourceMetadata", dataSource, context)
    end
end
DatasourceMetadata(;dataSource, context=nothing) = DatasourceMetadata(dataSource; context)

function execute(client::Client, query::DatasourceMetadata; pretty=false)
    res = JSON.parse(execute_native_query(client, query; pretty))
    d = Dict{String, String}()
    d["timestamp"] = res[1]["timestamp"]
    d["maxIngestedEventTime"] = res[1]["result"]["maxIngestedEventTime"]
    d
end


text:
using Luna
import Luna.PhysData: wlfreq
import FFTW
import Hankel
import NumericalIntegration: integrate, SimpsonEven

a = 15e-6
gas = :Ar
pres = 5.0

τ = 30e-15
λ0 = 800e-9

w0 = 0.64*a
energy = 1e-6
L = 0.15

R = a
N = 32

grid = Grid.RealGrid(L, λ0, (200e-9, 3000e-9), 0.6e-12)
q = Hankel.QDHT(R, N, dim=2)

energyfun, energyfun_ω = Fields.energyfuncs(grid, q)

densityfun = let dens0=PhysData.density(gas, pres)
    z -> dens0
end

ionpot = PhysData.ionisation_potential(gas)
ionrate = Ionisation.ionrate_fun!_PPTcached(gas, λ0)

responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
             #Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot))

linop = LinearOps.make_const_linop(grid, q, PhysData.ref_index_fun(gas, pres))

normfun = NonlinearRHS.const_norm_radial(grid, q, PhysData.ref_index_fun(gas, pres))

inputs = Fields.GaussGaussField(λ0=λ0, τfwhm=τ, energy=energy, w0=w0)

Eω, transform, FT = Luna.setup(grid, q, densityfun, normfun, responses, inputs)

# statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid))
#output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω), length(q.r)))
output = Output.MemoryOutput(0, grid.zmax, 201)
Luna.run(Eω, grid, linop, transform, FT, output)

ω = grid.ω
t = grid.t

zout = output.data["z"]
Eout = output.data["Eω"]

Erout = (q \ Eout)
Iωr = abs2.(Erout)
# Iω0 = Iωr[:, 1, :]
Er0 = dropdims(Hankel.onaxis(Eout, q), dims=2);
Iω0 = abs2.(Er0);
Iω0log = log10.(Maths.normbymax(Iω0));
Etout = FFTW.irfft(Erout, length(grid.t), 1)

Ilog = log10.(Maths.normbymax(abs2.(Eout)))

It = PhysData.c * PhysData.ε_0/2 * abs2.(Maths.hilbert(Etout));
Itlog = log10.(Maths.normbymax(It))

Ir = zeros(Float64, (length(q.r), length(zout)))

Et = Maths.hilbert(Etout)
energy = zeros(length(zout))
for ii = 1:size(Etout, 3)
    energy[ii] = energyfun(Etout[:, :, ii]);
    Ir[:, ii] = integrate(grid.ω, Iωr[:, :, ii], SimpsonEven());
end

ω0idx = argmin(abs.(grid.ω .- 2π*PhysData.c/λ0))

zr = π*w0^2/λ0
points = L/2 .+ [-15, 3, 21].*zr
idcs = [argmin(abs.(zout .- point)) for point in points]

Epoints = [Hankel.symmetric(Et[:, :, idxi], q) for idxi in idcs]
rsym = Hankel.Rsymmetric(q);

import PyPlot:pygui, plt
pygui(true)

plt.figure()
plt.pcolormesh(zout*1e2, q.r*1e3, Iωr[ω0idx, :, :])
plt.colorbar()
plt.ylim(0, 4)
plt.xlabel("z (cm)")
plt.ylabel("r (mm)")
plt.title("I(r, ω=ω0)")

plt.figure()
plt.pcolormesh(zout*1e2, q.r*1e3, It[length(grid.t)÷2, :, :])
plt.colorbar()
plt.ylim(0, 4)
plt.xlabel("z (cm)")
plt.ylabel("r (mm)")
plt.title("I(r, t=0)")

plt.figure()
plt.pcolormesh(zout*1e2, q.r*1e3, Ir)
# plt.pcolormesh(zout*1e2, q.r*1e3, log10.(Maths.normbymax(Ir)))
plt.colorbar()
plt.ylim(0, R*1e3)
# plt.clim(-6, 0)
plt.xlabel("z (cm)")
plt.ylabel("r (mm)")
plt.title("\$\\int I(r, \\omega) d\\omega\$")

plt.figure()
plt.pcolormesh(zout*1e2, grid.ω*1e-15/2π, Iω0log)
plt.colorbar()
plt.clim(0, -6)
plt.xlabel("z (cm)")
plt.ylabel("f (PHz)")
plt.title("I(r=0, ω)")

plt.figure()
plt.plot(zout.*1e2, energy.*1e6)
plt.xlabel("Distance [cm]")
plt.ylabel("Energy [μJ]")

jw = Plotting.cmap_white("jet", 512, 10)
fig = plt.figure()
fig.set_size_inches(12, 4)
for ii in 1:3
    plt.subplot(1, 3, ii)
    plt.pcolormesh(grid.t*1e15, rsym*1e3, abs2.(Epoints[ii]'), cmap=jw)
    plt.xlim(-42, 42)
    plt.ylim(-1.8, 1.8)
end


text:
function ITensors.OpSum(::Model"heisenberg", n1, n2)
  opsum = OpSum()
  opsum += 0.5, "S+", n1, "S-", n2
  opsum += 0.5, "S-", n1, "S+", n2
  opsum += "Sz", n1, "Sz", n2
  return opsum
end

# H = Σⱼ (½ S⁺ⱼS⁻ⱼ₊₁ + ½ S⁻ⱼS⁺ⱼ₊₁ + SᶻⱼSᶻⱼ₊₁)
function ITensors.MPO(::Model"heisenberg", s)
  N = length(s)
  os = OpSum()
  for j in 1:(N - 1)
    os .+= 0.5, "S+", j, "S-", j + 1
    os .+= 0.5, "S-", j, "S+", j + 1
    os .+= "Sz", j, "Sz", j + 1
  end
  return MPO(os, s)
end

"""
    reference(::Model"heisenberg", ::Observable"energy"; N)

Compute the analytic isotropic heisenberg chain ground energy per site for length `N`.
Assumes the heisenberg model is defined with spin
operators not pauli matrices (overall factor of 2 smaller). Taken from [1].

[1] Nickel, Bernie. "Scaling corrections to the ground state energy
of the spin-½ isotropic anti-ferromagnetic Heisenberg chain." Journal of
Physics Communications 1.5 (2017): 055021
"""
function reference(::Model"heisenberg", ::Observable"energy"; N=∞)
  isinf(N) && return (0.5 - 2 * log(2)) / 2
  E∞ = (0.5 - 2 * log(2)) * N
  Eᶠⁱⁿⁱᵗᵉ = π^2 / (6N)
  correction = 1 + 0.375 / log(N)^3
  return (E∞ - Eᶠⁱⁿⁱᵗᵉ * correction) / (2N)
end


text:
"""
## formnf!

Returns nodal freedom numbering array nf

### Function
```julia
formnf!(nodof, nn, nf)
```

### Arguments
```julia
* nodof::Int       : Number of degrees of freedom for each node
* nn::Int          : Number of nodes in mesh
* nf::Array{Int,2} : Nodal freedom matrix (updated)
```
"""
function formnf!(nodof::Int, nn::Int, nf::Matrix{Int})
  local m = 0
  for j in 1:nn
    for i in 1:nodof
      if nf[i,j] != 0
        m += 1
        nf[i,j] = m
      end
    end
  end
end


text:
# Particularly useful or beautiful grammar of graphics invocations.

## Return a DataFrame with x, y column suitable for plotting a function.
#
# Args:
#  f: Function/Expression to be evaluated.
#  a: Lower bound.
#  b: Upper bound.
#  n: Number of points to evaluate the function at.
#
# Returns:
#  A data frame with "x" and "f(x)" columns.
#
function evalfunc(f::Function, a, b, n)
    xs = [x for x in a:(b - a)/n:b]
    df = DataFrame(xs, map(f, xs))
    colnames!(df, ["x", "f(x)"])
    df
end

evalfunc(f::Expr, a, b, n) = evalfunc(eval(:(x -> $f)), a, b, n)

# A convenience plot function for quickly plotting functions are expressions.
#
# Args:
#   fs: An array in which each object is either a single argument function or an
#       expression computing some value on x.
#   a: Lower bound on x.
#   b: Upper bound on x.
#   elements: One ore more grammar elements.
#
# Returns:
#   A plot objects.
#
function plot(fs::Array, a, b, elements::Element...)
    df = DataFrame()
    for (i, f) in enumerate(fs)
        df_i = evalfunc(f, a, b, 250)
        name = typeof(f) == Expr ? string(f) : @sprintf("f<sub>%d</sub>", i)
        df_i = cbind(df_i, [name for _ in 1:size(df_i)[1]])
        colnames!(df_i, ["x", "f(x)", "f"])
        df = rbind(df, df_i)
    end

    mapping = {:x => "x", :y => "f(x)"}
    if length(fs) > 1
        mapping[:color] = "f"
    end

    plot(df, mapping, Geom.line, elements...)
end

# Plot a single function.
function plot(f::Function, a, b, elements::Element...)
    plot([f], a, b, elements...)
end

# Plot a single expression.
function plot(f::Expr, a, b, elements::Element...)
    plot([f], a, b, elements...)
end

# Plot an expression from a to b.
macro plot(expr, a, b)
    quote
        plot(x -> $(expr), $(a), $(b))
    end
end

# Simple heatmap plots of matrices.
#
# Args:
#   M: A matrix.
#
# Returns:
#   A plot object.
#
function spy(M)
    is, js, values = findn_nzs(M)
    df = DataFrame({"i" => is, "j" => js, "value" => values})
    plot(df, x="j", y="i", color="value",
         Scale.color_gradient,
         Geom.rectbin, Stat.identity)
end


text:
using Bio.Seq
using Bio.Align
using BioBenchmarks

seq1 = first(open(datafile("Q8WZ42.fasta"), FASTA)).seq
seq2 = first(open(datafile("A2ASS6.fasta"), FASTA)).seq
affinegap = AffineGapScoreModel(BLOSUM62, -10, -1)
@benchmark pairalign(GlobalAlignment(), seq1, seq2, affinegap)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -37.92418043260588 <= q <= 19.457860872935036)

                     add_NL_constraint(m, :(swish(-0.5948345021955626 + 0.7909419235855366*swish(0.16572281076834727 + 0.7701313689304001*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.3960198610694281*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.006655112331988011*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.5604852483433853*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.5611932928423773*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.9350053029495258*swish(0.6292157071886879 + 0.7447123643459479*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.8071073342022652*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.7192221993944883*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1911667575027427*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.30728485701247665*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.6696024550400765*swish(-0.5341808090129341 + 0.4032770914369115*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.6061337759397625*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.5534397384266061*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.16208042832895408*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.7376214833764658*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.9733717899660395*swish(0.43846326345255404 + -0.02862602361460853*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.5029397949622734*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.20686891235628035*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.9777195341190259*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.2597965376833278*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.07261217388131413*swish(0.02496868545483233 + -0.04485855107115411*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.47774406442351136*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.5903487696420049*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1719483424118926*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.8023455576856944*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2])))) + swish(-0.7044745597894595 + -0.85619817917634*swish(0.16572281076834727 + 0.7701313689304001*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.3960198610694281*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.006655112331988011*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.5604852483433853*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.5611932928423773*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.37830231556965854*swish(0.6292157071886879 + 0.7447123643459479*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.8071073342022652*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.7192221993944883*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1911667575027427*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.30728485701247665*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.12156626196118481*swish(-0.5341808090129341 + 0.4032770914369115*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.6061337759397625*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.5534397384266061*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.16208042832895408*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.7376214833764658*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.3744591738356289*swish(0.43846326345255404 + -0.02862602361460853*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.5029397949622734*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.20686891235628035*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.9777195341190259*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.2597965376833278*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.4985780257939494*swish(0.02496868545483233 + -0.04485855107115411*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.47774406442351136*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.5903487696420049*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1719483424118926*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.8023455576856944*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2])))) + swish(-0.9683275825894269 + -0.7637278723892034*swish(0.16572281076834727 + 0.7701313689304001*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.3960198610694281*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.006655112331988011*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.5604852483433853*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.5611932928423773*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.9143661573707975*swish(0.6292157071886879 + 0.7447123643459479*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.8071073342022652*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.7192221993944883*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1911667575027427*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.30728485701247665*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.7824888723979511*swish(-0.5341808090129341 + 0.4032770914369115*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.6061337759397625*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.5534397384266061*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.16208042832895408*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.7376214833764658*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.6866950641423624*swish(0.43846326345255404 + -0.02862602361460853*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.5029397949622734*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.20686891235628035*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.9777195341190259*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.2597965376833278*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.09633196432882363*swish(0.02496868545483233 + -0.04485855107115411*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.47774406442351136*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.5903487696420049*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1719483424118926*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.8023455576856944*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2])))) + swish(-0.5143076265809325 + -0.11312975790997548*swish(0.16572281076834727 + 0.7701313689304001*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.3960198610694281*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.006655112331988011*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.5604852483433853*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.5611932928423773*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.7215302240627146*swish(0.6292157071886879 + 0.7447123643459479*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.8071073342022652*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.7192221993944883*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1911667575027427*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.30728485701247665*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.8436266175918532*swish(-0.5341808090129341 + 0.4032770914369115*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.6061337759397625*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.5534397384266061*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.16208042832895408*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.7376214833764658*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.8029468027318871*swish(0.43846326345255404 + -0.02862602361460853*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.5029397949622734*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.20686891235628035*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.9777195341190259*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.2597965376833278*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.07363029133667531*swish(0.02496868545483233 + -0.04485855107115411*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.47774406442351136*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.5903487696420049*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1719483424118926*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.8023455576856944*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2])))) + swish(-0.8626839052591899 + -0.6726505975051698*swish(0.16572281076834727 + 0.7701313689304001*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.3960198610694281*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.006655112331988011*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.5604852483433853*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.5611932928423773*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.7431517048783043*swish(0.6292157071886879 + 0.7447123643459479*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.8071073342022652*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.7192221993944883*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1911667575027427*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.30728485701247665*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.03587478002549904*swish(-0.5341808090129341 + 0.4032770914369115*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.6061337759397625*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.5534397384266061*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.16208042832895408*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.7376214833764658*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.7126331133171049*swish(0.43846326345255404 + -0.02862602361460853*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.5029397949622734*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.20686891235628035*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.9777195341190259*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.2597965376833278*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.7945114361960401*swish(0.02496868545483233 + -0.04485855107115411*swish(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.47774406442351136*swish(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.5903487696420049*swish(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1719483424118926*swish(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.8023455576856944*swish(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module Catlab

include("core/Core.jl")
include("theories/Theories.jl")

include("categorical_algebra/CSetDataStructures.jl")
include("categorical_algebra/Permutations.jl")
include("graphs/Graphs.jl")
include("categorical_algebra/CategoricalAlgebra.jl")

include("wiring_diagrams/WiringDiagrams.jl")
include("graphics/Graphics.jl")
include("programs/Programs.jl")
include("linear_algebra/LinearAlgebra.jl")

end


text:
### Computing fields and other physical quantities ###

"""
    temperature(sol::ODESolution,sys::HeatConduction,t)

Return the vorticity field associated with solution vector `sol` on the grid in `sys`,
at time(s) `t`.
"""
@inline temperature(T::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY},t::Real) where {NX,NY} = T

function temperature_linesource(T::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY},t::Real;index=1) where {NX,NY}
    sys.qline[index].E*T
end

"""
    temperature(integrator)

Return the temperature field associated with `integrator` at its current state.
""" function temperature end

"""
    temperature_linesource(integrator,[index=1])

Using the solution associated with `integrator` at its current state, evaluate
the temperature field along line source `index` (which defaults to 1).
""" function temperature_linesource end

for fcn in (:temperature,:temperature_linesource)

  @eval $fcn(s::ConstrainedSystems.ArrayPartition,sys::HeatConduction,t;kwargs...) = $fcn(state(s),sys,t;kwargs...)

end

for fcn in (:temperature,:temperature_linesource)
  @eval $fcn(integ::ConstrainedSystems.OrdinaryDiffEq.ODEIntegrator;kwargs...) = $fcn(integ.u,integ.p,integ.t;kwargs...)

  @eval $fcn(sol::ConstrainedSystems.OrdinaryDiffEq.ODESolution,sys::HeatConduction,t;kwargs...) = $fcn(sol(t),sys,t;kwargs...)

  @eval $fcn(sol::ConstrainedSystems.OrdinaryDiffEq.ODESolution,sys::HeatConduction,t::AbstractArray;kwargs...) =
      map(ti -> $fcn(sol(ti),sys,ti;kwargs...),t)

end

# Surface field quantities

#=
"""
    traction(integrator)

Return the traction at all of the surface points associated with `integrator` at
its current state.
""" function traction end

"""
    pressurejump(integrator)

Return the pressure jump at all of the surface points associated with `integrator` at
its current state.
""" function pressurejump end

for fcn in (:traction,:pressurejump)
  @eval $fcn(s::ConstrainedSystems.ArrayPartition,a...;kwargs...) = $fcn(constraint(s),a...;kwargs...)
  @eval $fcn(integ::ConstrainedSystems.OrdinaryDiffEq.ODEIntegrator) = $fcn(integ.u,integ.p,integ.t)
  @eval $fcn(sol::ODESolution,sys::HeatConduction,t) = $fcn(sol(t),sys,t)

  @eval $fcn(sol::ODESolution,sys::HeatConduction,t::AbstractArray) =
        map(ti -> $fcn(sol(ti),sys,ti),t)

end

"""
    traction(sol::ODESolution,sys::HeatConduction,t)

Return the traction at all of the surface points associated with solution `sol` on
grid in `sys` at time(s) t.
"""
traction(τ,sys::HeatConduction{NX,NY,0,MT,FS,SD},t) where {NX,NY,MT,FS,SD} = τ

traction(τ::VectorData{N},sys::HeatConduction{NX,NY,N,MT,FS,ExternalInternalFlow},t) where {NX,NY,N,MT,FS} = τ

function traction(τ::VectorData{N},sys::HeatConduction{NX,NY,N,MT,FS,SD},t) where {NX,NY,N,MT,FS,SD}
    nrm = normals(sys.bodies)
    relative_surface_velocity!(sys.Δus,sys,t)
    pointwise_dot!(sys.Sb,nrm,sys.Δus)
    surface_velocity_jump!(sys.Δus,sys,t)
    product!(sys.Vb,sys.Δus,sys.Sb)
    return τ + sys.Vb
end

"""
    pressurejump(sol::ODESolution,sys::HeatConduction,t)

Return the pressure jump at all of the surface points associated with solution `sol` on
grid in `sys` at time(s) t.
"""
function pressurejump(τ::VectorData{N},sys::HeatConduction{NX,NY,N},t) where {NX,NY,N}

    #_hasfilter(sys) ? (τf = similar(τ); τf .= sys.Cf*force(τ,sys)) : τf = deepcopy(force(τ,sys))
    #τf = deepcopy(τ)
    #return nrm.u∘τf.u + nrm.v∘τf.v # This might need some modification

    press = zero(sys.Sb)
    nrm = normals(sys.bodies)
    pointwise_dot!(press,nrm,traction(τ,sys,t))
    return -press
end

## Total quantities

force(τ,sys::HeatConduction{NX,NY,0},t,bodyi::Int) where {NX,NY} = Vector{Float64}(), Vector{Float64}()

function force(τ::VectorData{N},sys::HeatConduction,t,bodyi::Int) where {N}
    product!(sys.τ,traction(τ,sys,t),areas(sys.bodies))
    fxvec = sum(sys.τ.u,sys.bodies,bodyi)
    fyvec = sum(sys.τ.v,sys.bodies,bodyi)
    return fxvec, fyvec
end

force(s::ConstrainedSystems.ArrayPartition,sys,t,bodyi) = force(constraint(s),sys,t,bodyi)

"""
    force(integ,bodyindex)

Given the state of the system in integrator `integ`, return the current force
on the body with index `bodyindex`.
"""
force(integ::ConstrainedSystems.OrdinaryDiffEq.ODEIntegrator,bodyi) = force(integ.u,integ.p,integ.t,bodyi)

"""
    force(sol,sys::HeatConduction,bodyindex) -> Tuple{Vector,Vector}

Given the solution history vector `sol` and the system `sys`, return the force
history on the body with index `bodyindex` as a tuple of vectors, one for
each component.
"""
function force(sol::ConstrainedSystems.OrdinaryDiffEq.ODESolution,sys::HeatConduction,bodyi::Int)
    fx = map((u,t) -> force(u,sys,t,bodyi)[1],sol.u,sol.t)
    fy = map((u,t) -> force(u,sys,t,bodyi)[2],sol.u,sol.t)
    fx, fy
end

function moment(τ::VectorData{N},bodies::BodyList,sys::HeatConduction,t,bodyi::Int;center=(0.0,0.0)) where {N}
    xc, yc = center
    x, y = collect(bodies)
    dX = VectorData(x .- xc,  y .- yc)
    product!(sys.τ,traction(τ,sys,t),areas(bodies))
    return sum(cross(dX,sys.τ),bodies,bodyi)
end

function moment(s::ConstrainedSystems.ArrayPartition,sys::HeatConduction{NX,NY,N,MovingPoints},t,bodyi;center=(0.0,0.0)) where {NX,NY,N}
    bodies = deepcopy(sys.bodies)
    tl! = RigidTransformList(aux_state(s))
    tl!(bodies)
    moment(constraint(s),bodies,sys,t,bodyi;center=center)
end

moment(s::ConstrainedSystems.ArrayPartition,sys::HeatConduction{NX,NY,N,StaticPoints},t,bodyi;center=(0.0,0.0)) where {NX,NY,N} =
        moment(constraint(s),sys.bodies,sys,t,bodyi;center=center)

"""
    moment(integ,bodyindex[,center=(0,0))

Given the state of the system in integrator `integ`, return the current moment
on the body with index `bodyindex`.  The center of the
moment can be passed as an optional tuple argument.
"""
moment(integ::ConstrainedSystems.OrdinaryDiffEq.ODEIntegrator,bodyi;center=(0.0,0.0)) =
      moment(integ.u,integ.p,integ.t,bodyi,center=center)

"""
    moment(sol,sys::HeatConduction,bodyindex[,center=(0,0)]) -> Vector

Given the solution history vector `sol` and the system `sys`, return the moment
history on the body with index `bodyindex` as a vector. The center of the
moment can be passed as an optional tuple argument.
"""
function moment(sol::ConstrainedSystems.OrdinaryDiffEq.ODESolution,sys::HeatConduction,bodyi::Int;center=(0.0,0.0))
    mom = map((u,t) -> moment(u,sys,t,bodyi,center=center),sol.u,sol.t)
end
=#


text:
export iauUt1utc
"""
Time scale transformation:  Universal Time, UT1, to Coordinated
Universal Time, UTC.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  canonical.

Given:
   ut11,ut12  double   UT1 as a 2-part Julian Date (Note 1)
   dut1       double   Delta UT1: UT1-UTC in seconds (Note 2)

Returned:
   utc1,utc2  double   UTC as a 2-part quasi Julian Date (Notes 3,4)

Returned (function value):
               int      status: +1 = dubious year (Note 5)
                                 0 = OK
                              -1 = unacceptable date

Notes:

   1. ut11+ut12 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where ut11 is the Julian
   Day Number and ut12 is the fraction of a day.  The returned utc1
   and utc2 form an analogous pair, except that a special convention
   is used, to deal with the problem of leap seconds - see Note 3.

   2. Delta UT1 can be obtained from tabulations provided by the
   International Earth Rotation and Reference Systems Service.  The
   value changes abruptly by 1s at a leap second;  however, close to
   a leap second the algorithm used here is tolerant of the "wrong"
   choice of value being made.

   3. JD cannot unambiguously represent UTC during a leap second unless
   special measures are taken.  The convention in the present
   function is that the returned quasi JD day UTC1+UTC2 represents
   UTC days whether the length is 86399, 86400 or 86401 SI seconds.

   4. The function iauD2dtf can be used to transform the UTC quasi-JD
   into calendar date and clock time, including UTC leap second
   handling.

   5. The warning status "dubious year" flags UTCs that predate the
   introduction of the time scale or that are too far in the future
   to be trusted.  See iauDat for further details.

Called:
   iauJd2cal    JD to Gregorian calendar
   iauDat       delta(AT) = TAI-UTC
   iauCal2jd    Gregorian calendar to JD

References:

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

   Explanatory Supplement to the Astronomical Almanac,
   P. Kenneth Seidelmann (ed), University Science Books (1992)

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauUt1utc(double ut11, double ut12, double dut1,
#               double *utc1, double *utc2)

function iauUt1utc(ut11::Real, ut12::Real, dta::Real)
   ref_tt1 = Ref{Float64}(0.0)
   ref_tt2 = Ref{Float64}(0.0)

   status = ccall((:iauUt1utc, libsofa_c), Cint,
                   (Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                   convert(Float64, ut11), convert(Float64, ut12), convert(Float64, dta),
                   ref_tt1, ref_tt2)

   return status, ref_tt1[], ref_tt2[]
end


text:
function jacobian(m::HetDSGELag)
    truncate_distribution!(m)

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

    # Load in parameters, steady-state parameters, and grids
    r::Float64     = m[:r].value
    α::Float64     = m[:α].value
    H::Float64     = m[:H].value
    δ::Float64     = m[:δ].value
    μ_sp::Float64  = m[:μ_sp].value
    ρ_sp::Float64  = m[:ρ_sp].value
    σ_sp::Float64  = m[:σ_sp].value
    γ::Float64     = m[:γ].value
    g::Float64     = m[:g].value
    η::Float64     = m[:η].value
    ρ_B::Float64    = m[:ρ_B].value
    ρ_G::Float64    = m[:ρ_G].value
    ρ_z::Float64    = m[:ρ_z].value
    ρ_μ::Float64    = m[:ρ_μ].value
    ρ_lamw::Float64 = m[:ρ_lamw].value
    ρ_lamf::Float64 = m[:ρ_lamf].value
    ρ_mon::Float64  = m[:ρ_mon].value
    spp::Float64   = m[:spp].value
    lamw::Float64  = m[:lamw].value
    ϕh::Float64    = m[:ϕh].value
    Φw::Float64    = m[:Φw].value
    lamf::Float64  = m[:lamf].value
    Φp::Float64    = m[:Φp].value
    ρ_R::Float64    = m[:ρR].value
    ψπ::Float64    = m[:ψπ].value
    ψy::Float64    = m[:ψy].value

    R = 1 + r

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

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

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

    zlo::Float64 = get_setting(m, :zlo)
    zhi::Float64 = get_setting(m, :zhi)

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

    qp(z) = dmollifier_hetdsge_lag(z, zhi, zlo)
    qfunction_hetdsge_lag(x) = mollifier_hetdsge_lag(x, zhi, zlo) #/sumz

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

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

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

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

    setup_indices!(m)
    normalize_model_state_indices!(m)

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

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

    # KF eqn
    JJ[eq[:eq_kolmogorov_fwd],endo[:μ_t1]]   = dF2_dM
    JJ[eq[:eq_kolmogorov_fwd],endo[:R_t1]]  = dF2_dRZ
    JJ[eq[:eq_kolmogorov_fwd],endo[:l_t1]] = dF2_dELL
    JJ[eq[:eq_kolmogorov_fwd],endo[:z_t]]    = -dF2_dRZ
    JJ[eq[:eq_kolmogorov_fwd],endo[:μ_t]]    = -eye(nxns)
    JJ[eq[:eq_kolmogorov_fwd],endo[:w_t]]    = dF2_dWH
    JJ[eq[:eq_kolmogorov_fwd],endo[:L_t]]   = dF2_dWH
    JJ[eq[:eq_kolmogorov_fwd],endo[:t_t]]   = dF2_dTT

    # update lagged ELL
    JJ[eq[:eq_lag_ell],endo[:l′_t1]] = Matrix{Float64}(I, nxns, nxns)
    JJ[eq[:eq_lag_ell],endo[:l_t]]   = -Matrix{Float64}(I, nxns, nxns)

    # update lagged M
    JJ[eq[:eq_lag_wealth],endo[:μ′_t1]] = Matrix{Float64}(I, nxns, nxns)
    JJ[eq[:eq_lag_wealth],endo[:μ_t]]   = -Matrix{Float64}(I, nxns, nxns)

    # mkt clearing
    JJ[first(eq[:eq_market_clearing]),first(endo[:y_t])]   = ystar/g
    JJ[first(eq[:eq_market_clearing]),first(endo[:G])]   = -ystar/g
    JJ[first(eq[:eq_market_clearing]),first(endo[:I_t])]   = -xstar
    JJ[eq[:eq_market_clearing],endo[:l_t]] = (μ.*unc.*xswts.*c)'
    JJ[eq[:eq_market_clearing],endo[:μ_t]]   = -(xswts.*c)' # note, now we linearize

    # lambda = average marginal utility
    JJ[first(eq[:eq_lambda]),first(endo[:mu_t])] = lam
    JJ[first(eq[:eq_lambda]),endo[:μ_t]]   = -(xswts./c)' # note, now we linearize
    JJ[first(eq[:eq_lambda]),endo[:l_t]] = -(xswts.*unc.*μ./c)'

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

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

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

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

    # wage phillips curve
    JJ[first(eq[:eq_wage_phillips]),first(endo[:wageinflation_t])]  = -1.
    JJ[first(eq[:eq_wage_phillips]),first(endo[:LAMW])] = (ϕ*H^ϕh)/Φw
    JJ[first(eq[:eq_wage_phillips]),first(endo[:L_t])]   = (ϕ*(H^ϕh)*(1+lamw)/lamw*Φw)*ϕh
    JJ[first(eq[:eq_wage_phillips]),first(endo[:mu_t])]  = -(ϕ*(H^ϕh)*(1+lamw)/lamw*Φw)
    JJ[first(eq[:eq_wage_phillips]),first(endo[:w_t])]    = -(ϕ*(H^ϕh)*(1+lamw)/lamw*Φw)
    JJ[first(eq[:eq_wage_phillips]),first(endo[:wageinflation′_t])]  = β

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

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

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

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

    # taylor rule
    JJ[first(eq[:eq_taylor]),first(endo[:i_t])]   = -1.
    JJ[first(eq[:eq_taylor]),first(endo[:i_t1])]  = ρ_R
    JJ[first(eq[:eq_taylor]),first(endo[:π_t])]   = (1-ρ_R)*ψπ
    JJ[first(eq[:eq_taylor]),first(endo[:π_t1])]  = (1-ρ_R)*ψπ
    JJ[first(eq[:eq_taylor]),first(endo[:π_t2])] = (1-ρ_R)*ψπ
    JJ[first(eq[:eq_taylor]),first(endo[:π_t3])] = (1-ρ_R)*ψπ
    JJ[first(eq[:eq_taylor]),first(endo[:y_t])]    = (1-ρ_R)*ψy
    JJ[first(eq[:eq_taylor]),first(endo[:y_t4])]  = -(1-ρ_R)*ψy
    JJ[first(eq[:eq_taylor]),first(endo[:z_t])]    = (1-ρ_R)*ψy
    JJ[first(eq[:eq_taylor]),first(endo[:z_t1])]   = (1-ρ_R)*ψy
    JJ[first(eq[:eq_taylor]),first(endo[:z_t2])]  = (1-ρ_R)*ψy
    JJ[first(eq[:eq_taylor]),first(endo[:z_t3])]  = (1-ρ_R)*ψy
    JJ[first(eq[:eq_taylor]),first(endo[:MON])]  = 1.

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

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

    # update lagged variables
    JJ[first(eq[:LR]),first(endo[:R′_t1])] = 1.
    JJ[first(eq[:LR]),first(endo[:R_t])]   = -1.

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

    JJ[first(eq[:LPI]),first(endo[:π′_t1])] = 1.
    JJ[first(eq[:LPI]),first(endo[:π_t])]   = -1.

    JJ[first(eq[:L2PI]),first(endo[:π′_t2])] = 1.
    JJ[first(eq[:L2PI]),first(endo[:π_t1])]   = -1.

    JJ[first(eq[:L3PI]),first(endo[:π′_t3])] = 1.
    JJ[first(eq[:L3PI]),first(endo[:π_t2])]  = -1.

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

    JJ[first(eq[:L2Y]),first(endo[:y′_t2])] = 1.
    JJ[first(eq[:L2Y]),first(endo[:y_t1])]   = -1.

    JJ[first(eq[:L3Y]),first(endo[:y′_t3])] = 1.
    JJ[first(eq[:L3Y]),first(endo[:y_t2])]  = -1.

    JJ[first(eq[:L4Y]),first(endo[:y′_t4])] = 1.
    JJ[first(eq[:L4Y]),first(endo[:y_t3])]  = -1.

    JJ[first(eq[:LZ]),first(endo[:z′_t1])] = 1.
    JJ[first(eq[:LZ]),first(endo[:z_t])]   = -1.

    JJ[first(eq[:L2Z]),first(endo[:z′_t2])] = 1.
    JJ[first(eq[:L2Z]),first(endo[:z_t1])]   = -1.

    JJ[first(eq[:L3Z]),first(endo[:z′_t3])] = 1.
    JJ[first(eq[:L3Z]),first(endo[:z_t2])]  = -1.

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

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

    # discount factor shock
    JJ[first(eq[:F33]),first(endo[:B′])] = 1.
    JJ[first(eq[:F33]),first(endo[:B])]  = -ρ_B

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

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

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

    # wage mkup shock
    JJ[first(eq[:F37]),first(endo[:LAMW′])] = 1.
    JJ[first(eq[:F37]),first(endo[:LAMW])]  = -ρ_lamw

    # price mkup shock
    JJ[first(eq[:F38]),first(endo[:LAMF′])] = 1.
    JJ[first(eq[:F38]),first(endo[:LAMF])]  = -ρ_lamf

    # monetary policy shock
    JJ[first(eq[:F39]),first(endo[:MON′])] = 1.
    JJ[first(eq[:F39]),first(endo[:MON])]  = -ρ_mon

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

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

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

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

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

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

    m <= Setting(:n_predetermined_variables, size(Qx, 1))

	Jac1 = Qleft*sparse(JJ)*Qright

    return Jac1
end

function compose_normalization_matrices(m::HetDSGELag)
    nx::Int = get_setting(m, :nx)
    ns::Int = get_setting(m, :ns)
    nscalars::Int = get_setting(m, :nscalars)
    nyscalars::Int = get_setting(m, :nyscalars)
    nxscalars::Int = get_setting(m, :nxscalars)

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

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

    nxns = nx*ns

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

    return Qx, Qy, Qleft, Qright
end

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

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

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

        setup_indices!(m)

        endo = m.endogenous_states_unnormalized
        state_indices = [endo[:l′_t1]; endo[:μ′_t1];  endo[:k′_t]; endo[:R′_t1];endo[:i′_t1];endo[:π′_t1];
                         endo[:π′_t2];endo[:π′_t3];endo[:y′_t1];endo[:y′_t2];endo[:y′_t3];endo[:y′_t4];
                         endo[:z′_t1];endo[:z′_t2];endo[:z′_t3];endo[:w′_t1];endo[:I′_t1];endo[:B′];
                         endo[:G′];endo[:z′_t];endo[:MU′];endo[:LAMW′]; endo[:LAMF′];endo[:MON′]]
                         #=endo[:l_t1];
                         endo[:μ_t1];endo[:k_t];endo[:R_t1];endo[:i_t1]; endo[:π_t1];endo[:π_t2];
                         endo[:π_t3]; endo[:y_t1];endo[:y_t2];endo[:y_t3]; endo[:y_t4];endo[:z_t1];
                         endo[:z_t2]; endo[:z_t3];endo[:w_t1];endo[:I_t1]; endo[:B];endo[:G];endo[:Z];
                         endo[:MU]; endo[:LAMW];endo[:LAMF];endo[:MON]]=#

        jump_indices = [endo[:l′_t]; endo[:μ′_t];endo[:R′_t];endo[:i′_t];endo[:t′_t];endo[:w′_t];
                        endo[:L′_t];endo[:π′_t];endo[:wageinflation′_t];endo[:mu′_t];endo[:y′_t];
                        endo[:I′_t];endo[:mc′_t]; endo[:Q′_t];endo[:capreturn′_t]]
                        #=endo[:ELL];
                        endo[:μ_t];endo[:R_t];endo[:i_t];endo[:t_t];endo[:w_t]; endo[:L_t];endo[:π_t];
                        endo[:wageinflation_t];endo[:mu_t];endo[:y_t];endo[:I_t];endo[:mc_t];endo[:Q_t];
                        endo[:capreturn_t]]=#

        m <= Setting(:state_indices, state_indices, "Which indices of m.endogenous_states correspond to
                     backward looking state variables")
        m <= Setting(:jump_indices, jump_indices, "Which indices of m.endogenous_states correspond to jump
                     variables")

        n_backward_looking_vars = length(get_setting(m, :state_indices))
        n_backward_looking_function_valued_vars = get_setting(m, :n_function_valued_backward_looking_states)
        n_backward_looking_scalar_vars = n_backward_looking_vars - 2*nxns #n_backward_looking_function_valued_vars
        n_backward_looking_distr_vars = get_setting(m, :n_backward_looking_distributional_vars)
        m <= Setting(:n_backward_looking_states, nxns*n_backward_looking_distr_vars +
                     n_backward_looking_scalar_vars - get_setting(m, :backward_looking_states_normalization_factor))
        n_jump_vars = length(get_setting(m, :jump_indices))
        n_jump_function_valued_vars = get_setting(m, :n_function_valued_jumps)
        n_jump_scalar_vars = n_jump_vars - n_jump_function_valued_vars*nxns
        m <= Setting(:n_jumps, nxns*n_jump_function_valued_vars + n_jump_scalar_vars - get_setting(m, :jumps_normalization_factor))
        m <= Setting(:n_model_states, get_setting(m, :n_backward_looking_states) + get_setting(m, :n_jumps))
    end
end


text:
#
# fix-fitsio.jl --
#
# Some fixes and additional routines which should be part of FITSIO.jl package.
#
#------------------------------------------------------------------------------
#
# This file is part of OIFITS.jl which is licensed under the MIT "Expat"
# License:
#
# Copyright (C) 2015-2019: Éric Thiébaut.
#
#------------------------------------------------------------------------------

using FITSIO
using FITSIO.Libcfitsio

import FITSIO: libcfitsio, fits_assert_ok, fits_assert_open
import FITSIO.Libcfitsio: fits_get_errstatus

# The exported functions cfitsio_datatype and fits_bitpix deal with conversion
# between CFITSIO type code or BITPIX value and actual Julia data types.
# They can be used as follows (assuming `T` is a Julia data type, while
# `code` and `bitpix` are integers):
#
#     cfitsio_datatype(T) --------> code (e.g., TBYTE, TFLOAT, etc.)
#     cfitsio_datatype(code) -----> T
#
#     fits_bitpix(T) -------------> bitpix (e.g., BYTE_IMG, FLOAT_IMG, etc.)
#     fits_bitpix(bitpix) --------> T
#
export cfitsio_datatype, fits_bitpix

# The following table gives the correspondances between CFITSIO "types",
# the BITPIX keyword and Julia types.
#
#     -------------------------------------------------
#     BITPIX   CFISTIO         Julia     Comments
#     -------------------------------------------------
#              int             Cint
#              long            Clong
#              LONGLONG        Int64     64-bit integer
#     -------------------------------------------------
#        8     BYTE_IMG        UInt8
#       16     SHORT_IMG       Int16
#       32     LONG_IMG        Int32
#       64     LONGLONG_IMG    Int64
#      -32     FLOAT_IMG       Float32
#      -64     DOUBLE_IMG      Float64
#     -------------------------------------------------
#              TBIT
#              TBYTE           Cuchar = UInt8
#              TSBYTE          Cchar = Int8
#              TLOGICAL        Bool
#              TSHORT          Cshort
#              TUSHORT         Cushort
#              TINT            Cint
#              TUINT           Cuint
#              TLONG           Clong
#              TLONGLONG       Int64
#              TULONG          Culong
#              TFLOAT          Cfloat
#              TDOUBLE         Cdouble
#              TCOMPLEX        Complex{Cfloat}
#              TDBLCOMPLEX     Complex{Cdouble}
#     -------------------------------------------------

# BITPIX routines and table.
const _BITPIX = Dict{Cint, DataType}()
for (sym, val, T) in ((:BYTE_IMG,        8,       UInt8),
                      (:SHORT_IMG,      16,       Int16),
                      (:LONG_IMG,       32,       Int32),
                      (:LONGLONG_IMG,   64,       Int64),
                      (:FLOAT_IMG,     -32,       Float32),
                      (:DOUBLE_IMG,    -64,       Float64))
    val = Cint(val)
    _BITPIX[val] = T
    @eval begin
        fits_bitpix(::Type{$T}) = $val
    end
end
fits_bitpix(code::Integer) = get(_BITPIX, Cint(code), Nothing)

# Data type routines and table.
const _DATATYPE = Dict{Cint, DataType}()
const _REVERSE_DATATYPE = Dict{DataType, Cint}()
for (sym, val, T) in ((:TBIT       , Cint(  1), Nothing),
                      (:TBYTE      , Cint( 11), UInt8),
                      (:TSBYTE     , Cint( 12), Int8),
                      (:TLOGICAL   , Cint( 14), Bool),
                      (:TSTRING    , Cint( 16), Name),
                      (:TUSHORT    , Cint( 20), Cushort),          # Uint16
                      (:TSHORT     , Cint( 21), Cshort),           # Int16
                      (:TUINT      , Cint( 30), Cuint),            # Uint32
                      (:TINT       , Cint( 31), Cint),             # Int32
                      (:TULONG     , Cint( 40), Culong),
                      (:TLONG      , Cint( 41), Clong),
                      (:TFLOAT     , Cint( 42), Cfloat),           # Float32
                      (:TLONGLONG  , Cint( 81), Int64),
                      (:TDOUBLE    , Cint( 82), Cdouble),          # Float64
                      (:TCOMPLEX   , Cint( 83), Complex{Cfloat}),  # Complex64
                      (:TDBLCOMPLEX, Cint(163), Complex{Cdouble})) # Complex128
    _DATATYPE[val] = T
    if ! haskey(_REVERSE_DATATYPE, T)
        _REVERSE_DATATYPE[T] = val
        if T == AbstractString
            @eval cfitsio_datatype{S<:AbstractString}(::Type{S}) = $val
        elseif T != Nothing
            @eval cfitsio_datatype(::Type{$T}) = $val
        end
    end
end
cfitsio_datatype(code::Integer) = get(_DATATYPE, code, Nothing)


text:
export evaluate

function validate_size(recommender::Recommender, truth_data::DataAccessor)
    n_rec_user, n_rec_item = size(recommender.data.R)
    n_truth_user, n_truth_item = size(truth_data.R)

    if n_rec_user != n_truth_user
        error("number of users is mismatched: (recommenre, truth) = ($(n_rec_user), $(n_truth_user)")
    elseif n_rec_item != n_truth_item
        error("number of items is mismatched: (recommenre, truth) = ($(n_rec_item), $(n_truth_item)")
    end

    n_truth_user, n_truth_item
end

function evaluate(recommender::Recommender, truth_data::DataAccessor,
                  metric::AccuracyMetric)
    check_build_status(recommender)
    n_user, n_item = validate_size(recommender, truth_data)

    accum = 0.0

    for u in 1:n_user
        pred = zeros(n_item)
        for i in 1:n_item
            pred[i] = predict(recommender, u, i)
        end
        accum += measure(metric, truth_data.R[u, :], pred)
    end

    # return average accuracy over the all target users
    accum / n_user
end

function evaluate(recommender::Recommender, truth_data::DataAccessor,
                  metric::RankingMetric, k::Int=0)
    check_build_status(recommender)
    n_user, n_item = validate_size(recommender, truth_data)

    accum = 0.0

    candidates = Array(1:n_item)
    for u in 1:n_user
        truth = [first(t) for t in sort(collect(zip(candidates, truth_data.R[u, :])), by=t->last(t), rev=true)]
        recos = recommend(recommender, u, k, candidates)
        pred = [first(t) for t in sort(recos, by=t->last(t), rev=true)]
        accum += measure(metric, truth, pred, k)
    end

    # return average accuracy over the all target users
    accum / n_user
end


text:
module CalendarPuzzle

# ref: https://twitter.com/ataboh/status/1390710476800073733

using Dates

const BOARD_H = 8
const BOARD_W = 8
const PADDING_B = 1
const PADDING_R = 1

const PIECES = [:O, :L, :P, :S, :U, :V, :Y, :Z]
const PIECE_DATA = let dy=1, dx=BOARD_H
    Dict(
        :O => [
            cumsum([0, dy, dy, dx-2dy, dy, dy]),
            cumsum([0, dy, dx-dy, dy, dx-dy, dy]),
        ],
        :L => [
            cumsum([0, dy, dy, dy, dx]),
            cumsum([0, dy, dy, dy, dx-3dy]),
            cumsum([0, dx, dy, dy, dy]),
            cumsum([0, dx-3dy, dy, dy, dy]),
            cumsum([0, dx, dx, dx, dy]),
            cumsum([0, dy, dx-dy, dx, dx]),
            cumsum([0, dy, dx, dx, dx]),
            cumsum([0, dx, dx, dx-dy, dy]),
        ],
        :P => [
            cumsum([0, dy, dy, dx-2dy, dy]),
            cumsum([0, dy, dy, dx-dy, dy]),
            cumsum([0, dy, dx-dy, dy, dy]),
            cumsum([0, dy, dx-2dy, dy, dy]),
            cumsum([0, dy, dx-dy, dy, dx-dy]),
            cumsum([0, dx, dy, dx-dy, dy]),
            cumsum([0, dy, dx-dy, dy, dx]),
            cumsum([0, dx-dy, dy, dx-dy, dy]),
        ],
        :S => [
            cumsum([0, dy, dx, dy, dy]),
            cumsum([0, dy, dy, dx, dy]),
            cumsum([0, dy, dx-3dy, dy, dy]),
            cumsum([0, dy, dy, dx-3dy, dy]),
            cumsum([0, dx, dy, dx, dx]),
            cumsum([0, dx, dx, dy, dx]),
            cumsum([0, dx-dy, dy, dx-dy, dx]),
            cumsum([0, dx, dx-dy, dy, dx-dy]),
        ],
        :U => [
            cumsum([0, dy, dy, dx-2dy, 2dy]),
            cumsum([0, 2dy, dx-2dy, dy, dy]),
            cumsum([0, dy, dx-dy, dx, dy]),
            cumsum([0, dy, dx, dx-dy, dy]),
        ],
        :V => [
            cumsum([0, dy, dy, dx, dx]),
            cumsum([0, dy, dy, dx-2dy, dx]),
            cumsum([0, dx, dx, dy, dy]),
            cumsum([0, dx, dx-2dy, dy, dy]),
        ],
        :Y => [
            cumsum([0, dy, dy, dy, dx-2dy]),
            cumsum([0, dy, dy, dy, dx-dy]),
            cumsum([0, dx-dy, dy, dy, dy]),
            cumsum([0, dx-2dy, dy, dy, dy]),
            cumsum([0, dx, dy, dx-dy, dx]),
            cumsum([0, dx, dx, dy, dx-dy]),
            cumsum([0, dx-dy, dy, dx, dx]),
            cumsum([0, dx, dx-dy, dy, dx]),
        ],
        :Z => [
            cumsum([0, dy, dx, dx, dy]),
            cumsum([0, dy, dx-dy, dx-dy, dy]),
            cumsum([0, dx, dy, dy, dx]),
            cumsum([0, dx-2dy, dy, dy, dx-2dy]),
        ],
    )
end

function makeboard(
    board_size=(BOARD_H, BOARD_W),
    padding_bottom=PADDING_B,
    padding_right=PADDING_R
)
    board = fill(-1, board_size)
    field_h, field_w = board_size .- (padding_bottom, padding_right)
    board[1:field_h, 1:field_w] .= 0
    for (x, y) in fldmod1.([28, 35, 42, 43, 44, 49], 7)
        board[y, x] = -1
    end
    board
end

function makeboard(
    dt::Date,
    board_size=(BOARD_H, BOARD_W),
    padding_bottom=PADDING_B,
    padding_right=PADDING_R
)
    m, d = Dates.monthday(dt)
    board = makeboard(board_size, padding_bottom, padding_right)
    my, mx = fldmod1(m, 6)
    dy, dx = fldmod1(d + 14, 7)
    board[my, mx] = -1
    board[dy, dx] = -1
    board
end

function printboard(io::IO, board::AbstractMatrix, padding_bottom=PADDING_B, padding_right=PADDING_R)
    h, w = size(board) .- (padding_bottom, padding_right)
    function calcrow(y, row)
        w = length(row)
        strrow = if y < 3
            _r = lpad.((1:w) .+ (y - 1) .* 6, 3)
            _r[7:w] .= "   "
            _r
        else
            _ds = (1:w) .+ (y - 3) .* 7
            _r = lpad.(_ds, 3)
            _r[8:w] .= "   "
            _r[_ds .> 31] .= "   "
            _r
        end
        strrow[row .> 0] .= lpad.(PIECES[row[row .> 0]], 3)
        strrow
    end
    for y = 1:h
        println(io, calcrow(y, board[y, 1:w])...)
    end
end

availablepiece(board, pos, piece) = all(board[pos + idx] == 0 for idx in piece)

function putpiece(board, pos, piece, pidx)
    newboard = copy(board)
    newboard[piece .+ pos] .= pidx
    newboard
end

function find_nextpos(board, pos)
    findnext(eachindex(IndexLinear(), board), pos) do idx
        board[idx] == 0
    end
end

function solve_rec(fn, board, pos, pflags)
    if !any(pflags)
        # solved
        fn(board)
        return
    end
    pids = findall(pflags)
    for pid in pids
        for piece in PIECE_DATA[PIECES[pid]]
            if availablepiece(board, pos, piece)
                # recursive call
                newboard = putpiece(board, pos, piece, pid)
                newflags = copy(pflags)
                newflags[pid] = false
                nextpos = find_nextpos(newboard, pos)
                solve_rec(fn, newboard, nextpos, newflags)
            end
        end
    end
end

args2date() = Dates.today()
args2date(dt::Date) = dt
args2date(m::Int, d::Int) = Date(2020, m, d)  # Must be leap year!

function printresultfn(io::IO)
    function (board)
        printboard(io, board)
        println(io)
    end
end

function solve_all(fn, board)
    solve_rec(fn, board, find_nextpos(board, 1), trues(length(PIECES)))
end

function solve_one(fn, board)
    chnl = Channel{Matrix{Int}}(1) do chnl
        first_pos = find_nextpos(board, 1)
        solve_rec(board, first_pos, trues(length(PIECES))) do result
            put!(chnl, result)
        end
    end
    #
    try
        for result in chnl
            fn(result)
            return
        end
        println(stderr, "No Answers!")
    finally
        close(chnl)
    end
end

function hasanswer(args...)
    dt = args2date(args...)
    board = makeboard(dt)
    chnl = Channel{Matrix{Int}}(1) do chnl
        first_pos = find_nextpos(board, 1)
        solve_rec(board, first_pos, trues(length(PIECES))) do result
            put!(chnl, result)
        end
    end
    # 1つでも解が見つかれば終了
    try
        return !isnothing(iterate(chnl))
    finally
        close(chnl)
    end
end

function countanswers(args...)
    dt = args2date(args...)
    board = makeboard(dt)
    cnt = Ref(0)
    first_pos = find_nextpos(board, 1)
    solve_rec(board, first_pos, trues(length(PIECES))) do _
        cnt[] += 1
    end
    return cnt[]
end

function solve(args...; kwds...)
    dt = args2date(args...)
    fn = printresultfn(get(kwds, :io, stdout))
    board = makeboard(dt)
    if get(kwds, :all, false)
        solve_all(fn, board)
    else
        solve_one(fn, board)
    end
end

# solve()

function countallanswers()
    for dt in Date(2020,1,1):Day(1):Date(2020,12,31)
        print("#answer of $(dt): ")
        println(countanswers(dt), ".")
    end
end

function outputallanswers(destdir=pwd())
    for dt in Date(2020,1,1):Day(1):Date(2020,12,31)
        print("solving $(dt)...")
        destpath = joinpath(destdir, Dates.format(dt, dateformat"yyyymmdd.txt"))
        open(destpath, "w") do f
            solve(dt; io=f, all=true)
        end
        println("Done.")
    end
end

end # module


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# Reverse McCormick.jl
# A Reverse McCormick operator library in Julia
# See https://github.com/PSORLab/ReverseMcCormick.jl
#############################################################################
# src/reverse_operators/extrema.jl
# Contains definitions of reverse min and max.
#############################################################################

"""
$(SIGNATURES)

Reverse McCormick operator for `max(::MC, ::MC)`. Note that this function is not
invertible. If `b` dominates `c`, the  `b = b ∩ a`. If `c` dominates `b`,  `c = c ∩ a`.
No update otherwise. TODO: See if constructing an reverse from the identity
`max(x,y) = 0.5*(x + y + abs(x - y))` yields a tigher relaxation in this case.
"""
function max_rev(a::MC, b::MC, c::MC)
    isempty(a) && (return a, a, a)
    if hi(b) <= lo(c)
        c = c ∩ a
    elseif lo(b) >= hi(c)
        b = b ∩ a
    end
    return a, b, c
end

function max_rev(a::MC{N,T}, b::MC{N,T}, c::C) where {N, T<:RelaxTag, C<:NumberNotRelax}
    isempty(a) && (return a, a, a)
    if c < lo(a)
        b = b ∩ a
    end
    a, b, MC{N,T}(c)
end

function max_rev(a::MC{N,T}, c::C, b::MC{N,T}) where {N, T<:RelaxTag, C<:NumberNotRelax}
    isempty(a) && (return a, a, a)
    if c < lo(a)
        b = b ∩ a
    end
    a, MC{N,T}(c), b
end

"""
$(SIGNATURES)

Reverse McCormick operator for `min`. Note that this function is not
invertible. If `b` dominates `c`, the  `b = b ∩ a`. If `c` dominates `b`,  `c = c ∩ a`.
No update otherwise.  TODO: See if constructing an reverse from the identity
`min(x, y) = 0.5*(x + y - abs(x - y))` yields a tigher relaxation in this case.
"""
function min_rev(a::MC, b::MC, c::MC)
    isempty(a) && (return a, a, a)
    if hi(b) <= lo(c)
        b = b ∩ a
    elseif lo(b) >= hi(c)
        c = c ∩ a
    end
    a, b, c
end

function min_rev(a::MC{N,T}, b::MC{N,T}, c::C) where {N, T<:RelaxTag, C<:NumberNotRelax}
    isempty(a) && (return a, a, a)
    if c > hi(a)
        b = b ∩ a
    end
    a, b, MC{N,T}(c)
end

function min_rev(a::MC{N,T}, c::C, b::MC{N,T}) where {N, T<:RelaxTag, C<:NumberNotRelax}
    isempty(a) && (return a, a, a)
    if c > hi(a)
        b = b ∩ a
    end
    a, MC{N,T}(c), b
end