text:
"""
    test_bqgql_rossbywave(; kwargs...)

Evolvesa a Rossby wave and compares with the analytic solution.
"""
function test_bqgql_rossbywave(stepper, dt, nsteps)
    nx = 64
  beta = 2.0
    Lx = 2π
    mu = 0.0
    nu = 0.0

  # the following if statement is called so that all the cases of
  # Problem() fuction are tested
  if stepper=="ForwardEuler"
    eta = zeros(nx, nx)
  else
    eta(x, y) = 0*x
  end

  prob = BarotropicQGQL.InitialValueProblem(nx=nx, Lx=Lx, eta=eta, beta=beta, mu=mu, nu=nu, stepper=stepper, dt=dt)
  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

  x, y = gridpoints(g)

  # the Rossby wave initial condition
   ampl = 1e-2
  kwave = 3.0*2π/g.Lx
  lwave = 2.0*2π/g.Ly
      ω = -p.beta*kwave/(kwave^2.0 + lwave^2.0)
     ζ0 = @. ampl*cos(kwave*x)*cos(lwave*y)
    ζ0h = rfft(ζ0)

  BarotropicQGQL.set_zeta!(prob, ζ0)

  stepforward!(prob, nsteps)
  dealias!(sol, g)
  BarotropicQGQL.updatevars!(prob)

  ζ_theory = @. ampl*cos(kwave*(x - ω/kwave*cl.t))*cos(lwave*y)

  isapprox(ζ_theory, v.zeta, rtol=g.nx*g.ny*nsteps*1e-12)
end

"""
    test_stochasticforcingbudgets(; kwargs...)

Tests if the energy budgets are closed for BarotropicQG with stochastic forcing.
"""
function test_bqgql_stochasticforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7, nnu=2, mu=1e-1)
  n, L  = 256, 2π
  nu, nnu = 1e-7, 2
  mu = 1e-1
  dt, tf = 0.005, 0.1/mu
  nt = round(Int, tf/dt)
  ns = 1

  # Forcing
  kf, dkf = 12.0, 2.0
  ε = 0.1
  gr  = TwoDGrid(n, L)
  x, y = gridpoints(gr)

  Kr = [ gr.kr[i] for i=1:gr.nkr, j=1:gr.nl]

  force2k = zero(gr.Krsq)
  @. force2k = exp.(-(sqrt(gr.Krsq)-kf)^2/(2*dkf^2))
  @. force2k[gr.Krsq .< 2.0^2 ] = 0
  @. force2k[gr.Krsq .> 20.0^2 ] = 0
  @. force2k[Kr .< 2π/L] = 0
  ε0 = parsevalsum(force2k.*gr.invKrsq/2.0, gr)/(gr.Lx*gr.Ly)
  force2k .= ε/ε0 * force2k

  Random.seed!(1234)

  function calcF!(F, sol, t, cl, v, p, g)
    eta = exp.(2π*im*rand(Float64, size(sol)))/sqrt(cl.dt)
    eta[1, 1] = 0
    @. F = eta*sqrt(force2k)
    nothing
  end

  prob = BarotropicQGQL.ForcedProblem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, dt=dt,
   stepper="RK4", calcF=calcF!, stochastic=true)

  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

  BarotropicQGQL.set_zeta!(prob, 0*x)
  E = Diagnostic(BarotropicQGQL.energy,      prob, nsteps=nt)
  D = Diagnostic(BarotropicQGQL.dissipation, prob, nsteps=nt)
  R = Diagnostic(BarotropicQGQL.drag,        prob, nsteps=nt)
  W = Diagnostic(BarotropicQGQL.work,        prob, nsteps=nt)
  diags = [E, D, W, R]

  # Step forward

  stepforward!(prob, diags, round(Int, nt))

  BarotropicQGQL.updatevars!(prob)

  cfl = cl.dt*maximum([maximum(v.v)/g.dx, maximum(v.u)/g.dy])

  E, D, W, R = diags

  t = round(mu*cl.t, digits=2)

  i₀ = 1
  dEdt = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt
  ii = (i₀):E.i-1
  ii2 = (i₀+1):E.i

  # dEdt = W - D - R?
  # If the Ito interpretation was used for the work
  # then we need to add the drift term
  # total = W[ii2]+ε - D[ii] - R[ii]      # Ito
  total = W[ii2] - D[ii] - R[ii]        # Stratonovich

  residual = dEdt - total

  # println(mean(abs.(residual)))
  isapprox(mean(abs.(residual)), 0, atol=1e-4)
end

"""
    test_stochasticforcingbudgets(; kwargs...)

Tests if the energy budgets are closed for BarotropicQG with stochastic forcing.
"""
function test_bqgql_deterministicforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7, nnu=2, mu=1e-1)
  n, L  = 256, 2π
  nu, nnu = 1e-7, 2
  mu = 1e-1
  dt, tf = 0.005, 0.1/mu
  nt = round(Int, tf/dt)
  ns = 1

  # Forcing = 0.01cos(4x)cos(5y)cos(2t)
  gr  = TwoDGrid(n, L)
  x, y = gridpoints(gr)
  f = @. 0.01*cos(4*gr.kr[2]*x)*cos(5*gr.l[2]*y)
  fh = rfft(f)

  function calcF!(Fh, sol, t, cl, v, p, g)
    @. Fh = fh*cos(2*t)
    nothing
  end

  prob = BarotropicQGQL.ForcedProblem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, dt=dt,
   stepper="RK4", calcF=calcF!, stochastic=false)

  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

  BarotropicQGQL.set_zeta!(prob, 0*x)
  E = Diagnostic(BarotropicQGQL.energy,      prob, nsteps=nt)
  D = Diagnostic(BarotropicQGQL.dissipation, prob, nsteps=nt)
  R = Diagnostic(BarotropicQGQL.drag,        prob, nsteps=nt)
  W = Diagnostic(BarotropicQGQL.work,        prob, nsteps=nt)
  diags = [E, D, W, R]

  # Step forward

  stepforward!(prob, diags, round(Int, nt))

  BarotropicQGQL.updatevars!(prob)

  cfl = cl.dt*maximum([maximum(v.v)/g.dx, maximum(v.u)/g.dy])

  E, D, W, R = diags

  t = round(mu*cl.t, digits=2)

  i₀ = 1
  dEdt = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt
  ii = (i₀):E.i-1
  ii2 = (i₀+1):E.i

  # dEdt = W - D - R?
  total = W[ii2] - D[ii] - R[ii]

  residual = dEdt - total

  # println(mean(abs.(residual)))
  isapprox(mean(abs.(residual)), 0, atol=1e-8)
end

"""
    test_bqgql_nonlinearadvection(dt, stepper; kwargs...)

Tests the advection term in the twodturb module by timestepping a
test problem with timestep dt and timestepper identified by the string stepper.
The test problem is derived by picking a solution ζf (with associated
streamfunction ψf) for which the advection term J(ψf, ζf) is non-zero. Next, a
forcing Ff is derived according to Ff = ∂ζf/∂t + J(ψf, ζf) - nuΔζf. One solution
to the vorticity equation forced by this Ff is then ζf. (This solution may not
be realized, at least at long times, if it is unstable.)
"""
function test_bqgql_advection(dt, stepper; n=128, L=2π, nu=1e-2, nnu=1, mu=0.0)
  n, L  = 128, 2π
  nu, nnu = 1e-2, 1
  mu = 0.0
  tf = 1.0
  nt = round(Int, tf/dt)

  gr  = TwoDGrid(n, L)
  x, y = gridpoints(gr)

  psif = @.    cos(3y) +  sin(2x)*cos(2y) +  2sin(x)*cos(3y)
    qf = @. - 9cos(3y) - 8sin(2x)*cos(2y) - 20sin(x)*cos(3y)

  Ff = @. nu*( 81cos(3y) + 200cos(3y)*sin(x) + 64cos(2y)*sin(2x) ) -
    3sin(3y)*(-16cos(2x)*cos(2y) - 20cos(x)*cos(3y)) -
 27sin(3y)*(2cos(2x)*cos(2y) + 2cos(x)*cos(3y)) + 0*(-8cos(x)*cos(3y)*sin(2x)*sin(2y) +
 24*cos(2x)*cos(2y)*sin(x)*sin(3y))

  Ffh = -rfft(Ff)

  # Forcing
  function calcF!(Fh, sol, t, cl, v, p, g)
    Fh .= Ffh
    nothing
  end

  prob = BarotropicQGQL.ForcedProblem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, dt=dt, stepper=stepper, calcF=calcF!)
  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
  BarotropicQGQL.set_zeta!(prob, qf)

  # Step forward
  stepforward!(prob, round(Int, nt))
  BarotropicQGQL.updatevars!(prob)
  isapprox(v.zeta+v.Zeta, qf, rtol=1e-13)
end

function test_bqgql_energyenstrophy()
  nx, Lx  = 64, 2π
  ny, Ly  = 64, 3π
  g  = TwoDGrid(nx, Lx, ny, Ly)
  k0, l0 = g.k[2], g.l[2] # fundamental wavenumbers
  x, y = gridpoints(g)

  energy_calc = 29/9
  enstrophy_calc = 2701/162

    eta = @. cos(10k0*x)*cos(10l0*y)
   psi0 = @. sin(2k0*x)*cos(2l0*y) + 2sin(k0*x)*cos(3l0*y)
  zeta0 = @. -((2k0)^2+(2l0)^2)*sin(2k0*x)*cos(2l0*y) - (k0^2+(3l0)^2)*2sin(k0*x)*cos(3l0*y)

  prob = BarotropicQGQL.InitialValueProblem(nx=nx, Lx=Lx, ny=ny, Ly=Ly, eta=eta, stepper="ForwardEuler")
  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

  BarotropicQGQL.set_zeta!(prob, zeta0)
  BarotropicQGQL.updatevars!(prob)

  energyzeta0 = BarotropicQGQL.energy(prob)
  enstrophyzeta0 = BarotropicQGQL.enstrophy(prob)

  isapprox(energyzeta0, energy_calc, rtol=1e-13) && isapprox(enstrophyzeta0, enstrophy_calc, rtol=1e-13) && BarotropicQGQL.addforcing!(prob.timestepper.N, sol, cl.t, cl, v, p, g)==nothing
end

function test_bqgql_problemtype(T=Float32)
  prob = BarotropicQGQL.Problem(T=T)

  (typeof(prob.sol)==Array{Complex{T},2} && typeof(prob.grid.Lx)==T && typeof(prob.grid.x)==Array{T,2} && typeof(prob.vars.u)==Array{T,2})
end


text:
@symbol_func function T_rho(cur_reactor::AbstractReactor, cur_rho::AbstractSymbol)
  cur_T_bar = cur_reactor.T_bar

  cur_nu_T = cur_reactor.nu_T

  cur_T_rho = 1 - cur_rho ^ 2

  cur_T_rho ^= cur_nu_T

  cur_T_rho *= cur_T_bar

  cur_T_rho *= 1 + cur_nu_T

  cur_T_rho
end

ext: jl


text:
function satisfy(interpreter::AbstractInterpreter,
                 domain::Domain, state::GenericState,
                 terms::AbstractVector{<:Term})
     # Quick check that avoids SLD resolution unless necessary
     sat = all(check(interpreter, domain, state, to_nnf(t)) for t in terms)
     if sat !== missing return sat end
     # Call SLD resolution only if there are free variables or axioms
     return !isempty(satisfiers(interpreter, domain, state, terms))
end

function satisfy(interpreter::AbstractInterpreter,
                 domain::Domain, state::GenericState, term::Term)
    # Convert to NNF
    term = to_nnf(term)
    # Quick check that avoids SLD resolution unless necessary
    sat = check(interpreter, domain, state, term)
    if sat !== missing return sat end
    # Call SLD resolution only if there are free variables or axioms
    return !isempty(satisfiers(interpreter, domain, state, term))
end

function satisfiers(interpreter::AbstractInterpreter,
                    domain::Domain, state::GenericState,
                    terms::AbstractVector{<:Term})
    # Reify negations
    terms = [reify_negations(to_nnf(t), domain) for t in terms]
    # Initialize Julog knowledge base
    clauses = Clause[get_clauses(domain); get_negation_clauses(domain);
                     collect(state.types); collect(state.facts)]
    # Pass in fluents and function definitions as a dictionary of functions
    funcs = merge(GLOBAL_PREDICATES, state.values, get_funcdefs(domain))
    return resolve(collect(terms), clauses; funcs=funcs, mode=:all)[2]
end

function satisfiers(interpreter::AbstractInterpreter,
                    domain::Domain, state::GenericState, term::Term)
    return satisfiers(domain, state, [term])
end

function check(interpreter::AbstractInterpreter,
               domain::Domain, state::GenericState, term::Compound)
    sat = if term.name == :and
        all(check(interpreter, domain, state, a) for a in term.args)
    elseif term.name == :or
        any(check(inte, domain, state, a) for a in term.args)
    elseif term.name == :imply
        !check(interpreter, domain, state, term.args[1]) |
        check(interpreter, domain, state, term.args[2])
    elseif term.name == :not
        !check(interpreter, domain, state, term.args[1]) |
        check(interpreter, domain, state, negate(term.args[1]))
    elseif term.name == :forall
        missing
    elseif term.name == :exists
        missing
    elseif !is_ground(term)
        missing
    elseif is_derived(term, domain)
        missing
    elseif is_type(term, domain)
        if has_subtypes(term, domain)
            missing
        elseif term in state.types
            true
        elseif !isempty(get_constants(domain))
            domain.constypes[term.args[1]] == term.name
        else
            false
        end
    elseif is_global_pred(term)
        evaluate(interpreter, domain, state, term)::Bool
    elseif is_func(term, domain)
        evaluate(interpreter, domain, state, term)::Bool
    else
        term = partialeval(domain, state, term)
        term in state.facts
    end
    return sat
end

function check(interpreter::AbstractInterpreter,
               domain::Domain, state::GenericState, term::Const)
    if term in state.facts || term in state.types || term.name == true
        return true
    elseif is_func(term, domain) || is_derived(term, domain)
        return missing
    else
        return false
    end
