text:
using RobustAndOptimalControl: ControlSystems
using RobustAndOptimalControl, ControlSystems
hinorm(x) = hinfnorm2(x)[1]

sys = let
    _A = [-2.467030078096955 -1.7579086436230489 0.2642212618287617 -1.8139967744112642 0.9581654870157803 -1.1573282803258818 -2.409554294298775 1.2006959533714947 0.2661208318398341 0.9999177692959037 -0.3558865668962797 0.9416230396873052 -0.04719968633569187 1.9048861053182367 0.04105479772871481; 1.7913628319734 -2.762572481036251 1.1195768845996525 -0.5685580615144388 0.6726614103668934 1.0083399722116977 0.5812702064616616 -0.7829110624242633 -0.5724789585031453 0.2739159931739683 0.2757716078817949 -1.1772430269658822 -0.6381809965413997 -0.5567081736321586 0.3625746395689792; -0.005143392441108892 -1.7086991517875587 -1.717585106651333 1.5009119207685369 1.8738081845828842 -1.0939749683914888 -0.5264010960315267 0.22497823448076093 1.246044372509282 -0.5303445465251139 -0.5548425351067188 0.33884707949495035 0.5101563763162115 0.20539506245514677 0.26057781689650356; 0.604556580615853 1.2138308974005823 -2.059380841591183 -1.191665990395876 1.150735331516496 0.8110103414198591 -0.154085845809975 1.1119072467025846 0.11430616747450531 1.2819431451075372 0.921519933726758 0.4697140300022479 0.005096755657399703 -0.6758568389946026 0.23258764043708394; -0.3279280635059737 -0.3899618604207743 0.12657097470583745 0.8299061544706142 -2.6661102102915484 -0.33139258711550823 -0.012925000635212677 -0.42098301162071433 1.5144416025925511 -0.796229814982033 -0.0657636569172887 1.2193694935643702 -1.700730366738557 -0.480614226127818 0.5646449944635483; 0.2450793035368965 -1.8805339894110746 0.30667982435893143 0.9479179427233335 -1.2301806768500736 -3.136669168734015 -0.2481694826668052 0.380557579742639 -1.416193374954648 0.6987044103237035 -0.8853533466891377 -1.030250851383474 1.4525483702515944 0.36288849878688056 0.385057654662975; -2.0022001785391392 -0.2835451807794339 -0.26028704624052934 0.07217513465675776 0.08414471225474451 0.024542390758292842 -2.6763138550617667 0.011591697680752557 1.7133530919443034 0.17602929220838714 -0.6726695861811227 1.5441809892864893 -0.00015597605070589597 1.301000782058375 0.03551676106571688; -0.16249013289239514 0.13198122357867412 1.3236770560172366 0.3301510019655505 0.45179991512285006 0.07486588236983925 -1.6288494827036406 -1.5392724756838398 -0.6217318765063355 0.2858777881749128 -1.4166215418299406 -0.35825226416886896 1.8697114631330003 -0.3059163232549132 1.934355556146854; -0.3097986537718366 1.9026121572125902 -0.6714075159339921 -0.2504951318869757 -1.0960812903814519 0.551662724714743 -0.5206166054513874 -0.10057595046577751 -4.293450503414131 0.24751394898682735 0.5853194271595069 1.2157066833567851 -0.5437491323644809 0.03154589360470189 -1.3601989491821653; -1.8758967387335015 -0.08090769024414653 1.703446780809511e-5 -0.15226980080851651 -0.4474269369804122 -1.0627549839203736 1.1462306659249069 0.4825852860273927 1.1227319809099574 -2.5566746127781825 0.2744599619933535 1.514041120509779 0.8781499781073471 -0.7005275854271636 0.8540288603061493; 1.3102799390114126 2.0268412603985713 -0.26445369495117343 0.1673277757620624 0.33995168379334895 1.3424350975481887 -0.5262869773852599 1.772445455267548 0.7427203416097469 -0.6894548775913741 -4.648846973081957 0.6160480706048977 0.27693416312465485 -0.5843844597444985 -0.05046410580231136; -0.2667927151895562 -0.8394567537703295 0.7356804516400021 0.021352734365649183 1.6248465860949257 1.0375631106350516 -0.44557649819583367 -1.318671094466728 0.1198337174336245 0.028769936034846715 1.1651225420369067 -2.5802262184585376 -0.000961942820959872 0.3531297512430618 0.8801838685743447; -0.2340108920478918 -0.47020830048023965 0.0905085276240309 -1.159187923078835 -0.32263237957704327 -0.20398562485896923 -1.3882045935474685 0.6799676640930128 0.34600863452034686 -0.8950627077509554 0.37192323444628994 -0.17983996993004395 -3.254198004959321 -0.12498721862068969 -0.7674453121406706; -0.48890022859302995 0.7856693831160491 0.5530727968750053 0.28679439003867857 -0.353796629578632 1.3881939385941933 0.04935762134562992 -0.6880042763129538 1.3414828049193512 -1.6887109143288497 -0.846493928570519 -2.8372692994687916 -0.8703006934200097 -2.0491003542049913 -0.04580815825737827; 0.14776326405901308 -1.0181952707299127 0.17409898161459436 0.5045755253054325 1.2756549793646714 0.0513146326521474 1.9399237561334002 -1.3740316671021375 -1.6331290401813572 -0.4569927755409647 -0.6655757376846796 -0.11273958657543114 0.17478370341836297 0.9495603386771154 -2.069076828262977]
    _B = [-0.9890109448948233; 0.94489034588656; -0.6141178465083617; 0.16165201706201165; 0.3159060979974231; 0.48695887811075705; 1.4986210850273796; 0.20368746793706913; 0.30197970526128687; -0.8327748143726169; 0.5412418165171653; -0.6001356417667288; 0.9563775528650693; 3.0913194636365113; -0.34466316244242773]
    _C = [0.3132381087641251 -0.2983781324325293 -0.27755322605976035 0.722187317484002 0.10756083023764713 -0.8432926363233019 -1.2942110080251952 0.08989673121535467 -1.6845175332374567 -2.535783822800354 -0.6469091057919979 0.6817108268640213 -1.4395236987685638 -0.3003320643423357 0.49320084051761887]
    _D = [0.0]
    ss(_A, _B, _C, _D)
end
sysr = frequency_weighted_reduction(sys, 1, 1, 10, residual=false)
sysr2 = baltrunc(sys, n=10)[1]
@test sysr.nx == 10
@test hinorm(sys-sysr) < 1e-5
# bodeplot([sys, sysr])

# sys = ssrand(1,1,10, proper=true)
sys = let
    _A = [-3.2457568884577848 -0.9709881648628496 2.5986830066501487 0.10929577556477561 0.05027607854258857 0.27953443726273464 0.37869575253657956 -0.7415670467124911 -0.8723306000440615 0.9186219277892705; -0.67951689956761 -2.853540754527736 0.6121903826231613 2.167022775942329 -1.3362518751231087 0.27052292738707584 -0.41205997059440236 0.9675328597726843 -0.31902340467214135 -0.2624694054518293; -0.9693825264182581 0.2706161544633993 -1.433312733210014 -0.4195076769611032 0.4703049252849061 1.47290014682162 0.09970803641137543 -0.7189595131691185 -1.7372665057580614 0.6496704162176036; -0.7315165984188576 -3.438587829345989 -0.2646150665230416 -1.2369179161091948 -0.5714931604453887 0.2655752856563833 0.2244428073746226 1.3200417901754988 -0.5139418130660682 0.8869070876329531; 0.6678034518746887 -0.28242407386536494 -0.01306571590718952 0.8180780644546597 -1.5636655476823766 -1.432041330305812 -0.3186555690625707 0.4047899977465902 0.7766543740761176 1.1332916952645569; 1.0439462961677213 0.5603761347480476 -1.0673946370803042 2.065297667030172 2.0622680121683077 -1.5407850641702399 2.562320986628486 0.4291456610040989 -0.5853999459761943 0.17886564635348326; -0.13285626489018923 0.47503952348791406 1.105921739516141 1.0933062321124096 1.2703952887202579 -1.0825698120819383 -1.4243942957100129 -1.3994571653345012 0.6890508640402956 0.3244453993527589; -0.27988553587589793 -0.2865641537928437 -0.020379326141595905 -0.8113230885198695 -0.29383504852361897 -0.7637611675988282 0.29682298030397475 -4.065031083750159 0.5691551157574606 -1.8639256107234736; 0.4537930380958493 0.6614234998446303 -0.4277870018921912 1.7342323590344288 0.2542882624030592 0.002642962031894451 -1.1148591940907127 -0.9394549359866202 -2.262544410098742 -0.9420467603687506; -1.088408576121221 -0.13854618864389867 -0.9251283170605236 1.5905768356892043 0.011581927582191199 -2.4924299911211043 -0.4921574784104917 -0.10267947834304443 -0.5622726271236202 -2.499827020287407]
    _B = [2.109834811943383; 0.7122212935939083; 0.9477960944077879; -1.123537196617033; -0.5570048263462403; -0.6970324371311861; -0.6094269119334326; -0.5234946661946995; -0.18596966094535153; -0.4420098228227759]
    _C = [0.3210661400839473 -1.4299373597711122 -2.193683276878009 -1.7203675776227303 -0.8359692872531289 -0.438485194671709 -0.6987593761582938 -0.7233564542540655 -0.9744869837106208 -0.31839951662879357]
    _D = [0.0]
    ss(_A, _B, _C, _D)
end
sysi = fudge_inv(sys)
sysr = frequency_weighted_reduction(sys, sysi, 1, 5, residual=false)
sysr2 = baltrunc(sys, n=5)[1]
@test sysr.nx == 5
@test norm(sys-sysr) < 0.1
@test norm(sysi*(sys-sysr)) < 0.05
@test norm(sysi*(sys-sysr)) <= norm(sysi*(sys-sysr2))
@test hinorm(minreal(sysi*(sys-sysr))) <= hinorm(minreal(sysi*(sys-sysr2)))

Wi = tf(1, [1, 1])
sysr = frequency_weighted_reduction(sys, 1, Wi, 5, residual=false)
sysr2 = baltrunc(sys, n=5)[1]
@test sysr.nx == 5
@test norm(sys-sysr) < 0.1
@test norm((sys-sysr)*Wi) < 0.05
@test norm((sys-sysr)*Wi) <= norm((sys-sysr2)*Wi)
@test hinorm((sys-sysr)*Wi) <= hinorm((sys-sysr2)*Wi)

Wo = sysi
Wi = tf(1, [1, 1])
sysr = frequency_weighted_reduction(sys, Wo, Wi, 5, residual=false)
sysr2 = baltrunc(sys, n=5)[1]
@test sysr.nx == 5
@test_broken norm(sys-sysr) < 0.1
@test norm(Wo*(sys-sysr)*Wi) < 0.05
# @test_broken norm(Wo*(sys-sysr)*Wi) <= norm(Wo*(sys-sysr2)*Wi)
@test_broken hinorm(Wo*(sys-sysr)*Wi) <= hinorm(Wo*(sys-sysr2)*Wi)

# residual
sysi = fudge_inv(sys)
sysr = frequency_weighted_reduction(sys, sysi, 1, 3, residual=true)
sysr2 = baltrunc(sys, n=3, residual=true)[1]
@test sysr.nx == 3
@test norm(sys-sysr, Inf) < 3
@test norm(sysi*(sys-sysr), Inf) < 0.4
@test hinorm(minreal(sysi*(sys-sysr))) <= hinorm(minreal(sysi*(sys-sysr2)))
@test dcgain(sys)[] ≈ dcgain(sysr)[] rtol=1e-10 # test the residual property

## Test from paper https://elib.dlr.de/11746/1/varga_cdc01p2.pdf
G = tf([8, 6, 2], [1, 4, 5, 2])
Wi = tf(1, [1, 3])
Wo = tf(1, [1, 4])
sysr = frequency_weighted_reduction(ss(G), Wo, Wi, 1, residual=true)
@test isstable(sysr)

@test tf(sysr) ≈ tf([2.398, 1.739], [1, 1.739]) rtol=1e-1
@test norm(Wo*(G-sysr)*Wi, Inf) ≈ 0.0855 rtol=0.02

# Example 3
A = diagm(-1:-1:-4)
B = [
    0 5
    1/2 -3/2
    1 -5
    -1/2 1/6
]
C = [
    1 0 1 0
    4/15 1 0 1
]
G = ss(A,B,C,0)
Aw = -4.5I(2)
Bw =    3I(2)
Cw =  1.5I(2)
Dw =     I(2)
Wo = Wi = ss(Aw,Dw,Cw,Dw)

errors = map(1:3) do r
    Gr = frequency_weighted_reduction(G, Wo, Wi, r, residual=false)
    norm(Wo*(G-Gr)*Wi, Inf)
end

@test errors[1] <= 1.1 * 2.12
@test_broken errors[2] <= 1.1 * 0.265
@test errors[3] <= 1.1 * 0.112

errors = map(1:3) do r
    Gr = frequency_weighted_reduction(G, Wo, Wi, r, residual=true)
    norm(Wo*(G-Gr)*Wi, Inf)
end

@test errors[1] <= 1.1 * 1.405
@test_broken errors[2] <= 1.1 * 0.25
@test errors[3] <= 1.1 * 0.065

## Controller reduction

# A = [
#     -0.015948101119319734 0.0 0.04185849126302005 0.0
#     0.0 -0.011069870050970518 0.0 0.03334113387614947
#     0.0 0.0 -0.04185849126302005 0.0
#     0.0 0.0 0.0 -0.03334113387614947
# ]
# B = [0.08325 0.0; 0.0 0.0628125; 0.0 0.04785714285714286; 0.031218750000000007 0.0]
# C = [0.5 0.0 0.0 0.0; 0.0 0.5 0.0 0.0]
# D = [0.0 0.0; 0.0 0.0]
# G = ss(A, B, C, D)

# A = [-0.001 0.0; 0.0 -0.001]
# B = [1.0 0.0; 0.0 1.0]
# C = [0.0999 0.0; 0.0 0.0999]
# D = [0.1 0.0; 0.0 0.1]
# WS = ss(A, B, C, D)

# D = [0.01 0.0; 0.0 0.01]
# WU = ss(D)

# A = [-1.0 0.0; 0.0 -1.0]
# B = [1.0 0.0; 0.0 1.0]
# C = [-9.9 0.0; 0.0 -9.9]
# D = [10.0 0.0; 0.0 10.0]
# WT = ss(A, B, C, D)
# P = hinfpartition(G, WS, WU, WT)

# flag, C, γ = hinfsynthesize(P, γrel=1)

# Pcl, S, CS, T = hinfsignals(P, G, C)

# Cr = controller_reduction(P,C,7, false)
# Cr2 = baltrunc(C,n=7)[1]
# Pclr, Sr, CSr, Tr = hinfsignals(P, G, Cr)
# Pclr2, Sr2, CSr2, Tr2 = hinfsignals(P, G, Cr2)
# bodeplot([Pcl, Pclr], plotphase=false, size=(1900,920))
# bodeplot([C, Cr, Cr2])

# hinfn = ControlSystems._infnorm_two_steps_ct(minreal(Pcl-Pclr), :hinf, 1e-9, 1000, 1e-6)[1]
# @test hinfn ≈ 14.88609988 rtol=1e-3
# @test isstable(Pclr)

# IRKA
# for n = 4:10
#     r = n÷2
#     sys = ssrand(rand(1:2), rand(1:2), n, proper=true)
#     sysr0 = baltrunc(sys, n=r)[1]
#     sysr = RobustAndOptimalControl.model_reduction_irka(sys, r; tol = 1e-6)
#     @test abs(norm(sys-sysr)) < norm(sys-sysr0)
# end


text:
@reexport module MultiArmBandits

export MultiArmBanditsEnv

using ReinforcementLearningBase

mutable struct MultiArmBanditsEnv <: AbstractEnv
    truevalues::Vector{Float64}
    truereward::Float64
    bestaction::Int
    isbest::Bool
    reward::Float64
    isdone::Bool
    observation_space::DiscreteSpace
    action_space::DiscreteSpace
    function MultiArmBanditsEnv(; truereward::Float64 = 0.0, k::Int = 10)
        truevalues = truereward .+ randn(k)
        new(
            truevalues,
            truereward,
            findmax(truevalues)[2],
            false,
            0.0,
            false,
            DiscreteSpace(1),
            DiscreteSpace(length(truevalues)),
        )
    end
end

## interfaces

RLBase.get_observation_space(env::MultiArmBanditsEnv) = env.observation_space
RLBase.get_action_space(env::MultiArmBanditsEnv) = env.action_space

function RLBase.reset!(env::MultiArmBanditsEnv)
    env.isbest = false
    nothing
end

function (env::MultiArmBanditsEnv)(a::Int)
    env.isbest = a == env.bestaction
    env.reward = randn() + env.truevalues[a]
    env.isdone = true
    nothing
end

RLBase.observe(env::MultiArmBanditsEnv) =
    (state = 1, terminal = env.isdone, reward = env.reward)

end

ext: jl


text:
#
# patches.jl -
#
# Methods that should be part of FITSIO or CFITSIO.
#

# Should be in FITSIO.jl/src/header.jl
"""
    get(hdr, key, def)

yields value of keyword `key` in FITS header `hdr` or `def` if not found.

"""
Base.get(hdr::FITSHeader, key::String, def) =
    ((i = get(hdr.map, key, -1)) > 0 ? hdr.values[i] : def)

# Should be in CFITSIO.jl/src/CFITSIO.jl
Base.isopen(f::FITSFile) = (f.ptr != C_NULL)


text:
function do_hello()
    comm = MPI.COMM_WORLD
    println("Hello world, I am $(MPI.Comm_rank(comm)) of $(MPI.Comm_size(comm))")
    MPI.Barrier(comm)
end


text:
export σx, σy, σz, 𝟙, up, down, ⊗, speye, identity_op, single_spin_op, correlator, op_list, symmetrize_state, symmetrize_op, nspins

const σx = sparse([0 1; 1 0])
const σy = sparse([0 -im; im 0])
const σz = sparse([1 0; 0 -1])
const 𝟙 = spdiagm([1, 1])

const up = [1,0]
const down = [0,1]

const ⊗ = kron

speye(k::Integer) = spdiagm([1 for _ in 1:k])
identity_op(k::Integer) = speye(2^k)

"""
    nspins(obj)

Give number of spins in the system this object belongs to.
"""
function nspins(system_size::Int)
    N = Int(ceil(log2(system_size)))
    2^N == system_size || throw(ArgumentError("$(system_size) is not a power of 2!"))
    N
end
nspins(state::AbstractVector) = nspins(length(state))
nspins(op::AbstractMatrix) = size(op,1) == size(op, 2) ? nspins(size(op, 1)) : throw(ArgumentError("The operator is not a square-matrix!"))

single_spin_op(op, k::Integer, N::Integer) = identity_op(k-1) ⊗ op ⊗ identity_op(N-k)
correlator(op, i::Integer, j::Integer, N::Integer) = i > j ? correlator(op, j, i, N) : identity_op(i-1) ⊗ op ⊗ identity_op(j-i-1) ⊗ op ⊗ identity_op(N-j)

op_list(op, N::Integer) = [single_spin_op(op, k, N) for k in 1:N]


text:
#
# Declare abstract types
#
abstract type IEntity end
abstract type IModification <: IEntity end


text:
@testcase "issue 62 and 98" begin
    function test_substring_match(x::AbstractString)
        @match x begin
            "1" => 1
            _   => 0
        end
    end
    @test test_substring_match("1") == 1
    @test test_substring_match("11") == 0

    @test test_substring_match(SubString("11", 2)) == 1
    @test test_substring_match(SubString("11", 1)) == 0
end


text:
function prepcomp(;companame=COMP_A_NAME, compqname = COMP_Q_NAME, compname = COMP_NAME,
   datapath=DATA_PATH, refreshcomp=true,
   incsvstream::Function = IN_CSV_STREAM,
   injlsstream::Function = IN_JLS_STREAM,
   outjlsstream::Function = OUT_JLS_STREAM,
   usequarterlydates::Bool=USE_QUARTERLY_DATES,
   parallel::NBool=nothing,
   retainedcolscomp::Vector{Symbol} = RETAINED_COLS_COMP,
   archive::Bool = true)::DataFrame

   (!isnothing(parallel)) && (PARALLEL[] = parallel)

  local compa::DataFrame
  local compq::DataFrame
  local comp::DataFrame

  #makes a serialization object so we don't have to keep parsing the CSV
  if refreshcomp || (!isfile("$datapath\\$compname.jls.lz4"))
    compa = incsvstream("$datapath\\$companame.csv.lz4") |> CSV.read |> DataFrame
    compq = incsvstream("$datapath\\$compqname.csv.lz4") |> CSV.read |> DataFrame

    comp = preprocesscomp(compa, compq, usequarterlydates=usequarterlydates)
    comp = mergecompccm!(comp)
    comp = select!(comp,retainedcolscomp)
    outjlsstream("$datapath\\$compname.jls.lz4") do s
      serialize(s, comp) #NOTE: DEBUG ONLY
    end

    archive && archivefile("$compname.jls.lz4") #never hurts to have a backup

  else
    try
      injlsstream("$datapath\\$compname.jls.lz4") do s
        comp = deserialize(s) #NOTE: DEBUG ONLY
      end
    catch err
      @warn "Could not load $compname. Error: $err"
      @info "Attempting to restore back-up and load from archive..."
      unarchivefile("$compname.jls.lz4")
      injlsstream("$datapath\\$compname.jls.lz4") do s
        comp = deserialize(s) #NOTE: DEBUG ONLY
      end
      @info "Successfully restored from archive."
    end

  end

  @info "Compustat and CCM linking data loaded and saved into file $compname.jls.lz4"

  return comp
end

#helper function to clean up names in the quarterly file
function propername(s::Symbol)
  str::String = string(s)

  if str[end] == 'q' || (length(str) ≥ 3 && str[(end-2):end] == "qtr")
    return s
  end

  s = Symbol(s,:q)
  return s
end

#this processes both the q and the a comp frames, with as much commonly used code as possible
function preprocesscomp(compa::DataFrame, compq::DataFrame;
  compdateformat::DateFormat = COMP_DATE_FORMAT,
  testoutput::Bool = TEST_OUTPUT,
  usequarterlydates::Bool = true)

  local comp::DataFrame
  #local usedcompfields::Vector{Symbol}

  #do operations that are common to both the annual and the quarterly data
  compa = preprocesscomp!(compa) #methods common to compa and compq
  compa = preprocesscompa!(compa) #methods specific to compa

  if usequarterlydates #need to merge this in if we are using it
    compq = preprocesscomp!(compq)
    compq = preprocesscompq!(compq)

    #(if we were using quarterly data it would be joined in here
    comp = join(compa, compq, on=[(:gvkey, :gvkeyq),
      (:fyear, :fyearq), (:indfmt,:indfmtq)], kind=:left)

    #the below performs a lookup and makes sure each row has an
    #actual quarterly announcement date if available
    gdf::GroupedDataFrame = groupby(comp, [:gvkey, :fyear])
    @mpar for i ∈ 1:length(gdf)
      scomp::SubDataFrame = gdf[i]

      for r ∈ eachrow(scomp)
        local adateq::MDate = r.adateq
        (!ismissing(adateq)) && (r.adate = adateq)
        if (size(scomp,1) > 1) && (!ismissing(adateq))
          scomp.adate .= ((ad::MDate,adq::MDate)->
            (ismissing(adq)) ? adateq : ad).(scomp.adate, scomp.adateq)
        end
      end
    end

    comp.adate = ((ad::MDate,adq::MDate)->
      ismissing(adq) ? ad : adq).(comp.adate, comp.adateq)

      #usedcompfields = [USED_COMPA_VALUE_FIELDS; USED_COMPQ_VALUE_FIELDS]
  else
    comp = compa
    #usedcompfields = USED_COMPA_VALUE_FIELDS
  end

  println("rows: $(size(comp,1))")
  filter!(comp) do r::DataFrameRow
      r[:adate] ≥ eom(r[:fyrmonth])
  end

  println("rows adate≥fdate: $(size(comp,1))")
  validatefadates!(comp; Fgroup=:gvkey, Ffdate=:fyear, Fadate=:adate)
  lagwithin!(comp, :fdate, :gvkey, :fyear)

  #form the initial interval candidates
  #these are accurate if the link doesn't changes (adjust for this case later)
  comp = formannouncementintervals!(comp)

  #println(describe(comp))
  if testoutput
    d = Date(2015,12,31)
    testoutput && CSV.write("output\\comp_$d.csv", comp[((bd,ed)->
      (d > bd) && (d ≤ ed)).(comp.begindate, comp.enddate), :])
  end

  return comp
end

#operations that are common to both the annual and the quarterly data
function preprocesscomp!(comp::DataFrame)
  #make all column names lower case
  rename!(comp, (s::Symbol->Symbol(lowercase(string(s)))).(names(comp)))

  #format dates
  comp.ddate = (s->parsecomp(Date, s)).(comp.datadate)
  comp.dyear = (d::MDate-> ismissing(d) ? missing : Dates.year(d)).(comp.ddate)

  comp = comp[(!ismissing).(comp.dyear), :] #require the calendar date

  rename!(comp, :indfmt=>:indfmtold)
  comp.indfmt = (s::MString->parsecomp(Symbol, s)).(comp.indfmtold)
  select(comp, Not(:indfmtold))

  return comp
end

function preprocesscompa!(compa::DataFrame;
  usedcompafields::Vector{Symbol} = USED_COMPA_VALUE_FIELDS,
  yearrange::UnitRange{Int} = YEAR_RANGE[])

  #try to use the more efficient year-month format when possible
  compa = compa[(!ismissing).(compa.fyear), :]
  compa.fyrmonth = ((y::Integer,m::Integer)->
    YearMonth(y,m)).(compa.fyear, compa.fyr)
  compa.fdate = (eom).(compa.fyrmonth)

  filter!(r::DataFrameRow->year(r.fdate) ∈ yearrange, compa) #restrict to desired years

  #sets the fixed lead of each entry
  @inline function ayrmonth(ym::YearMonth)::Date
    Δ::Int = 4 + (ym.y < 2004)
    return eom(ym + Δ)
  end

  compa.adate = (ayrmonth).(compa.fyrmonth)

  compa = dedupcomp!(compa, [:gvkey, :fyear], usedcompafields, :indfmt, :adate)

  #make the cash field
  compa.cash = ((che::MFloat64, ch::MFloat64)-> coalesce(che, ch)).(compa.che, compa.ch)
  compa = compa[(!ismissing).(compa.cash), :] #this is an essential field
  compa.bknegcash = 1.0 .- (compa.cash ./ compa.at)

  #compute #flev field
  #no negative values
  @inline function maxval0!(df::AbstractDataFrame, s::Symbol)
    df[!,s] .= (f::MFloat64->max(f,0.0)).(compa[!,s])
  end

  #=@mpar for s ∈ [:dltt, :dlc, :ceq, :txditc, :txdb, :pstk]
    maxval0!(compa, s)
  end=#

  compa.netdebt = compa.dltt .+ compa.dlc .- compa.cash
  compa.bkequity = compa.ceq .- compa.pstk .- ((txditc::MFloat64, txdb::MFloat64)->
    coalesce(txditc, txdb, 0.0)).(compa.txditc,compa.txdb)
  @mpar for s ∈ [:netdebt, :bkequity] #WARNING TESTING THIS
    maxval0!(compa, s)
  end
  compa.bkflev = compa.netdebt ./ (compa.netdebt .+ compa.bkequity)

  #compute liab
  compa.netliab = compa.lt .- compa.cash
  compa.bkliab = compa.netliab ./ compa.at

  #gross profit is proxy for operating profits as in Novy Marx 2013
  compa.op = (r::DataFrameRow->ismissing(r.gp) ? r.revt - r.cogs : r.gp).(eachrow(compa))
  #maxval0!(compa, :op)  #WARNING TESTING THIS
  compa.nop = compa.op ./ compa.bkequity #this is used for the operating profit ratio

  #compa = compa[completecases(compa[:, [:bknegcash, :bkflev, :bkliab]]),:]

  #make finite
  @mpar for s ∈ [:bknegcash, :bkflev, :bkliab, :nop]
    compa[!,s] = (f::MFloat64->coalesce(isfinite(f),false) ? f : missing).(compa[!,s])
  end

  #get the market value if available
  compa.mkequity = ((mkvalt::MFloat64, pso::MFloat64)->
    coalesce(mkvalt, pso)).(compa.mkvalt, compa.prcc_f .* compa.csho)

  #NOTE: we won't translate these into our measures until after we do the primary merge

  #sort the data
  sort!(compa, [:gvkey, :adate])
  issorted(compa, [:gvkey, :fyear]) || error("compa not sorted on fyear")
  lagwithin!(compa, [:at], :gvkey, :fyear, sorted=true)

  #Now make the BRS metrics
  compa.Aat = (r::DataFrameRow->coalesce((r.at + r.Lat)/2, r.at, r.Lat)).(eachrow(compa))
  compa.brse = (compa.sstk .- compa.prstkc .- compa.dv) ./ compa.Aat
  compa.brsd = (r::DataFrameRow->(r.dltis - r.dltr - coalesce(r.dlcch,0.))/r.Aat).(eachrow(compa))
  compa.brsx = compa.brse .+ compa.brsd

  (f::Symbol->finiteormissing!(compa, f)).([:brse, :brsd, :brsx])

  #println(describe(compa))

  return compa