end

function check(interpreter::AbstractInterpreter,
               domain::Domain, state::GenericState, term::Var)
    return missing
end


text:
mutable struct FactorLogPredictiveDensity <: MyXMLElement
    fac_el::XMLOrNothing
    like_el::XMLOrNothing
    cl_el::XMLOrNothing # compound likelihood
    intfac::IntegratedFactorsXMLElement
    like::TraitLikelihoodXMLElement
    trait_ind::Int

    function FactorLogPredictiveDensity(
                intfac::IntegratedFactorsXMLElement,
                like::TraitLikelihoodXMLElement;
                trait_name::String = bn.DEFAULT_TRAIT_NAME)

        name_ind = findfirst(isequal(trait_name), intfac.treeModel.node_traits)

        intfac = copy(intfac)
        intfac.tree_trait_ind = name_ind
        like = copy(like)
        like.attrs[bn.ID] = "$(trait_name).treeLikelihood"
        intfac.id = "$(trait_name).factorLikelihood"
        intfac.standardize_traits = false

        like.extension_el = intfac

        return new(nothing, nothing, nothing, intfac, like, name_ind)
    end
end

function make_xml(lpd::FactorLogPredictiveDensity)

    @unpack intfac, like = lpd
    fac_el = make_xml(intfac, reference_precision = true)
    lpd.fac_el = fac_el
    trait_name = get_trait_name(intfac)

    like.extension_el = intfac
    lpd.like_el = make_xml(like)

    cl_el = new_element(bn.LIKELIHOOD)
    set_id!(cl_el, "$trait_name.likelihood")
    add_ref_el(cl_el, lpd.fac_el)
    add_ref_el(cl_el, lpd.like_el)
    lpd.cl_el = cl_el

end

function get_loggables(lpd::FactorLogPredictiveDensity)
    make_xml(lpd)
    return lpd.cl_el
end


text:
module ModelClockSystem

    export ModelAlarm, ModelClock, advanceClock!, setClock!, addAlarm!, clock2str, dt2str, dt2tuple, dropRungAlarm!

    using CFTime, Formatting

    using Dates
    using Dates: CompoundPeriod

    mutable struct ModelAlarm
        name      :: String
        time      :: AbstractCFDateTime
        priority  :: Integer   # There might be multiple alarms ring at the same time. And if order is important, then the higher the number the higher the priority.
        callbacks :: Array{Function, 1}
        recurring :: Union{Nothing, Period, CompoundPeriod, Function}
        done      :: Bool
    end

    function isEarlier(x :: ModelAlarm, y :: ModelAlarm)
        return (x.time < y.time) || ( (x.time == y.time) && (x.priority > y.priority))
    end

    function isEqual(x :: ModelAlarm, y :: ModelAlarm)
        return (x.time == y.time) && (x.priority == y.priority)
    end

    function isLater(x :: ModelAlarm, y :: ModelAlarm)
        return ! ( isEarlier(x, y) ||  isEqual(x, y) )
    end

    mutable struct ModelClock

        name         :: String
        time         :: AbstractCFDateTime
        alarms       :: Array{ModelAlarm, 1}
        alarms_dict  :: Dict
        alarm_ptr    :: Integer

        function ModelClock(
            name :: String,
            time :: AbstractCFDateTime,
        )
            alarms = Array{ModelAlarm}(undef, 0)
            alarms_dict = Dict()
            return new(
                name,
                time + Dates.Second(0),
                alarms,
                alarms_dict,
                0,
            )
        end

    end

    function checkType(time1 :: AbstractCFDateTime, time2 :: AbstractCFDateTime)
        if typeof(time1) != typeof(time2)
            throw(ErrorException("Time type does not match!"))
        end
    end

    function setClock!(
        mc :: ModelClock,
        time :: AbstractCFDateTime,
    )

        checkType(mc.time, time)
        mc.time = time + Second(0)

    end

    function advanceClock!(
        mc :: ModelClock,
        t :: Union{Second, AbstractCFDateTime},
    )
        if typeof(t) <: Second  # inteprete as time interval
            mc.time += t
        elseif typeof(Δt) <: AbstractCFDateTime
            setClock!(mc, t)
        end

        # Only check alarms when there is any
        if length(mc.alarms) > 0
            if mc.alarm_ptr == 0 && mc.time >= mc.alarms[1].time
                mc.alarm_ptr = 1
            end

            while (0 < mc.alarm_ptr <= length(mc.alarms)) && (mc.time >= mc.alarms[mc.alarm_ptr].time)
            #    println("Before ring alarm: mc.alarm_ptr  = $(mc.alarm_ptr)")
                # In case we are at the end of alarm array

                alarm_ptr_tmp = mc.alarm_ptr
                if mc.alarm_ptr < length(mc.alarms)
                    mc.alarm_ptr += 1
                end
                ringAlarm!(mc, mc.alarms[alarm_ptr_tmp])
             #   println("xxx alarm_ptr = $(mc.alarm_ptr)")
                if mc.alarm_ptr == length(mc.alarms)
                    break
                end
            end
        end

    end

    function ringAlarm!(mc :: ModelClock, alarm :: ModelAlarm)
        #println("alarm.time = $(alarm.time). done = $(alarm.done)")
        if ! alarm.done
            println(format("Alarm '{:s}' rings at {:s}", alarm.name, dt2str(alarm.time)))

            for callback in alarm.callbacks
                callback(mc, alarm)
            end

            # IMPORTANT: alarm.done has to be here.
            # If it is after addAlarm!, there is a potential recursive loop because
            # advanceClock! might be called when calling addAlarm! (i.e. ring immediately)
            alarm.done = true

            if alarm.recurring != nothing

                if isa(alarm.recurring, Function)
                    next_time = alarm.recurring(alarm.time)
                else
                    next_time = alarm.time + alarm.recurring
                end
                #println("Current time: $(string(alarm.time))")
                #println("Next alarm: " * string(next_time))
                addAlarm!(
                    mc,
                    alarm.name,
                    next_time,
                    alarm.priority;
                    callback = alarm.callbacks,
                    recurring = alarm.recurring,
                )
            end
        end
    end

    function clock2str(mc :: ModelClock)
        return dt2str(mc.time)
    end

    function dt2str(dt)
        return format("{:04d}-{:02d}-{:02d} {:02d}:{:02d}:{:02d}", dt2tuple(dt)...)
    end

    function dt2tuple(dt)
        return (Dates.year(dt), Dates.month(dt), Dates.day(dt), Dates.hour(dt), Dates.minute(dt), Dates.second(dt))
    end

    function addAlarm!(
        mc   :: ModelClock,
        name :: String,
        time :: AbstractCFDateTime,
        priority :: Integer;
        callback :: Union{Array{Function, 1}, Function, Nothing} = nothing,
        recurring :: Union{ Nothing, Period, CompoundPeriod, Function } = nothing,
    )

        checkType(time, mc.time)

        if callback == nothing
            callbacks = Array{Function, 1}(undef, 0)
        elseif typeof(callback) <: Function
            callbacks = [ callback ]
        else # it is already an array
            callbacks = callback
        end

        alarm = ModelAlarm(
            name,
            time,
            priority,
            callbacks,
            recurring,
            false,
        )

        if mc.alarm_ptr == 0
            mc.alarm_ptr = 1
        end

        if alarm.time < mc.time
            println("alarm.time = ", dt2str(alarm.time))
            throw(ErrorException("Alarm needs to be set in the future."))
        end

        push!(mc.alarms, alarm)
        if ! haskey(mc.alarms_dict, name)
            mc.alarms_dict[name] = []
        end

        if (mc.alarm_ptr > 1) && (isEarlier(alarm, mc.alarms[mc.alarm_ptr]) || isEqual(alarm, mc.alarms[mc.alarm_ptr]))
            mc.alarm_ptr -= 1
        end

        push!(mc.alarms_dict[name], alarm)
        # lt = less than = earlier and higher priority
        sort!(mc.alarms, lt = isEarlier)

        if alarm.time == mc.time
            advanceClock!(mc, Second(0))
        end

    end

    function dropRungAlarm!(
        mc   :: ModelClock,
    )
        while length(mc.alarms) > 0 && mc.alarms[1].done
            alarm = popfirst!(mc.alarms)
            delete!(mc.alarms_dict, alarm.name)
            mc.alarm_ptr -= 1
            if mc.alarm_ptr < 0
                throw(ErrorException("Error happens in alarm counts. Please check."))
            end
        end
    end

    function listAlarms(
        mc :: ModelClock,
    )
        for (i, alarm) in enumerate(mc.alarms)
            println("Alarm[$(i)] = $(alarm.time). [$( (alarm.done) ? "v" : " " )] $( (i==mc.alarm_ptr) ? "*" : "" )")
        end
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
n = rand(2: 10)

# test Hessenberg

A = matrixdepot("rohess", n)

B = tril(A) - diagm(diag(A)) - diagm(diag(A, -1), -1)

@test_approx_eq_eps B zeros(n,n) 1.0e-7

# test orithogonality

@test_approx_eq_eps A'*A eye(n) 1.0e-7

println("'rohess' passed test...")


text:
function add_2!(a::Int8, b::Int8, c::Int8, i::Int8)

    c = c + a + b
    c = [c]

    writedlm("./src/reprexs/data_$i.csv", c, ",")
end

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
@testset "DoublyStrictlyLexicographicallyGreaterThan2LexicographicallyGreaterThan: $(fct_type), $(x_size) × $(y_size), $(T)" for fct_type in ["vector of variables", "vector affine function"], x_size in [2, 4], y_size in [2, 4], T in [Int, Float64]
    base_model = if T == Int
        IntStrictlyLexicographicallyGreaterThanModel{T}()
    elseif T == Float64
        FloatStrictlyLexicographicallyGreaterThanModel{T}()
    else
        @assert false
    end
    mock = MOIU.MockOptimizer(base_model)
    model = COIB.DoublyStrictlyLexicographicallyGreaterThan2LexicographicallyGreaterThan{T}(mock)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Strictly{CP.LexicographicallyGreaterThan, T},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Strictly{CP.DoublyLexicographicallyGreaterThan, T},
    )

    if T == Int
        x_array, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:(x_size * y_size)])
    elseif T == Float64
        x_array = MOI.add_variables(model, x_size * y_size)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables(x_array)
    elseif fct_type == "vector affine function"
        MOIU.vectorize(x_array)
    else
        @assert false
    end
    c = MOI.add_constraint(model, fct, CP.Strictly{CP.DoublyLexicographicallyGreaterThan, T}(CP.DoublyLexicographicallyGreaterThan(x_size, y_size)))

    for i in 1:(x_size * y_size)
        @test MOI.is_valid(model, x_array[i])
    end
    @test MOI.is_valid(model, c)

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

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.Strictly{CP.DoublyLexicographicallyGreaterThan, T}) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == Tuple{Type}[]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.VectorAffineFunction{T}, CP.Strictly{CP.LexicographicallyGreaterThan, T}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 0
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, CP.Strictly{CP.LexicographicallyGreaterThan, T}}()) == 2

        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, CP.Strictly{CP.LexicographicallyGreaterThan, T}}()) == [bridge.con, bridge.con_transposed]
    end

    @testset "Column constraint" begin
        @test MOI.is_valid(model, bridge.con)
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con) == CP.Strictly{CP.LexicographicallyGreaterThan, T}(CP.LexicographicallyGreaterThan(x_size, y_size))

        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con)
        @test length(f.terms) == x_size * y_size
        @test f.constants == zeros(T, x_size * y_size)

        idx = LinearIndices((x_size, y_size))
        for i in 1:x_size
            for j in 1:y_size
                t = f.terms[idx[i, j]]
                @test t.output_index == idx[i, j]
                @test t.scalar_term.coefficient === one(T)
                @test t.scalar_term.variable === x_array[idx[i, j]]
            end
        end
    end

    @testset "Row constraint (transposed)" begin
        @test MOI.is_valid(model, bridge.con_transposed)
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_transposed) == CP.Strictly{CP.LexicographicallyGreaterThan, T}(CP.LexicographicallyGreaterThan(y_size, x_size))

        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_transposed)
        @test length(f.terms) == x_size * y_size
        @test f.constants == zeros(T, x_size * y_size)

        idx = LinearIndices((y_size, x_size))
        for i in 1:x_size
            for j in 1:y_size
                t = f.terms[idx[j, i]]
                @test t.output_index == idx[j, i]
                @test t.scalar_term.coefficient === one(T)
                @test t.scalar_term.variable === x_array[idx[j, i]]
            end
        end
    end
end


text:
using InelasticScattering
using Plots
λ = 440*1.e-7 #cm
ν̃ = 1/λ
T = 250 #k
RRS = true
VS_0to1 = false
VS_1to0 = false

o2 = InelasticScattering.getMolecularConstants(InelasticScattering.O₂(), (0.2));
compute_effective_coefficents!(ν̃, T, o2);
compute_energy_levels!(o2);
compute_σ_Rayl_coeff!(o2);
compute_σ_Rayl_VibRaman_coeff_hires!(T, o2);
compute_σ_VibRaman_coeff!(T, o2);
compute_σ_RoVibRaman_coeff!(T, o2);
n2 = InelasticScattering.getMolecularConstants(InelasticScattering.N₂(), (0.8));
compute_effective_coefficents!(ν̃, T, n2);
compute_energy_levels!(n2);
compute_σ_Rayl_coeff!(n2);
compute_σ_Rayl_VibRaman_coeff_hires!(T, n2);
compute_σ_VibRaman_coeff!(T, n2);
compute_σ_RoVibRaman_coeff!(T, n2);

# Create grid:
grid_out = -250:0.002:250
σ_out = similar(collect(grid_out));
σ_out .= 0;

apply_lineshape!(0, n2.effCoeff.σ_Rayl_coeff, λ, collect(grid_out), σ_out, 1, 300.0, 28);
atmo_σ_Rayl = n2.vmr * σ_out #cross section in cm^2