end

function preprocesscompq!(compq::DataFrame;
  usedcompqfields::Vector{Symbol} = USED_COMPQ_VALUE_FIELDS,
  month2stale::Int = MONTHS_2_STALE_SHORT,
  yearrange::UnitRange{Int} = YEAR_RANGE[])
  #use a consistent naming convention
  rename!(compq, (propername).(names(compq)))

  #helper function to convert a fiscal year-quarter to a datw
  @inline function fyearqtr2date(yq::YearQuarter, fyr::Int)
    ym::YearMonth = YearMonth(yq.y, fyr)
    offsetmonths::Int = -(4 - yq.q)*3
    ym += offsetmonths
    return eom(ym)
  end

  #pre-proess date info
  compq = compq[((fq)->(!ismissing(fq))).(compq.datafqtr),:]
  compq.fyearqtr =  (s->parsecomp(YearQuarter, s)).(compq.datafqtr)
  #compq.fyearqtr =  (s->parsecomp(Symbol, s)).(compq.indfmtq)
  compq.fdateq = ((yq::YearQuarter,fyr::Int)->
    fyearqtr2date(yq,fyr)).(compq.fyearqtr, compq.fyrq)

  filter!(r::DataFrameRow->year(r.fdateq) ∈ yearrange, compq) #restrict to desired years

  #WARNING: We drop all quarters except Q4 here since we only want annnual data
  compq = compq[(yq::YearQuarter->yq.q==4).(compq.fyearqtr),:]

  compq.adateq = (s->parsecomp(Date, s)).(compq.rdq)
  compq = compq[(!ismissing).(compq.adateq),:]
  compq = compq[((yq::YearQuarter, fyearq::Int)->yq.y == fyearq).(compq.fyearqtr, compq.fyearq),:]

  #NOTE: REMINDER- convert to first day of month for date math
  #Date math works better using the first day of the month,
  #even though the last day of the month is technically correct

  #don't let the fiscal date be too far from the announcement date
  period2stale::Month = Month(month2stale)
  compq = compq[(r->(r.adateq - period2stale ≤ r.fdateq)).(eachrow(compq)),:]

  validatefadates!(compq; Fgroup=:gvkeyq, Ffdate=:fyearqtr, Fadate=:adateq)

  compq = dedupcomp!(compq, [:gvkeyq, :fyearq], usedcompqfields, :indfmt, :adateq)
  compq = dedupcomp!(compq, [:gvkeyq, :adateq], usedcompqfields, :indfmt, :fyearq)
  #now process cases where firms have duplicate date-firm key pairs
  #the date field is used to affirm the latest record
  #(for annual records, that is the annual date,
  #for quarterly its the announcement date)
  #compq = dedupcomp!(compq, [:gvkeyq, :adateq], usedcompqfields, :indfmtq, :fyearq)

  return compq
end

  #under no circumstances allow restatements past the subsequent quarter's announcement date
function validatefadates!(comp::DataFrame; Fgroup::NSymbol=nothing, Ffdate::NSymbol=nothing,
  Fadate::NSymbol=nothing, sort::Bool = true)::Nothing

  (isnothing(Fgroup)) && error("Fgroup missing from validatefadates")
  (isnothing(Ffdate)) && error("Ffdate missing from validatefadates")
  (isnothing(Fadate)) && error("Fadate missing from validatefadates")

  sort!(comp, [Fgroup, Fadate])
  T::Type = eltype(comp[!,Ffdate])
  comp.maxfdate_at_adate = Vector{Union{T,Missing}}(undef, size(comp,1))

  compgrp::GroupedDataFrame = groupby(comp, Fgroup)
  @mpar for i ∈ 1:length(compgrp)
    subdf::SubDataFrame = compgrp[i]
    maxfdate::T = subdf[1,Ffdate]
    for r ∈ eachrow(subdf)
      (r[Ffdate] > maxfdate) &&  (maxfdate = r[Ffdate])
      r.maxfdate_at_adate = maxfdate
    end
  end

  filter!(comp) do r::DataFrameRow
    r.maxfdate_at_adate ≤ r[Ffdate]
  end

  select!(comp, Not(:maxfdate_at_adate))

  return nothing
end

function dedupcomp!(comp::DataFrame,
  groupbyfields::Vector{Symbol}, usedcompfields::Vector{Symbol},
  formatfield::Symbol, datefield::Symbol)

  #keepflags
  comp.keep = trues(size(comp,1))

  #use this for deduping
  comp.valuesmissings = Vector{MFloat64}(undef, size(comp,1))

  ncols::Int = size(comp,2)
  scomps::GroupedDataFrame = groupby(comp, groupbyfields)
  @mpar for i ∈ 1:length(scomps)
    scomp::SubDataFrame = scomps[i]

    if size(scomp,1) > 1

      #keep the most recent version
      #keep the most recent version
      maxdate::Union{Int, YearQuarter, Date} = maximum(scomp[!,datefield])
      mindate::Union{Int, YearQuarter, Date} = minimum(scomp[!,datefield])
      #println("got here")

      #why could this happen? not clear to me, but it does happen a bit
      if maxdate ≠ mindate
        scomp.keep .= false
      end
      #println("got here")
      #scomp.keep .= (yd->isequal(yd,maxdate)).(scomp[!,datefield])

      #special case where we have two records of two different formats (indl and fs)
      #in this case, keep the format with the least missing values
      if sum(scomp.keep) > 1
        sscomp::SubDataFrame = view(scomp, scomp.keep, :)
        Nsscomp::Int = size(sscomp, 1)

        for k ∈ 1:Nsscomp
          sscomp.valuesmissings[k] = sum((c->ismissing(sscomp[k,c])).(usedcompfields))
        end

        minmissing::Int = minimum(sscomp.valuesmissings)
        sscomp.keep .= (sscomp.valuesmissings .== minmissing)

        #tie breaker is to keep the INDL record (industrial since we are looking at non-financials)
        if sum(sscomp.keep) > 1
          sscomp.keep .= ((keep::MBool, form::MSymbol)->
            (!ismissing(form)) && keep && (form == :INDL)).(sscomp.keep, sscomp[!, formatfield])
        end

      end

      #if we still have a dup, delete all the records
      if sum(scomp.keep) > 1
        @warn "dedup failed group fields:
          $groupbyfields with values $((f->scomp[1,f]).(groupbyfields)) Records removed.\n
          $scomp"
        scomp.keep .= false
      end
    end
  end

  comp = comp[comp.keep, :] #bandaid until filter is fixed

  select!(comp, Not([:keep, :valuesmissings]))

  return comp
end

#NOTE: messy function
#basically to merge in crsp, the crsp day or month must occur
#in the intersection of the link dates, the period six months following
#the announcement date, and the next announcement date
#this handles everything but the intersection with the link dates
function formannouncementintervals!(comp::DataFrame;
    lastdate::Date = LAST_DATE, months2stale::Int =  MONTHS_2_STALE_LONG)

  sort!(comp, [:gvkey, :adate])

  comp.Nadate = Vector{MDate}(undef, size(comp,1))

  scomps = groupby(comp, :gvkey)
  @mpar for i ∈ 1:length(scomps)
    scomp::SubDataFrame = scomps[i]
    Nsub::Int = size(scomp,1)

    for i ∈ 1:(Nsub-1)
      scomp[i,:Nadate] = scomp[i+1,:adate]

      #check some things that shouldn't happen
      (scomp[i,:adate] ≥ scomp[i,:Nadate]) && @warn (
        "announcement dates out of order.\n $(scomp[i:(i+1),:])")

      (scomp[i,:fyear] ≥ scomp[i+1,:fyear]) && @warn (
        "fiscal year dates out of order.\n $(scomp[i:(i+1),:])")
    end

    #final vliad interval extends to the present
    scomp[Nsub, :Nadate] = lastdate
  end

  comp.fdatemax = (eom).(comp.fyrmonth .+ months2stale)

  comp.begindate = comp.adate
  comp.enddate = (min).(comp.Nadate,comp.fdatemax)

  return comp
end

parsecomp(::Type{Date}, s::String;
    compdateformat::DateFormat = COMP_DATE_FORMAT) = Dates.Date(s, compdateformat)

#helper methods and special cases
parsecomp(::Type{<:Any}, v::Missing) = missing
parsecomp(::Type{Date}, i::Int) = parsecomp(Date, "$i")
parsecomp(::Type{Symbol}, s::String) = Symbol(s)

#parse a value where the type is known: Generic case
function parsecomp(::Type{T}, s::String;
    ##NOTE: uncomment when we have the missings parsedmissings = COMP_PARSED_MISSINGS
    ) where T

  #if it doesn't parse to the right type, set to missing
  v::Union{T,Missing} = something(tryparse(T, s), missing)

  #check against the list of missing codes
  #NOTE: uncomment when we have the missings
  #(!ismissing(v)) && (v ∈ parsedmissings) && (v=missing)

  return v
end

function parsecomp(::Type{YearQuarter}, s::String)
  local yq::MYearQuarter

  if length(s) ≠ 6
    yq = missing
  else
    y::MInt = parsecomp(Int, "$(s[1:4])")
    q::MInt = parsecomp(Int, "$(s[6])")
    (q==1) || (q==2) || (q==3) || (q==4) || (q=missing)

    yq = (!ismissing(y)) && (!ismissing(q)) ? YearQuarter(y,q) : missing
  end

  return yq
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
import LinearAlgebra: ldiv!, mul!
using LinearAlgebra: det, UniformScaling, isdiag
import Base: length, getindex, setindex!, deleteat!

import NeighbourLists: cutoff

import ACEbase: fltype, fltype_intersect, rfltype,
                alloc_temp, alloc_temp_d, alloc_temp_dd,
                evaluate, evaluate_d, evaluate_dd, evaluate_ed,
                evaluate!, evaluate_d!, evaluate_dd!, evaluate_ed!,
                precon!

# function defined primarily on AbstractAtoms
export positions, get_positions, set_positions!,
       momenta, get_momenta, set_momenta!,
       masses, get_masses, set_masses,
       cell, get_cell, set_cell!, is_cubic, cell_vecs,
       pbc, get_pbc, set_pbc!,
       set_data!, get_data, has_data,
       set_calculator!, calculator, get_calculator,
       neighbourlist, cutoff,
       apply_defm!,
       energy, potential_energy, forces, gradient, hessian,
       dipole,
       site_energies,
       site_energy, partial_energy,
       site_energy_d, partial_energy_d,
       stress, virial, site_virials,
       position_dofs, set_position_dofs!,
       momentum_dofs, set_momentum_dofs!,
       dofs, set_dofs!,
       preconditioner,
       volume,
       chemical_symbols,
       atomic_numbers,
       AbstractAtoms, AbstractCalculator

# temporary prototypes while rewriting

# todo: prototype and document these here
function chemical_symbols end
function atomic_numbers end
function site_energy_d end
function partial_energy_d end

"""
`function dipole end ` : prototype for a global dipole function,
not currently implemented in JuLIP.
"""
function dipole end

# -----

# here we define and document the prototypes that are implemented

"""
`AbstractAtoms{T}`

the abstract supertype for storing atomistic configurations. A basic
implementation might simply store a list of positions, see e.g. `Atoms`.

Some Base functions that should be overloaded for atoms objects:
- length
- deleteat
- getindex
- setindex!
"""
abstract type AbstractAtoms{T} end

fltype(at::AbstractAtoms{T}) where {T} = T

abstract type AbstractConstraint end

"""
`AbstractCalculator`: theabstractsupertype of calculators. These
store model information, and are linked to the implementation of energy,
forces, and so forth.
"""
abstract type AbstractCalculator end

fltype(V::AbstractCalculator) = Float64  # a sensible default!

"""
`Dofs{T}`: dof vector; simply an alias for `Vector{T}`
"""
const Dofs{T} = Vector{T}

"""
`positions(at)` : Return copy of positions of all atoms as a `Vector{JVec}`
"""
function positions end

"alias for `positions`"
get_positions = positions

"`set_positions!(at, X) -> at` : Set positions of all atoms"
function set_positions! end