if (RRS)
    # N2
    apply_lineshape!(n2.effCoeff.Δν̃_RoRaman_coeff_JtoJp2, n2.effCoeff.σ_RoRaman_coeff_JtoJp2,  λ, collect(grid_out), σ_out, 1, 300.0, 28);
    atmo_σ_RRS_JtoJp2 = n2.vmr * σ_out #cross section in cm^2

    apply_lineshape!(n2.effCoeff.Δν̃_RoRaman_coeff_JtoJm2, n2.effCoeff.σ_RoRaman_coeff_JtoJm2, λ, collect(grid_out), σ_out, 1, 300.0, 40);
    atmo_σ_RRS_JtoJm2 = n2.vmr * σ_out #cross section in cm^2

    # O2
    apply_lineshape!(o2.effCoeff.Δν̃_RoRaman_coeff_JtoJp2, o2.effCoeff.σ_RoRaman_coeff_JtoJp2, λ, collect(grid_out), σ_out, 1, 300.0, 28);
    atmo_σ_RRS_JtoJp2 += o2.vmr * σ_out #cross section in cm^2

    apply_lineshape!(o2.effCoeff.Δν̃_RoRaman_coeff_JtoJm2, o2.effCoeff.σ_RoRaman_coeff_JtoJm2, λ, collect(grid_out), σ_out, 1, 300.0, 40);
    atmo_σ_RRS_JtoJm2 += o2.vmr * σ_out #cross section in cm^2

elseif (VS_0to1)
    # N2
    xin = [n2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJp2 n2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJm2]
    yin = [n2.effCoeff.σ_RoVibRaman_coeff_0to1_JtoJp2 n2.effCoeff.σ_RoVibRaman_coeff_0to1_JtoJm2];
    apply_lineshape!(xin, yin, λ, collect(grid_out), σ_out, 1, 300.0, 40);
    atmo_σ_RVRS_0to1 = n2.vmr * σ_out #cross section in cm^2

    xin = n2.effCoeff.Δν̃_VibRaman_coeff_0to1_hires
    yin = n2.effCoeff.σ_VibRaman_coeff_0to1_hires
    apply_lineshape!(xin, yin, λ, collect(grid_out), σ_out, 1, 300.0, 40);
    atmo_σ_VRS_0to1 = n2.vmr * σ_out #cross section in cm^2

    # O2
    xin = [o2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJp2 o2.effCoeff.Δν̃_RoVibRaman_coeff_0to1_JtoJm2]
    yin = [o2.effCoeff.σ_RoVibRaman_coeff_0to1_JtoJp2 o2.effCoeff.σ_RoVibRaman_coeff_0to1_JtoJm2];
    apply_lineshape!(xin, yin, λ, collect(grid_out), σ_out, 1, 300.0, 40);
    atmo_σ_RVRS_0to1 += o2.vmr * σ_out #cross section in cm^2

    xin = o2.effCoeff.Δν̃_VibRaman_coeff_0to1_hires
    yin = o2.effCoeff.σ_VibRaman_coeff_0to1_hires
    apply_lineshape!(xin, yin, λ, collect(grid_out), σ_out, 1, 300.0, 40);
    atmo_σ_VRS_0to1 += o2.vmr * σ_out #cross section in cm^2
elseif (VS_1to0)
    # N2
    xin = [n2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJp2 n2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJm2]
    yin = [n2.effCoeff.σ_RoVibRaman_coeff_1to0_JtoJp2 n2.effCoeff.σ_RoVibRaman_coeff_1to0_JtoJm2];
    apply_lineshape!(xin, yin, λ, collect(grid_out), σ_out, 1, 300.0, 40);
    atmo_σ_RVRS_1to0 = n2.vmr * σ_out #cross section in cm^2

    xin = n2.effCoeff.Δν̃_VibRaman_coeff_1to0_hires
    yin = n2.effCoeff.σ_VibRaman_coeff_1to0_hires
    apply_lineshape!(xin, yin, λ, collect(grid_out), σ_out, 1, 300.0, 40);
    atmo_σ_VRS_1to0 = n2.vmr * σ_out #cross section in cm^2

    # O2
    xin = [o2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJp2 o2.effCoeff.Δν̃_RoVibRaman_coeff_1to0_JtoJm2]
    yin = [o2.effCoeff.σ_RoVibRaman_coeff_1to0_JtoJp2 o2.effCoeff.σ_RoVibRaman_coeff_1to0_JtoJm2];
    apply_lineshape!(xin, yin, λ, collect(grid_out), σ_out, 1, 300.0, 40);
    atmo_σ_RVRS_1to0 += o2.vmr * σ_out #cross section in cm^2

    xin = o2.effCoeff.Δν̃_VibRaman_coeff_1to0_hires
    yin = o2.effCoeff.σ_VibRaman_coeff_1to0_hires
    apply_lineshape!(xin, yin, λ, collect(grid_out), σ_out, 1, 300.0, 40);
    atmo_σ_VRS_1to0 += o2.vmr * σ_out #cross section in cm^2
end
plotly()
plot(grid_out, σ_out*1e50, yscale=:log10)

max_stars_repo_name: RadiativeTransfer/InelasticScattering

max_issues_repo_name: RadiativeTransfer/InelasticScattering

max_forks_repo_name: RadiativeTransfer/InelasticScattering


text:
using Millboard
using Base.Test

str = "The quick brown fox jumps over the lazy dog α,β,γ"
regexz = [r"b[\w]*n", r"[\w]{4,}"]
functions = [match, matchall]
board = [func(regex, str) for regex=regexz, func=functions]
@test """
+=============+=====================+==========+
|             |               match | matchall |
+=============+=====================+==========+
|  r"b[\\w]*n" | RegexMatch("brown") |    brown |
+-------------+---------------------+----------+
|             | RegexMatch("quick") |    quick |
|             |                     |    brown |
|             |                     |    jumps |
|             |                     |     over |
| r"[\\w]{4,}" |                     |     lazy |
+-------------+---------------------+----------+""" == table(board, :colnames=>functions, :rownames=>regexz) |> string


text:
## MatrixBeta

function Distributions.logpdf(d::MatrixBeta, X::AbstractArray{<:TrackedMatrix{<:Real}})
    return mapvcat(x -> logpdf(d, x), X)
end
@adjoint function Distributions.logpdf(d::MatrixBeta, X::AbstractArray{<:Matrix{<:Real}})
    f(d, X) = map(x -> logpdf(d, x), X)
    return pullback(f, d, X)
end

# Adapted from Distributions.jl

## Wishart

using StatsFuns: logtwo, logmvgamma

struct TuringWishart{T<:Real, ST <: Cholesky} <: ContinuousMatrixDistribution
    df::T     # degree of freedom
    chol::ST  # the Cholesky of scale matrix
    c0::T     # the logarithm of normalizing constant in pdf
end

#### Constructors

function TuringWishart(d::Wishart)
    return TuringWishart(d.df, getchol(d.S), d.c0)
end
getchol(p::PDMat) = p.chol
getchol(p::PDiagMat) = Diagonal(map(sqrt, p.diag))
getchol(p::ScalMat) = Diagonal(fill(sqrt(p.value), p.dim))

function TuringWishart(df::T, S::AbstractMatrix) where {T <: Real}
    p = size(S, 1)
    df > p - 1 || error("dpf should be greater than dim - 1.")
    C = cholesky(S)
    return TuringWishart(df, C)
end
function TuringWishart(df::T, C::Cholesky) where {T <: Real}
    c0 = _wishart_c0(df, C)
    R = Base.promote_eltype(T, c0)
    return TuringWishart(R(df), C, R(c0))
end

function _wishart_c0(df::Real, C::Cholesky)
    h_df = df / 2
    p = size(C, 1)
    h_df * (logdet(C) + p * float(typeof(df))(logtwo)) + logmvgamma(p, h_df)
end

#### Properties

Distributions.insupport(::Type{TuringWishart}, X::Matrix) = isposdef(X)
Distributions.insupport(d::TuringWishart, X::Matrix) = size(X) == size(d) && isposdef(X)

Distributions.dim(d::TuringWishart) = size(d.chol, 1)
Base.size(d::TuringWishart) = (p = Distributions.dim(d); (p, p))
Base.size(d::TuringWishart, i) = size(d)[i]
LinearAlgebra.rank(d::TuringWishart) = Distributions.dim(d)

#### Statistics

Distributions.mean(d::TuringWishart) = d.df * Matrix(d.chol)

function Distributions.mode(d::TuringWishart)
    r = d.df - Distributions.dim(d) - 1.0
    if r > 0.0
        return Matrix(d.chol) * r
    else
        error("mode is only defined when df > p + 1")
    end
end

function Distributions.meanlogdet(d::TuringWishart)
    p = Distributions.dim(d)
    df = d.df
    v = logdet(d.chol) + p * logtwo
    for i = 1:p
        v += digamma(0.5 * (df - (i - 1)))
    end
    return v
end

function Distributions.entropy(d::TuringWishart)
    p = Distributions.dim(d)
    df = d.df
    d.c0 - 0.5 * (df - p - 1) * Distributions.meanlogdet(d) + 0.5 * df * p
end

#  Gupta/Nagar (1999) Theorem 3.3.15.i
function Distributions.cov(d::TuringWishart, i::Integer, j::Integer, k::Integer, l::Integer)
    S = Matrix(d.chol)
    d.df * (S[i, k] * S[j, l] + S[i, l] * S[j, k])
end

function Distributions.var(d::TuringWishart, i::Integer, j::Integer)
    S = Matrix(d.chol)
    d.df * (S[i, i] * S[j, j] + S[i, j] ^ 2)
end

#### Evaluation

function Distributions.logpdf(d::Wishart, X::TrackedMatrix)
    return logpdf(TuringWishart(d), X)
end
function Distributions.logpdf(d::Wishart, X::AbstractArray{<:TrackedMatrix})
    return logpdf(TuringWishart(d), X)
end
function Distributions.logpdf(d::TuringWishart, X::AbstractMatrix{<:Real})
    df = d.df
    p = Distributions.dim(d)
    Xcf = cholesky(X)
    return 0.5 * ((df - (p + 1)) * logdet(Xcf) - tr(d.chol \ X)) - d.c0
end
function Distributions.logpdf(d::TuringWishart, X::AbstractArray{<:AbstractMatrix{<:Real}})
    return mapvcat(x -> logpdf(d, x), X)
end
function Distributions.logpdf(d::TuringWishart, X::AbstractArray{<:Matrix{<:Real}})
    return mapvcat(x -> logpdf(d, x), X)
end

#### Sampling
function Distributions._rand!(rng::AbstractRNG, d::TuringWishart, A::AbstractMatrix)
    _wishart_genA!(rng, Distributions.dim(d), d.df, A)
    unwhiten!(d.chol, A)
    A .= A * A'
end

function _wishart_genA!(rng::AbstractRNG, p::Int, df::Real, A::AbstractMatrix)
    # Generate the matrix A in the Bartlett decomposition
    #
    #   A is a lower triangular matrix, with
    #
    #       A(i, j) ~ sqrt of Chisq(df - i + 1) when i == j
    #               ~ Normal()                  when i > j
    #
    A .= zero(eltype(A))
    for i = 1:p
        @inbounds A[i,i] = rand(rng, Chi(df - i + 1.0))
    end
    for j in 1:p-1, i in j+1:p
        @inbounds A[i,j] = randn(rng)
    end
end

function unwhiten!(C::Cholesky, x::StridedVecOrMat)
    cf = C.U
    lmul!(transpose(cf), x)
end

## InverseWishart

struct TuringInverseWishart{T<:Real, ST<:AbstractMatrix} <: ContinuousMatrixDistribution
    df::T     # degree of freedom
    S::ST     # Scale matrix
    c0::T     # log of normalizing constant
end

#### Constructors

function TuringInverseWishart(d::InverseWishart)
    d = TuringInverseWishart(d.df, getmatrix(d.Ψ), d.c0)
end
getmatrix(p::PDMat) = p.mat
getmatrix(p::PDiagMat) = Diagonal(p.diag)
getmatrix(p::ScalMat) = Diagonal(fill(p.value, p.dim))

function TuringInverseWishart(df::T, Ψ::AbstractMatrix) where T<:Real
    p = size(Ψ, 1)
    df > p - 1 || error("df should be greater than dim - 1.")
    C = cholesky(Ψ)
    c0 = _invwishart_c0(df, C)
    R = Base.promote_eltype(T, c0)
    return TuringInverseWishart(R(df), Ψ, R(c0))
end
function _invwishart_c0(df::Real, C::Cholesky)
    h_df = df / 2
    p = size(C, 1)
    h_df * (p * float(typeof(df))(logtwo) - logdet(C)) + logmvgamma(p, h_df)
end

#### Properties

Distributions.insupport(::Type{TuringInverseWishart}, X::Matrix) = isposdef(X)
Distributions.insupport(d::TuringInverseWishart, X::Matrix) = size(X) == size(d) && isposdef(X)

Distributions.dim(d::TuringInverseWishart) = size(d.S, 1)
Base.size(d::TuringInverseWishart) = (p = Distributions.dim(d); (p, p))
Base.size(d::TuringInverseWishart, i) = size(d)[i]
LinearAlgebra.rank(d::TuringInverseWishart) = Distributions.dim(d)

#### Statistics

function Distributions.mean(d::TuringInverseWishart)
    df = d.df
    p = Distributions.dim(d)
    r = df - (p + 1)
    if r > 0.0
        return d.S * (1.0 / r)
    else
        error("mean only defined for df > p + 1")
    end
end

Distributions.mode(d::TuringInverseWishart) = d.S * inv(d.df + Distributions.dim(d) + 1.0)

#  https://en.wikipedia.org/wiki/Inverse-Wishart_distribution#Moments
function Distributions.cov(d::TuringInverseWishart, i::Integer, j::Integer, k::Integer, l::Integer)
    p, ν, Ψ = (Distributions.dim(d), d.df, d.S)
    ν > p + 3 || throw(ArgumentError("cov only defined for df > dim + 3"))
    inv((ν - p)*(ν - p - 3)*(ν - p - 1)^2)*(2Ψ[i,j]*Ψ[k,l] + (ν-p-1)*(Ψ[i,k]*Ψ[j,l] + Ψ[i,l]*Ψ[k,j]))
end

function Distributions.var(d::TuringInverseWishart, i::Integer, j::Integer)
    p, ν, Ψ = (Distributions.dim(d), d.df, d.S)
    ν > p + 3 || throw(ArgumentError("var only defined for df > dim + 3"))
    inv((ν - p)*(ν - p - 3)*(ν - p - 1)^2)*(ν - p + 1)*Ψ[i,j]^2 + (ν - p - 1)*Ψ[i,i]*Ψ[j,j]
end

#### Evaluation

function Distributions.logpdf(d::InverseWishart, X::TrackedMatrix)
    return logpdf(TuringInverseWishart(d), X)
end
function Distributions.logpdf(d::InverseWishart, X::AbstractArray{<:TrackedMatrix})
    return logpdf(TuringInverseWishart(d), X)
end
function Distributions.logpdf(d::TuringInverseWishart, X::AbstractMatrix{<:Real})
    p = Distributions.dim(d)
    df = d.df
    Xcf = cholesky(X)
    # we use the fact: tr(Ψ * inv(X)) = tr(inv(X) * Ψ) = tr(X \ Ψ)
    Ψ = d.S
    -0.5 * ((df + p + 1) * logdet(Xcf) + tr(Xcf \ Ψ)) - d.c0
end
function Distributions.logpdf(d::TuringInverseWishart, X::AbstractArray{<:AbstractMatrix{<:Real}})
    return mapvcat(x -> logpdf(d, x), X)
end
function Distributions.logpdf(d::TuringInverseWishart, X::AbstractArray{<:Matrix{<:Real}})
    return mapvcat(x -> logpdf(d, x), X)
end

#### Sampling

function Distributions._rand!(rng::AbstractRNG, d::TuringInverseWishart, A::AbstractMatrix)
    X = Distributions._rand!(rng, TuringWishart(d.df, inv(cholesky(d.S))), A)
    A .= inv(cholesky!(X))
end

## Adjoints

@adjoint function Distributions.Wishart(df::Real, S::AbstractMatrix{<:Real})
    return pullback(TuringWishart, df, S)
end
@adjoint function Distributions.InverseWishart(df::Real, S::AbstractMatrix{<:Real})
    return pullback(TuringInverseWishart, df, S)
end

Distributions.Wishart(df::TrackedReal, S::Matrix{<:Real}) = TuringWishart(df, S)
Distributions.Wishart(df::TrackedReal, S::AbstractMatrix{<:Real}) = TuringWishart(df, S)
Distributions.Wishart(df::Real, S::TrackedMatrix) = TuringWishart(df, S)
Distributions.Wishart(df::TrackedReal, S::TrackedMatrix) = TuringWishart(df, S)

Distributions.InverseWishart(df::TrackedReal, S::Matrix{<:Real}) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::TrackedReal, S::AbstractMatrix{<:Real}) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::Real, S::TrackedMatrix) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::TrackedReal, S::TrackedMatrix) = TuringInverseWishart(df, S)


text:
@testset "H-W test" begin
    # <https://en.wikipedia.org/wiki/Hardy–Weinberg_principle#Example_.7F.27.22.60UNIQ--postMath-0000001E-QINU.60.22.27.7F_test_for_deviation>
    #@code_warntype hwe(1469, 138, 5)
    @test VCFTools.hwe(1469, 138, 5) ≈ 0.406002899
end

@testset "openvcf" begin
    # download and extract test file if not exist
    isfile("test.08Jun17.d8b.vcf.gz") || download("http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz",
        joinpath(dirname(pathof(VCFTools)), "..", "test/test.08Jun17.d8b.vcf.gz"))
    @test_throws ArgumentError openvcf("test.08Jun17.d8b.vcf.gz", "r+")
    @test_throws ArgumentError openvcf("test.08Jun17.d8b.vcf.gz", "w+")
    @test_throws ArgumentError openvcf("test.08Jun17.d8b.vcf.gz", "c")
    # VCF file name has to end with vcf or vcf.gz
    @test_throws ArgumentError openvcf("test.08Jun17.d8b")
end

@testset "nrecords(vcf)" begin
    # download and extract test file if not exist
    isfile("test.08Jun17.d8b.vcf.gz") ||   download("http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz",
        joinpath(dirname(pathof(VCFTools)), "..", "test/test.08Jun17.d8b.vcf.gz"))
    @test nrecords("test.08Jun17.d8b.vcf.gz") == 1356
end

@testset "nsamples(vcf)" begin
    # download test file if not exist
    isfile("test.08Jun17.d8b.vcf.gz") ||  download("http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz",
        joinpath(dirname(pathof(VCFTools)), "..", "test/test.08Jun17.d8b.vcf.gz"))
    @test nsamples("test.08Jun17.d8b.vcf.gz") == 191
end

@testset "gtstats(record)" begin
    record = VCF.Record("20\t14370\trs6054257\tG\tA\t29\tPASS\tNS=3;DP=14;AF=0.5;DB;H2\tGT:GQ:DP:HQ\t0|0:48:1:51,51\t1|0:48:8:51,51")
    #@code_warntype gtstats(record)
    @test gtstats(record)[1:end-1] == (1, 1, 0, 3, 1, 0.25, 0.75, 0, false, 0.25)
    @test gtstats(record)[end] ≈ 0.82414089574
    # with missing data
    record = VCF.Record("20\t14370\trs6054257\tG\tA\t29\tPASS\tNS=3;DP=14;AF=0.5;DB;H2\tGT:GQ:DP:HQ\t.\t1|0:48:8:51,51")
    @test gtstats(record)[1:end-1] == (0, 1, 0, 1, 1, 0.5, 0.5, 1, false, 0.5)
    @test gtstats(record)[end] ≈ 0.31731050786291415
end

@testset "gtstats(vcf)" begin
    #@code_warntype gtstats("test.08Jun17.d8b.vcf.gz", STDOUT)
    #@time records, samples, lines = gtstats("test.08Jun17.d8b.vcf.gz")
    #@time gtstats("test.08Jun17.d8b.vcf")
    # output tuple: (records, samples, lines, missings_by_sample, missings_by_record, maf_by_record, minorallele_by_record)

    # download and extract test file and unzip
    isfile("test.08Jun17.d8b.vcf.gz") || download("http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz",
        joinpath(dirname(pathof(VCFTools)), "..", "test/test.08Jun17.d8b.vcf.gz"))
    write("test.08Jun17.d8b.vcf", read(GzipDecompressorStream(open("test.08Jun17.d8b.vcf.gz", "r"))))

    @testset "input: text file, output: text file" begin
        @time out = gtstats("test.08Jun17.d8b.vcf", "gtstats.out.txt")
        @test out[1:3] == (1356, 191, 1356)
        outtxt = readdlm("gtstats.out.txt")
        @test all(outtxt[:, 9]  .== out[5]) # missings_by_record
        @test all(outtxt[:, 14] .== out[6]) # maf_by_record
    end

    @testset "input: text file, output: gz file" begin
        @time out = gtstats("test.08Jun17.d8b.vcf", "gtstats.out.gz")
        @test out[1:3] == (1356, 191, 1356)
        outtxt = readdlm(GzipDecompressorStream(open("gtstats.out.gz", "r")))
        @test all(outtxt[:, 9]  .== out[5]) # missings_by_record
        @test all(outtxt[:, 14] .== out[6]) # maf_by_record
    end

    @testset "input: gz file, output: text file" begin
        @time out = gtstats("test.08Jun17.d8b.vcf.gz", "gtstats.out.txt")
        @test out[1:3] == (1356, 191, 1356)
        outtxt = readdlm("gtstats.out.txt")
        @test all(outtxt[:, 9]  .== out[5]) # missings_by_record
        @test all(outtxt[:, 14] .== out[6]) # maf_by_record
    end

    @testset "input: gz file, output: gz file" begin
        @time out = gtstats("test.08Jun17.d8b.vcf.gz", "gtstats.out.gz")
        @test out[1:3] == (1356, 191, 1356)
        outtxt = readdlm(GzipDecompressorStream(open("gtstats.out.gz", "r")))
        @test all(outtxt[:, 9]  .== out[5]) # missings_by_record
        @test all(outtxt[:, 14] .== out[6]) # maf_by_record
    end

    #@time records, samples, lines = gtstats("test.08Jun17.d8b.vcf.gz", "gtstats.out.gz")
    #@time records, samples, lines = gtstats("chr22.1kg.ref.phase1_release_v3.20101123.vcf.gz", "gtstats.out.gz") # about 180 seconds
    #@show maf_by_record
end

@testset "get sampleID" begin
    # download and extract test file if not exist
    isfile("test.08Jun17.d8b.vcf.gz") || download("http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz",
        joinpath(dirname(pathof(VCFTools)), "..", "test/test.08Jun17.d8b.vcf.gz"))
    ids = sampleID("test.08Jun17.d8b.vcf.gz")
    @test length(ids) == 191
    @test typeof(ids) == Array{String, 1}
    @test ids[1] == "HG00096"
    @test ids[end] == "HG00428"
end

@testset "Fishers exact HWE test" begin
    # See p20-21 in http://courses.washington.edu/b516/lectures_2009/HWE_Lecture.pdf
    n0 = 179
    n1 = 21
    N = (n0 + n1) >> 1
    n01 = 1
    @test VCFTools.fisher_exact(n01, N, n0) < .000001
    n01 = 3
    @test VCFTools.fisher_exact(n01, N, n0) < .000001
    n01 = 5
    @test VCFTools.fisher_exact(n01, N, n0) < .000001
    n01 = 7
    @test round(VCFTools.fisher_exact(n01, N, n0), digits=6) == .000001
    n01 = 9
    @test round(VCFTools.fisher_exact(n01, N, n0), digits=6) == .000047
    n01 = 11
    @test round(VCFTools.fisher_exact(n01, N, n0), digits=6) == .000870
    n01 = 13
    @test round(VCFTools.fisher_exact(n01, N, n0), digits=6) == .009375
    n01 = 15
    @test round(VCFTools.fisher_exact(n01, N, n0), digits=6) == .059283
    n01 = 17
    @test round(VCFTools.fisher_exact(n01, N, n0), digits=6) == .214465
    n01 = 19
    @test round(VCFTools.fisher_exact(n01, N, n0), digits=6) == .406355
    n01 = 21
    @test round(VCFTools.fisher_exact(n01, N, n0), digits=6) == .309604

    # Picked a few values in Table 1 of
    # AJHG Volume 76, Issue 5, May 2005, Pages 887-893
    # https://www.sciencedirect.com/science/article/pii/S0002929707607356
    N  = 100
    n0 = 21
    o01 = 11
    hwepval, ts_low, ts_high = VCFTools.hwe_fisher(o01, N, n0)
    @test round(hwepval, digits=6) == .000919
    @test round(ts_low, digits=6)  == .000919
    @test round(ts_high, digits=6) == .999952
    o01 = 17
    hwepval, ts_low, ts_high = VCFTools.hwe_fisher(o01, N, n0)
    @test round(hwepval, digits=6) == .284042
    @test round(ts_low, digits=6)  == .284042
    @test round(ts_high, digits=6) == .930424
    o01 = 21
    hwepval, ts_low, ts_high = VCFTools.hwe_fisher(o01, N, n0)
    @test round(hwepval, digits=6) == .593645
    @test round(ts_low, digits=6)  == 1.0
    @test round(ts_high, digits=6) == .309604
end


text:
module Execution

using ..HSA

debug_print = HSA.debug_print

include("arguments.jl")
include("kernel.jl")
include("packet.jl")
include("hsa_kernel.jl")
end # module Execution

using .Execution

export @hsa
macro hsa(range, call)
    # support plain integer for 1D ranges
    if isa(range, Integer) || isa(range, Symbol)
        range = :(($range,))
    end

    if range.head != :tuple
        error("The first argument must be the range for the kernel as an integer (1D) or tuple")
    end

    if call.head != :call
        error("The second argument must be a simple call to a kernel function")
    end

    kernel = call.args[1]
    args = call.args[2:end]

    return quote
        cfg = get_or_init_defaults()
        agent = cfg.agent
        queue = cfg.queue

        call_args = Any[$([esc(a) for a in args]...)]
        kernel_args = prepare_args(call_args)
        arg_types = [ typeof(x) for x in kernel_args ]

        kernel_info = build_kernel(agent, $(esc(kernel)), arg_types)

        karg_memory = allocate_args(agent, kernel_info, kernel_args)

        signal = Signal(value = 1)
        packet = build_dispatch($(esc(range)), kernel_info, karg_memory, signal)

        # Enqueue Dispatch
        #idx = HSA.add_write_index!(queue, UInt64(1))
        idx = HSA.load_write_index(queue)

        queue[idx] = packet

        HSA.store_write_index!(queue, UInt64(idx + 1))

        HSA.store!(queue.doorbell_signal, Int64(idx))

        wait(signal, :(<), 1, wait_state_hint = HSA.HSA_WAIT_STATE_ACTIVE)

        cleanup_args(call_args)

        finalize(karg_memory)
        finalize(signal)
    end
end


text:
using LinearAlgebra, GPUifyLoops, StaticArrays, Random, Test

@static if Base.find_package("CuArrays") !== nothing
    using CuArrays, CUDAnative
end

Random.seed!(4)

"""
Implementation of the 2 x 2 eignvalue and eigenvector solver found in:
  @misc{Borges2017,
    title = {Numerically sane solution of the 2 x 2 real symmetric eigenvalue problem},
    author = {Carlos F. Borges},
    year = {2017},
    url = {https://calhoun.nps.edu/handle/10945/51506},
  }
"""
@inline function cfbeigen(A::SHermitianCompact{2,T,3}) where {T <: Real}
    @inbounds α₁ = A[1, 1]
    @inbounds α₂ = A[2, 2]
    @inbounds β = A[2, 1]

    δ = (α₁ - α₂)/2
    h = hypot(δ, β)

    if h == 0
        c = one(T)
        s = zero(T)
    elseif δ ≥ 0
        c = δ + h
        s = β
    else
        c = β
        s = h - δ
    end

    ρ = 1/hypot(s, c)

    c = ρ*c
    s = ρ*s

    V = @SMatrix [c -s;
                  s  c]
    λ = @SVector [c^2*α₁ + s*(s*α₂ + 2c*β),
                  c^2*α₂ + s*(s*α₁ - 2c*β)]

    return Eigen(λ, V)