set_positions!(at::AbstractAtoms, p::AbstractMatrix) = set_positions!(at, vecs(p))
set_positions!(at::AbstractAtoms, x, y, z) = set_positions!(at, [x'; y'; z'])

# this is already documented at its first definition in `arrayconversions`?
xyz(at::AbstractAtoms) = xyz(positions(at))

"`momenta(at)` : Return copy of momenta of all atoms as a `Vector{<:JVec}`."
function momenta end

"alias for `momenta`"
get_momenta = momenta

"`set_momenta!(at, P) -> at` : Set momenta of all atoms"
function set_momenta! end

set_momenta!(at::AbstractAtoms, p::AbstractMatrix) = set_momenta!(at, vecs(p))

"`masses` : return vector of all masses"
function masses end

"`get_masses` : alias for `masses`"
get_masses = masses

"`set_masses!(at, M) -> at` : set the atom masses"
function set_masses! end

"`cell(at) -> JMat` : get computational cell (the rows are the lattice vectors)"
function cell end

"alias for `cell`"
get_cell = cell

"`set_cell!(at, C) -> at` : set computational cell; cf. `?cell`"
function set_cell! end

"""
`cell_vecs(at) -> (SVec, SVec, SVec)` : return the three cell vectors
"""
function cell_vecs(at::AbstractAtoms)
   C = cell(at)
   @assert size(C) == (3,3)
   return C[1,:], C[2,:], C[3,:]
end

"`volume(at)` : return volume of computational cell"
volume(at) = abs(det(cell(at)))

"""
`apply_defm!(at, F, t = zero(JVec)) -> at` : for a 3 x 3 matrix `F` and
`at::AbstractAtoms` the affine deformation `x -> Fx + t` is applied
to the configuration, modifying `at` inplace and returning it. This
modifies both the cell and the positions.
"""
function apply_defm!(at::AbstractAtoms{T}, F::AbstractMatrix,
                     t::AbstractVector = zero(JVec{T})) where {T}
   @assert size(F) == (3,3)
   @assert length(t) == 3
   C = cell(at)
   X = positions(at)
   Cnew = C * F'
   for n = 1:length(X)
      # TODO: project back into the cell?
      X[n] = F * X[n] + t
   end
   set_cell!(at, Cnew)
   set_positions!(at, X)
   return at
end

"`set_pbc!(at, p) -> at` : set periodic boundary conditions"
function set_pbc! end
set_pbc!(at::AbstractAtoms, p::Bool) = set_pbc!(at, (p,p,p))

"`pbc(at)` : get array or tuple determining which directions are periodic"
function pbc end

"`get_pbc(at)` : alias for `pbc`"
get_pbc = pbc

"""
`is_cubic(at) -> Bool` : determines whether a cubic cell is used
(i.e. cell is a diagonal matrix)
"""
is_cubic(a::AbstractAtoms) = isdiag(cell(a))

"""
`set_data!(at, name, value)`:
associate some data with `at`; to be stored in a Dict within `at`

if `name::Union{Symbol, AbstractString}`, then `setindex!` is an alias
for `set_data!`, i.e., we may write `at[:id] = val`
"""
function set_data! end
setindex!(at::AbstractAtoms, value,
          name::Union{Symbol, AbstractString}) = set_data!(at, name, value)

"""
`get_data(at, name)`:
obtain some data stored with `set_data!`

if `name::Union{Symbol, AbstractString}`, then `getindex` is an alias
for `get_data`.
"""
function get_data end

getindex(at::AbstractAtoms,
         name::Union{Symbol, AbstractString}) = get_data(at, name)

"`has_data(at, key)` : check whether some data with id `key` is already stored"
function has_data end

"alias for `deleteat!`"
delete_atom! = deleteat!

"`calculator(at)` : return an attached calculator"
function calculator end

"`set_calculator!(at, calc) -> at` : attach a calculator"
function set_calculator! end

"""
`neighbourlist(at, rcut)`

construct a suitable neighbourlist for this atoms object. `rcut` should
be allowed to be either a scalar or a vector.
"""
function neighbourlist end

"""
`static_neighbourlist(at::AbstractAtoms, cutoff; key = :staticnlist)`

This function first checks whether a static neighbourlist already exists
with cutoff `cutoff` and if it does then it returns the existing list.
If it does not, then it computes a new neighbour list with the current
configuration, stores it for later use and returns it.
"""
function static_neighbourlist(at::AbstractAtoms, cutoff; key=:staticnlist)
   recompute = false
   if has_data(at, key)
      nlist = get_data(at, key)
      if cutoff(nlist) != cutoff
         recompute = true
      end
   else
      recompute = true
   end
   if recompute
      set_data!( at, key, neighbourlist(at, cutoff) )
   end
   return get_data(at, key)
end

static_neighbourlist(at::AbstractAtoms) = static_neighbourlist(at, cutoff(at))

#######################################################################
#     CALCULATOR
#######################################################################

"""
`cutoff(calc)` : Returns the cut-off radius of the attached potential.

`cutoff(at)` : returns the cutoff of the attached calculator
"""
function cutoff end

cutoff(at::AbstractAtoms) = cutoff(calculator(at))

"""
`energy(calc, at)`: Return the total potential energy of a configuration of
atoms `at`, using the calculator `calc`. In addition, `energy` can be
called in the following ways:
*  `energy(calc, at)`: base definition; this is normally the only method that
     needs to be overloaded when a new calculator is implemented.
* `energy(at)`: if a calculator is attached to `at`
* `energy(calc, at, const, dof)`
* `energy(calc, at, dof)`
* `energy(at, dofs)`
"""
function energy end
potential_energy = energy

"""
`site_energy(calc, at, n)` : return site energy at atom idx `n`
"""
function site_energy end
site_energy(at::AbstractAtoms, n::Integer) = site_energy(calculator(at), at, n)

"""
`site_energies(calc, at)` : return vector of all site energies
"""
function site_energies end
site_energies(a::AbstractAtoms) = site_energies(calculator(a), a)

"""
`partial_energy(calc, at, subset)` : return energy contained in a `subset` of the
configuration
"""
function partial_energy end
partial_energy(at::AbstractAtoms, subset) =
      partial_energy(calculator(at), at, subset)

"""
`forces(calc, at)` : forces as `Vector{<:JVec}`  (negative gradient w.r.t. atom positions only)
"""
function forces end
forces(at::AbstractAtoms) = forces(calculator(at), at)

"""
`hessian_pos(calc, at):` block-hessian with respect to all atom positions
"""
function hessian_pos end
hessian_pos(at::AbstractAtoms) = hessian_pos(calculator(at), at)

"""
* `virial(c::AbstractCalculator, a::AbstractAtoms) -> JMat`
* `virial(a::AbstractAtoms) -> JMat`

returns virial, (- ∂E / ∂F) where `F = cell(a)'`
"""
function virial end
virial(at::AbstractAtoms) = virial(calculator(at), at)

"""
* `stress(c::AbstractCalculator, a::AbstractAtoms) -> JMatF`
* `stress(a::AbstractAtoms) -> JMatF`

stress = - virial / volume; this function should *not* be overloaded;
instead overload virial
"""
stress(calc::AbstractCalculator, at::AbstractAtoms) =
      - virial(calc, at) / volume(at)

stress(at::AbstractAtoms) = stress(calculator(at), at)

#######################################################################
#  Constraints and DoF Handling
#######################################################################

"""
`position_dofs(at::AbstractAtoms) -> Dofs`
`dofs(at::AbstractAtoms) -> Dofs`

Take an atoms object `at` and return a Dof-vector that fully describes the
state given the potential constraints placed on it.
"""
function position_dofs end

"""
`dofs` is an alias for `position_dofs`
"""
dofs = position_dofs

"""
* `set_position_dofs!(at::AbstractAtoms, cons::AbstractConstraint, x::Dofs) -> at`
* `set_position_dofs!(at::AbstractAtoms, x::Dofs) -> at`
* `set_dofs!(at::AbstractAtoms, cons::AbstractConstraint, x::Dofs) -> at`
* `set_dofs!(at::AbstractAtoms, cons::AbstractConstraint) -> at`

change configuration stored in `at` according to `cons` and `x`.
"""
function set_position_dofs! end

"""
`set_dofs!` is an alias for `set_position_dofs!`
"""
set_dofs! = set_position_dofs!

"""
`momentum_dofs(at::AbstractAtoms) -> Dofs`

Take an atoms object `at` and return a Dof-vector that fully describes the
momenta (given any constraints placed on `at`)
"""
function momentum_dofs end

"""
* `set_momentum_dofs!(at::AbstractAtoms, cons::AbstractConstraint, p::Dofs) -> at`
* `set_momentum_dofs!(at::AbstractAtoms, p::Dofs) -> at`

change configuration stored in `at` according to `cons` and `p`.
"""
function set_momentum_dofs! end

"""
* `project!(at::AbstractAtoms) -> at`

project the `at` onto the constraint manifold.
"""
function project! end

# converting calculator functionality
# The following lines define variants of the `energy` function, which
# get around the issue that energy may or may not include an external
# potential; the rule is as follows:
#  - a new calculator must implement `energy(calc, at)`
#  - a new constraint must implement `energy(calc, at, cons)`
#    and at the end *must* use `energy(calc, at)` but may not call any
#    other variants of `energy`.

"""
* `energy(at, x::Dofs) -> Float64`

`potential_energy`; with potentially added cell terms, e.g., if there is
applied stress or applied pressure.
"""
energy(at::AbstractAtoms, x::Dofs) =
      energy(set_dofs!(at, x))
energy(calc::AbstractCalculator, at::AbstractAtoms, x::Dofs) =
      energy(calc, set_dofs!(at, x))
energy(at::AbstractAtoms) =
      energy(calculator(at), at)

"""
* `gradient(calc, at, cons::AbstractConstraint, c::Dofs) -> Dofs`
* `gradient(at) -> Dofs`
* `gradient(at, x::Dofs) -> Dofs`

gradient of `potential_energy` in dof-format;
depending on the constraint this could include e.g. a gradient w.r.t. cell
shape; normally only the form `gradient(calc, at, cons)` should be
overloaded by an `AbstractConstraints` object.
"""
function gradient end
gradient(at::AbstractAtoms, x::Dofs) =
      gradient(set_dofs!(at, x))
gradient(calc::AbstractCalculator, at::AbstractAtoms, x::Dofs) =
      gradient(calc, set_dofs!(at, x))
gradient(at::AbstractAtoms) =
      gradient(calculator(at), at)

"""
`hessian`: compute hessian of total energy with respect to DOFs!

same as gradient, just returns the hessian; new constraints should only
overload `hessian(calc, at, cons)`
"""
function hessian end
hessian(at::AbstractAtoms, x::Dofs) =
      hessian(set_dofs!(at, x))
hessian(calc::AbstractCalculator, at::AbstractAtoms, x::Dofs) =
      hessian(calc, set_dofs!(at, x))
hessian(at::AbstractAtoms) =
      hessian(calculator(at), at)

#######################################################################
#                     PRECONDITIONER
# we don't create an abstract type, but somewhere we should
# document that a preconditioner must implement the following
# operations:
#   - update!(precond, at)
#   - ldiv!(out, P, f)
#   - mul!)(out, P, x)
#######################################################################

"""
`update!(precond, at)`

Update the preconditioner with the new geometry information.
"""
function update! end
update!(precond, at::AbstractAtoms, x::Dofs) =
            update!(precond, set_dofs!(at, x))

# TODO: make a Julia PR? do this in-place
ldiv!(out::Dofs, P::UniformScaling, x::Dofs) =
      copyto!(out, P \ x)

# mul! is already correctly defined -> no need to define it here.
update!(P::UniformScaling, at::AbstractAtoms) = P

"""
construct a preconditioner suitable for this atoms object

default is `I`
"""
preconditioner(at::AbstractAtoms) =
            preconditioner(at, calculator(at))   # should be preconditioner(calc,  at) ???

preconditioner(at::AbstractAtoms, calc::AbstractCalculator) = I

#######################################################################
## Experimental Prototypes for in-place versions
#######################################################################

"""
`energy!`: non-allocating version of `energy`
"""
energy!(calc::AbstractCalculator, at::AbstractAtoms, temp::Nothing) =
      energy(calc, at)

"""
`energy!`: non-allocating version of `forces`
"""
forces!(F::AbstractVector{JVec{T}}, calc::AbstractCalculator,
        at::AbstractAtoms{T}, temp::Nothing) where {T} =
      copy!(F, forces(calc, at))

"""
`energy!`: non-allocating version of `virial`
"""
virial!(calc::AbstractCalculator, at::AbstractAtoms, temp::Nothing) =
      virial(calc, at)


text:
using Test
using ManyBody
using LinearAlgebra

@testset "Term" begin

@testset "pauli" begin
    @test iszero(zero(Term(1.0)))
    @test iszero(creator(3) * creator(3))
    @test iszero(creator(1) * annihilator(4) * annihilator(4))
end

@testset "convert" begin
    @test Term(2.0) == Term(2)
    @test creator(3) == Term{Int64}(creator(3))
end

@testset "normal ordering" begin
    c = annihilator
    cdag = creator

    t = cdag(7) * cdag(4) * c(2) * c(6)
    @test -t == cdag(4) * cdag(7) * c(2) * c(6)
    @test t == cdag(4) * cdag(7) * c(6) * c(2)

    t2 = c(9) * c(5) * c(3)
    @test t * t2 == +cdag(4) * cdag(7) * c(9) * c(6) * c(5)  * c(3) * c(2)

    # Sign not sure
    @test t2 * t == +cdag(4) * cdag(7) * c(9) * c(6) * c(5)  * c(3) * c(2)
end

@testset "vacancy" begin
    c = annihilator
    cdag = creator
    n = occupation
    v = vacancy

    t1 = 10 * n(2)
    @test t1 * t1 == 10 * t1

    t1 = 4 * v(4)
    @test t1 * t1 == 4 * t1

    t1 = 3 * c(3) * v(2) * cdag(1)
    @test t1 * v(2) == t1
end

@testset "quadraticterm" begin
    t = 2.0 * creator(1) * annihilator(3)
    tprime = 2.0 * creator(3) * annihilator(1)

    @test t' == tprime
    @test transpose(t) == tprime
end

end


text:
@testset "ReificationDifferentFrom2MILP: $(fct_type), type $(T)" for fct_type in ["vector of variables", "vector affine function"], T in [Int, Float64]
    mock = MOIU.MockOptimizer(AbsoluteValueModel{T}())
    model = COIB.ReificationDifferentFrom2MILP{T}(mock)

    if T == Int
        @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    end
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{T},
        MOI.LessThan{T},
    )
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.AbsoluteValue,
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Reification{CP.DifferentFrom{T}},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorOfVariables,
        CP.Reification{CP.DifferentFrom{T}},
    )

    x, _ = MOI.add_constrained_variable(model, MOI.ZeroOne())
    if T == Int
        y, _ = MOI.add_constrained_variable(model, MOI.Integer())
    elseif T == Float64
        y = MOI.add_variable(model)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables([x, y])
    elseif fct_type == "vector affine function"
        MOIU.vectorize([x, y])
    else
        @assert false
    end

    @test_throws AssertionError MOI.add_constraint(model, fct, CP.Reification(CP.DifferentFrom(zero(T))))

    MOI.add_constraint(model, y, MOI.LessThan(5 * one(T)))
    MOI.add_constraint(model, y, MOI.GreaterThan(2 * one(T)))

    c = MOI.add_constraint(model, fct, CP.Reification(CP.DifferentFrom(zero(T))))

    @test MOI.is_valid(model, x)
    @test MOI.is_valid(model, y)
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.Reification{CP.DifferentFrom{T}}}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.Reification{CP.DifferentFrom{T}}) == typeof(bridge)
        if T == Int
            @test MOIB.added_constrained_variable_types(typeof(bridge)) == [(MOI.Integer,)]
        elseif T == Float64
            @test MOIB.added_constrained_variable_types(typeof(bridge)) == Tuple{Type}[]
        else
            @assert false
        end
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.VectorAffineFunction{T}, CP.AbsoluteValue),
            (MOI.VectorAffineFunction{T}, MOI.LessThan{T}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 1
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VariableIndex, MOI.Integer}()) == ((T == Int) ? 1 : 0)
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, CP.AbsoluteValue}()) == 1
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}()) == 2

        @test MOI.get(bridge, MOI.ListOfVariableIndices()) == [bridge.var_abs]
        if T == Int
            @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VariableIndex, MOI.Integer}()) == [bridge.var_abs_int]
        end
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, CP.AbsoluteValue}()) == [bridge.con_abs]
        @test Set(MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}())) == Set([bridge.con_bigm, bridge.con_smallm])
    end

    @testset "Constraint: absolute value" begin
        @test MOI.is_valid(model, bridge.con_abs)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_abs)
        @test length(f.terms) == 2
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_abs) == CP.AbsoluteValue()

        t1 = f.terms[1]
        @test t1.output_index == 1
        @test t1.scalar_term.coefficient === one(T)
        @test t1.scalar_term.variable == bridge.var_abs

        t2 = f.terms[2]
        @test t2.output_index == 2
        @test t2.scalar_term.coefficient === one(T)
        @test t2.scalar_term.variable == y
    end

    @testset "Constraint: big-M" begin
        @test MOI.is_valid(model, bridge.con_bigm)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_bigm)
        @test length(f.terms) == 2
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_bigm) == MOI.LessThan(zero(T))

        t1 = f.terms[1]
        @test t1.coefficient === -5 * one(T)
        @test t1.variable == x

        t2 = f.terms[2]
        @test t2.coefficient === one(T)
        @test t2.variable == bridge.var_abs
    end

    @testset "Constraint: small-M" begin
        @test MOI.is_valid(model, bridge.con_smallm)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_smallm)
        @test length(f.terms) == 2
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_smallm) == MOI.LessThan(zero(T))

        t1 = f.terms[1]
        @test t1.coefficient === -one(T)
        @test t1.variable == x

        smallm = if T == Int
            -1
        elseif T == Float64
            -100_000.0
        else
            @assert false
        end

        t2 = f.terms[2]
        @test t2.coefficient ≈ smallm
        @test t2.variable == bridge.var_abs
    end