end

# initialize with random symmetric matrices
N = 2^25
T = Float64

As = (rand(T, N), rand(T, N), rand(T, N))
Vs = ntuple(_->similar(As[1]), 4)
λs = ntuple(_->similar(As[1]), 2)

function kernel(As, Vs, λs)
    @inbounds @loop for i in (1:size(As[1],1);
                              (blockIdx().x-1)*blockDim().x + threadIdx().x)
        A = SHermitianCompact(@SVector([As[1][i], As[2][i], As[3][i]]))

        F = cfbeigen(A)

        @unroll for j = 1:length(F.vectors)
            Vs[j][i] = F.vectors[j]
        end

        @unroll for j = 1:length(F.values)
            λs[j][i] = F.values[j]
        end
    end
    nothing
end

@launch CPU() kernel(As, Vs, λs)

@static if Base.find_package("CuArrays") !== nothing
    cAs = CuArray.(As)
    cVs = ntuple(_->similar(cAs[1]), length(Vs))
    cλs = ntuple(_->similar(cAs[1]), length(λs))

    threads = 1024
    blocks = cld(size(cAs[1],1), threads)
    @launch(CUDA(), threads=threads, blocks=blocks, kernel(cAs, cVs, cλs))

    for j = 1:length(As)
      @test As[j] ≈ Array(cAs[j])
    end

    for j = 1:length(Vs)
      @test Vs[j] ≈ Array(cVs[j])
    end

    for j = 1:length(λs)
      @test λs[j] ≈ Array(cλs[j])
    end
end


text:
module segment
import Term
import Term: remove_markup, remove_ansi
import ..style: apply_style, MarkupStyle
import ..measure: Measure

export Segment

# ---------------------------------------------------------------------------- #
#                                    SEGMENT                                   #
# ---------------------------------------------------------------------------- #

"""
    Segment

stores one piece of text with all the styles applied to it.
"""
struct Segment
    text::AbstractString   # text with ANSI codes injected
    plain::AbstractString  # plain text with no style
    measure::Measure       # measure of plain text
end

# ------------------------------- constructors ------------------------------- #
"""
    Segment(text::AbstractString)

Construct a Segment out of a string with markup.
"""
function Segment(text::Union{Segment,AbstractString})
    if typeof(text) == Segment
        return text
    end
    plain = remove_ansi(remove_markup(text))

    # len(x) = (length ∘ remove_ansi ∘ remove_markup)(x)
    # @info "Creating segment" len(text) len(apply_style(text)) len(remove_markup(apply_style(text)))
    return Segment(remove_markup(apply_style(text)), plain, Measure(plain))
end

"""
    Segment(text::Union{Segment, AbstractString}, markup::AbstractString)

Construct a Segment out of a plain string and a markup string with style info
"""
function Segment(text::Union{Segment,AbstractString}, markup::Union{Nothing,AbstractString})
    return isnothing(markup) ? Segment(text) : Segment("[$markup]" * text)
end

"""
    Segment(text::Union{Segment, AbstractString}, style::MarkupStyle)

Construct a Segment out of a plain string and a MarkupStyle object.
"""
function Segment(text::Union{Segment,AbstractString}, style::Union{Nothing,MarkupStyle})
    if isnothing(style)
        Segment(text, text, Measure(text))
    else
        Segment(apply_style(text, style), text, Measure(text))
    end
end

Segment(text::Union{AbstractString,Segment}, null::Nothing) = Segment(text)

# --------------------------------- printing --------------------------------- #
"""print styled in stdout, info otherwise"""
function Base.show(io::IO, seg::Segment)
    if io == stdout
        print(io, seg.text)

    else
        print(io, "Segment \e[2m(size: $(seg.measure))\e[0m")
    end
end

# ---------------------------------------------------------------------------- #
#                                    LAYOUT                                    #
# ---------------------------------------------------------------------------- #

"""
    Term.fillin(segments::Vector{Segment})::Vector{Segment}

Ensure that for each segment the text has the same width
"""
function Term.fillin(segments::Vector{Segment})::Vector{Segment}
    widths = [seg.measure.w for seg in segments]
    w = max(widths...)

    filled::Vector{Segment} = []
    for seg in segments
        push!(filled, Segment(seg.text * " "^(w-seg.measure.w)))
    end
    return filled
end

# -------------------------------- concatenate ------------------------------- #
"""
concatenate strings and segments
"""
Base.:*(seg::Segment, str::AbstractString) = Segment(seg.text * str)
Base.:*(str::AbstractString, seg::Segment) = Segment(str * seg.text)
Base.:*(seg1::Segment, seg2::Segment) = Segment(seg1.text * seg2.text)

end


text:
input = joinpath(@__DIR__, "input")
entries = readlines(input)

req = ["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"]

validators = Dict(
  "byr" => v -> begin byr = parse(Int, v); 1920 <= byr && byr <= 2002 end,
  "iyr" => v -> begin iyr = parse(Int, v); 2010 <= iyr && iyr <= 2020 end,
  "eyr" => v -> begin eyr = parse(Int, v); 2020 <= eyr && eyr <= 2030 end,
  "hgt" => v -> occursin(r"^(1[5-8][0-9]|19[0-3])(cm)$|^(59|6[0-9]|7[0-6])(in)$", v),
  "hcl" => v -> occursin(r"^#([0-9]|[a-f]){6}$", v),
  "ecl" => v -> (v in ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]),
  "pid" => v -> occursin(r"^[0-9]{9}$", v),
  "cid" => v -> true
)

function completed(passport)
    ks = keys(passport)
    issubset(req, ks)
end

function validated(passport)
    ks = collect(keys(passport))
    !(false in map(k -> validators[k](passport[k]), ks))
end

function col(entries)
    acc = []
    curr = []
    for entry in entries
        if entry == ""
            push!(acc, Dict(curr))
            curr = []
        else
            kvs = split(entry, " ")
            attrs = (split(ks, ":") for ks in kvs)
            append!(curr, attrs)
        end
    end

    if !isempty(curr)
        push!(acc, Dict(curr))
    end

    return acc
end

passports = col(entries)
presentable = filter(p -> completed(p), passports)
valid = filter(p -> validated(p), presentable)

p1 = length(presentable)
p2 = length(valid)

@assert(p1 == 213)
@assert(p2 == 147)

println("-----------------------------------------------------------------------")
println("passport processing -- part one :: $p1")
println("passport processing -- part two :: $p2")
println("-----------------------------------------------------------------------")


text:
abstract type Spec end

import Base: @kwdef

# Dublicate present in utlis.jl
function _hex2bytes(x::String)

    normalized = join(split(x, " "), "")

    N = length(normalized)

    if mod(N, 2) != 0
        normalized = "0" * normalized
    end

    return hex2bytes(normalized)
end

function bytes2bits(x::Vector{UInt8})
    bv = BitVector(u << -i % Bool for u in x for i in 7:-1:0)
    return bv
end

hex2bits(x::String) = bytes2bits(_hex2bytes(x))

_parse_bits(x::String, N::Int) = hex2bits(x)[end - N + 1:end]
_parse_bits(x::BitVector, m::Int) = x
_parse_bits(x::Nothing, m::Int) = x

_parse_int(x::String) = parse(BigInt, x, base=16)
_parse_int(x::Integer) = BigInt(x)
_parse_int(::Nothing) = nothing
_parse_int(x::Vector{UInt8}) = octet2int(x)

struct ECP <: Spec
    p::BigInt
    n::Union{BigInt, Nothing} # = nothing
    a::BigInt # = -3
    b::BigInt
    Gx::Union{BigInt, Nothing} # = nothing
    Gy::Union{BigInt, Nothing} # = nothing

    function ECP(p, n, a, b, Gx, Gy)

        _a = mod(_parse_int(a), p) # taking mod as conventually a=-3

        _b = _parse_int(b)
        _Gx = _parse_int(Gx)
        _Gy = _parse_int(Gy)

        return new(p, n, _a, _b, _Gx, _Gy)
    end

end

function ECP(; p, n::Union{BigInt, Nothing} = nothing, a = -3, b, h = 1, G=nothing, Gx=nothing, Gy=nothing)

    if !isnothing(G)
        ecp = ECP(p, n, a, b, nothing, nothing)
        (Gx, Gy) = point(G, ecp) # Note that the method is defined later in conversions.jl
    end

    return ECP(p, n, a, b, Gx, Gy)
end

order(curve::ECP) = curve.n
generator(curve::ECP) = (curve.Gx, curve.Gy)

modulus(curve::ECP) = curve.p

bitlength(curve::ECP) = bitlength(modulus(curve))

# I could always add a field for equation to be used

Base.:(==)(x::ECP, y::ECP) = x.p == y.p && x.n == y.n && x.a == y.a && x.b == y.b && x.Gx == y.Gx && x.Gy == y.Gy

a(curve::ECP) = curve.a
b(curve::ECP) = curve.b

abstract type BinaryBasis end

struct PB <: BinaryBasis
    f::Vector{Int}
    PB(f::Vector{Int}) = new(f)
    PB(f::BitVector) = PB(((f.len - 1):-1:0)[f])
end

PB(x::Vector{UInt8}, N::Int) = PB(octet2bits(x, N + 1))

bitlength(x::PB) = maximum(x.f) #- 1

struct GNB <: BinaryBasis
    m::Int
    T::Int
end

bitlength(x::GNB) = x.m

struct EC2N{B<:BinaryBasis} <: Spec
    basis::B
    n::Union{BigInt, Nothing}
    a::BitVector
    b::BitVector
    Gx::Union{BitVector, Nothing}
    Gy::Union{BitVector, Nothing}

    function EC2N(basis::B, n::BigInt, a::BitVector, b::BitVector, Gx::BitVector, Gy::BitVector) where B <: BinaryBasis
        @assert bitlength(basis) == length(a) == length(b) == length(Gx) == length(Gy)

        new{B}(basis, n, a, b, Gx, Gy)
    end

    function EC2N(basis::B, n::Union{BigInt, Nothing}, a::BitVector, b::BitVector) where B <: BinaryBasis
        @assert bitlength(basis) == length(a) == length(b)
        new{B}(basis, n, a, b, nothing, nothing)
    end

    function EC2N(basis::B, a::BitVector, b::BitVector) where B <: BinaryBasis
        @assert bitlength(basis) == length(a) == length(b)
        new{B}(basis, nothing, a, b, nothing, nothing)
    end
end

_parse_bits(x::BitVector, basis::BinaryBasis) = x
_parse_bits(x::String, basis::BinaryBasis) = _parse_bits(x, bitlength(basis))

function _int2bits_gnb(x::Int, m::Int)
    if x == 0
        return BitVector(0 for i in 1:m)
    elseif x == 1
        return BitVector(1 for i in 1:m)
    else
        error("Conversion of $x not possible")
    end
end

_parse_bits(x::Int, basis::GNB) = _int2bits_gnb(x, bitlength(basis))

bitlength(x::EC2N) = bitlength(x.basis)

function _int2bits_pb(x::Int, m::Int)
    if x == 0
        return BitVector(0 for i in 1:m)
    elseif x == 1
        return BitVector(((0 for i in 1:m-1)..., 1))
    else
        error("Conversion of $x not possible")
    end
end

_parse_bits(x::Int, basis::PB) = _int2bits_pb(x, bitlength(basis))
_parse_bits(x::Vector{UInt8}, basis::BinaryBasis) = octet2bits(x, bitlength(basis))

# Basis could be nonoptional argument here
function EC2N(basis::BinaryBasis; n=nothing, a, b, h=nothing, G=nothing, Gx=nothing, Gy=nothing)

    #_n = convert(BigInt, n)
    _n = _parse_int(n)
    _a = _parse_bits(a, basis)
    _b = _parse_bits(b, basis)

    if !isnothing(G)
        sp = EC2N(basis, _a, _b)
        _Gx, _Gy = point(G, sp) # The method is defined in conversions.jl
    else
        _Gx = _parse_bits(Gx, bitlength(basis))
        _Gy = _parse_bits(Gy, bitlength(basis))
    end

    if isnothing(_Gx) && isnothing(_Gy)
        return EC2N(basis, _n, _a, _b)
    elseif !isnothing(_Gx) && !isnothing(_Gy)
        return EC2N(basis, _n, _a, _b, _Gx, _Gy)
    else
        error("Incompatible input.")
    end
end

function EC2N{PB}(; f, kwargs...)
    basis = PB(f)
    return EC2N(basis; kwargs...)
end

function EC2N{GNB}(; m, T, kwargs...)
    basis = GNB(m, T)
    return EC2N(basis; kwargs...)
end

function EC2N(; kwargs...)

    if :f in kwargs
        return EC2N{PB}(kwargs...)
    elseif :m in kwargs && :T in kwargs
        return EC2N{GNB}(kwargs...)
    else
        error("Dispatch can not be infered from passed keyword arguments")
    end
end

order(curve::EC2N) = curve.n
generator(curve::EC2N) = (curve.Gx, curve.Gy)

a(curve::EC2N) = curve.a
b(curve::EC2N) = curve.b

struct Koblitz{B<:BinaryBasis} <: Spec
    bec::EC2N{B}

    function Koblitz{B}(; a, kwargs...) where B <: BinaryBasis

        b = 1
        @assert a in [0, 1]

        bec =  EC2N{B}(; a, b, kwargs...)

        return new{B}(bec)
    end

    function Koblitz(; a, kwargs...)

        b = 1
        @assert a in [0, 1]

        bec =  EC2N(; a, b, kwargs...)

        B = typeof(bec.basis)

        return new{B}(bec)
    end