end


text:
int_rules_7_5_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 7.5.2*Miscellaneous*inverse*hyperbolic*secant =#
    @apply_utils Antiderivative(asech(~c + ~(d') * ~x), ~x) => ((~c + ~d * ~x) * asech(~c + ~d * ~x)) / ~d + Antiderivative(sqrt(((1 - ~c) - ~d * ~x) / (1 + ~c + ~d * ~x)) / ((1 - ~c) - ~d * ~x), ~x) <-- FreeQ([~c, ~d], ~x)
    @apply_utils Antiderivative(acsch(~c + ~(d') * ~x), ~x) => ((~c + ~d * ~x) * acsch(~c + ~d * ~x)) / ~d + Antiderivative(1 / ((~c + ~d * ~x) * sqrt(1 + 1 / (~c + ~d * ~x) ^ 2)), ~x) <-- FreeQ([~c, ~d], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * asech(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (1 / ~d) * Subst(Antiderivative((~a + ~b * asech(~x)) ^ ~p, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * acsch(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (1 / ~d) * Subst(Antiderivative((~a + ~b * acsch(~x)) ^ ~p, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * asech(~c + ~(d') * ~x)) ^ ~p, ~x) => Unintegrable((~a + ~b * asech(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && Not(IGtQ(~p, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * acsch(~c + ~(d') * ~x)) ^ ~p, ~x) => Unintegrable((~a + ~b * acsch(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && Not(IGtQ(~p, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * asech(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (1 / ~d) * Subst(Antiderivative(((~f * ~x) / ~d) ^ ~m * (~a + ~b * asech(~x)) ^ ~p, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (EqQ(~d * ~e - ~c * ~f, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * acsch(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (1 / ~d) * Subst(Antiderivative(((~f * ~x) / ~d) ^ ~m * (~a + ~b * acsch(~x)) ^ ~p, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (EqQ(~d * ~e - ~c * ~f, 0) && IGtQ(~p, 0))
    #= Antiderivative((~x)^(~m')*asech((~a)+(~b')*(~x)),~x) :=  -((-a)^(m+1)-b^(m+1)*x^(m+1))*asech(a+b*x)/(b^(m+1)*(m+1)) + 1/(b^(m+1)*(m+1))*Subst(Antiderivative(((-a*x)^(m+1)-(1-a*x)^(m+1))/(x^(m+1)* sqrt(-1+x)*sqrt(1+x)),x),x,1/(a+b*x)) <-- FreeQ([a,b],x) &&  IntegerQ(m) && NeQ(m,-1) =#
    #= Antiderivative((~x)^(~m')*acsch((~a)+(~b')*(~x)),~x) :=  -((-a)^(m+1)-b^(m+1)*x^(m+1))*acsch(a+b*x)/(b^(m+1)*(m+1)) + 1/(b^(m+1)*(m+1))*Subst(Antiderivative(((-a*x)^(m+1)-(1-a*x)^(m+1))/(x^(m+1)* sqrt(1+x^2)),x),x,1/(a+b*x)) <-- FreeQ([a,b],x) && IntegerQ(m) &&  NeQ(m,-1) =#
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * asech(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (-1 / (~d) ^ (~m + 1)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~p * sech(~x) * tanh(~x) * ((~d * ~e - ~c * ~f) + ~f * sech(~x)) ^ ~m, ~x), ~x, asech(~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 0) && IntegerQ(~m))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * acsch(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (-1 / (~d) ^ (~m + 1)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~p * csch(~x) * coth(~x) * ((~d * ~e - ~c * ~f) + ~f * csch(~x)) ^ ~m, ~x), ~x, acsch(~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 0) && IntegerQ(~m))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * asech(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (1 / ~d) * Subst(Antiderivative(((~d * ~e - ~c * ~f) / ~d + (~f * ~x) / ~d) ^ ~m * (~a + ~b * asech(~x)) ^ ~p, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * acsch(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (1 / ~d) * Subst(Antiderivative(((~d * ~e - ~c * ~f) / ~d + (~f * ~x) / ~d) ^ ~m * (~a + ~b * acsch(~x)) ^ ~p, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * asech(~c + ~(d') * ~x)) ^ ~p, ~x) => Unintegrable((~e + ~f * ~x) ^ ~m * (~a + ~b * asech(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~p], ~x) && Not(IGtQ(~p, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * acsch(~c + ~(d') * ~x)) ^ ~p, ~x) => Unintegrable((~e + ~f * ~x) ^ ~m * (~a + ~b * acsch(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~p], ~x) && Not(IGtQ(~p, 0))
    @apply_utils Antiderivative(~(u') * asech(~(c') / (~(a') + ~(b') * (~x) ^ ~(n'))) ^ ~(m'), ~x) => Antiderivative(~u * acosh(~a / ~c + (~b * (~x) ^ ~n) / ~c) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~m], ~x)
    @apply_utils Antiderivative(~(u') * acsch(~(c') / (~(a') + ~(b') * (~x) ^ ~(n'))) ^ ~(m'), ~x) => Antiderivative(~u * asinh(~a / ~c + (~b * (~x) ^ ~n) / ~c) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~m], ~x)
    @apply_utils Antiderivative(ℯ ^ asech(~(a') * ~x), ~x) => ~x * ℯ ^ asech(~a * ~x) + log(~x) / ~a + (1 / ~a) * Antiderivative((1 / (~x * (1 - ~a * ~x))) * sqrt((1 - ~a * ~x) / (1 + ~a * ~x)), ~x) <-- FreeQ(~a, ~x)
    @apply_utils Antiderivative(ℯ ^ asech(~(a') * (~x) ^ ~p), ~x) => ~x * ℯ ^ asech(~a * (~x) ^ ~p) + (~p / ~a) * Antiderivative(1 / (~x) ^ ~p, ~x) + ((~p * sqrt(1 + ~a * (~x) ^ ~p)) / ~a) * sqrt(1 / (1 + ~a * (~x) ^ ~p)) * Antiderivative(1 / ((~x) ^ ~p * sqrt(1 + ~a * (~x) ^ ~p) * sqrt(1 - ~a * (~x) ^ ~p)), ~x) <-- FreeQ([~a, ~p], ~x)
    @apply_utils Antiderivative(ℯ ^ acsch(~(a') * (~x) ^ ~(p')), ~x) => (1 / ~a) * Antiderivative(1 / (~x) ^ ~p, ~x) + Antiderivative(sqrt(1 + 1 / ((~a) ^ 2 * (~x) ^ (2 * ~p))), ~x) <-- FreeQ([~a, ~p], ~x)
    @apply_utils Antiderivative(ℯ ^ (~(n') * asech(~u)), ~x) => Antiderivative((1 / ~u + sqrt((1 - ~u) / (1 + ~u)) + (1 / ~u) * sqrt((1 - ~u) / (1 + ~u))) ^ ~n, ~x) <-- IntegerQ(~n)
    @apply_utils Antiderivative(ℯ ^ (~(n') * acsch(~u)), ~x) => Antiderivative((1 / ~u + sqrt(1 + 1 / (~u) ^ 2)) ^ ~n, ~x) <-- IntegerQ(~n)
    @apply_utils Antiderivative(ℯ ^ asech(~(a') * (~x) ^ ~(p')) / ~x, ~x) => -1 / (~a * ~p * (~x) ^ ~p) + (sqrt(1 + ~a * (~x) ^ ~p) / ~a) * sqrt(1 / (1 + ~a * (~x) ^ ~p)) * Antiderivative((sqrt(1 + ~a * (~x) ^ ~p) * sqrt(1 - ~a * (~x) ^ ~p)) / (~x) ^ (~p + 1), ~x) <-- FreeQ([~a, ~p], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * ℯ ^ asech(~(a') * (~x) ^ ~(p')), ~x) => ((~x) ^ (~m + 1) * ℯ ^ asech(~a * (~x) ^ ~p)) / (~m + 1) + (~p / (~a * (~m + 1))) * Antiderivative((~x) ^ (~m - ~p), ~x) + ((~p * sqrt(1 + ~a * (~x) ^ ~p)) / (~a * (~m + 1))) * sqrt(1 / (1 + ~a * (~x) ^ ~p)) * Antiderivative((~x) ^ (~m - ~p) / (sqrt(1 + ~a * (~x) ^ ~p) * sqrt(1 - ~a * (~x) ^ ~p)), ~x) <-- FreeQ([~a, ~m, ~p], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~x) ^ ~(m') * ℯ ^ acsch(~(a') * (~x) ^ ~(p')), ~x) => (1 / ~a) * Antiderivative((~x) ^ (~m - ~p), ~x) + Antiderivative((~x) ^ ~m * sqrt(1 + 1 / ((~a) ^ 2 * (~x) ^ (2 * ~p))), ~x) <-- FreeQ([~a, ~m, ~p], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * ℯ ^ (~(n') * asech(~u)), ~x) => Antiderivative((~x) ^ ~m * (1 / ~u + sqrt((1 - ~u) / (1 + ~u)) + (1 / ~u) * sqrt((1 - ~u) / (1 + ~u))) ^ ~n, ~x) <-- FreeQ(~m, ~x) && IntegerQ(~n)
    @apply_utils Antiderivative((~x) ^ ~(m') * ℯ ^ (~(n') * acsch(~u)), ~x) => Antiderivative((~x) ^ ~m * (1 / ~u + sqrt(1 + 1 / (~u) ^ 2)) ^ ~n, ~x) <-- FreeQ(~m, ~x) && IntegerQ(~n)
    @apply_utils Antiderivative(ℯ ^ asech(~(c') * ~x) / (~a + ~(b') * (~x) ^ 2), ~x) => (1 / (~a * ~c)) * Antiderivative(sqrt(1 / (1 + ~c * ~x)) / (~x * sqrt(1 - ~c * ~x)), ~x) + (1 / ~c) * Antiderivative(1 / (~x * (~a + ~b * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && EqQ(~b + ~a * (~c) ^ 2, 0)
    @apply_utils Antiderivative(ℯ ^ acsch(~(c') * ~x) / (~a + ~(b') * (~x) ^ 2), ~x) => (1 / (~a * (~c) ^ 2)) * Antiderivative(1 / ((~x) ^ 2 * sqrt(1 + 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x) + (1 / ~c) * Antiderivative(1 / (~x * (~a + ~b * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && EqQ(~b - ~a * (~c) ^ 2, 0)
    @apply_utils Antiderivative(((~(d') * ~x) ^ ~(m') * ℯ ^ asech(~(c') * ~x)) / (~a + ~(b') * (~x) ^ 2), ~x) => (~d / (~a * ~c)) * Antiderivative(((~d * ~x) ^ (~m - 1) * sqrt(1 / (1 + ~c * ~x))) / sqrt(1 - ~c * ~x), ~x) + (~d / ~c) * Antiderivative((~d * ~x) ^ (~m - 1) / (~a + ~b * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && EqQ(~b + ~a * (~c) ^ 2, 0)
    @apply_utils Antiderivative(((~(d') * ~x) ^ ~(m') * ℯ ^ acsch(~(c') * ~x)) / (~a + ~(b') * (~x) ^ 2), ~x) => ((~d) ^ 2 / (~a * (~c) ^ 2)) * Antiderivative((~d * ~x) ^ (~m - 2) / sqrt(1 + 1 / ((~c) ^ 2 * (~x) ^ 2)), ~x) + (~d / ~c) * Antiderivative((~d * ~x) ^ (~m - 1) / (~a + ~b * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && EqQ(~b - ~a * (~c) ^ 2, 0)
    @apply_utils Antiderivative(asech(~u), ~x) => ~x * asech(~u) + (sqrt(1 - (~u) ^ 2) / (~u * sqrt(-1 + 1 / ~u) * sqrt(1 + 1 / ~u))) * Antiderivative(SimplifyIntegrand((~x * D(~u, ~x)) / (~u * sqrt(1 - (~u) ^ 2)), ~x), ~x) <-- InverseFunctionFreeQ(~u, ~x) && Not(FunctionOfExponentialQ(~u, ~x))
    @apply_utils Antiderivative(acsch(~u), ~x) => ~x * acsch(~u) - (~u / sqrt(-((~u) ^ 2))) * Antiderivative(SimplifyIntegrand((~x * D(~u, ~x)) / (~u * sqrt(-1 - (~u) ^ 2)), ~x), ~x) <-- InverseFunctionFreeQ(~u, ~x) && Not(FunctionOfExponentialQ(~u, ~x))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * asech(~u)), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * (~a + ~b * asech(~u))) / (~d * (~m + 1)) + ((~b * sqrt(1 - (~u) ^ 2)) / (~d * (~m + 1) * ~u * sqrt(-1 + 1 / ~u) * sqrt(1 + 1 / ~u))) * Antiderivative(SimplifyIntegrand(((~c + ~d * ~x) ^ (~m + 1) * D(~u, ~x)) / (~u * sqrt(1 - (~u) ^ 2)), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (NeQ(~m, -1) && (InverseFunctionFreeQ(~u, ~x) && (Not(FunctionOfQ((~c + ~d * ~x) ^ (~m + 1), ~u, ~x)) && Not(FunctionOfExponentialQ(~u, ~x)))))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * acsch(~u)), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * (~a + ~b * acsch(~u))) / (~d * (~m + 1)) - ((~b * ~u) / (~d * (~m + 1) * sqrt(-((~u) ^ 2)))) * Antiderivative(SimplifyIntegrand(((~c + ~d * ~x) ^ (~m + 1) * D(~u, ~x)) / (~u * sqrt(-1 - (~u) ^ 2)), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (NeQ(~m, -1) && (InverseFunctionFreeQ(~u, ~x) && (Not(FunctionOfQ((~c + ~d * ~x) ^ (~m + 1), ~u, ~x)) && Not(FunctionOfExponentialQ(~u, ~x)))))
    @apply_utils Antiderivative(~v * (~(a') + ~(b') * asech(~u)), ~x) => With([w = IntHide(~v, ~x)], Dist(~a + ~b * asech(~u), w, ~x) + ((~b * sqrt(1 - (~u) ^ 2)) / (~u * sqrt(-1 + 1 / ~u) * sqrt(1 + 1 / ~u))) * Antiderivative(SimplifyIntegrand((w * D(~u, ~x)) / (~u * sqrt(1 - (~u) ^ 2)), ~x), ~x) <-- InverseFunctionFreeQ(w, ~x)) <-- FreeQ([~a, ~b], ~x) && (InverseFunctionFreeQ(~u, ~x) && Not(MatchQ(~v, (~(c') + ~(d') * ~x) ^ ~(m') <-- FreeQ([c, d, m], ~x))))
    @apply_utils Antiderivative(~v * (~(a') + ~(b') * acsch(~u)), ~x) => With([w = IntHide(~v, ~x)], Dist(~a + ~b * acsch(~u), w, ~x) - ((~b * ~u) / sqrt(-((~u) ^ 2))) * Antiderivative(SimplifyIntegrand((w * D(~u, ~x)) / (~u * sqrt(-1 - (~u) ^ 2)), ~x), ~x) <-- InverseFunctionFreeQ(w, ~x)) <-- FreeQ([~a, ~b], ~x) && (InverseFunctionFreeQ(~u, ~x) && Not(MatchQ(~v, (~(c') + ~(d') * ~x) ^ ~(m') <-- FreeQ([c, d, m], ~x))))
end

max_stars_repo_path: src/rules/7 Inverse hyperbolic functions/7.5 Inverse hyperbolic secant/7.5.2 Miscellaneous inverse hyperbolic secant.jl


text:
using Pkg; for p in ("Knet","ArgParse"); haskey(Pkg.installed(),p) || Pkg.add(p); end

"""

This example learns to classify hand-written digits from the
[MNIST](http://yann.lecun.com/exdb/mnist) dataset.  There are 60000
training and 10000 test examples. Each input x consists of 784 pixels
representing a 28x28 image.  The pixel values are normalized to
[0,1]. Each output y is converted to a ten-dimensional one-hot vector
(a vector that has a single non-zero component) indicating the correct
class (0-9) for a given image.  10 is used to represent 0.

You can run the demo using `julia mlp.jl` on the command line or
`julia> MLP.main()` at the Julia prompt.  Options can be used like
`julia mlp.jl --epochs 3` or `julia> MLP.main("--epochs 3")`.  Use
`julia mlp.jl --help` for a list of options.  The dataset will be
automatically downloaded.  By default a softmax model will be trained
for 10 epochs.  You can also train a multi-layer perceptron by
specifying one or more --hidden sizes.  The accuracy for the training
and test sets will be printed at every epoch and optimized parameters
will be returned.

"""
module MLP
using Knet,ArgParse
include(Knet.dir("data","mnist.jl"))

function predict(w,x)
    for i=1:2:length(w)
        x = w[i]*mat(x) .+ w[i+1]
        if i<length(w)-1
            x = relu.(x) # max(0,x)
        end
    end
    return x
end

loss(w,x,ygold) = nll(predict(w,x),ygold)

lossgradient = grad(loss)

function train(w, dtrn; lr=.5, epochs=10)
    for epoch=1:epochs
        for (x,y) in dtrn
            g = lossgradient(w, x, y)
            update!(w,g;lr=lr)
        end
    end
    return w
end

function weights(h...; atype=Array{Float32}, winit=0.1)
    w = Any[]
    x = 28*28
    for y in [h..., 10]
        push!(w, convert(atype, winit*randn(y,x)))
        push!(w, convert(atype, zeros(y, 1)))
        x = y
    end
    return w
end

function main(args="")
    s = ArgParseSettings()
    s.description="mlp.jl (c) Deniz Yuret, 2016. Multi-layer perceptron model on the MNIST handwritten digit recognition problem from http://yann.lecun.com/exdb/mnist."
    s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        ("--seed"; arg_type=Int; default=-1; help="random number seed: use a nonnegative int for repeatable results")
        ("--batchsize"; arg_type=Int; default=100; help="minibatch size")
        ("--epochs"; arg_type=Int; default=10; help="number of epochs for training")
        ("--hidden"; nargs='*'; arg_type=Int; help="sizes of hidden layers, e.g. --hidden 128 64 for a net with two hidden layers")
        ("--lr"; arg_type=Float64; default=0.5; help="learning rate")
        ("--winit"; arg_type=Float64; default=0.1; help="w initialized with winit*randn()")
        ("--fast"; action=:store_true; help="skip loss printing for faster run")
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}" : "Array{Float32}"); help="array type: Array for cpu, KnetArray for gpu")
        ("--gcheck"; arg_type=Int; default=0; help="check N random gradients per parameter")
        # These are to experiment with sparse arrays
        # ("--xtype"; help="input array type: defaults to atype")
        # ("--ytype"; help="output array type: defaults to atype")
    end
    isa(args, AbstractString) && (args=split(args))
    if in("--help", args) || in("-h", args)
        ArgParse.show_help(s; exit_when_done=false)
        return
    end
    o = parse_args(args, s; as_symbols=true)
    if !o[:fast]
        println(s.description)
        println("opts=",[(k,v) for (k,v) in o]...)
    end
    o[:seed] > 0 && Knet.seed!(o[:seed])
    atype = eval(Meta.parse(o[:atype]))
    w = weights(o[:hidden]...; atype=atype, winit=o[:winit])
    xtrn,ytrn,xtst,ytst = mnist()
    global dtrn = minibatch(xtrn, ytrn, o[:batchsize]; xtype=atype)
    global dtst = minibatch(xtst, ytst, o[:batchsize]; xtype=atype)
    report(epoch)=println((:epoch,epoch,:trn,accuracy(w,dtrn,predict),:tst,accuracy(w,dtst,predict)))
    if o[:fast]
        (train(w, dtrn; lr=o[:lr], epochs=o[:epochs]); gpu()>=0 && Knet.cudaDeviceSynchronize())
    else
        report(0)
        @time for epoch=1:o[:epochs]
            train(w, dtrn; lr=o[:lr], epochs=1)
            report(epoch)
            if o[:gcheck] > 0
                gradcheck(loss, w, first(dtrn)...; gcheck=o[:gcheck], verbose=true)
            end
        end
    end
    return w
end

# This allows both non-interactive (shell command) and interactive calls like:
# $ julia mlp.jl --epochs 10
# julia> MLP.main("--epochs 10")
PROGRAM_FILE == "mlp.jl" && main(ARGS)

end # module

# SAMPLE RUN 65f57ff+ Wed Sep 14 10:02:30 EEST 2016
#
# mnist2d.jl (c) Deniz Yuret, 2016. Multi-layer perceptron model on the MNIST handwritten digit recognition problem from http://yann.lecun.com/exdb/mnist.
# opts=(:seed,-1)(:batchsize,100)(:hidden,Int64[])(:epochs,10)(:lr,0.5)(:atype,"KnetArray{Float32}")(:gcheck,0)(:winit,0.1)(:fast,true)
# (:epoch,0,:trn,0.079066664f0,:tst,0.0842f0)
#   2.168927 seconds (2.95 M allocations: 115.993 MB, 1.84% gc time)
# (:epoch,10,:trn,0.9195333f0,:tst,0.9158f0)


text:
"""
    QNumber

Abstract type representing any expression involving operators.
"""
abstract type QNumber end

"""
    QSym <: QNumber

Abstract type representing fundamental operator types.
"""
abstract type QSym <: QNumber end

# Generic hash fallback for interface -- this will be slow
function Base.hash(op::T, h::UInt) where T<:QSym
    n = fieldcount(T)
    if n == 3
        # These three fields need to be defined for any QSym
        return hash(T, hash(op.hilbert, hash(op.name, hash(op.aon, h))))
    else
        # If there are more we'll need to iterate through
        h_ = copy(h)
        for k = n:-1:4
            if fieldname(typeof(op), k) !== :metadata
                h_ = hash(getfield(op, k), h_)
            end
        end
        return hash(T, hash(op.hilbert, hash(op.name, hash(op.aon, h_))))
    end
end

"""
    QTerm <: QNumber

Abstract type representing noncommutative expressions.
"""
abstract type QTerm <: QNumber end

Base.isless(a::QSym, b::QSym) = a.name < b.name

## Interface for SymbolicUtils

TermInterface.exprhead(::QNumber) = :call
TermInterface.istree(::QSym) = false
TermInterface.istree(::QTerm) = true
TermInterface.istree(::Type{T}) where {T<:QTerm} = true

# Symbolic type promotion
SymbolicUtils.promote_symtype(f, Ts::Type{<:QNumber}...) = promote_type(Ts...)
SymbolicUtils.promote_symtype(f, T::Type{<:QNumber}, Ts...) = T
SymbolicUtils.promote_symtype(f,T::Type{<:QNumber},S::Type{<:Number}) = T
SymbolicUtils.promote_symtype(f,T::Type{<:Number},S::Type{<:QNumber}) = S
SymbolicUtils.promote_symtype(f,T::Type{<:QNumber},S::Type{<:QNumber}) = promote_type(T,S)

SymbolicUtils.symtype(x::T) where T<:QNumber = T

# Standard simplify
function SymbolicUtils.simplify(x::QNumber;kwargs...)
    avg = average(x)
    avg_ = SymbolicUtils.simplify(avg;kwargs...)
    return undo_average(avg_)
end

## End of interface

## Methods
import Base: *, +, -
const SNuN = Union{<:SymbolicUtils.Symbolic{<:Number}, <:Number}

Base.:~(a::QNumber, b::QNumber) = Symbolics.Equation(a, b)

## Multiplication
"""
    QMul <: QTerm

Represent a multiplication involving [`QSym`](@ref) types.

Fields:
======

* arg_c: The commutative prefactor.
* args_nc: A vector containing all [`QSym`](@ref) types.
"""
struct QMul{M} <: QTerm
    arg_c
    args_nc::Vector{Any}
    metadata::M
    function QMul{M}(arg_c, args_nc, metadata) where {M}
        if SymbolicUtils._isone(arg_c) && length(args_nc)==1
            return args_nc[1]
        else
            return new(arg_c, args_nc, metadata)
        end
    end
end
QMul(arg_c, args_nc; metadata::M=NO_METADATA) where {M} = QMul{M}(arg_c, args_nc, metadata)
Base.hash(q::QMul, h::UInt) = hash(QMul, hash(q.arg_c, SymbolicUtils.hashvec(q.args_nc, h)))
Base.isless(a::QMul, b::QMul) = isless(a.h, b.h)

SymbolicUtils.operation(::QMul) = (*)
SymbolicUtils.arguments(a::QMul) = vcat(a.arg_c, a.args_nc)

function SymbolicUtils.similarterm(::QMul, ::typeof(*), args; metadata=NO_METADATA, exprhead=nothing)
    args_c = filter(x->!(x isa QNumber), args)
    args_nc = filter(x->x isa QNumber, args)
    arg_c = *(args_c...)
    return QMul(arg_c, args_nc; metadata)
end

SymbolicUtils.metadata(a::QMul) = a.metadata

function Base.adjoint(q::QMul)
    args_nc = map(adjoint, q.args_nc)
    reverse!(args_nc)
    sort!(args_nc, by=acts_on)
    return QMul(conj(q.arg_c), args_nc; q.metadata)
end

function Base.isequal(a::QMul, b::QMul)
    isequal(a.arg_c, b.arg_c) || return false
    length(a.args_nc)==length(b.args_nc) || return false
    for (arg_a, arg_b) ∈ zip(a.args_nc, b.args_nc)
        isequal(arg_a,arg_b) || return false
    end
    return true
end

function *(a::QSym,b::QSym)
    check_hilbert(a, b)
    args = [a,b]
    sort!(args, by=acts_on)
    QMul(1,args)
end

function *(a::QSym, b::SNuN)
    SymbolicUtils._iszero(b) && return b
    SymbolicUtils._isone(b) && return a
    return QMul(b,[a])
end
*(b::SNuN, a::QNumber) = a*b

function *(a::QMul, b::SNuN)
    SymbolicUtils._iszero(b) && return b
    SymbolicUtils._isone(b) && return a
    arg_c = a.arg_c * b
    return QMul(arg_c, a.args_nc)
end

function *(a::QSym, b::QMul)
    check_hilbert(a, b)
    args_nc = vcat(a,b.args_nc)
    sort!(args_nc, by=acts_on)
    return merge_commutators(b.arg_c,args_nc)
end
function *(a::QMul, b::QSym)
    check_hilbert(a, b)
    args_nc = vcat(a.args_nc, b)
    sort!(args_nc, by=acts_on)
    return merge_commutators(a.arg_c,args_nc)
end

function *(a::QMul, b::QMul)
    check_hilbert(a, b)
    args_nc = vcat(a.args_nc, b.args_nc)
    sort!(args_nc, by=acts_on)
    arg_c = a.arg_c*b.arg_c
    return merge_commutators(arg_c,args_nc)
end

Base.:/(a::QNumber, b::SNuN) = (1/b) * a

function merge_commutators(arg_c,args_nc)
    i = 1
    was_merged = false
    while i<length(args_nc)
        if _ismergeable(args_nc[i], args_nc[i+1])
            args_nc[i] = *(args_nc[i], args_nc[i+1])
            iszero(args_nc[i]) && return 0
            deleteat!(args_nc, i+1)
            was_merged = true
        end
        i += 1
    end
    if was_merged
        return *(arg_c, args_nc...)
    else
        return QMul(arg_c, args_nc)
    end
end

_ismergeable(a,b) = isequal(acts_on(a),acts_on(b)) && ismergeable(a,b) && isequal(hilbert(a),hilbert(b))
ismergeable(a,b) = false

## Powers
function Base.:^(a::QNumber, n::Integer)
    iszero(n) && return 1
    isone(n) && return a
    return *((a for i=1:n)...)
end

## Addition
"""
    QAdd <: QTerm

Represent an addition involving [`QNumber`](@ref) and other types.
"""
struct QAdd <: QTerm
    arguments::Vector{Any}
end

Base.hash(q::T, h::UInt) where T<:QAdd = hash(T, SymbolicUtils.hashvec(q.arguments, h))
function Base.isequal(a::QAdd,b::QAdd)
    length(a.arguments)==length(b.arguments) || return false
    for (arg_a,arg_b) ∈ zip(a.arguments, b.arguments)
        isequal(arg_a, arg_b) || return false
    end
    return true
end

SymbolicUtils.operation(::QAdd) = (+)
SymbolicUtils.arguments(a::QAdd) = a.arguments
SymbolicUtils.similarterm(::QAdd, ::typeof(+), args; metadata=NO_METADATA, exprhead=nothing) = QAdd(args; metadata)

SymbolicUtils.metadata(q::QAdd) = q.metadata

Base.adjoint(q::QAdd) = QAdd(map(adjoint, q.arguments))

-(a::QNumber) = -1*a
-(a,b::QNumber) = a + (-b)
-(a::QNumber,b) = a + (-b)
-(a::QNumber,b::QNumber) = a + (-b)

+(a::QNumber, b::SNuN) = QAdd([a,b])
+(a::SNuN,b::QNumber) = +(b,a)
function +(a::QAdd,b::SNuN)
    SymbolicUtils._iszero(b) && return a
    args = vcat(a.arguments, b)
    return QAdd(args)
end

function +(a::QNumber, b::QNumber)
    check_hilbert(a, b)
    args = [a,b]
    return QAdd(args)
end

function +(a::QAdd,b::QNumber)
    check_hilbert(a, b)
    args = vcat(a.arguments, b)
    return QAdd(args)
end
function +(b::QNumber,a::QAdd)
    check_hilbert(a, b)
    args = vcat(a.arguments, b)
    return QAdd(args)
end
function +(a::QAdd,b::QAdd)
    check_hilbert(a, b)
    args = vcat(a.arguments, b.arguments)
    return QAdd(args)
end

function *(a::QAdd, b)
    check_hilbert(a, b)
    args = Any[a_ * b for a_ ∈ a.arguments]
    flatten_adds!(args)
    q = QAdd(args)
    return q
end
function *(a::QNumber, b::QAdd)
    check_hilbert(a, b)
    args = Any[a * b_ for b_ ∈ b.arguments]
    flatten_adds!(args)
    q = QAdd(args)
    return q
end

function *(a::QAdd, b::QAdd)
    check_hilbert(a, b)
    args = []
    for a_ ∈ a.arguments, b_ ∈ b.arguments
        push!(args, a_ * b_)
    end
    flatten_adds!(args)
    q = QAdd(args)
    return q
end

function flatten_adds!(args)
    i = 1
    while i <= length(args)
        if args[i] isa QAdd
            append!(args,args[i].arguments)
            deleteat!(args, i)
        end
        i += 1
    end
    return args
end

## Hilbert space checks
check_hilbert(a::QNumber,b::QNumber) = (hilbert(a) == hilbert(b)) || error("Incompatible Hilbert spaces $(a.hilbert) and $(b.hilbert)!")
check_hilbert(x,y) = nothing

hilbert(a::QSym) = a.hilbert
hilbert(a::QMul) = hilbert(a.args_nc[1])
function hilbert(a::QAdd)
    idx = findfirst(x->x isa QNumber, a.arguments)
    hilbert(a.arguments[idx])
end

const AonType = Union{<:Int,<:ClusterAon}
"""
    acts_on(op)

Shows on which Hilbert space `op` acts. For [`QSym`](@ref) types, this
returns an Integer, whereas for a `Term` it returns a `Vector{Int}`
whose entries specify all subspaces on which the expression acts.
"""
acts_on(op::QSym) = op.aon
function acts_on(q::QMul)
    aon = AonType[]
    for arg ∈ q.args_nc
        aon_ = acts_on(arg)
        aon_ ∈ aon || push!(aon, aon_)
    end
    return aon
end
function acts_on(q::QAdd)
    aon = AonType[]
    for arg ∈ q.arguments
        append!(aon, acts_on(arg))
    end
    unique!(aon)
    sort!(aon)
    return aon
end
acts_on(x) = Int[]

Base.one(::T) where T<:QNumber = one(T)
Base.one(::Type{<:QNumber}) = 1
Base.isone(::QNumber) = false
Base.zero(::T) where T<:QNumber = zero(T)
Base.zero(::Type{<:QNumber}) = 0
Base.iszero(::QNumber) = false

"""
    @qnumbers

Convenience macro for the construction of operators.

Examples
========
```
julia> h = FockSpace(:fock)
ℋ(fock)

julia> @qnumbers a::Destroy(h)
(a,)

julia> h = FockSpace(:one) ⊗ FockSpace(:two)
ℋ(one) ⊗ ℋ(two)

julia> @qnumbers b::Destroy(h,2)
(b,)
```
"""
macro qnumbers(qs...)
    ex = Expr(:block)
    qnames = []
    for q in qs
        @assert q isa Expr && q.head==:(::)
        q_ = q.args[1]
        @assert q_ isa Symbol
        push!(qnames, q_)
        f = q.args[2]
        @assert f isa Expr && f.head==:call
        op = _make_operator(q_, f.args...)
        ex_ = Expr(:(=), esc(q_), op)
        push!(ex.args, ex_)
    end
    push!(ex.args, Expr(:tuple, map(esc, qnames)...))
    return ex
end

function _make_operator(name, T, h, args...)
    name_ = Expr(:quote, name)
    d = source_metadata(:qnumbers, name)
    return Expr(:call, T, esc(h), name_, args..., Expr(:kw, :metadata, Expr(:quote, d)))
end


text:
import Distances
using Parameters

# TODO: very inefficient for sparse connections,
# upgrade to SparseArrays when performance becomes limiting
mutable struct Contraption
    # Columns are points, row 1 is x-coordinates, row 2 is y-coordinates
    position::Array{Float64, 2}
    velocity::Array{Float64, 2}

    # Mass of point i
    mass::Array{Float64, 1}

    # Zero value implies lack of spring
    # Spring strength between point i and j
    springs::Array{Float64, 2}

    # Length of springs at which no force is applied
    rest_length::Array{Float64, 2}
end

function Contraption(position::Array{Float64, 2},
                     velocity::Array{Float64, 2},
                     mass::Array{Float64, 1},
                     springs::Array{Float64, 2})

    # Basic shape checking
    n = length(mass)

    @assert size(position) ==  (2, n)

    @assert size(velocity) == (2, n)

    @assert size(mass) == (n,)

    @assert size(springs) == (n, n)

    # Spring field checking
    for i in 1:n, j in (i+1):n
        # Symmetry checking
        @assert springs[i, j] == springs[j, i]

        # Positivity checking
        @assert springs[i, j] >= 0.
    end

    # Self-connection checking
    for i in 1:n
        @assert springs[i, i] == 0.
    end

    # Initialization
    rest_length = zeros(n, n)
    for i in 1:n, j in (i+1):n
        if springs[i, j] != 0.
            initial_length = Distances.norm(position[:, i] - position[:, j])
            @assert initial_length > 0.

            rest_length[i, j], rest_length[j, i] = initial_length, initial_length
        end
    end
    @assert size(rest_length) == (n, n)

    return Contraption(position, velocity, mass,
                       springs, rest_length)

end

struct Bounds
    x::Tuple{Float64, Float64}
    y::Tuple{Float64, Float64}
    width::Float64
    height::Float64

    function Bounds(x, y)
        @assert x[1] < x[2]
        @assert y[1] < y[2]

        width = x[2] - x[1]
        height = y[2] - y[1]

        return new(x, y, width, height)
    end
end

@with_kw struct PhysicsSettings
    g::Float64  = 10.
    drag::Float64 = 0.1
    elasticity::Float64 = 0.4
    bounds::Bounds = Bounds((-5., 5.),
                            (-5., 5.))

    # Check that constants are sane
    function PhysicsSettings(g, drag, elasticity, bounds)
        if !(drag >= 0)
            error("drag must be at least 0")
        end

        if !(0 <= elasticity <= 1)
            error("elasticity must be between 0 and 1")
        end

        return new(g, drag, elasticity, bounds)
    end
end


text:
# Helper functions used in test routines

"Test helper functions"
module testhelpers

"Closure to check if something is between `a` and `b`"
function between(a, b)
    @assert a < b
    x -> a < x < b
end

"Placeholder function so that we can test that something doesn't error"
noerror(x) = true

end # module testhelpers


text:
passes = open("boarding_passes.txt") do f
    readlines(f)
end

seatids = sort(map(bits -> parse(Int, bits[1:7], base = 2) * 8 + parse(Int, bits[8:10], base = 2), map(pass -> replace(replace(pass, r"B|R" => "1"), r"F|L" => "0"), passes)))

println(maximum(seatids))

println(argmin(map(i -> seatids[i+1]-seatids[i] == 1, 1:(length(seatids)-1))) + minimum(seatids))


text:
using Test
using Pkg
using Random

const jot_path = abspath(joinpath(pwd(), ".."))
Pkg.develop(PackageSpec(path=jot_path))
using Jot

Pkg.develop(PackageSpec(path=joinpath(jot_path, "test", "JotTest1")))
using JotTest1

const aws_config = AWSConfig(account_id="513118378795", region="ap-northeast-1")
const test_suffix = randstring("abcdefghijklmnopqrstuvwxyz1234567890", 12)

function reset_response_suffix(test_path::String)::String
  response_suffix = randstring(12)
  open(joinpath(jot_path, test_path), "w") do rsfile
    write(rsfile, response_suffix)
  end
  response_suffix
end

function get_response_suffix(test_path::String)::String
  open(joinpath(jot_path, test_path), "r") do rsfile
    readchomp(rsfile)
  end
end

struct ExpectedLabels
  RESPONDER_PACKAGE_NAME::Union{Nothing, String}
  RESPONDER_FUNCTION_NAME::String
  RESPONDER_PKG_SOURCE::String
  user_defined_labels::Dict{String, String}
end

function test_actual_labels_against_expected(
    actual::Jot.Labels,
    expected::ExpectedLabels,
  )::Bool
  @info [getfield(actual, fn) for fn in fieldnames(ExpectedLabels) if !isnothing(getfield(expected, fn))]
  @info [getfield(expected, fn) for fn in fieldnames(ExpectedLabels) if !isnothing(getfield(expected, fn))]
  all([getfield(actual, fn) == getfield(expected, fn) for fn in fieldnames(ExpectedLabels)])
end

function run_tests(;
    clean_up::Bool=true,
    example_simple::Bool=false,
    example_components::Bool=false,
    quartet::Bool=false,
    quartet_partial::Bool=true,
    quartet_to::AbstractString="lambda_function",
  )
  ENV["JOT_TEST_RUNNING"] = "true"
  if all([example_simple, example_components, quartet] .== false)
    example_simple = true; example_components = true; quartet = true
  end
  example_simple && run_example_simple_test(clean_up)
  example_components && run_example_components_test(clean_up)
  quartet && run_quartet_test(quartet_partial, quartet_to, clean_up)
  ENV["JOT_TEST_RUNNING"] = "false"
end

function run_example_components_test(clean_up::Bool)
  clean_up_example_test()
  @testset "Example components" begin
    # Create a simple script to use as a lambda function
    open("increment_vector.jl", "w") do f
      write(f, "increment_vector(v::Vector{Int}) = map(x -> x + 1, v)")
    end

    # Turn it into a Responder; this specifies the function we will create a Lambda from
    @info "Creating responder"
    increment_responder = get_responder("./increment_vector.jl", :increment_vector, Vector{Int})

    # Create a LambdaComponents instance from the responder
    @info "Creating LambdaComponents"
    lambda_components = create_lambda_components(increment_responder; image_suffix="increment-vector")

    lambda_components |> with_remote_image! |> with_lambda_function! |> run_test

    # Clean up
    if clean_up
      delete!(lambda_components)
      rm("./increment_vector.jl")

      @test isnothing(get_local_image("increment-vector"))
      @test isnothing(get_aws_role(lambda_components.lambda_function.Role))
      @test isnothing(get_ecr_repo("increment-vector"))
      @test isnothing(get_remote_image("increment-vector"))
    end
  end
end

function run_example_simple_test(clean_up::Bool)
  clean_up_example_test()
  @testset "Example simple" begin
    # Create a simple script to use as a lambda function
    open("increment_vector.jl", "w") do f
      write(f, "increment_vector(v::Vector{Int}) = map(x -> x + 1, v)")
    end

    # Turn it into a Responder; this specifies the function we will create a Lambda from
    @info "Creating responder"
    increment_responder = get_responder("./increment_vector.jl", :increment_vector, Vector{Int})

    # Create a local docker image from the responder
    @info "Creating local image"
    local_image = create_local_image(increment_responder; image_suffix="increment-vector")
    @test get_local_image("increment-vector") == local_image
    @test get_lambda_name(local_image) == "increment-vector"

    # Push this local docker image to AWS ECR
    @info "Creating remote image"
    remote_image = push_to_ecr!(local_image)
    @test get_remote_image("increment-vector") == remote_image
    @test get_lambda_name(remote_image) == "increment-vector"

    # Create a lambda function from this remote_image...
    @info "Creating lambda function"
    increment_vector_lambda = create_lambda_function(remote_image)
    @test get_lambda_function("increment-vector") == increment_vector_lambda
    @test get_lambda_name(increment_vector_lambda) == "increment-vector"

    # ... and test it to see if it's working OK
    @info "Testing lambda function"
    @test run_test(increment_vector_lambda, [2,3,4], [3,4,5]; check_function_state=true) |> first

    # Clean up
    if clean_up
      delete!(increment_vector_lambda)
      delete!(remote_image)
      delete!(local_image)
      rm("./increment_vector.jl")
      @test isnothing(get_aws_role(increment_vector_lambda.Role))
      @test isnothing(get_ecr_repo("increment-vector"))
    end
    # Check that this has also cleaned up the ECR Repo and the AWS Role
  end
end

function clean_up_example_test()
  # Before the test, delete anything which might be left over
  existing_lf = get_lambda_function("increment-vector")
  !isnothing(existing_lf) && delete!(existing_lf)

  existing_ri = get_remote_image("increment-vector")
  !isnothing(existing_ri) && delete!(existing_ri)

  existing_ecr = get_ecr_repo("increment-vector")
  !isnothing(existing_ecr) && delete!(existing_ecr)

  existing_li = get_local_image("increment-vector")
  !isnothing(existing_li) && delete!(existing_li)
end

function run_quartet_test(partial::Bool, to::AbstractString, clean_up::Bool)
  reset_response_suffix("test/JotTest1/response_suffix")
  reset_response_suffix("test/JotTest2/response_suffix")

  responder_inputs = [
    ((JotTest1, :response_func, Dict), Dict()),
    ((PackageSpec(path=joinpath(jot_path, "test/JotTest2")), :response_func, Dict), Dict()),
    (("https://github.com/harris-chris/JotTest3", :response_func, Vector{Float64}), Dict()),
    ((joinpath(jot_path, "test/JotTest4/jot-test-4.jl"), :map_log_gamma, Vector{Float64}), Dict(:dependencies => ["SpecialFunctions"])),
  ]

  responders = Vector{AbstractResponder}()
  @testset "Test Responder" begin
    foreach(responder_inputs) do input
      args = first(input)
      kwargs = last(input)
      push!(responders, test_responder(args...; kwargs))
    end
  end
  if to == "responder"
    clean_up && quartet_clean_up()
    return
  end

  test_data = [ # Actual, expected, bad input
    (Dict("double" => 4.5), 9.0, [1,2]),
    (Dict("add suffix" => "test-"), "test-"*get_response_suffix("test/JotTest2/response_suffix"), [1,2]),
    ([1, 2, 3, 4], Vector{Float64}([1.0, 1.0, 2.0, 6.0]), "string arg"),
    ([1, 2, 3, 4], Vector{Float64}([0.0, 0.0, 0.6931471805599453, 1.791759469228055]), Dict("this" => "that")),
  ]

  local_image_inputs = [ # Responder, number, use_config, package_compile
    (responders[1], 1, false, true),
    (responders[2], 2, true, false),
    (responders[3], 3, false, false),
    (responders[4], 4, false, false),
  ]
  user_labels = [
                 Dict("TEST"=>"1"),
                 Dict("TEST"=>"2"),
                 Dict("TEST"=>"3"),
                 Dict("TEST"=>"4"),
                ]

  expected_labels = [
    ExpectedLabels("JotTest1", "response_func", joinpath(jot_path, "test/JotTest1"), user_labels[1]),
    ExpectedLabels("JotTest2", "response_func", joinpath(jot_path, "test/JotTest2"), user_labels[2]),
    ExpectedLabels("JotTest3", "response_func", "https://github.com/harris-chris/JotTest3", user_labels[3]),
    ExpectedLabels(
      Jot.get_package_name_from_script_name("jot-test-4.jl"),
      "map_log_gamma",
      joinpath(jot_path, "test/JotTest4/jot-test-4.jl"),
      user_labels[4]),
  ]

  if !(length(responders) == length(test_data) == length(local_image_inputs))
    error("Input lengths do not match")
  end

  for res in responders
    @show res.package_name
  end

  for li in local_image_inputs
    @show li[1].package_name
  end

  local_images = Vector{LocalImage}()
  @testset "Local Images" begin
    foreach(zip(local_image_inputs, expected_labels, test_data)) do (li_inputs, labels, test_datum)
      (test_input, expected_result) = (test_datum[1], test_datum[2])
      this_li = test_local_image(li_inputs..., test_input, expected_result, labels)
      push!(local_images, this_li)
    end
  end

  if to == "local_image"
    clean_up && quartet_clean_up()
    return
  end

  @testset "Package compiler" begin
    test_package_compile(;
      compiled_image=local_images[1],
      uncompiled_image=local_images[2],
      compiled_test_data=test_data[1][1:2],
      uncompiled_test_data=test_data[2][1:2],
    )
  end
  if to == "package_compiler"
    clean_up && quartet_clean_up()
    return
  end

  # If partial, randomly select one of our lambdas
  test_list = partial ? [i == rand(1:4) ? true : false for i in 1:4] : fill(true, length(responders))

  repos = Vector{Union{Nothing, ECRRepo}}()
  remote_images = Vector{Union{Nothing, RemoteImage}}()
  @testset "ECR Repo" begin
    foreach(enumerate(test_list)) do (num, use_bl)
      if use_bl
        (this_repo, this_remote_image) = test_ecr_repo(responders[num], local_images[num], expected_labels[num])
      else
        (this_repo, this_remote_image) = (nothing, nothing)
      end
      push!(repos, this_repo)
      push!(remote_images, this_remote_image)
    end
  end

  if to == "ecr_repo"
    clean_up && quartet_clean_up()
    return
  end

  aws_role = test_aws_role()
  if to == "aws_role"
    clean_up && quartet_clean_up()
    return
  end

  lambda_functions = Vector{Union{Nothing, LambdaFunction}}()
  @testset "Lambda Function" begin
    foreach(enumerate(test_list)) do (num, use_bl)
      if use_bl
        this_lambda_function = test_lambda_function(repos[num], remote_images[num], aws_role, test_data[num]...)
      else
        this_lambda_function = nothing
      end
      push!(lambda_functions, this_lambda_function)
    end
  end

  if to == "lambda_function"
    clean_up && quartet_clean_up()
    return
  end
  clean_up && clean_up()
end

function test_responder(
    res_obj::Any,
    res_func::Symbol,
    res_type::Type{IT};
    kwargs::Dict = Dict(),
  )::AbstractResponder{IT} where {IT}
  this_res = get_responder(res_obj, res_func, IT; kwargs...)
  @test isa(Jot.get_tree_hash(this_res), String)
  @test isa(Jot.get_commit(this_res), String)
  this_res
end

function test_local_image(
    res::AbstractResponder,
    num::Int64,
    use_config::Bool,
    package_compile::Bool,
    test_request::Any,
    expected_test_result::Any,
    expected_labels::ExpectedLabels,
  )::LocalImage
  local_image = create_local_image(res;
                                   aws_config = use_config ? aws_config : nothing,
                                   package_compile = package_compile,
                                   user_defined_labels = expected_labels.user_defined_labels,
                                  )
  @test Jot.matches(res, local_image)
  @test Jot.is_jot_generated(local_image)
  @test test_actual_labels_against_expected(get_labels(local_image), expected_labels)
  # Test that container runs
  cont = run_image_locally(local_image)
  @test is_container_running(cont)
  conts = get_all_containers(local_image)
  @test length(conts) == 1
  foreach(stop_container, conts)
  # Check containers have stopped
  @test length(get_all_containers(local_image)) == 0
  # Run local test of container, without value
  @test run_test(local_image) |> first
  # Run local test of container, with expected response
  @show test_request
  @test run_test(local_image, test_request, expected_test_result; then_stop=true) |> first
  sleep(1)
  return local_image
end

function test_package_compile(;
    uncompiled_image::LocalImage,
    compiled_image::LocalImage,
    uncompiled_test_data::Tuple{Any, Any},
    compiled_test_data::Tuple{Any, Any},
  )
  @show "test_package_compile"
  sleep(2)
  @show "running test on compiled"
  @show compiled_test_data
  (_, compiled_time) = run_test(compiled_image, compiled_test_data...; then_stop=true)
  @show "first test ran"
  @show (readchomp(`docker container ls`))
  sleep(2)
  (_, uncompiled_time) = run_test(uncompiled_image, uncompiled_test_data...; then_stop=true)
  @show "second test ran"
  @test compiled_time < (uncompiled_time / 2)
end

function test_ecr_repo(
    res::AbstractResponder,
    local_image::LocalImage,
    expected_labels::ExpectedLabels,
  )::Tuple{ECRRepo, RemoteImage}
  remote_image = push_to_ecr!(local_image)
  ecr_repo = remote_image.ecr_repo
  @testset "Test remote image" begin
    @test Jot.matches(local_image, ecr_repo)
    @test Jot.is_jot_generated(remote_image)
    @test test_actual_labels_against_expected(get_labels(ecr_repo), expected_labels)
    # Check we can find the repo
    @test !isnothing(Jot.get_ecr_repo(local_image))
    # Check that we can find the remote image which matches our local image
    ri_check = Jot.get_remote_image(local_image)
    @test !isnothing(ri_check)
    @test Jot.matches(local_image, remote_image)
    @test Jot.matches(res, remote_image)
  end
  (ecr_repo, remote_image)
end

function test_aws_role()::AWSRole
  aws_role =  create_aws_role("jot-test-role-"*test_suffix)
  @testset "Test AWS role" begin
    @test aws_role in get_all_aws_roles()
  end
  aws_role
end

function test_lambda_function(
    ecr_repo::ECRRepo,
    remote_image::RemoteImage,
    aws_role::AWSRole,
    test_request::Any,
    expected::Any,
    exception_request::Any,
  )::LambdaFunction
  lambda_function = create_lambda_function(ecr_repo; role = aws_role)
  @testset "Lambda Function test" begin
    @test Jot.is_jot_generated(lambda_function)
    @test Jot.matches(remote_image, lambda_function)
    # Check that we can find it
    @test lambda_function in Jot.get_all_lambda_functions()
    # Invoke it
    response = invoke_function(test_request, lambda_function; check_state=true)
    @test response == expected
    # Create the same thing using a remote image
    @test lambda_function == create_lambda_function(
      remote_image; role=aws_role, function_name="addl"*test_suffix
    )
    @test_throws LambdaException invoke_function(exception_request, lambda_function; check_state=true)
  end
  lambda_function
end

function quartet_clean_up()
  # Clean up
  # TODO clean up based on lambdas, eventually
  @show "running clean up"
  @testset "Clean up" begin
    test_lfs = [x for x in Jot.get_all_lambda_functions() if occursin(test_suffix, x.FunctionName)]
    test_repos = [x for x in Jot.get_all_ecr_repos() if occursin(test_suffix, x.repositoryName)]
    test_roles = [x for x in Jot.get_all_aws_roles() if occursin(test_suffix, x.RoleName)]
    test_local_images = [x for x in Jot.get_all_local_images() if occursin(test_suffix, x.Repository)]
    test_containers = [x for img in test_local_images for x in get_all_containers(img)]

    foreach(delete!, test_lfs)
    foreach(delete!, test_repos)
    foreach(delete!, test_roles)
    foreach(x -> delete!(x, force=true), test_local_images)
    foreach(delete!, test_containers)

    @test all([!x.exists for x in test_lfs])
    @test all([!x.exists for x in test_repos])
    @test all([!x.exists for x in test_roles])
    @test all([!x.exists for x in test_local_images])
    @test all([!x.exists for x in test_containers])
  end
end

@testset "All Tests" begin
  test_args = Dict(key => val for (key, val) in map(x -> split(x, '='), ARGS))
  test_args = Dict(Symbol(key) => (val in ["true", "false"] ? Meta.parse(val) : val) for (key, val) in test_args)
  run_tests(;test_args...)
end


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

JLLWrappers.@generate_main_file_header("SwiftShader")
JLLWrappers.@generate_main_file("SwiftShader", UUID("c404ff5a-e271-5628-8eaa-423d39d64c15"))
end  # module SwiftShader_jll


text:
using PyCall
@pyimport matplotlib.pyplot as plt
include("Gaussian.jl")

# Make a gaussian
gaussian = Gaussian([1.0, 0.5, 0.0, 0.0, 1.0, pi/6])

# Define a grid
xmin, xmax, ymin, ymax = -5.0, 5.0, -5.0, 5.0
ni, nj = 201, 201
x = Array(linspace(xmin, xmax, nj))
y = Array(linspace(ymax, ymin, ni))
dx, dy = (xmax - xmin)/(nj - 1), (ymax - ymin)/(ni - 1)
dA = dx*dy

# Evaluate the gaussian on the grid
f = Array(Float64, (ni, nj))
for(j in 1:nj)
	for(i in 1:ni)
		f[i, j] = evaluate(gaussian, x[j], y[i])
	end
end

println(sum(f)*dA)

plt.imshow(f, interpolation="nearest")
plt.show()

max_stars_repo_name: eggplantbren/MogTrack

max_issues_repo_name: eggplantbren/MogTrack

max_forks_repo_name: eggplantbren/MogTrack