end

order(curve::Koblitz) = order(curve.bec)
generator(curve::Koblitz) = generator(curve.bec)

a(curve::Koblitz) = a(curve.bec)
b(curve::Koblitz) = b(curve.bec)

### MODP spec

tobint(x::String) = parse(BigInt, x, base=16)
tobint(x::BigInt) = x
tobint(x::Integer) = BigInt(x)
tobint(x::Nothing) = nothing

struct MODP <: Spec
    p::BigInt
    g::Union{BigInt, Nothing}
    q::Union{BigInt, Nothing}

    MODP(p, g, q) = new(tobint(p), tobint(g), tobint(q))
    MODP(p, g) = MODP(p, g, nothing)

    MODP(p::BigInt; g=nothing, q=nothing) = new(p, g, q)
    MODP(;p::BigInt, g=nothing, q=nothing) = new(p, g, q)
end

generator(spec::MODP) = spec.g
modulus(spec::MODP) = spec.p
order(spec::MODP) = spec.q


text:
mutable struct SRC0 <: SRCType
end
function (obj::SRC0)(
              xy::AbstractVector,
              f::AbstractArray)
  f[:] = 0.0
  return nothing
end

mutable struct SRC1 <: SRCType
end
function (obj::SRC1,
              xy::AbstractVector,
              f::AbstractArray)
  f[:] = 1.0
  return nothing
end

#
# a scalar diffusion coefficient, 10
#
mutable struct SrcExpTrigPoly0thDiffn <: SRCType
end
function (obj::SrcExpTrigPoly0thDiffn)(
              xy::AbstractVector,
              src::AbstractArray)
  k = 2.0
  ss = sin(2*k*pi*xy[1])*sin(2*k*pi*xy[2])
  cs = cos(2*k*pi*xy[1])*sin(2*k*pi*xy[2])
  sc = sin(2*k*pi*xy[1])*cos(2*k*pi*xy[2])
  ex = exp(xy[1] + xy[2])
  src[:] = (2. - 8*k*k*pi*pi) * ss + 4*k*pi*(cs + sc)
  src[:] *= -10.*ex
  return nothing
end
#
# a scalar diffusion coefficient, 10
#
mutable struct SrcTrigPoly0thDiffn <: SRCType
end
function (obj::SrcTrigPoly0thDiffn)(
              xy::AbstractVector,
              src::AbstractArray)
  n = 2.0
  factor = 8.0*n*n*pi*pi
  ss = sin(2*n*pi*xy[1])*sin(2*n*pi*xy[2])
  src[:] = 0.0
  for i = 1:length(src)
    src[i] = factor*ss * 1.0e1
  end
  return nothing
end
#
# a 2nd order diffusion coefficient tensor
#
mutable struct SrcTrigPoly2ndDiffn <: SRCType
end
function (obj::SrcTrigPoly2ndDiffn)(
              xy::AbstractVector,
              src::AbstractArray)
  n = 2.0
  sc = sin(2*n*pi*xy[1])*cos(2*n*pi*xy[2])
  cs = cos(2*n*pi*xy[1])*sin(2*n*pi*xy[2])
  ss = sin(2*n*pi*xy[1])*sin(2*n*pi*xy[2])
  cc = cos(2*n*pi*xy[1])*cos(2*n*pi*xy[2])
  src[:] = 0.0
  for i = 1:length(src)
    src[i] += 6*n*pi*xy[1]*cs
    src[i] += 6*n*pi*xy[2]*sc
    src[i] += 4*n*n*pi*pi*(-2.0 - xy[1]*xy[1] - xy[2]*xy[2])*ss
    src[i] += 8*n*n*pi*pi*xy[1]*xy[2]*cc
    src[i] = -src[i]
  end
  return nothing
end

#
# a 6th order diffusion coefficient tensor
#
mutable struct SrcTrigPoly6thDiffn <: SRCType
end
function (obj::SrcTrigPoly6thDiffn)(
              xy::AbstractVector,
              src::AbstractArray)
  n = 2.0
  sc = sin(2*n*pi*xy[1])*cos(2*n*pi*xy[2])
  cs = cos(2*n*pi*xy[1])*sin(2*n*pi*xy[2])
  ss = sin(2*n*pi*xy[1])*sin(2*n*pi*xy[2])
  cc = cos(2*n*pi*xy[1])*cos(2*n*pi*xy[2])
  src[:] = 0.0
  for i = 1:length(src)
    src[i] += 2*n*pi*(6*xy[1]^5 + 3*xy[1]^3*xy[2]^2)*cs
    src[i] += 2*n*pi*(6*xy[2]^5 + 3*xy[1]^2*xy[2]^3)*sc
    src[i] += 4*n*n*pi*pi*(-2.0 - xy[1]^6 - xy[2]^6)*ss
    src[i] += 8*n*n*pi*pi*xy[1]^3*xy[2]^3*cc
    src[i] = -src[i]
  end
  return nothing
end

global const SRCDict = Dict{String, SRCType}(
 "SRC0" => SRC0(),
 "SRC1" => SRC1(),
 "SrcTrigPoly0thDiffn" => SrcTrigPoly0thDiffn(),
 "SrcExpTrigPoly0thDiffn" => SrcExpTrigPoly0thDiffn(),
 "SrcTrigPoly2ndDiffn" => SrcTrigPoly2ndDiffn(),
 "SrcTrigPoly6thDiffn" => SrcTrigPoly6thDiffn(),
)

function getSrcFuntors(mesh::AbstractMesh,
                       sbp::AbstractOperator,
                       eqn::EllipticData,
                       opts)
  eqn.src_func = SRCDict[opts["SRC_name"]]
  return nothing
end

function calcSource(mesh::AbstractMesh{Tmsh},
                    sbp::AbstractOperator,
                    eqn::EllipticData{Tsol, Tres, Tdim},
                    srcFunc::SRCType,
                    src::AbstractArray{Tres, 3}) where {Tmsh, Tsol, Tres, Tdim}

  numElems = mesh.numEl
  numNodesPerElement = mesh.numNodesPerElement
  # xy = Array{Tmsh}(Tdim)
  src_tmp = Array{eltype(src)}(size(src, 1))
  for elem = 1:numElems
    for n = 1:numNodesPerElement
      xy = sview(mesh.coords, :, n, elem)

      srcFunc(xy, sview(src, :, n, elem))
    end
  end
end


text:
module BufferedArrays

import Base: size, getindex, setindex!, unsafe_convert, @pure
using Base.Libc

export Buffer, ArrayOffsetBuffer, BArray, BVector, BMatrix

"""
    Buffer(bytes)

Allocate `bytes` bytes of memory, and return a `Buffer` object, which will be
garbage collected when (or after) no references to the `Buffer` remain.
"""
type Buffer
    ptr::Ptr{Void} # Ideally this would be const

    function Buffer(bytes::Int)
        out = new(malloc(bytes))
        finalizer(out, _free_buffer)
        return out
    end
end

"""
`_free_buffer(buffer)`

An internal function for freeing buffers. Users should never call this function -
it will be called automatically by the GC (and will result in errors if called
twice).
"""
_free_buffer(buf::Buffer) = free(buf.ptr)

"""
    ArrayOffsetBuffer{T,N}

A specialized buffer suitable for implementing the `BArray` type. See also `Buffer`.
"""
type ArrayOffsetBuffer{T,N}
    ptr::Ptr{Void} # Ideally this would be const

    function ArrayOffsetBuffer(bytes::Int)
        out = new(malloc(bytes+offset_bytes(T,N)) + offset_size(T,N))
        finalizer(out, _free_buffer)
        return out
    end
end

_free_buffer{T,N}(buf::ArrayOffsetBuffer{T,N}) = free(buf.ptr - offset_size(T,N))

# Calculate an offset which is in multiples of 4 Ints. This *should* result in
# good memory alignment for the first piece of data in the array (assuming
# malloc returns an well-aligned address)
# TODO do something faster w.r.t. one-based indexing
@pure function offset_size{T}(::Type{T}, N::Int)
    return offset_bytes(T,N) # - sizeof(T) # future correction for 1-based indexing
end

# The total size of the offset padding (no 1-based indexing correction, if exists)
@pure function offset_bytes{T}(::Type{T}, N::Int)
    @assert N >= 0
    tmp = N
    while tmp % 4 != 0
        tmp += 1
    end
    return sizeof(Int)*tmp
end

"""
    BArray{T}(dims)
    BArray{T}(dims...)

Constructs a N-dimensional array backed by a native Julia memory buffer, where
`dims` is an `N`-tuple of integers giving the array size and `T` is the element
type. Note: the returned array has unintialized data.
"""
immutable BArray{T,N} <: DenseArray{T,N}
    buf::ArrayOffsetBuffer{T,N}

    function BArray(sizes::NTuple{N,Int})
        bytes = sizeof(T)*prod(sizes)
        out = new(ArrayOffsetBuffer{T,N}(bytes))
        unsafe_store!(unsafe_convert(Ptr{NTuple{N,Int}}, out.buf.ptr - offset_size(T,N)), sizes)
        return out
    end
end

@inline (::Type{BA}){BA<:BArray}(x::Int...) = BA(x)
@inline (::Type{BArray{T}}){T,N}(x::NTuple{N,Int}) = BArray{T,N}(x)

Base.linearindexing{BA<:BArray}(::Union{BA,Type{BA}}) = Base.LinearFast()

size{T,N}(a::BArray{T,N}) = unsafe_load(unsafe_convert(Ptr{NTuple{N,Int}}, a.buf.ptr - offset_size(T,N)))

@inline function getindex{T,N}(a::BArray{T,N}, i::Int)
    # @boundscheck bla
    unsafe_load(unsafe_convert(Ptr{T}, a.buf.ptr), i) # oh dear! This is broken since unsafe_load assumes 1-based indexing
end

@inline function setindex!{T,N}(a::BArray{T,N}, x::T, i::Int)
    # @boundscheck bla
    unsafe_store!(unsafe_convert(Ptr{T}, a.buf.ptr), x, i)  # oh dear! This is broken since unsafe_store! assumes 1-based indexing
end

@inline unsafe_convert{T}(::Type{Ptr{T}}, a::BArray{T}) = unsafe_convert(Ptr{T}, a.buf.ptr)

"""
   BVector{T}(len)

Create a vector of length `len` and element type `T` backed by a native Julia
memory buffer. Note: the returned array has unintialized data.
"""
typealias BVector{T} BArray{T,1}

"""
   BMatrix{T}(size)
   BMatrix{T}(m, n)

Create a matrix of dimensions `size = (m,n)` and element type `T` backed by a
native Julia  memory buffer. Note: the returned array has unintialized data.
"""
typealias BMatrix{T} BArray{T,2}

end # module


text:
struct IPNewton{F,Tμ<:Union{Symbol,Number}} <: IPOptimizer{F}
    linesearch!::F
    μ0::Tμ      # Initial value for the barrier penalty coefficient μ
    show_linesearch::Bool
    # TODO: μ0, and show_linesearch were originally in options
end

Base.summary(::IPNewton) = "Interior Point Newton"

promote_objtype(method::IPNewton, x, autodiff::Symbol, inplace::Bool, f) = TwiceDifferentiable(f, x, real(zero(eltype(x))); autodiff = autodiff)
promote_objtype(method::IPNewton, x, autodiff::Symbol, inplace::Bool, f, g) = TwiceDifferentiable(f, g, x, real(zero(eltype(x))); inplace = inplace, autodiff = autodiff)
promote_objtype(method::IPNewton, x, autodiff::Symbol, inplace::Bool, f, g, h) = TwiceDifferentiable(f, g, h, x, real(zero(eltype(x))); inplace = inplace)

# TODO: Add support for InitialGuess from LineSearches
"""
# Interior-point Newton
## Constructor
```jl
IPNewton(; linesearch::Function = Optim.backtrack_constrained_grad,
         μ0::Union{Symbol,Number} = :auto,
         show_linesearch::Bool = false)
```

The initial barrier penalty coefficient `μ0` can be chosen as a number, or set
to `:auto` to let the algorithm decide its value, see `initialize_μ_λ!`.

*Note*: For constrained optimization problems, we recommend
always enabling `allow_f_increases` and `successive_f_tol` in the options passed to `optimize`.
The default is set to `Optim.Options(allow_f_increases = true, successive_f_tol = 2)`.

As of February 2018, the line search algorithm is specialised for constrained
interior-point methods. In future we hope to support more algorithms from
`LineSearches.jl`.

## Description
The `IPNewton` method implements an interior-point primal-dual Newton algorithm for solving
nonlinear, constrained optimization problems. See Nocedal and Wright (Ch. 19, 2006) for a discussion of
interior-point methods for constrained optimization.

## References
The algorithm was [originally written by Tim Holy](https://github.com/JuliaNLSolvers/Optim.jl/pull/303) (@timholy, tim.holy@gmail.com).

 - J Nocedal, SJ Wright (2006), Numerical optimization, second edition. Springer.
 - A Wächter, LT Biegler (2006), On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Mathematical Programming 106 (1), 25-57.
"""
IPNewton(; linesearch::Function = backtrack_constrained_grad,
         μ0::Union{Symbol,Number} = :auto,
         show_linesearch::Bool = false) =
  IPNewton(linesearch, μ0, show_linesearch)

mutable struct IPNewtonState{T,Tx} <: AbstractBarrierState
    x::Tx
    f_x::T
    x_previous::Tx
    g::Tx
    f_x_previous::T
    H::Matrix{T}          # Hessian of the Lagrangian?
    HP # TODO: remove HP? It's not used
    Hd::Vector{Int8}  # TODO: remove Hd? It's not used
    s::Tx  # step for x
    # Barrier penalty fields
    μ::T                  # coefficient of the barrier penalty
    μnext::T              # μ for the next iteration
    L::T                  # value of the Lagrangian (objective + barrier + equality)
    L_previous::T
    bstate::BarrierStateVars{T}   # value of slack and λ variables (current "position")
    bgrad::BarrierStateVars{T}    # gradient of slack and λ variables at current "position"
    bstep::BarrierStateVars{T}    # search direction for slack and λ
    constr_c::Vector{T}   # value of the user-supplied constraints at x
    constr_J::Matrix{T}   # value of the user-supplied Jacobian at x
    ev::T                 # equality violation, ∑_i λ_Ei (c*_i - c_i)
    Optim.@add_linesearch_fields() # x_ls and alpha
    b_ls::BarrierLineSearchGrad{T}
    gtilde::Tx
    Htilde               # Positive Cholesky factorization of H from PositiveFactorizations.jl
end

# TODO: Do we need this convert thing? (It seems to be used with `show(IPNewtonState)`)
function Base.convert(::Type{IPNewtonState{T,Tx}}, state::IPNewtonState{S, Sx}) where {T,Tx,S,Sx}
    IPNewtonState(convert(Tx, state.x),
                  T(state.f_x),
                  convert(Tx, state.x_previous),
                  convert(Tx, state.g),
                  T(state.f_x_previous),
                  convert(Matrix{T}, state.H),
                  state.HP,
                  state.Hd,
                  convert(Tx, state.s),
                  T(state.μ),
                  T(state.μnext),
                  T(state.L),
                  T(state.L_previous),
                  convert(BarrierStateVars{T}, state.bstate),
                  convert(BarrierStateVars{T}, state.bgrad),
                  convert(BarrierStateVars{T}, state.bstep),
                  convert(Vector{T}, state.constr_c),
                  convert(Matrix{T}, state.constr_J),
                  T(state.ev),
                  convert(Tx, state.x_ls),
                  T(state.alpha),
                  convert(BarrierLineSearchGrad{T}, state.b_ls),
                  convert(Tx, state.gtilde),
                  state.Htilde,
                  )
end

function initial_state(method::IPNewton, options, d::TwiceDifferentiable, constraints::TwiceDifferentiableConstraints, initial_x::AbstractArray{T}) where T
    # Check feasibility of the initial state
    mc = nconstraints(constraints)
    constr_c = Array{T}(undef, mc)
    # TODO: When we change to `value!` from NLSolversBase instead of c!
    # we can also update `initial_convergence` for ConstrainedOptimizer in interior.jl
    constraints.c!(constr_c, initial_x)
    if !isinterior(constraints, initial_x, constr_c)
        @warn("Initial guess is not an interior point")
        Base.show_backtrace(stderr, backtrace())
        println(stderr)
    end
    # Allocate fields for the objective function
    n = length(initial_x)
    g = similar(initial_x)
    s = similar(initial_x)
    f_x_previous = NaN
    f_x, g_x = value_gradient!(d, initial_x)
    g .= g_x # needs to be a separate copy of g_x
    H = Matrix{T}(undef, n, n)
    Hd = Vector{Int8}(undef, n)
    hessian!(d, initial_x)
    copyto!(H, hessian(d))

    # More constraints
    constr_J = Array{T}(undef, mc, n)
    gtilde = copy(g)
    constraints.jacobian!(constr_J, initial_x)
    μ = T(1)
    bstate = BarrierStateVars(constraints.bounds, initial_x, constr_c)
    bgrad = copy(bstate)
    bstep = copy(bstate)
    # b_ls = BarrierLineSearch(similar(constr_c), similar(bstate))
    b_ls = BarrierLineSearchGrad(copy(constr_c), copy(constr_J), copy(bstate), copy(bstate))

    state = IPNewtonState(
        copy(initial_x), # Maintain current state in state.x
        f_x, # Store current f in state.f_x
        copy(initial_x), # Maintain previous state in state.x_previous
        g, # Store current gradient in state.g (TODO: includes Lagrangian calculation?)
        T(NaN), # Store previous f in state.f_x_previous
        H,
        0,    # will be replaced
        Hd,
        similar(initial_x), # Maintain current x-search direction in state.s
        μ,
        μ,
        T(NaN),
        T(NaN),
        bstate,
        bgrad,
        bstep,
        constr_c,
        constr_J,
        T(NaN),
        Optim.@initial_linesearch()..., # Maintain a cache for line search results in state.lsr
        b_ls,
        gtilde,
        0)

    Hinfo = (state.H, hessianI(initial_x, constraints, 1 ./ bstate.slack_c, 1))
    initialize_μ_λ!(state, constraints.bounds, Hinfo, method.μ0)
    update_fg!(d, constraints, state, method)
    update_h!(d, constraints, state, method)

    state
end

function update_fg!(d, constraints::TwiceDifferentiableConstraints, state, method::IPNewton)
    state.f_x, state.L, state.ev = lagrangian_fg!(state.g, state.bgrad, d, constraints.bounds, state.x, state.constr_c, state.constr_J, state.bstate, state.μ)
    update_gtilde!(d, constraints, state, method)
end

function update_gtilde!(d, constraints::TwiceDifferentiableConstraints, state, method::IPNewton)
    # Calculate the modified x-gradient for the block-eliminated problem
    # gtilde is the gradient for the affine-scaling problem, i.e.,
    # with μ=0, used in the adaptive setting of μ. Once we calculate μ we'll correct it
    gtilde, bstate, bgrad = state.gtilde, state.bstate, state.bgrad
    bounds = constraints.bounds
    copyto!(gtilde, state.g)
    JIc = view(state.constr_J, bounds.ineqc, :)
    if !isempty(JIc)
        Hssc = Diagonal(bstate.λc./bstate.slack_c)
        # TODO: Can we use broadcasting / dot-notation here and eliminate gc?
        gc = JIc'*(Diagonal(bounds.σc) * (bstate.λc - Hssc*bgrad.λc))  # NOT bgrad.slack_c
        gtilde .+= gc
    end
    for (i,j) in enumerate(bounds.ineqx)
        gxi = bounds.σx[i]*(bstate.λx[i] -  bgrad.λx[i]*bstate.λx[i]/bstate.slack_x[i])
        gtilde[j] += gxi
    end
    state
end

function update_h!(d, constraints::TwiceDifferentiableConstraints, state, method::IPNewton)
    x, μ, Hxx, J = state.x, state.μ, state.H, state.constr_J
    bstate, bgrad, bounds = state.bstate, state.bgrad, constraints.bounds
    m, n = size(J, 1), size(J, 2)

    hessian!(d, state.x)  # objective's Hessian
    copyto!(Hxx, hessian(d))  # objective's Hessian
    # accumulate the constraint second derivatives
    λ = userλ(bstate.λc, constraints)
    λ[bounds.eqc] = -bstate.λcE  # the negative sign is from the Hessian
    # Important! We are assuming that constraints.h! adds the hessian of the
    # non-objective Lagrangian terms to the existing objective Hessian Hxx.
    # This follows the approach by the CUTEst interface
    constraints.h!(Hxx, x, λ)

    # Add the Jacobian terms (JI'*Hss*JI)
    JIc = view(J, bounds.ineqc, :)
    Hssc = Diagonal(bstate.λc./bstate.slack_c)
    HJ = JIc'*Hssc*JIc
    for j = 1:n, i = 1:n
        Hxx[i,j] += HJ[i,j]
    end
    # Add the variable inequalities portions of J'*Hssx*J
    for (i,j) in enumerate(bounds.ineqx)
        Hxx[j,j] += bstate.λx[i]/bstate.slack_x[i]
    end
    state.Htilde = cholesky(Positive, Hxx, Val{true})

    state
end

function update_state!(d, constraints::TwiceDifferentiableConstraints, state::IPNewtonState{T}, method::IPNewton, options) where T
    state.f_x_previous, state.L_previous = state.f_x, state.L
    bstate, bstep, bounds = state.bstate, state.bstep, constraints.bounds
    qp = solve_step!(state, constraints, options, method.show_linesearch)
    # If a step α=1 will not change any of the parameters, we can quit now.
    # This prevents a futile linesearch.
    if is_smaller_eps(state.x, state.s) &&
        is_smaller_eps(bstate.slack_x, bstep.slack_x) &&
        is_smaller_eps(bstate.slack_c, bstep.slack_c) &&
        is_smaller_eps(bstate.λx, bstep.λx) &&
        is_smaller_eps(bstate.λc, bstep.λc)
        return false
    end

    # Estimate αmax, the upper bound on distance of movement along the search line
    αmax = convert(eltype(bstate), Inf)
    αmax = estimate_maxstep(αmax, bstate.slack_x, bstep.slack_x)
    αmax = estimate_maxstep(αmax, bstate.slack_c, bstep.slack_c)
    αmax = estimate_maxstep(αmax, bstate.λx, bstep.λx)
    αmax = estimate_maxstep(αmax, bstate.λc, bstep.λc)

    # Determine the actual distance of movement along the search line
    ϕ = linesearch_anon(d, constraints, state, method)
    # TODO: This only works for method.linesearch = backtracking_constrained_grad
    # TODO: How are we meant to implement backtracking_constrained?.
    #       It requires both an alpha and an alphaI (αmax and αImax) ...
    state.alpha =
        method.linesearch!(ϕ, T(1), αmax, qp; show_linesearch=method.show_linesearch)

    # Maintain a record of previous position
    copyto!(state.x_previous, state.x)

    # Update current position # x = x + alpha * s
    ls_update!(state.x, state.x, state.s, state.alpha)
    ls_update!(bstate, bstate, bstep, state.alpha)

    # Ensure that the primal-dual approach does not deviate too much from primal
    # (See Waechter & Biegler 2006, eq. 16)
    μ = state.μ
    for i = 1:length(bstate.slack_x)
        p = μ / bstate.slack_x[i]
        bstate.λx[i] = max(min(bstate.λx[i], 1e10*p), 1e-10*p)
    end
    for i = 1:length(bstate.slack_c)
        p = μ / bstate.slack_c[i]
        bstate.λc[i] = max(min(bstate.λc[i], 1e10*p), 1e-10*p)
    end
    state.μ = state.μnext

    # Evaluate the constraints at the new position
    constraints.c!(state.constr_c, state.x)
    constraints.jacobian!(state.constr_J, state.x)
    state.ev == equality_violation(constraints, state)

    false
end

function solve_step!(state::IPNewtonState, constraints, options, show_linesearch::Bool = false)
    x, s, μ, bounds = state.x, state.s, state.μ, constraints.bounds
    bstate, bstep, bgrad = state.bstate, state.bstep, state.bgrad
    J, Htilde = state.constr_J, state.Htilde
    # Solve the Newton step
    JE = jacobianE(state, bounds)
    gE = [bgrad.λxE;
          bgrad.λcE]
    M = JE*(Htilde \ JE')
    MF = cholesky(Positive, M, Val{true})
    # These are a solution to the affine-scaling problem (with μ=0)
    ΔλE0 = MF \ (gE + JE * (Htilde \ state.gtilde))
    Δx0 = Htilde \ (JE'*ΔλE0 - state.gtilde)
    # Check that the solution to the linear equations represents an improvement
    Hpstepx, HstepλE = Matrix(Htilde)*Δx0 - JE'*ΔλE0, -JE*Δx0  # TODO: don't use full here
    # TODO: How to handle show_linesearch?
    # This was originally in options.show_linesearch, but I removed it as none of the other Optim algorithms have it there.
    # We should move show_linesearch back to options when we refactor
    # LineSearches to work on the function ϕ(α)
    if show_linesearch
        println("|gx| = $(norm(state.gtilde)), |Hstepx + gx| = $(norm(Hpstepx+state.gtilde))")
        println("|gE| = $(norm(gE)), |HstepλE + gE| = $(norm(HstepλE+gE))")
    end
    if norm(gE) + norm(state.gtilde) < max(norm(HstepλE + gE),
                                           norm(Hpstepx  + state.gtilde))
        # Precision problems gave us a worse solution than the one we started with, abort
        fill!(s, 0)
        fill!(bstep, 0)
        return state
    end
    # Set μ (see the predictor strategy in Nodecal & Wright, 2nd ed., section 19.3)
    solve_slack!(bstep, Δx0, bounds, bstate, bgrad, J, zero(state.μ)) # store temporarily in bstep
    αs = convert(eltype(bstate), 1.0)
    αs = estimate_maxstep(αs, bstate.slack_x, bstep.slack_x)
    αs = estimate_maxstep(αs, bstate.slack_c, bstep.slack_c)
    αλ = convert(eltype(bstate), 1.0)
    αλ = estimate_maxstep(αλ, bstate.λx, bstep.λx)
    αλ = estimate_maxstep(αλ, bstate.λc, bstep.λc)
    m = max(1, length(bstate.slack_x) + length(bstate.slack_c))
    μaff = (dot(bstate.slack_x + αs*bstep.slack_x, bstate.λx + αλ*bstep.λx) +
            dot(bstate.slack_c + αs*bstep.slack_c, bstate.λc + αλ*bstep.λc))/m
    μmean = (dot(bstate.slack_x, bstate.λx) + dot(bstate.slack_c, bstate.λc))/m
    # When there's only one constraint, μaff can be exactly zero. So limit the decrease.
    state.μnext = NaNMath.max((μaff/μmean)^3 * μmean, μmean/10)
    μ = state.μ
    # Solve for the *real* step (including μ)
    μsinv = μ * [bounds.σx./bstate.slack_x; bounds.σc./bstate.slack_c]
    gtildeμ = state.gtilde  - jacobianI(state, bounds)' * μsinv
    ΔλE = MF \ (gE + JE * (Htilde \ gtildeμ))
    Δx = Htilde \ (JE'*ΔλE - gtildeμ)
    copyto!(s, Δx)
    k = unpack_vec!(bstep.λxE, ΔλE, 0)
    k = unpack_vec!(bstep.λcE, ΔλE, k)
    k == length(ΔλE) || error("exhausted targets before ΔλE")
    solve_slack!(bstep, Δx, bounds, bstate, bgrad, J, μ)
    # Solve for the quadratic parameters (use the real H, not the posdef H)
    Hstepx, HstepλE  = state.H*Δx - JE'*ΔλE, -JE*Δx
    qp = state.L, slopealpha(state.s, state.g, bstep, bgrad), dot(Δx, Hstepx) + dot(ΔλE, HstepλE)
    qp
end

function solve_slack!(bstep, s, bounds, bstate, bgrad, J, μ)
    # Solve for the slack variable and λI updates
    for (i, j) in enumerate(bounds.ineqx)
        bstep.slack_x[i] = -bgrad.λx[i] + bounds.σx[i]*s[j]
        # bstep.λx[i] = -bgrad.slack_x[i] - μ*bstep.slack_x[i]/bstate.slack_x[i]^2
        # bstep.λx[i] = -bgrad.slack_x[i] - bstate.λx[i]*bstep.slack_x[i]/bstate.slack_x[i]
        bstep.λx[i] = -(-μ/bstate.slack_x[i] + bstate.λx[i]) - bstate.λx[i]*bstep.slack_x[i]/bstate.slack_x[i]
    end
    JIc = view(J, bounds.ineqc, :)
    SigmaJIΔx = Diagonal(bounds.σc)*(JIc*s)
    for i = 1:length(bstep.λc)
        bstep.slack_c[i] = -bgrad.λc[i] + SigmaJIΔx[i]
        # bstep.λc[i] = -bgrad.slack_c[i] - μ*bstep.slack_c[i]/bstate.slack_c[i]^2
        # bstep.λc[i] = -bgrad.slack_c[i] - bstate.λc[i]*bstep.slack_c[i]/bstate.slack_c[i]
        bstep.λc[i] = -(-μ/bstate.slack_c[i] + bstate.λc[i]) - bstate.λc[i]*bstep.slack_c[i]/bstate.slack_c[i]
    end
    bstep
end

function is_smaller_eps(ref, step)
    ise = true
    for (r, s) in zip(ref, step)
        ise &= (s == 0) | (abs(s) < eps(r))
    end
    ise
end

function default_options(method::ConstrainedOptimizer)
    Dict(:allow_f_increases => true, :successive_f_tol => 2)
end

# Utility functions that assist in testing: they return the "full
# Hessian" and "full gradient" for the equation with the slack and λI
# eliminated.
# TODO: should we put these elsewhere?
function Hf(bounds::ConstraintBounds, state)
    JE = jacobianE(state, bounds)
    Hf = [Matrix(state.Htilde) -JE';
          -JE zeros(eltype(JE), size(JE, 1), size(JE, 1))]
end
Hf(constraints, state) = Hf(constraints.bounds, state)
function gf(bounds::ConstraintBounds, state)
    bstate, μ = state.bstate, state.μ
    μsinv = μ * [bounds.σx./bstate.slack_x; bounds.σc./bstate.slack_c]
    gtildeμ = state.gtilde  - jacobianI(state, bounds)' * μsinv
    [gtildeμ; state.bgrad.λxE; state.bgrad.λcE]
end
gf(constraints, state) = gf(constraints.bounds, state)


text:
using LinearAlgebra, FastSolversForWeightedTV, CUDA, Flux, Test, TestImages
CUDA.allowscalar(false)
include("./test_utils.jl")

## 2-D

# flag_gpu = true
flag_gpu = false

# Random data
T = Float64
CT = Complex{T}
step = T(1e-5)
rtol = T(1e-4)
n = (8,8)
y = CT.(T.(TestImages.shepp_logan(n...))); flag_gpu && (y = y |> gpu)
h = (T(1),T(1))

g = gradient_norm(2, 1, n, h; T=CT); flag_gpu && (g = g |> gpu)

# Optimization options
opt = opt_fista(T(1/8); niter=100000, tol_x=nothing, nesterov=true)

# Proxy
λ = T(0.5)*norm(y)^2/g(y)
x = proxy(y, λ, g, opt)

# Gradient test
fun = proxy_objfun(λ, g, opt)
test_grad(fun, y; step=step, rtol=rtol)

# Projection test
ε = T(0.1)*g(y)
x = project(y, ε, g, opt)
@test g(x) ≈ ε rtol=rtol

# Gradient test
fun = proj_objfun(ε, g, opt)
test_grad(fun, y; step=step, rtol=rtol)

## 3-D

# Random data
n = (8,8,8)
h = (T(1),T(1),T(1))
y = T.(TestImages.shepp_logan(n[1:2]...)); flag_gpu && (y = y |> gpu)
y = CT.(repeat(y; outer=(1,1,8)))

g = gradient_norm(2, 1, n, h; T=CT); flag_gpu && (g = g |> gpu)

# Optimization options
opt = opt_fista(T(1/12); niter=100000, tol_x=nothing, nesterov=true)

# Proxy
λ = T(0.5)*norm(y)^2/g(y)
x = proxy(y, λ, g, opt)

# Gradient test
fun = proxy_objfun(λ, g, opt)
test_grad(fun, y; step=step, rtol=rtol)

# Projection test
ε = T(0.1)*g(y)
x = project(y, ε, g, opt)
@test g(x) ≈ ε rtol=rtol

# Gradient test
fun = proj_objfun(ε, g, opt)
test_grad(fun, y; step=step, rtol=rtol)

max_stars_repo_name: grizzuti/FastSolversForWeightedTV

max_issues_repo_name: grizzuti/FastSolversForWeightedTV

max_forks_repo_name: grizzuti/FastSolversForWeightedTV


text:
import CUDA.CUDNN: cudnnSoftmaxForwardAD
using CUDA.CUDNN: cudnnSoftmaxBackward, handle
using AutoGrad: AutoGrad, @primitive1, value

@primitive1((cudnnSoftmaxForwardAD(x; algo, mode, alpha, xDesc, beta, yDesc, y),
             _dy,_y),
            ((x,y,dy,dx) = (value(x),value(_y),value(_dy),similar(x));
             if alpha[] != 1; y = y ./ alpha[]; end;
             cudnnSoftmaxBackward(handle(), algo, mode, alpha, yDesc, y, yDesc, dy, beta, xDesc, dx);
             dx))


text:
# https://www.sciencedirect.com/science/article/pii/S0730725X07001701?via%3Dihub
T2s = Dict( :B7T    => Dict(:gm => 33.2, :wm => 26.8, :ven => 7.4, :csf => 500, :ms => 50, :iron_ring => 21), # csf and ms estimated from 7T scan)
            :B3T    => Dict(:gm => 66.0, :wm => 53.2, :ven => 25), # ven visually estimated from swi_3T scan
            :B1p5T  => Dict(:gm => 84.0, :wm => 66.2))
#https://www.ncbi.nlm.nih.gov/pubmed/18259791
T1 = Dict(  :B7T    => Dict(:gm => 1939, :wm => 1126, :putamen => 1644, :caudate_head => 1684, #=estimated=#:ven => 2000, :ms => 1200),
            :B3T    => Dict(:gm => 1607, :wm => 838, :putamen => 1332, :caudate_head => 1395, #=estimated=#:ven => 1600, :ms => 1000),
            :B1p5T  => Dict(:gm => 1197, :wm => 646, :putamen => 1084, :caudate_head => 1109, #=estimated=#:ven => 1400, :ms => 800))
#https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3375320/
PD = Dict(:gm => 0.71, :wm => 0.63, :ms => 0.68, :ven => 0.9 #=estimated=#)
factor = Dict(  :B7T => Dict(:gm => 1.0, :wm => 0.95, :ven => 0.95, :csf => 0.7, :ms => 0.8, :iron_ring => 0.96), # estimated from 7T scan
                :B3T => Dict(:gm => 1.0, :wm => 0.95, :ven => 0.95)) # estimated from 7T scan

function gettissue(field)
    T2s[field], T1[field], PD
end

function gettissue_all(field)
    T2s[field], T1[field], PD, factor[field]
end

function gettissue_easy(field)
    T2s[field], factor[field]
end


text:
using Documenter
using YaoArrayRegister, YaoBase, YaoBlocks

makedocs(
    modules = [YaoBlocks],
    doctest = true,
    clean = false,
    sitename = "YaoBlocks.jl",
    pages = ["index.md"]
)


text:
module DSP

using FFTW
using LinearAlgebra: mul!, rmul!
using IterTools: subsets

export conv, deconv, filt, filt!, xcorr

include("dspbase.jl")

include("util.jl")
include("unwrap.jl")
include("windows.jl")
include("periodograms.jl")
include("Filters/Filters.jl")
include("lpc.jl")
include("estimation.jl")

using Reexport
@reexport using .Util, .Windows, .Periodograms, .Filters, .LPC, .Unwrap, .Estimation

include("deprecated.jl")
end


text:
using YAML

const info = Dict{String, Any}([
    "title"   => "Doge to the moon"
    "version" => "1.0.5"
])

@swagger """
/doge:
  get:
    description: doge
    responses:
      '200':
        description: Returns a doge.
"""

@swagger """
/bad:
  get:
    description: doge
    responses:
      '200':
        description: Returns a doge.
        bad: I'm causing trouble!!
"""

const all_paths = SwaggerMarkdown.parse_spec()

function get_paths(ps)
    res = Dict{String, Any}()
    for p in ps
        res[p] = all_paths[p]
    end
    return res
end

paths = get_paths(["/doge"])
invalid_paths = get_paths(["/bad"])

function get_optional_fields(version)
    if version == "2.0"
        return Dict{String, Any}(["host" => "localhost", "basePath" => "/doge"])
    end
    if startswith(version, "3.0")
        return Dict{String, Any}(["servers" => [Dict(["url" => "localhost"])]])
    end
end

function get_invalid_optional_fields(version)
    optional_fields = get_optional_fields(version)
    optional_fields["bad"] = "bad"
    return optional_fields
end

function valid_basic(openApi, version)
    spec = build(openApi)
    if version == "2.0"
        @test spec["swagger"] == version
    elseif version == "3.0"
        @test spec["openapi"] == "3.0.0"
    else
        @test spec["openapi"] == version
    end
    @test spec["info"] == info
    @test spec["paths"] == paths
end

function throws(openApi, _)
    @test_throws SwaggerMarkdown.InvalidSwaggerSpecificationException build(openApi)
end

versions = ["2.0", "3.0", "3.0.1"]
cases = Dict([
    "Required Fields excluding 'paths' tests" => Dict([
        "Fields missing" => Dict([
            "info" => Dict{String, Any}(["title" => "Doge to the moon"]),
            "paths" => paths,
            "test" => throws
        ]),
        "Invalid fields" => Dict([
            "info" => Dict{String, Any}(["title" => "Doge to the moon"]),
            "paths" => paths,
            "test" => throws
        ]),
        "All valid" => Dict([
            "info" => info,
            "paths" => paths,
            "test" => valid_basic
        ])
    ]),
    "Optional Fields tests" => Dict([
        "Invalid fields" => Dict([
            "info" => info,
            "paths" => paths,
            "optional_fields" => get_invalid_optional_fields,
            "test" => throws
        ]),
        "All valid" => Dict([
            "info" => info,
            "paths" => paths,
            "optional_fields" => get_optional_fields,
            "test" => valid_basic
        ])
    ]),
    "paths' field tests" => Dict([
        "Invalid fields" => Dict([
            "info" => info,
            "paths" => invalid_paths,
            "test" => throws
        ])
    ])
])

for version in versions
    for (case, val) in cases
        @testset "Version $(version) - $(case)" begin
            for (subcase_name, subcase) in val
                @testset "$(subcase_name)" begin
                    openApi = OpenAPI(version, subcase["info"])
                    openApi.paths = subcase["paths"]
                    if haskey(subcase, "optional_fields")
                        openApi.optional_fields = subcase["optional_fields"](version)
                    end
                    subcase["test"](openApi, version)
                end
            end
        end
    end
end

const valid_specs = ["spec1_valid_v3.yml", "spec1_valid_v2.yml", "spec1_valid_v3.json", "spec1_valid_v2.json"]
const invalid_specs = ["spec1_invalid_v3.yml", "spec1_invalid_v2.yml", "spec1_invalid_v3.json", "spec1_invalid_v2.json"]

# More complex tests
@testset "Validate more complex specifications" begin
    for spec in valid_specs
        @test SwaggerMarkdown.validate_spec(joinpath(SPECS_PATH, spec)) === nothing
    end
    for spec in invalid_specs
        @test_throws SwaggerMarkdown.InvalidSwaggerSpecificationException SwaggerMarkdown.validate_spec(joinpath(SPECS_PATH, spec))
    end
end


text:
export Diff
"""
    Diff{GT, N} <: TagBlock{GT, N}
    Diff(block) -> Diff

Mark a block as quantum differentiable.
"""
struct Diff{GT, N} <: TagBlock{GT, N}
    content::GT
    function Diff(content::AbstractBlock{N}) where {N}
        @warn "Diff block has been deprecated, please use `Yao.AD.NoParams` to block non-differential parameters."
        new{typeof(content), N}(content)
    end
end
content(cb::Diff) = cb.content
chcontent(cb::Diff, blk::AbstractBlock) = Diff(blk)
YaoBlocks.PropertyTrait(::Diff) = YaoBlocks.PreserveAll()

apply!(reg::AbstractRegister, db::Diff) = apply!(reg, content(db))
mat(::Type{T}, df::Diff) where T = mat(T, df.content)
Base.adjoint(df::Diff) = chcontent(df, content(df)')

function YaoBlocks.print_annotation(io::IO, df::Diff)
    printstyled(io, "[∂] "; bold=true, color=:yellow)
end

#### interface #####
export markdiff

"""
    markdiff(mode::Symbol, block::AbstractBlock) -> AbstractBlock
    markdiff(mode::Symbol) -> Function

automatically mark differentiable items in a block tree as differentiable.
"""
function markdiff end

# for QC
markdiff(block::Union{RotationGate, CPhaseGate}) = Diff(block)
# escape control blocks.
markdiff(block::ControlBlock) = block

function markdiff(blk::AbstractBlock)
    blks = subblocks(blk)
    isempty(blks) ? blk : chsubblocks(blk, markdiff.(blks))
end

YaoBlocks.AD.mat_back!(::Type{T}, db::Diff, adjm::AbstractMatrix, collector) where T = AD.mat_back!(T, content(db), adjm, collector)