text:
# Compute a few eigenvalues of the dielectric matrix at (q=0, ω=0) iteratively
using DFTK
using Plots
import KrylovKit: ArnoldiIterator, Orthogonalizer, OrthonormalBasis, KrylovDefaults, orthogonalize!
using KrylovKit
using Dates
using HDF5
using WriteVTK
using Printf

#
# For next time: Consider residual error in stopping criterion and
# stop only if estimated error in κ is less than 5%
#

# KrylovKit custom orthogonaliser
struct OrthogonalizeAndProject{F, O <: Orthogonalizer} <: Orthogonalizer
    projector!::F
    orth::O
end
OrthogonalizeAndProject(projector) = OrthogonalizeAndProject(projector, KrylovDefaults.orth)
function KrylovKit.orthogonalize!(v::T, b::OrthonormalBasis{T}, x::AbstractVector, alg::OrthogonalizeAndProject) where {T}
    v, x = orthogonalize!(v, b, x, alg.orth)
    v .= alg.projector!(v)::T
    v, x
end
function KrylovKit.orthogonalize!(v::T, q::T, alg::OrthogonalizeAndProject) where {T}
    v, s = orthogonalize!(v, q, alg.orth)
    v .= alg.projector!(v)::T
    v, s
end

# A straightfoward Arnoldi eigensolver that diagonalizes the matrix at each step
# This is more efficient than Arpack when `f` is very expensive
# Convergence is flagged after a given number of iterations or if the condition number
# changes by less than a given relative tolerance
function arnoldi(f, x0; cond_rtol=1e-3, iterations=30, min_iterations=10, n_print=5,
                 projector=identity)
    println()
    println("Starting Arnoldi (now = $(Dates.now()) ...")
    println()
    flush(stdout)

    oldcond = NaN
    for (V, B, r, nr, b) in ArnoldiIterator(f, x0, OrthogonalizeAndProject(projector))
        @assert isreal(r)

        # A * V = V * B + r * b'
        V = hcat(V...)
        AV = V*B + r*b'

        ew, ev = eigen(B, sortby=real)
        Vr = V*ev
        AVr = AV*ev
        R = AVr - Vr * Diagonal(ew)

        N = size(V, 2)
        normr = [norm(r) for r in eachcol(R)]

        println("#--- $N ---#")
        println("idcs      evals   residnorms   |P(X)|")
        inds = unique(append!(collect(1:min(N, n_print)), max(1, N-n_print):N))
        for i in inds
            Xi = @view Vr[:, i]
            norm_asym = abs(imag(ew[i])) < 1e-6 ? norm(Xi - projector(Xi)) : NaN
            @printf "% 3i  %10.6g  %10.6g  %10.6g\n" i real(ew[i]) normr[i] norm_asym
        end
        any(imag.(ew[inds]) .> 1e-5) && println("Warn: Suppressed imaginary part.")
        println()
        flush(stdout)

        cond = real(ew[end] / ew[1])
        is_converged = abs(cond - oldcond) / oldcond < cond_rtol && N > min_iterations
        oldcond = cond
        if N ≥ iterations || is_converged
            return (λ=ew, X=Vr, AX=AVr, residual_norms=normr)
        end
    end
end

function dump_arnoldi_hdf5(scfres, arnoldires, mixing, filename)
    gridsize = size(scfres.ρ.real)
    n_epairs = length(arnoldires.λ)
    isfile(filename) && rm(filename)
    h5open(filename, "w") do h5f
        arnoldi = g_create(h5f, "arnoldi")
        arnoldi["λ"]               = arnoldires.λ
        arnoldi["residual_norms"]  = arnoldires.residual_norms
        arnoldi["X",  "shuffle", (), "compress", 3] = real(reshape(arnoldires.X,
                                                                   gridsize..., n_epairs))
        arnoldi["AX", "shuffle", (), "compress", 3] = real(reshape(arnoldires.AX,
                                                                   gridsize..., n_epairs))
    end
end

function dump_arnoldi_vtk(scfres, arnoldires, mixing, filename)
    gridsize = size(scfres.ρ.real)
    arrays = Array{Float64, 3}[]
    labels = String[]

    for i in 1:length(arnoldires.λ)
        push!(arrays, real(reshape(arnoldires.X[:, i], gridsize)))
        push!(labels, "X_$(@sprintf "%02d" i)")

        push!(arrays, real(reshape(arnoldires.AX[:, i], gridsize)))
        push!(labels, "AX_$(@sprintf "%02d" i)")
    end

    # Save the LDOS
    if scfres.ham.basis.model.temperature > 0
        ldos = LDOS(scfres.εF, scfres.ham.basis, scfres.eigenvalues, scfres.ψ)
        push!(arrays, ldos)
        push!(labels, "LDOS")
    end

    idx = findfirst(isa.(scfres.ham.basis.model.term_types, AtomicLocal))
    push!(arrays, scfres.ham.basis.terms[idx].potential)
    push!(labels, "pot_atomic_local")

    push!(arrays, DFTK.total_local_potential(scfres.ham))
    push!(labels, "pot_total_local")

    push!(arrays, scfres.ρ.real)
    push!(labels, "ρ")
    vtk_write_array(filename, Tuple(arrays), Tuple(labels))
end

function run_arnoldi_dielectric(scfres; mixing=nothing, iterations=30, cond_rtol=5e-4,
                                fileprefix="arnoldi", force_lda=false)
    if isfile(fileprefix * ".hdf5")
        arnoldifile = fileprefix * ".hdf5"
        return (λ=h5read(arnoldifile, "arnoldi/λ"),
                residual_norms=h5read(arnoldifile, "arnoldi/residual_norms"))
    end

    # self_consistent_field uses 3 extra bands, which are not converged by the eigensolver
    # For the χ0 application, bands need to be perfectly converged, so we discard them here
    ψ_cvg = [@view ψk[:, 1:end-3]  for ψk in scfres.ψ]
    eigenvalues_cvg = [εk[1:end-3] for εk in scfres.eigenvalues]

    ham = scfres.ham
    basis = ham.basis
    if force_lda
        origmodel = basis.model
        model = model_LDA(origmodel.lattice, origmodel.atoms, smearing=origmodel.smearing,
                          temperature=origmodel.temperature)
        kcoords = [kpt.coordinate for kpt in basis.kpoints]
        basis = PlaneWaveBasis(model, basis.Ecut, kcoords, basis.ksymops, basis.symops,
                               fft_size=basis.fft_size)
    end

    # Apply ε = 1 - χ0 (vc + fxc)
    function eps_fun(dρ)
        dρ = reshape(dρ, size(scfres.ρ.real))
        dρ = from_real(basis, dρ)
        dv = apply_kernel(basis, dρ; ρ=scfres.ρ)
        χdv = apply_χ0(ham, ψ_cvg, scfres.εF, eigenvalues_cvg, dv)
        vec((dρ - χdv).real)
    end
    function pinv_fun(dρ)
        dρ = reshape(dρ, size(scfres.ρ.real))
        dρ_dc = sum(dρ) / length(dρ)
        dρ .-= dρ_dc

        # Make dummy ρin and ρout to use mixing interface,
        # such that δρ = ρout - ρin is as required
        ρin  = scfres.ρ
        ρout = from_real(basis, dρ .+ scfres.ρ.real)
        ρnext = DFTK.mix(mixing, basis, ρin, ρout; scfres...)

        ρdiff = ρnext - ρin
        vec(ρdiff.real .+ dρ_dc)
    end

    # Projector into the subspace of interest
    # (i.e. modes symmetric wrt. the lattice symmetries of the problem)
    function projector(x)
        x_dc = sum(x) / length(x)
        x .-= x_dc
        x_real = reshape(real(x), size(scfres.ρ.real))
        vec(DFTK.symmetrize(from_real(basis, x_real)).real)
    end

    apply = projector ∘ eps_fun
    if !isnothing(mixing)
        apply = projector ∘ pinv_fun ∘ eps_fun
    end

    # Construct guess: Normalised, mostly random, but large mass on smallest G mode
    aguess = 0.1randn(size(scfres.ρ.real))
    Gmin = minimum(G -> iszero(G) ? Inf : norm(model.recip_lattice * G), G_vectors(basis))
    idx_Gmin = findfirst(G -> abs(norm(model.recip_lattice * G) - Gmin) < 1e-12,
                         collect(G_vectors(basis)))
    aguess[idx_Gmin] = 2.0
    aguess /= norm(aguess)
    aguess = projector(aguess)
    aguess = projector(aguess)

    arnoldires = arnoldi(apply, vec(aguess), cond_rtol=cond_rtol, iterations=iterations,
                         projector=projector)
    dump_arnoldi_hdf5(scfres, arnoldires, mixing, fileprefix * ".hdf5")
    dump_arnoldi_vtk(scfres, arnoldires, mixing, fileprefix)

    arnoldires
end


text:
#
# Methods for orbital systems.
#

"""
Returns the mass parameter of the R2BP system.
"""
massparameter(system::R2BPParameters) = OrbitalStates.get_μ(system) * massparamunit(system)
massparameter(orbit::R2BPOrbit) = massparameter(system(orbit))

"""
Returns the mass parameter of the CR3BP system.
"""
massparameter(system::CR3BPParameters) = OrbitalStates.get_μ(system)
massparameter(orbit::CR3BPOrbit) = massparameter(system(orbit))

"""
Returns the mass parameters of the CR3BP system.
"""
massparameters(system::CR3BPParameters) = (primary_massparameter(system), secondary_massparameter(system))
massparameters(orbit::CR3BPOrbit) = massparameters(system(orbit))

"""
Returns the primary mass parameter of the CR3BP system.
"""
function primary_massparameter(system::CR3BPParameters)
    DU  = lengthunit(system)
    TU  = timeunit(system)
    μ   = massparameter(system)
    ∑μᵢ = DU^3 / ((TU / 2π)^2)

    return ∑μᵢ - μ * ∑μᵢ
end
primary_massparameter(orbit::CR3BPOrbit) = primary_massparameter(system(orbit))

"""
Returns the secondary mass parameter of the CR3BP system.
"""
function secondary_massparameter(system::CR3BPParameters)
    DU  = lengthunit(system)
    TU  = timeunit(system)
    μ   = massparameter(system)
    ∑μᵢ = DU^3 / ((TU / 2π)^2)

    return μ * ∑μᵢ
end
secondary_massparameter(orbit::CR3BPOrbit) = secondary_massparameter(system(orbit))

ext: jl


text:
#Operadores
println(40 + 2, " Soma")
println(43 - 1," Subtração")
println(6 * 7," Multiplicação")
println(84 / 2," Divisão")
println(84.5 ÷ 2," Divisão inteira")
println(84.5 % 2," Resto divisão")
println(84 // 2," ", typeof(84 // 2))  #devem ser Int
println(2 \ 84," Divisão do contra")
println(6^2 + 6," Exponenciação")
println(1_000_000)

println("\n------ Cap 2 ------\n")
variavel = "valor" #convenção: tolower()
😎 = 👍 = 21
println("😎 + 👍 = ",😎 + 👍)
#Precedência PEMDAS, esq pra dir.
# Parenteses
# Exponenciação
# MD - Mult e Divisão
# AS - Adic e subt
println("Conca" * "tenar")
println("Spam" ^ 3)
println(4//3*π*5^3)

println("\n------ Cap 3: Funções ------\n")
println(parse(Int64, "42"))
println(trunc(Int64, 3.999))
println("dB = ", 10*log10(2))
println("radianos = ", 0.7,'\n',"altura = ",sin(0.7))
#angulos por padrão são em radianos
graus = 45
println("convertido em ", graus / 180 * π," rad")
println(sqrt(4))

function baskhara(a, b, c)
    Δ = delta(a,b,c)
    return [(-b+√Δ)/2a, (-b-√Δ)/2a]
end
function delta(a,b,c)
    return b^2 + 4a*c
end

roots = baskhara(2,2,0)
println(roots)
# definição de funções não mudam o fluxo de execução do programa.
#quando é chamada,interrompe o fluxo, executa seu bloco e retorna

#funções "prudutivas" retornam algo, as "nulas" retornam nothimg.

println("println") #imprime e pula linha
print("um print, ") #não pula linha
print("outro print")

println("\n------ Cap 4: Design de interface ------\n")
using ThinkJulia
🐢 = Turtle()

@svg begin
    forward(🐢, 100)
    turn(🐢, -90)
    forward(🐢, 100)
    turn(🐢, -90)
    forward(🐢, 100)
    turn(🐢, -90)
    forward(🐢, 100)
end

@svg begin
    for i in 1:4
        forward(🐢, 100)
        turn(🐢, -90)
    end
end

function quadrado(t, com)
    @svg begin
        for i in 1:4
            forward(t, com)
            turn(t, -90)
        end
    end
end
quadrado(🐢,50)

function polígono(t, com,n)
    angulo = 360/n
    polilinha(t, com, n, angulo)
end
polígono(🐢,50,6)

function circulo(t, r)
    arco(t, r, 360)
end
circulo(🐢,100)

function arco(t, r, θ)
    arco = 2π*r*θ / 360
    n = arco ÷ 3 + 1
    com = arco/n
    angulo = θ / n
    #fazendo uma leve curva para a esquerda antes de iniciar
    #reduz o erro causado pela aproximação linear do arco
    turn(t, -angulo/2)
    polilinha(t, com, n, angulo)
    turn(t, -angulo/2)
end

#Docstring
"""
polilinha(t, com, n, angulo)

Desenha n segmentos de linha dado o comprimento(com)
e o angulo(em graus) entre eles. t é uma tartaruga.
"""
function polilinha(t, com, n, angulo)
    @svg begin
        for i in 1:n
            forward(t, com)
            turn(t, -angulo)
        end
    end
end
arco(🐢,100,135)
arco(🐢,100,135)
@svg
"""🐱 = Turtle()
@svg begin
    turn(🐱, 135)
    forward(🐱, 200)
end"""

#exercicio 4
# 4.1

function flor()
    arco = 2π*r*θ / 360
    n = arco ÷ 3 + 1
    com = arco/n
    angulo = θ / n
    #fazendo uma leve curva para a esquerda antes de iniciar
    #reduz o erro causado pela aproximação linear do arco
    turn(t, -angulo/2)
    polilinha(t, com, n, angulo)
    turn(t, -angulo/2)

    @svg begin
        for i in 1:n
            forward(t, com)
            turn(t, -angulo)
        end
    end
end

println("\n------ Cap 5: Cond e Recursão ------\n")
println(true && false)  #AND
println(true || false)  #OR
println(!true)  #NOT

function contagem_regressiva(n)
    if n ≤ 0
        println("Feliz Ano Novo!")
    else
        print(n, " ")
        contagem_regressiva(n-1)
    end
end
contagem_regressiva(3)

function imprima_n(s, n)
    if n ≤ 0
        return
    end
    println(s)
    imprima_n(s, n-1)
end
imprima_n("spam", 4)

texto = readline()
print(texto)

println("Qual é a velocidade de voo de uma andorinha sem carga?");
velocidade = readline()
println(typeof(velocidade))
velocidade = parse(Int64, velocidade)
println(typeof(velocidade))

potência_do_sinal = 9
potência_do_ruido = 10
razão = potência_do_sinal / potência_do_ruido
decibéis = 10 * log10(razão)
println(decibéis)

seg = time()
dias = seg / 60 / 60 / 24
println(dias)
hora = dias - dias ÷ 1
hora *= 60
dias ÷= 1
println(dias)
println(hora)

function distância(x₁, y₁, x₂, y₂)
    dx = x₂-x₁
    dy = y₂-y₁
    d² = dx^2 + dy^2
    sqrt(d²)
end
distância(1, 2, 4, 6)

function é_divisivel(x, y)
    x%y == 0
end
é_divisivel(5,2)

if é_divisivel(21, 3)
    println("x é divisível por y")
end

function fat(n)
    espaço = " " ^ (4 * n)
    println(espaço, "fatorial ", n)
    if n==0
        return 1
    else
        recursão = fat(n-1)
        resultado = n * recursão
        println(espaço, "retornando ", resultado)
        return resultado
    end
end
fat(4)

function fib(n)
    if !(n isa Int64)
        error("Fatorial é definido somente para números inteiros.")
    elseif n<0
        error("Fatorial não é definido para números inteiros negativos.")
    elseif n==0
        return 0
    elseif n==1
        return 1
    else
        return fib(n-1)+fib(n-2)
    end
end
fib(1.5)

function A(m, n)
    if m==0
        return n+1
    elseif m>0 && n==0
        return A(m-1,1)
    else
        return A(m-1,A(m,n-1))
    end
end
A(3,4)

function bhaskara(a,b,c)
    Δ = b^2 - 4a*c
    x₁ = (-b + √Δ) / 2a
    x₂ = (-b - √Δ) / 2a
    return [x₁,x₂]
end
println(bhaskara(3,4,2))

#y=1/2*(x+a/x) #metodo de Newton sqrt
a, x = 4, 3
y = (x + a/x) / 2
x = y
y = (x + a/x) / 2
x = y
y = (x + a/x) / 2
x = y
y = (x + a/x) / 2
x = y
y = (x + a/x) / 2
x = y
y = (x + a/x) / 2
# resultado parou de mudar, y == x.
while true
    println(x)
    y = (x + a/x) / 2
    if abs(y-x) < eps()
        break
    end
    x = y
end

#Conjectura de Collatz
function collatz(n)
    while n != 1
        println(n)
        if n % 2 == 0
            n = n/2
        else
            n = 3n + 1
        end
    end
end
collatz(77031)

#exercicio 7-2
function minha_raiz(a)
    x = a>1 ? a-1 : a
    while true
        y = (x + a/x) / 2
        if abs(y-x) < eps()
            return y
            break
        end
        x = y
    end

end
function avalia_raiz_quadrada()
    println("a   mysqrt             sqrt               diff")
    println("-   ------             ----               ----")
    for i in 1.0:9.0
        raiz_f, raiz_my = sqrt(i), minha_raiz(i)
        println(i, raiz_my, raiz_f,abs(raiz_f-raiz_my))
        #print(raiz_my)
        #print(raiz_f)
        #println(abs(raiz_f-raiz_my))
    end
end
avalia_raiz_quadrada()

exp = Meta.parse("1+2*3")
eval(exp)

function estima_pi() #FAIL
    a = (2*√2)/9801
    k = 0
    um_sopi = 0
    while true
        um_sopi += (factorial(4big(k))*(1103+2690*k)) / (((factorial(k))^4) * (396^(4*k)))
        println(a*um_sopi)
        println(k)
        if abs(π-(a*um_sopi)) <= 1e-15
            println("Hm...")
            return a * um_sopi
        end
        k += 1
    end
end
println(estima_pi()) #something it's wrongs

max_stars_repo_name: Mutanne/hiworld

max_issues_repo_name: Mutanne/hiworld

max_forks_repo_name: Mutanne/hiworld


text:
using TMLEEpistasis
using ArgParse

function parse_commandline()
    s = ArgParseSettings(description="Compute the Sieve Variance Plateau estimate for each phenotype in the result file")

    @add_arg_table s begin
        "prefix"
            help = "Prefix to the .hdf5 files generated by TMLEEpistasis `ukbb.jl` script"
            arg_type = String
            required = true
        "grm-prefix"
            arg_type = String
            help = "Prefix of the aggregated GRM"
            required = true
        "out"
            arg_type = String
            help = "output filename"
            required = true
        "--nb-estimators", "-n"
            arg_type = Int
            help = "Number of variance estimators to compute"
            default = 10
        "--max-tau", "-m"
            arg_type = Float64
            help = "Maximum distance of individuals to take into account (maximum=2)"*
                   "It was witnessed that beyond 0.9, weird limit effects happen"
            default = 0.8
        "--pval", "-p"
            arg_type = Float64
            help = "Only traits/queries pairs that satisfy this pvalue threshold with an independence"*
                    " assumption between individuals will be considered for the sieve."
            default = 0.05

    end

    return parse_args(s)
end

parsed_args = parse_commandline()

sieve_variance_plateau(parsed_args)


text:
#--------------------------------------------------------------------------------------------------
# Functions on [a,b]
#--------------------------------------------------------------------------------------------------
using Test
using Hensel

#--------------------------------------------------------------------------------------------------

@test Hensel.hVector(11.5, 11, 12, 2, 3)    == [1, 0, 0]
@test Hensel.hVector(0.0, -1.0, 1.0, 2, 3)  == [1, 0, 0]
@test Hensel.hVector(-0.9, -1.0, 1.0, 2, 3) == [0, 0, 0]
@test Hensel.hVector(100.3, 100, 101, 2, 3) == [0, 1, 0]

#--------------------------------------------------------------------------------------------------

max_stars_repo_name: ayarodionov/Hensel

max_issues_repo_name: ayarodionov/Hensel

max_forks_repo_name: ayarodionov/Hensel


text:
# load comparison data
using JLD, LinearAlgebra, SparseArrays
(@isdefined randconf) || (global const randconf = load("data/O3.jld", "randconf"))

# function rewrite()
#     objs = Dict{String, Any}()

#     jldopen("data/O3.jld", "r") do f
#         objkeys = names(f)
#         for k in objkeys
#             objs[k] = read(f[k])
#         end
#     end

#     jldopen("O3new.jld", "w") do fout
#         for (k, v) in objs
#             fout[k] = v
#         end
#     end
#     nothing
# end

# set up minimal O3 simulation
mc = mc_from_inxml("parameters/O3_generic_small_system.in.xml", randconf)
mc_odd_L = mc_from_inxml("parameters/O3_odd_L_small_system.in.xml", load("data/O3.jld", "randconf_odd_L"))
mc_nob = mc_from_inxml("parameters/O3_no_bfield_small_system.in.xml", randconf)
mc_nob_nochkr = mc_from_inxml("parameters/O3_no_bfield_no_chkr_small_system.in.xml", randconf)
# set to nothing at the end of this file to trigger GC

@testset "O3 model" begin

    @testset "types" begin
        @test geltype(mc) == ComplexF64
        @test heltype(mc) == ComplexF64
        @test cbtype(mc) == CBAssaad

        @test geltype(mc_nob) == ComplexF64
        @test heltype(mc_nob) == Float64
        @test cbtype(mc_nob) == CBAssaad

        @test geltype(mc_nob_nochkr) == ComplexF64
        @test heltype(mc_nob_nochkr) == Float64
        @test cbtype(mc_nob_nochkr) == CBFalse
    end

    @testset "boson action" begin
        @test isapprox(calc_boson_action(mc, randconf), 75.57712964980982)
        @test isapprox(calc_boson_action_diff(mc, 3, 10, [0.370422, 0.797014, 0.956094]), 0.015453643701304112)

        # edrun test
        mc.p.edrun = true
        @test isapprox(calc_boson_action(mc, randconf), 16.348917129437076)
        @test isapprox(calc_boson_action_diff(mc, 3, 10, [0.370422, 0.797014, 0.956094]), 0.014491910524600016)
        mc.p.edrun = false
    end

    @testset "hopping matrices" begin
        # init_hopping_matrices(mc)
        @test isapprox(mc.l.hopping_matrix_exp, load("data/O3.jld", "hopping_matrix_exp"))
        @test isapprox(mc.l.hopping_matrix_exp_inv, load("data/O3.jld", "hopping_matrix_exp_inv"))
    end

    # test peirls phases
    # TODO: load a NNNN model for this (to also check higher order phases)
    include("tests_O3_peirls.jl")

    # check chkr_error ~ O(Δτ^2)
    function check_chkr_quadratic(mc)
        dt_before = mc.p.delta_tau
        for dt in (0.1, 0.01, 0.001, 0.0001, 0.00001)
            mc.p.delta_tau = dt
            init_hopping_matrices(mc)
            hop = mc.l.hopping_matrix_exp
            hop_chkr = foldl(*,mc.l.chkr_hop_half) * sqrt.(mc.l.chkr_mu)
            if maximum(absdiff(hop, hop_chkr)) > mc.p.delta_tau^2
                mc.p.delta_tau = dt_before
                init_hopping_matrices(mc)
                return false
            end
        end
        mc.p.delta_tau = dt_before
        init_hopping_matrices(mc)
        return true
    end

    @testset "checkerboard hopping matrices" begin
        # init_checkerboard_matrices(mc_nob)
        @test find_four_site_hopping_corners(mc_nob.l) == ([1, 3, 9, 11], [6, 8, 14, 16])

        @test isapprox(mc_nob.l.chkr_hop[1], load("data/O3.jld", "nob_chkr_hop1"))
        @test isapprox(mc_nob.l.chkr_hop[2], load("data/O3.jld", "nob_chkr_hop2"))
        @test isapprox(mc_nob.l.chkr_hop_inv[1], load("data/O3.jld", "nob_chkr_hop_inv1"))
        @test isapprox(mc_nob.l.chkr_hop_inv[2], load("data/O3.jld", "nob_chkr_hop_inv2"))

        @test isapprox(mc_nob.l.chkr_hop_half[1], load("data/O3.jld", "nob_chkr_hop_half1"))
        @test isapprox(mc_nob.l.chkr_hop_half[2], load("data/O3.jld", "nob_chkr_hop_half2"))
        @test isapprox(mc_nob.l.chkr_hop_half_inv[1], load("data/O3.jld", "nob_chkr_hop_half_inv1"))
        @test isapprox(mc_nob.l.chkr_hop_half_inv[2], load("data/O3.jld", "nob_chkr_hop_half_inv2"))

        @test isapprox(mc_nob.l.chkr_mu, load("data/O3.jld", "nob_chkr_mu"))
        @test isapprox(mc_nob.l.chkr_mu_half, load("data/O3.jld", "nob_chkr_mu_half"))
        @test isapprox(mc_nob.l.chkr_mu_inv, load("data/O3.jld", "nob_chkr_mu_inv"))
        @test isapprox(mc_nob.l.chkr_mu_half_inv, load("data/O3.jld", "nob_chkr_mu_half_inv"))

        @test check_chkr_quadratic(mc_nob)
    end

    @testset "checkerboard hopping matrices (Bfield)" begin
        @test isapprox(build_four_site_hopping_matrix_Bfield(mc, 3, 1, 2), load("data/O3.jld", "build_four_site_hopping_matrix_Bfield"))

        @test isapprox(mc.l.chkr_hop[1], load("data/O3.jld", "chkr_hop1"))
        @test isapprox(mc.l.chkr_hop[2], load("data/O3.jld", "chkr_hop2"))
        @test isapprox(mc.l.chkr_hop_inv[1], load("data/O3.jld", "chkr_hop_inv1"))
        @test isapprox(mc.l.chkr_hop_inv[2], load("data/O3.jld", "chkr_hop_inv2"))

        @test isapprox(mc.l.chkr_hop_half[1], load("data/O3.jld", "chkr_hop_half1"))
        @test isapprox(mc.l.chkr_hop_half[2], load("data/O3.jld", "chkr_hop_half2"))
        @test isapprox(mc.l.chkr_hop_half_inv[1], load("data/O3.jld", "chkr_hop_half_inv1"))
        @test isapprox(mc.l.chkr_hop_half_inv[2], load("data/O3.jld", "chkr_hop_half_inv2"))

        @test isapprox(mc.l.chkr_mu, load("data/O3.jld", "chkr_mu"))
        @test isapprox(mc.l.chkr_mu_half, load("data/O3.jld", "chkr_mu_half"))
        @test isapprox(mc.l.chkr_mu_inv, load("data/O3.jld", "chkr_mu_inv"))
        @test isapprox(mc.l.chkr_mu_half_inv, load("data/O3.jld", "chkr_mu_half_inv"))

        @test check_chkr_quadratic(mc)
    end

    @testset "generic checkerboard hopping matrices" begin
        @assert cbtype(mc_odd_L) == CBGeneric
        @test check_chkr_quadratic(mc_odd_L)
    end

    @testset "interactions" begin
        @testset "setblockdiag!" begin
            A = sprand(mc.l.sites*4, mc.l.sites*3, 0.1)
            B = copy(A)
            setblockdiag!(mc.l, B, 2,3, 1:mc.l.sites)
            @test diag(B[(mc.l.sites+1):(2*mc.l.sites),(2*mc.l.sites+1):(3*mc.l.sites)]) == collect(1.0:mc.l.sites)
            setblockdiag!(mc.l, B, 2,3, diag(A[(mc.l.sites+1):(2*mc.l.sites),(2*mc.l.sites+1):(3*mc.l.sites)]))
            @test A == B
        end

        @testset "interaction matrix exponentials" begin
            eV = similar(mc_nob.s.eV)
            interaction_matrix_exp!(mc_nob, 3, 1., eV)
            @test isapprox(eV, load("data/O3.jld", "nob_eVplus"))
            @test isapprox(eV, interaction_matrix_exp(mc_nob, 3, 1.))
            interaction_matrix_exp!(mc_nob, 3, -1., eV)
            @test isapprox(eV, load("data/O3.jld", "nob_eVminus"))

            @test isapprox(interaction_matrix_exp_op(mc_nob, [0.130018, 0.792039, 0.683411], 1.), load("data/O3.jld", "nob_eVexpop"))
        end

        @testset "interaction matrix exponentials (Bfield)" begin
            eV = similar(mc.s.eV)
            interaction_matrix_exp!(mc, 3, 1., eV)
            @test isapprox(eV, load("data/O3.jld", "eVplus"))
            @test isapprox(eV, interaction_matrix_exp(mc, 3, 1.))
            interaction_matrix_exp!(mc, 3, -1., eV)
            @test isapprox(eV, load("data/O3.jld", "eVminus"))

            @test isapprox(interaction_matrix_exp_op(mc, [0.130018, 0.792039, 0.683411], 1.), load("data/O3.jld", "eVexpop"))
        end
    end

    @testset "local updates" begin
        init!(mc)
        @test isapprox(calc_detratio(mc, 7, [0.488033, 0.0196912, 0.438309]), 1.000380293015979 + 1.0842021724855044e-18im)
        @test isapprox(mc.s.delta_i, load("data/O3.jld", "delta_i"))
        @test isapprox(mc.s.M, load("data/O3.jld", "M"))

        update_greens!(mc, 7)
        @test isapprox(mc.s.greens, load("data/O3.jld", "afterupdate_greens"))

        # TODO: test full local update
        Random.seed!(123456789) # set seed
        @test local_updates(mc) == 0.6875 # acc. rate
        @test isapprox(mc.s.greens, load("data/O3.jld", "afterlocal_greens"))
        @test isapprox(mc.p.hsfield, load("data/O3.jld", "afterlocal_hsfield"))
        @test isapprox(mc.p.boson_action, 75.18407422927604)
    end

    @testset "global updates" begin
        # TODO
        # Random.seed!(123456789)
        # hsfield = copy(mc.p.hsfield)
        # global_update_perform_shift!(mc)
        # randsite = rand(1:mc.l.sites)
        # randslice = rand(1:mc.p.slices)
        # p.hsfield[:,randsite,randslice]
        # isapprox( == )

        # TODO: test global_update_backup_swap!
    end

    @testset "stack" begin
        init!(mc)
        @test isapprox(mc.s.greens, calc_greens(mc, mc.s.current_slice)) # compare against "exact" greens

        @testset "greens (+ logdet) calculation" begin
            # load Ur, Dr, Tr, Ul, Dl, Tl of time slice 1, direction up
            mc.s.Ur = load("data/O3.jld", "Ur")
            mc.s.Dr = load("data/O3.jld", "Dr")
            mc.s.Tr = load("data/O3.jld", "Tr")
            mc.s.Ul = load("data/O3.jld", "Ul")
            mc.s.Dl = load("data/O3.jld", "Dl")
            mc.s.Tl = load("data/O3.jld", "Tl")
            greens = load("data/O3.jld", "greens")
            calculate_greens(mc)
            ld = calculate_logdet(mc)
            gfresh, ldfresh = calc_greens_and_logdet(mc, 1)
            @test isapprox(mc.s.greens, greens) # compare against dumped greens
            @test isapprox(mc.s.greens, gfresh) # compare against "exact" freshly calculated greens
            @test isapprox(ld, ldfresh) # compare against "exact" logdet
        end

        @testset "wrapping" begin
            g = calc_greens(mc,2)
            gup = wrap_greens(mc, g, 2, 1)
            gdwn = wrap_greens(mc, g, 2, -1)
            @test isapprox(gup, calc_greens(mc, 3))
            @test isapprox(gdwn, calc_greens(mc, 1))
        end

        init!(mc)
        @testset "propagation" begin
            @assert mc.s.current_slice == mc.p.slices
            while mc.s.current_slice != 1
                propagate(mc)
            end
            propagate(mc) # we want to go in direction up
            greens = load("data/O3.jld", "greens")
            @test isapprox(mc.s.greens, greens) # compare against dumped greens
        end

        init!(mc)
        @testset "propagation: greens always within bounds" begin
            @assert mc.s.current_slice == mc.p.slices
            @assert mc.s.direction == -1

            propagate(mc)
            while !(mc.s.current_slice == mc.p.slices && mc.s.direction == -1)
                mc.s.greens
                greens = calc_greens(mc, mc.s.current_slice)
                @test maximum(absdiff(greens, mc.s.greens)) < 1e-12
                propagate(mc)
            end
        end
    end

    @testset "lattice" begin
        @test mc.l.neighbors == [2 3 4 1 6 7 8 5 10 11 12 9 14 15 16 13;
                                5 6 7 8 9 10 11 12 13 14 15 16 1 2 3 4;
                                4 1 2 3 8 5 6 7 12 9 10 11 16 13 14 15;
                                13 14 15 16 1 2 3 4 5 6 7 8 9 10 11 12]
        @test mc.l.time_neighbors == [2 3 4 5 6 7 8 9 10 1; 10 1 2 3 4 5 6 7 8 9]
        @test mc.l.bonds == [1 2 0 0; 1 5 0 0; 2 3 0 0; 2 6 0 0; 3 4 0 0; 3 7 0 0;
                        4 1 0 0; 4 8 0 0; 5 6 0 0; 5 9 0 0; 6 7 0 0; 6 10 0 0;
                        7 8 0 0; 7 11 0 0; 8 5 0 0; 8 12 0 0; 9 10 0 0; 9 13 0 0;
                        10 11 0 0; 10 14 0 0; 11 12 0 0; 11 15 0 0; 12 9 0 0;
                        12 16 0 0; 13 14 0 0; 13 1 0 0; 14 15 0 0; 14 2 0 0;
                        15 16 0 0; 15 3 0 0; 16 13 0 0; 16 4 0 0]
        @test mc.l.dim == 2
    end

    @testset "slice matrices" begin
        init!(mc)
        init!(mc_nob_nochkr)
        init!(mc_odd_L)

        get_mc_and_suffix = (chkr) -> begin
            if chkr === CBAssaad
                return mc, ""
            elseif chkr === CBFalse
                return mc_nob_nochkr, "_CBFalse"
            else
                return mc_odd_L, "_CBGeneric"
            end
        end

        for chkr in [CBAssaad, CBFalse, CBGeneric]
            mc, suffix = get_mc_and_suffix(chkr)

            B = slice_matrix(mc, 3, 1.)
            Binv = slice_matrix(mc, 3, -1.)
            @test isapprox(B, load("data/O3.jld", "Bplus"*suffix))
            @test isapprox(Binv, load("data/O3.jld", "Bminus"*suffix))

            # multiply_daggered_B_left!
            A = rand!(similar(mc.s.greens))
            Ares = adjoint(B) * A
            multiply_daggered_B_left!(mc, 3, A)
            @test isapprox(A, Ares)

            # Binv * B * A == A and A * B * Binv == A
            A = rand!(similar(mc.s.greens))
            Aorig = copy(A)
            multiply_B_left!(mc, 3, A)
            multiply_B_inv_left!(mc, 3, A)
            @test isapprox(A, Aorig)
            multiply_B_right!(mc, 3, A)
            multiply_B_inv_right!(mc, 3, A)
            @test isapprox(A, Aorig)

            # B * A * Binv == A and Binv * A * B == A
            A = Matrix{geltype(mc)}(I, size(mc.s.greens)...)
            Aorig = copy(A)
            multiply_B_left!(mc, 3, A)
            multiply_B_inv_right!(mc, 3, A)
            @test isapprox(A, Aorig)
            multiply_B_right!(mc, 3, A)
            multiply_B_inv_left!(mc, 3, A)
            @test isapprox(A, Aorig)
        end
    end

    @testset "parameters constructors" begin
        # TODO: check if fields are set properly
    end

    @testset "xml read write" begin
        # TODO: xml read write
    end

    @testset "Measurements" begin
        include("tests_O3_measurements.jl")
    end

end # O3 model

mc = nothing
mc_odd_L = nothing
mc_nob = nothing
mc_nob_nochkr = nothing

nothing

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
flat2(arr) = (while any(a -> a isa Vector, arr) arr = collect(Iterators.flatten(arr)) end; arr)

arr = [[1], 2, [[3, 4], 5], [[[]]], [[[6]]], 7, 8, []]

@show flat(arr)
@show flat1(arr)
@show flat2(arr)

@btime flat(arr)
@btime flat1(arr)
@btime flat2(arr)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=
 The Computer Language Benchmarks Game
 https://salsa.debian.org/benchmarksgame-team/benchmarksgame/

 based on the C++ program of Alessandro Power
 with ideas borrowed from the C program of Mr Ledrug
 modified for Julia 1.0 by Luiz M. Faria
=#

# import a few functions from GMP wrapper to do in-place operations
using Base.GMP.MPZ: add_ui!, mul_ui!, add!, tdiv_q!
# and add a few more wrappers for GMP functions not currently in Base
const mpz_t = Ref{BigInt}
addmul_ui!(x::BigInt,a::BigInt,b) = ccall((:__gmpz_addmul_ui,"libgmp"), Cvoid, (mpz_t,mpz_t,Culong), x, a, b)
submul_ui!(x::BigInt,a::BigInt,b) = ccall((:__gmpz_submul_ui,"libgmp"), Cvoid, (mpz_t,mpz_t,Culong), x, a, b)
get_ui(x::BigInt) = ccall((:__gmpz_get_ui,"libgmp"), Culong, (mpz_t,), x)

function next!(q,r,t,k)
    mult  = 2*k + 1
    addmul_ui!(r,q,2)
    mul_ui!(r,mult)
    mul_ui!(t,mult)
    mul_ui!(q,k)
    return nothing
end

function extract!(q,r,t,tmp1,tmp2,n)
    mul_ui!(tmp1,q,n)
    add!(tmp1,r)
    tdiv_q!(tmp2,tmp1,t)
    return get_ui(tmp2)
end

function produce!(q,r,t,n)
    mul_ui!(q,10)
    submul_ui!(r,t,n)
    mul_ui!(r,10)
    return nothing
end

function pidigits(n::Int)
    q,r,t,tmp1,tmp2 = map(BigInt,(1,0,1,0,0))
    k = UInt(0)
    n_digits=0

    while n_digits<n
        i=0
        while (i<10 && n_digits<n)
            k +=1
            next!(q,r,t,k)
            q>r && continue

            digit = extract!(q,r,t,tmp1,tmp2,3)
            if digit == extract!(q,r,t,tmp1,tmp2,4)
                print(digit)
                produce!(q,r,t,digit)
                i +=1
                n_digits+=1
            end
        end
        while i < 10
            i+=1
            print(' ')
        end
        print("\t:$n_digits\n")
    end
end

n = parse(Int,ARGS[1])
pidigits(n)

max_stars_repo_name: bpecsek/Programming-Language-Benchmarks

max_issues_repo_name: bpecsek/Programming-Language-Benchmarks

max_forks_repo_name: bpecsek/Programming-Language-Benchmarks


text:
# Corrected by Emmanuel Lujan
using Printf
import Flux: σ
using ModelingToolkit
using GalacticOptim
using Optim
using DiffEqFlux
using NeuralPDE
using Quadrature, Cubature, Cuba

@parameters t,x
@variables c(..)
@derivatives Dt'~t
@derivatives Dxx''~x
@derivatives Dx'~x

# Parameters

v = 1
R = 0
D = 0.1 # diffusion
t_max = 1.0
x_min = -1.0
x_max = 1.0

# Equations, initial and boundary conditions
eqs = [ Dt(c(t, x)) ~ D * Dxx(c(t,x)) - Dx(c(t,x)) ]

bcs = [
         c(0, x) ~ cos(π*x) + 1.0,
         c(t, x_min) ~ c(t, x_max)
]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,t_max),
           x ∈ IntervalDomain(x_min,x_max)
]

# Discretization
nx = 32
dx = (x_max-x_min) / (nx - 1)
dt = t_max / 99. #

# Neural network
dim = length(domains)
output = length(eqs)
hidden = 8

chain = FastChain( FastDense(dim, hidden, σ),
                    FastDense(hidden, hidden, σ),
                    FastDense(hidden, 1))

strategy = GridTraining(dx=[dt,dx])

discretization = PhysicsInformedNN(chain, strategy=strategy)

pde_system = PDESystem(eqs, bcs, domains, [t,x], [c])
prob = discretize(pde_system,discretization)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

res = GalacticOptim.solve(prob,Optim.BFGS();cb=cb,maxiters=100)

# Plots

using Plots

phi = discretization.phi

initθ = discretization.initθ

acum =  [0;accumulate(+, length.(initθ))]
sep = [acum[i]+1 : acum[i+1] for i in 1:length(acum)-1]
minimizers = [res.minimizer[s] for s in sep]
ts,xs = [domain.domain.lower:dx:domain.domain.upper for domain in domains]

dt = t_max / 10.0
anim = @animate for (i, t) in enumerate(0:dt:t_max)
    @info "Animating frame $i..."
    c_predict = reshape([phi([t, x], res.minimizer)[1] for x in xs], length(xs))
    title = @sprintf("Advection-diffusion t = %.3f", t)
    plot(xs, c_predict, label="", title=title , ylims=(0., 2))
end

gif(anim, "advection_diffusion_pinn.gif", fps=15)

c_predict = reshape([ phi([0, x], res.minimizer)[1] for x in xs], length(xs))
plot(xs,c_predict)

max_stars_repo_name: drozzy/NeuralOceananigans

max_issues_repo_name: drozzy/NeuralOceananigans

max_forks_repo_name: drozzy/NeuralOceananigans


text:
module CompositeProblems

using LinearAlgebra
import StructuredProximalOperators:
    Manifold, g, prox_αg!, prox_αg, ∇M_g!, ∇M_g, ∇²M_g_ξ!, ∇²M_g_ξ
using StructuredProximalOperators
using Random
using Distributions
using DelimitedFiles
using JLD

using JuMP, OSQP

abstract type CompositeProblem end

const instances_dir = joinpath(dirname(pathof(CompositeProblems)), "..", "instances")

function __init__()
    !isdir(instances_dir) && mkdir(instances_dir)
end

##
function problem_dimension(pb::CompositeProblem)
    return error("problem_dimension not implemented for a problem $(typeof(pb)).")
end

function F(pb::CompositeProblem, x)
    return f(pb, x) + g(pb, x)
end

## Smooth side
# 0th order
function f(pb::CompositeProblem, x)
    return error("f not implemented for a problem $(typeof(pb)) and point $(typeof(x)).")
end

# 1st order
function ∇f!(pb::CompositeProblem, res, x)
    return error("∇f! not implemented for a problem $(typeof(pb)), result $(typeof(res)) and point $(typeof(x)).")
end
function ∇f(pb::CompositeProblem, x)
    res = zeros(size(x))
    ∇f!(pb, res, x)
    return res
end

# 2nd order
function ∇²f_h!(pb::CompositeProblem, res, x, h)
    return error("∇²f_h! not implemented for a problem $(typeof(pb)), result $(typeof(res)), point $(typeof(x)) and direction $(typeof(h)).")
end
function ∇²f_h(pb::CompositeProblem, x, h)
    res = zeros(size(h))
    ∇²f_h!(pb, res, x, h)
    return res
end

# conditioning
function get_gradlips(pb::CompositeProblem)
    return error("get_gradlips not implemented for a problem $(typeof(pb)).")
end
function get_μ_cvx(pb::CompositeProblem)
    return error("get_μ_cvx not implemented for a problem $(typeof(pb)).")
end

## Nonsmooth side
g(pb::CompositeProblem, x) = g(pb.regularizer, x)

prox_αg!(pb::CompositeProblem, res, x, α) = prox_αg!(pb.regularizer, res, x, α)
prox_αg(pb::CompositeProblem, x, α) = prox_αg(pb.regularizer, x, α)

∇M_g!(pb::CompositeProblem, M, ∇M_g, x) = ∇M_g!(pb.regularizer, M, ∇M_g, x)
∇M_g(pb::CompositeProblem, M, x) = ∇M_g(pb.regularizer, M, x)

∇²M_g_ξ(pb::CompositeProblem, M, x, ξ) = ∇²M_g_ξ(pb.regularizer, M, x, ξ)
∇²M_g_ξ!(pb::CompositeProblem, M, res, x, ξ) = ∇²M_g_ξ!(pb.regularizer, M, res, x, ξ)

export CompositeProblem
export F
export f, ∇f!, ∇f, ∇²f_h!, ∇²f_h
export get_gradlips, get_μ_cvx
export g, prox_αg!, prox_αg, ∇M_g!, ∇M_g, ∇²M_g_ξ!, ∇²M_g_ξ
export problem_dimension

include("problems/Leastsquares.jl")
include("problems/Logistic.jl")
include("problems/Logistic_instances.jl")
include("problems/Tracenorm.jl")

include("problems/Adrian's_testfunction.jl")

include("utils_instances.jl")
include("optimality_checkers.jl")

export LeastsquaresPb
export get_random_qualifiedleastsquares, get_random_qualifiedlasso, get_lasso_MLE

export LogisticPb
export get_logit_ionosphere, get_logit_gisette, get_random_logit, get_logit_MLE, get_logit_MLE_elastic_net

export TracenormPb
export get_tracenorm_MLE

export quadmaxquadAL

export firstorder_optimality_tangnorm
end


text:
@doc """
Calculates classic Dickey-Fuller (1979) statistic.

```math
  \\Delta y_t = \\alpha*y_{t-1} + \\epsilon_t, \\quad \\epsilon_t\\sim\\mathcal{N}(0,\\sigma^2)
```

Required Input Parameters
----------------------------------
`y::Vector{Float64}`: Observed time series in vector form.

Optional Input Parameters
----------------------------------
`trend::String`: Deterministic trend specification. It determines the asymptotic distribution of the t-statistic. Possible values are `"constant"`, which considers a constant deterministic term; `"linear"`, which also considers a linear time trend; and `"none"`, which considers no deterministic detrending.

Output
----------------------------------
`df_stat::Float64`: Dickey-Fuller (1979) test statistic.

`df_pval::Float64`: p-value obtained from 10000 simulations.
"""->
function df_test(y::Vector{Float64}; trend::String="constant", pval::Bool=true)
  #constante = false; trend = false;
  #Test regression: Δy_t = μ + αy_t-1 + ɛ_t, test if α = 0
  #                 Δy_t = μ + βt + αy_t-1 + ɛ_t
  if trend=="none"
    x = ar_lagmatrix(y, 1; trend="none")
    (pval==true) && (println("no deteriministic detrending."))
  elseif trend=="constant"
    (pval==true) && (println("constant deteriministic detrending."))
    x = ar_lagmatrix(y, 1; trend="constant")
  elseif trend=="linear"
    (pval==true) && (println("constant plus linear deteriministic detrending."))
    x = ar_lagmatrix(y, 1; trend="linear")
  else
    error(" trend must be either none, constant, or linear")
  end
  T = size(x, 1)
  Δy = diff(y)
  β, t_stat, bic = df_ols(Δy, x); # Test regression
  df_stat = t_stat[end]
  (pval==false) && return df_stat
  df_pval = mean(df_stat.>simulate_dfdist(T, trend=trend))
  return df_stat, df_pval
end
#
#
function simulate_dfdist(T; trend::String="constant")
 nsim = 10000
 t_dist = Vector{Float64}(nsim)
 for sim ∈ 1:nsim
   t_h0 = simul_arma(mu = 0.0, lin_trend = 0.0, ar = 1.0, ma = 0.0, T = T, y_0 = 0.0, stderr = 1.0) # simulate distribution of t under I(1) null hypothesis
   t_dist[sim] = df_test(t_h0[:]; pval=false) #here pval=true will break your computer
 end
 return t_dist
end
#
function df_critval(T; confidence_level::Float64=0.05, trend::String="constant")
  quantile(simulate_dfdist(T, trend=trend), confidence_level; sorted=false)
end
#
#
@doc """
Calculates the Augmented Dickey-Fuller (1979) ADF statistic.

```math
  \\Delta y_t = \\alpha*y_{t-1} + \\sum^{pmax}_{j=1} \\beta_j*\\Delta y_{t-j} + \\epsilon_t, \\quad \\epsilon_t\\sim\\mathcal{N}(0,\\sigma^2)
```

Required Input Parameters
----------------------------------
`y::Vector{Float64}`: Observed time series in vector form.

Optional Input Parameters
----------------------------------
`pmax::Int`: Maximum number of lags to be considered in the correction for serial correlation.

`trend::String`: Deterministic trend specification. It determines the asymptotic distribution of the t-statistic. Possible values are `"constant"`, which considers a constant deterministic term; `"linear"`, which also considers a linear time trend; and `"none"`, which considers no deterministic detrending.

Output
----------------------------------
`df_stat::Float64`: ADF test statistic.

`df_pval::Float64`: p-value obtained from 10000 simulations.
""" ->
function adf_test(y::Vector{Float64}; pmax::Integer=15, trend::String="constant")
  bics = Vector{Float64}(pmax)
  for p ∈ 1:pmax
    β_adf, adf_stat, bics[p,:] = adfreg(y, p; trend=trend)
  end
  bic_star, bic_index = findmin(bics)
  β_adf, adf_stat, bics = adfreg(y, bic_index; trend=trend)
  adf_pval = mean(adf_stat.>simulate_dfdist(length(y), trend=trend))
  return adf_stat, adf_pval
end
# http://www.eviews.com/help/helpintro.html#page/content/advtimeser-Unit_Root_Testing.html
function adfreg(y::Vector{Float64}, p;  trend::String="constant")
  x =  [ar_lagmatrix(y, 1, trend="none")[p+1:end] ar_lagmatrix(diff(y), p, trend=trend)] #ar_lagmatrix(y, 1, trend="none")
  dy = diff(y)[p+1:end]
  β, t_stat, bic = df_ols(dy, x)
  return β[1], t_stat[1], bic
end
#
#
function ar_lagmatrix(x::Vector{Float64}, p::Integer; trend::String="none")
  T = length(x)
  xlag = Array{Float64}(T-p, p)
  for ip ∈ 1:p
    xlag[:, ip] = x[(p-ip+1):(end-ip)]
  end
  (trend=="none") && (return xlag)
  if trend == "constant"
    xlag = hcat( ones(T-p), xlag)
    return xlag
  end
  if trend == "linear"
    xlag = hcat( ones(T-p), collect(linspace(1,T-p,T-p)), xlag )
    return xlag
  end
end

function df_ols(y::Vector{Float64}, x::Array{Float64,})
  T, k = size(x)
  β = inv(x'x)*(x'y)
  ϵ = y - x*β
  yhat = x*β
  sse = ϵ'ϵ
  σ2_y = sse/(T-k)
  Σₓ = inv(x'x)
  Σᵦ = σ2_y.*Σₓ
  t_stat = β./sqrt.(diag(Σᵦ))
  bic = (T-k).*log(sse./(T-k)) .+ k.*log(T-k)
  return β, t_stat, bic
end
#
@doc """
Simulates a stochastic process {yₜ} following an ARMA(1,1) with parameters φ and Θ for t = 1, ..., T periods.

```math
  y_t = \\mu_0 + \\mu_1 + \\phi*y_{t-1} + \\epsilon_t + \\theta*\\epsilon_{t-1}, \\quad \\epsilon_t \\sim \\mathbf{N}(0, \\sigma^2)
```

Required Input Parameters
----------------------------------
none

Optional Input Parameters
----------------------------------
`mu::Float64`: AR(1) process constant term. (default=0.0)

`lin_trend::Float64`: AR(1) process linear time trend. (default=0.0)

`ar::Float64`: AR(1) autorregressive parameter. (default=0.5)

`ma::Float64`: MA(1) autorregressive parameter. (default=0.5)

`T::Int64`: Length of simulated time series. (default=100)

`y_0::Float64`: Initial value of the series. (default=0.0)

`stderr::Float64` : Normally distributed innovation standard error. (default=1.0)

Output
----------------------------------
`y::Vector{Float64}`: Simulated ARMA(1,1) time series of length T.

""" ->
function simul_arma(;mu::Float64 = 0., lin_trend::Float64 = 0., ar::Float64 = .5, ma::Float64 = .5, T::Integer = 100, y_0::Float64 = 0.0, stderr::Float64 = 1.0)
  #y_t = μ + ϕy_t-1 + u_t/ u_t = ɛ_t + Θɛ_t-1
  y = zeros(T); ɛ = rand(Normal(0.0,stderr),T); u = zeros(T);
  for t in 2:T
    u[t, :] = ɛ[t, :] + ma*ɛ[t-1, :] #MA
    y[t, :] = mu + lin_trend*t + ar*y[t-1, :] + u[t, :] #AR
  end
  return y
end


text:
using Revise

using DataFrames, Statistics, StatsBase, LinearAlgebra, BenchmarkTools, Random, CUDA
CUDA.allowscalar(false)
Random.seed!(11)

function StatsBase.partialcor(X::AbstractArray{T,2}, Z::AbstractArray{T,2};
    zincludesintercept = false) where T
  N,K = size(X)
  @assert size(X,1) == size(Z,1)

  local Xe::Matrix
  if zincludesintercept
    Xₑ = X - Z * (cholesky!(Symmetric(Z'*Z))\(Z'*X))
  else
    Z1 = hcat(ones(T,N), Z)
    Xₑ = X - Z1 * (cholesky!(Symmetric(Z1'*Z1))\(Z1'*X))
  end
  return cor(Xₑ)
end

function StatsBase.pairwisepartialcor(X::AbstractArray{T,2}, Z::AbstractArray{T,2};
    zincludesintercept = false) where T
  N,K = size(X)
  @assert size(X,1) == size(Z,1)

  local Xe::Matrix
  if zincludesintercept
    Xₑ = X - Z * (cholesky!(Symmetric(Z'*Z))\(Z'*X))
  else
    Z1 = hcat(ones(T,N), Z)
    Xₑ = X - Z1 * (cholesky!(Symmetric(Z1'*Z1))\(Z1'*X))
  end
  return cor(Xₑ)
end

function testpartialcor(;N=252, K=500, R = 4)
  X = rand(N, K)
  Z = rand(N, R)
  #slow, manual version
  function partialcorver(X::AbstractArray{T,2}, Z::AbstractArray{T,2}) where T
    N,K = size(X)
    @assert size(X,1) == size(Z,1)

    ρ = ones(T, K,K)
    for c ∈ 1:K
      for r ∈ 1:(c-1)
        ρ[r,c] = partialcor(X[:,c], X[:,r],Z)
      end
    end

    ρ .= Symmetric(ρ)
    return ρ
  end

  #create an orthogonalized version for testing purposes
  Z1 = hcat(ones(N), Z)
  Xr = X .- Z1*(cholesky!(Z1'*Z1)\(Z1'*X))

  @info "testing partial correlations"
  @assert (partialcor(Xr, Z) .≈ cor(Xr)) |> all
  @assert (partialcor(X,Z) .≈ partialcor(X,Z1,zincludesintercept=true)) |> all
  #@btime partialcor($X, $Z1,zincludesintercept=true)
  @btime partialcor($X, $Z)

  if K ≤ 1000
    @assert (partialcorver(X,Z) .≈ cor(Xr)) |> all
    @assert (partialcor(X,Z) .≈ partialcorver(X,Z)) |> all
  end

end

@time testpartialcor(N=252, K=5000, R=4)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
@testset "bitcode" begin

Context() do ctx
    invalid_bitcode = "invalid"
    @test_throws LLVMException parse(LLVM.Module, unsafe_wrap(Vector{UInt8}, invalid_bitcode); ctx)
end

Context() do ctx
Builder(ctx) do builder
LLVM.Module("SomeModule"; ctx) do source_mod
    ft = LLVM.FunctionType(LLVM.VoidType(ctx))
    fn = LLVM.Function(source_mod, "SomeFunction", ft)

    entry = BasicBlock(fn, "entry"; ctx)
    position!(builder, entry)

    ret!(builder)

    verify(source_mod)

    bitcode_buf = convert(MemoryBuffer, source_mod)

    let
        mod = parse(LLVM.Module, bitcode_buf; ctx)
        verify(mod)
        @test haskey(functions(mod), "SomeFunction")
        dispose(mod)
    end

    bitcode = convert(Vector{UInt8}, source_mod)

    let
        mod = parse(LLVM.Module, bitcode; ctx)
        verify(mod)
        @test haskey(functions(mod), "SomeFunction")
        dispose(mod)
    end

    mktemp() do path, io
        mark(io)
        write(io, source_mod)
        flush(io)
        reset(io)

        @test read(io) == bitcode
    end
end
end
end

end


text:
using ASTInterpreter2
using Test
using REPL

using DebuggerFramework

include("evaling.jl")
include("stepping.jl")
include("interpret.jl")
include("ui.jl")
include("misc.jl")


text:
# Julia wrapper for header: /usr/include/scip/pub_dialog.h
# Automatically generated using Clang.jl wrap_c

function SCIPdialoghdlrGetRoot(dialoghdlr)
    ccall((:SCIPdialoghdlrGetRoot, libscip), Ptr{SCIP_DIALOG}, (Ptr{SCIP_DIALOGHDLR},), dialoghdlr)
end

function SCIPdialoghdlrClearBuffer(dialoghdlr)
    ccall((:SCIPdialoghdlrClearBuffer, libscip), Cvoid, (Ptr{SCIP_DIALOGHDLR},), dialoghdlr)
end

function SCIPdialoghdlrIsBufferEmpty(dialoghdlr)
    ccall((:SCIPdialoghdlrIsBufferEmpty, libscip), UInt32, (Ptr{SCIP_DIALOGHDLR},), dialoghdlr)
end

function SCIPdialoghdlrGetLine(dialoghdlr, dialog, prompt, inputline, endoffile)
    ccall((:SCIPdialoghdlrGetLine, libscip), SCIP_RETCODE, (Ptr{SCIP_DIALOGHDLR}, Ptr{SCIP_DIALOG}, Cstring, Ptr{Cstring}, Ptr{UInt32}), dialoghdlr, dialog, prompt, inputline, endoffile)
end

function SCIPdialoghdlrGetWord(dialoghdlr, dialog, prompt, inputword, endoffile)
    ccall((:SCIPdialoghdlrGetWord, libscip), SCIP_RETCODE, (Ptr{SCIP_DIALOGHDLR}, Ptr{SCIP_DIALOG}, Cstring, Ptr{Cstring}, Ptr{UInt32}), dialoghdlr, dialog, prompt, inputword, endoffile)
end

function SCIPdialoghdlrAddInputLine(dialoghdlr, inputline)
    ccall((:SCIPdialoghdlrAddInputLine, libscip), SCIP_RETCODE, (Ptr{SCIP_DIALOGHDLR}, Cstring), dialoghdlr, inputline)
end

function SCIPdialoghdlrAddHistory(dialoghdlr, dialog, command, escapecommand)
    ccall((:SCIPdialoghdlrAddHistory, libscip), SCIP_RETCODE, (Ptr{SCIP_DIALOGHDLR}, Ptr{SCIP_DIALOG}, Cstring, UInt32), dialoghdlr, dialog, command, escapecommand)
end

function SCIPdialogHasEntry(dialog, entryname)
    ccall((:SCIPdialogHasEntry, libscip), UInt32, (Ptr{SCIP_DIALOG}, Cstring), dialog, entryname)
end

function SCIPdialogFindEntry(dialog, entryname, subdialog)
    ccall((:SCIPdialogFindEntry, libscip), Cint, (Ptr{SCIP_DIALOG}, Cstring, Ptr{Ptr{SCIP_DIALOG}}), dialog, entryname, subdialog)
end

function SCIPdialogDisplayMenu(dialog, scip)
    ccall((:SCIPdialogDisplayMenu, libscip), SCIP_RETCODE, (Ptr{SCIP_DIALOG}, Ptr{SCIP_}), dialog, scip)
end

function SCIPdialogDisplayMenuEntry(dialog, scip)
    ccall((:SCIPdialogDisplayMenuEntry, libscip), SCIP_RETCODE, (Ptr{SCIP_DIALOG}, Ptr{SCIP_}), dialog, scip)
end

function SCIPdialogDisplayCompletions(dialog, scip, entryname)
    ccall((:SCIPdialogDisplayCompletions, libscip), SCIP_RETCODE, (Ptr{SCIP_DIALOG}, Ptr{SCIP_}, Cstring), dialog, scip, entryname)
end

function SCIPdialogGetPath(dialog, sepchar, path)
    ccall((:SCIPdialogGetPath, libscip), Cvoid, (Ptr{SCIP_DIALOG}, UInt8, Cstring), dialog, sepchar, path)
end

function SCIPdialogGetName(dialog)
    ccall((:SCIPdialogGetName, libscip), Cstring, (Ptr{SCIP_DIALOG},), dialog)
end

function SCIPdialogGetDesc(dialog)
    ccall((:SCIPdialogGetDesc, libscip), Cstring, (Ptr{SCIP_DIALOG},), dialog)
end

function SCIPdialogIsSubmenu(dialog)
    ccall((:SCIPdialogIsSubmenu, libscip), UInt32, (Ptr{SCIP_DIALOG},), dialog)
end

function SCIPdialogGetParent(dialog)
    ccall((:SCIPdialogGetParent, libscip), Ptr{SCIP_DIALOG}, (Ptr{SCIP_DIALOG},), dialog)
end

function SCIPdialogGetSubdialogs(dialog)
    ccall((:SCIPdialogGetSubdialogs, libscip), Ptr{Ptr{SCIP_DIALOG}}, (Ptr{SCIP_DIALOG},), dialog)
end

function SCIPdialogGetNSubdialogs(dialog)
    ccall((:SCIPdialogGetNSubdialogs, libscip), Cint, (Ptr{SCIP_DIALOG},), dialog)
end

function SCIPdialogGetData(dialog)
    ccall((:SCIPdialogGetData, libscip), Ptr{SCIP_DIALOGDATA}, (Ptr{SCIP_DIALOG},), dialog)
end

function SCIPdialogSetData(dialog, dialogdata)
    ccall((:SCIPdialogSetData, libscip), Cvoid, (Ptr{SCIP_DIALOG}, Ptr{SCIP_DIALOGDATA}), dialog, dialogdata)
end

function SCIPdialogWriteHistory(filename)
    ccall((:SCIPdialogWriteHistory, libscip), SCIP_RETCODE, (Cstring,), filename)
end


text:
module PSCOPFTest

using Test

using PSCOPF

include("test_wip//test_models/fixtures.jl")

println("Testing with solver : ", PSCOPF.OPTIMIZER)

#TODO : redirect logs
@testset verbose=true "PSCOPF_TESTS" begin
    include("test_wip/test_modes.jl")
    include("test_wip/test_target_ts.jl")
    include("test_wip/test_network.jl")
    include("test_wip/test_ech_generator.jl")
    include("test_wip/test_uncertainties_generator.jl")

    include("test_wip/test_coherence_checks/test_check_ptdf.jl")
    include("test_wip/test_coherence_checks/test_check_generator.jl")
    include("test_wip/test_coherence_checks/test_check_branch.jl")
    include("test_wip/test_coherence_checks/test_check_uncertainties.jl")
    include("test_wip/test_coherence_checks/test_check_initial_state.jl")
    include("test_wip/test_coherence_checks/test_check_timesteps.jl")

    include("test_wip/test_schedule.jl")
    include("test_wip/test_tso_actions.jl")
    include("test_wip/test_starts.jl")
    include("test_wip/test_flows.jl")
    include("test_wip/test_verify_firmness.jl")
    include("test_wip/test_verify_firmness_on_schedule.jl")
    include("test_wip/test_compute_firmness.jl")
    include("test_wip/test_init_firmness.jl")
    include("test_wip/test_custom_sequence.jl")
    include("test_wip/test_firmness_in_sequence.jl")
    include("test_wip/test_seq_generator.jl")

    include("test_wip/test_models/energy_market/test_constraints.jl")
    include("test_wip/test_models/energy_market/test_unit_priority.jl")
    include("test_wip/test_models/energy_market/test_start_cost.jl")
    include("test_wip/test_models/energy_market/test_dmo.jl")
    include("test_wip/test_models/energy_market/test_dp.jl")
    include("test_wip/test_models/energy_market/test_slacks.jl")
    include("test_wip/test_models/energy_market/test_capping.jl")

    include("test_wip/test_models/energy_market/test_energy_market_at_fo.jl")

    include("test_wip/test_models/energy_market/test_balance_market.jl")

    include("test_wip/test_models/tso_out_fo/test_dp_pilotable.jl")
    include("test_wip/test_models/tso_out_fo/test_dp_limitable.jl")
    include("test_wip/test_models/tso_out_fo/test_dmo.jl")
    include("test_wip/test_models/tso_out_fo/test_constraints.jl")
    include("test_wip/test_models/tso_out_fo/test_slacks.jl")
    include("test_wip/test_models/tso_out_fo/test_start_cost.jl")
    include("test_wip/test_models/tso_out_fo/test_unit_priority.jl")
    include("test_wip/test_models/tso_out_fo/test_limitation.jl")

    include("test_wip/test_models/tso_bilevel/test_bilevel_limitables.jl")
    include("test_wip/test_models/tso_bilevel/test_bilevel_limitables_scenarios.jl")
    include("test_wip/test_models/tso_bilevel/test_bilevel_limitables_timesteps.jl")
    include("test_wip/test_models/tso_bilevel/test_bilevel_pilotables_deltas.jl")
    include("test_wip/test_models/tso_bilevel/test_bilevel_pilotables.jl")
    include("test_wip/test_models/tso_bilevel/test_bilevel_pilotables_start.jl")
    include("test_wip/test_models/tso_bilevel/test_bilevel_pilotables_costing.jl")
    include("test_wip/test_models/tso_bilevel/test_bilevel_pilotables_scenarios.jl")

    include("test_wip/test_models/eod_assessment/test_simple_call.jl")
    include("test_wip/test_models/rso_assessment/test_simple_call.jl")

    # include("test_wip/test_seq_launch.jl")
    include("test_wip/example_usecase.jl")
    include("test_wip/example_usecase_from_folder.jl")

    include("test_wip/test_sequence/test_bileveltso_market.jl")

end

end #module PSCOPFTest

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using NFFT
using FFTW
import NFFT3
using DataFrames

FFTW.set_num_threads(Threads.nthreads())

### performance test ###
function nfft_performance_1()
  println("\n\n ##### nfft_performance_1 - simple ##### \n\n")

  timing = TimingStats()

  let m = 3, σ = 2.0
    @info "NFFT Performance Test 1D"
    let N = 2^19, M = N, x = rand(M) .- 0.5, fHat = rand(M)*1im
      for pre in [NFFT.LINEAR, NFFT.FULL]
        @info "* precomputation = $pre"
        p = plan_nfft(x, N; m, σ, precompute=pre, timing)
        fApprox = *(adjoint(p), fHat; timing)
        *(p, fApprox; timing)

        println(timing)
      end
    end

    @info "NFFT Performance Test 2D"
    let N = 1024, M = N*N, x2 = rand(2,M) .- 0.5, fHat = rand(M)*1im
      for pre in [NFFT.LINEAR, NFFT.FULL]
        @info "* precomputation = $pre"
        p = plan_nfft(x2, (N,N); m, σ, precompute=pre, timing)
        fApprox = *(adjoint(p), fHat; timing)
        *(p, fApprox; timing)

        println(timing)
      end
    end

    @info "NFFT Performance Test 3D"
    let N = 32, M = N*N*N, x3 = rand(3,M) .- 0.5, fHat = rand(M)*1im
      for pre in [NFFT.LINEAR, NFFT.FULL]
        @info "* precomputation = $pre"
        p = plan_nfft(x3, (N,N,N); m, σ, precompute=pre, timing)
        fApprox = *(adjoint(p), fHat; timing)
        *(p, fApprox; timing)

        println(timing)
      end
    end
  end
  return nothing
end

nfft_performance_1()

function nfft_performance_2(N = 64, M = N*N*N)
  println("\n\n ##### nfft_performance_2 - multithreading ##### \n\n")

  m = 3; σ = 2.0
  timing = TimingStats()

  let x = Float32.(rand(3,M) .- 0.5), fHat = ComplexF32.(rand(M)*1im)

    for pre in [NFFT.LINEAR, NFFT.FULL]
      for threading in [true, false]
        NFFT._use_threads[] = threading

        @info "* precomputation = $pre threading = $threading"
        p = plan_nfft(x, (N,N,N); m, σ, precompute=pre, timing, fftflags=FFTW.MEASURE)
        fApprox = *(adjoint(p), fHat; timing)
        *(p, fApprox; timing)

        println(timing)
      end
    end
  end
end

nfft_performance_2()
#nfft_performance_2(128,46_000)


text:
"""
    mutable struct Serialized{M <: MetabolicModel}
        m::Maybe{M}
        filename::String
    end

A meta-model that represents a model that is serialized on the disk. The
internal model will be loaded on-demand by using any accessor, or by calling
[`precache!`](@ref) directly.
"""
mutable struct Serialized{M} <: MetabolicModel where {M<:MetabolicModel}
    m::Maybe{M}
    filename::String

    Serialized{T}(filename::String) where {T} = new{T}(nothing, filename)
    Serialized(model::T, filename::String) where {T<:MetabolicModel} =
        new{T}(model, filename)
end

function _on_precached(m::Serialized, f, args...)
    precache!(m)
    f(m.m, args...)
end

reactions(m::Serialized) = _on_precached(m, reactions)
n_reactions(m::Serialized) = _on_precached(m, n_reactions)
metabolites(m::Serialized) = _on_precached(m, metabolites)
n_metabolites(m::Serialized) = _on_precached(m, n_metabolites)
stoichiometry(m::Serialized) = _on_precached(m, stoichiometry)
bounds(m::Serialized) = _on_precached(m, bounds)
balance(m::Serialized) = _on_precached(m, balance)
objective(m::Serialized) = _on_precached(m, objective)
coupling(m::Serialized) = _on_precached(m, coupling)
n_coupling_constraints(m::Serialized) = _on_precached(m, n_coupling_constraints)
coupling_bounds(m::Serialized) = _on_precached(m, coupling_bounds)
genes(m::Serialized) = _on_precached(m, genes)
n_genes(m::Serialized) = _on_precached(m, n_genes)
metabolite_formula(m::Serialized, id::String) = _on_precached(m, metabolite_formula, id)
metabolite_charge(m::Serialized, id::String) = _on_precached(m, metabolite_charge, id)
reaction_annotations(m::Serialized, id::String) = _on_precached(m, reaction_annotations, id)
metabolite_annotations(m::Serialized, id::String) =
    _on_precached(m, metabolite_annotations, id)
gene_annotations(m::Serialized, id::String) = _on_precached(m, gene_annotations, id)
reaction_notes(m::Serialized, id::String) = _on_precached(m, reaction_notes, id)
metabolite_notes(m::Serialized, id::String) = _on_precached(m, metabolite_notes, id)
gene_notes(m::Serialized, id::String) = _on_precached(m, gene_notes, id)
metabolite_compartment(m::Serialized, id::String) =
    _on_precached(m, metabolite_compartment, id)
reaction_subsystem(m::Serialized, id::String) = _on_precached(m, reaction_subsystem, id)

"""
    precache!(model::Serialized{MetabolicModel})::Nothing

Load the `Serialized` model from disk in case it's not alreadly loaded.
"""
function precache!(model::Serialized)::Nothing
    if isnothing(model.m)
        model.m = deserialize(model.filename)
    end
    nothing
end


text:
function main()
    game = readlines() .|> collect
    sy, sx = length(game), length(game[1])
    for iteration in 1:1000
        something_changed = false
        new_pos = Dict()
        for (dir, addx, addy) in [('>', 1, 0), ('v', 0, 1)]
            for y in 1:sy, x in 1:sx
                if game[y][x] == dir
                    ny, nx = mod1(y+addy, sy), mod1(x+addx, sx)
                    if (game[ny][nx] == '.' || dir == 'v' && game[ny][nx] == '>') && (ny, nx) ∉ keys(new_pos)
                        new_pos[(ny, nx)] = dir
                        something_changed = true
                    else
                        new_pos[(y, x)] = dir
                    end
                end
            end
        end
        game = [[get(new_pos, (y, x), '.') for x in 1:sx] for y in 1:sy]
        !something_changed && return println(iteration)
    end
end
main()

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
# augmented_dickey_fuller.jl
# (Augmented) Dickey-Fuller unit root test
#
# Copyright (C) 2017   Benjamin Born
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export ADFTest

struct ADFTest <: HypothesisTest
    n               ::Int               # number of observations
    deterministic   ::Symbol            # deterministic terms included in regression
    lag             ::Int               # number of lags in test statistic
    stat            ::Float64           # test statistic
    coef            ::Float64           # coefficient on lagged (non-differenced) variable
    cv              ::Vector{Float64}   # critical values
end

"""
    ADFTest(y, deterministic, lag)

Compute the augmented Dickey-Fuller unit root test.

`y` is the time series to be tested, `deterministic` determines the deterministic terms
(options: `none`, `constant`, `trend`, `squared_trend`) and `lag` the number of lagged
first-differences included in the test regression, respectively.

Critical values and asymptotic p-values are computed based on response surface regressions
following MacKinnon (2010) and MacKinnon (1994), respectively. These may differ slightly
from those reported in other regression packages as different algorithms might be used.

# References

  * James G. MacKinnon, 2010, "Critical values for cointegration tests,"
    QED Working Paper No. 1227, 2010, [link](http://ideas.repec.org/p/qed/wpaper/1227.html).
  * James G. MacKinnon, 1994, "Approximate Asymptotic Distribution Functions for
    Unit-Root and Cointegration Tests", Journal of Business & Economic Statistics,
    Vol. 12, No. 2, pp. 167-176, [link](http://www.jstor.org/stable/1391481).

# External links

  * [Augmented Dickey-Fuller test on Wikipedia](https://en.wikipedia.org/wiki/Augmented_Dickey–Fuller_test)
"""
function ADFTest(y::AbstractVector{T}, deterministic::Symbol, lag::Int) where T<:Real

    nobs = length(y)
    Δy = diff(y)

    if deterministic == :none
        Δyt   = Δy[lag+1:end]
        xt    = y[lag+1:end-1,1:1]
        dfexo = 0
    elseif deterministic == :constant
        Δyt   = Δy[lag+1:end]
        xt    = [ y[lag+1:end-1] ones(nobs-lag-1, 1) ]
        dfexo = 1
    elseif deterministic == :trend
        Δyt   = Δy[lag+1:end]
        xt    = [ y[lag+1:end-1] ones(nobs-lag-1, 1) (1:nobs-lag-1) ]
        dfexo = 2
    elseif deterministic == :squared_trend
        Δyt   = Δy[lag+1:end]
        xt    = [ y[lag+1:end-1] ones(nobs-lag-1, 1) (1:nobs-lag-1) [
                  i^2 for i = 1:nobs-lag-1] ]
        dfexo   = 3
    else
        throw(ArgumentError("deterministic = $(deterministic) is invalid"))
    end

    if lag > 0
        Δylag = hcat([Δy[lag-i+1:end-i] for i=1:lag]...)
        xt = [ xt Δylag ]
    end

    adf_coef = xt \ Δyt
    adf_res  = Δyt - xt * adf_coef
    n        = length(adf_res)
    k        = lag + 1 + dfexo
    sigma2   = dot(adf_res, adf_res) / (n - k)
    inv_xtxt = inv(xt' * xt)
    adf_stat = adf_coef[1] / sqrt(sigma2 * inv_xtxt[1,1])

    crit_val = adf_cv(nobs, deterministic)

    ADFTest(n, deterministic, lag, adf_stat, adf_coef[1], crit_val)

end

function adf_cv(nobs::Int, deterministic::Symbol)
    # computes critical values
    #
    # based on James G. MacKinnon, "Critical values for cointegration tests,"
    # QED Working Paper No. 1227, 2010. http://ideas.repec.org/p/qed/wpaper/1227.html

    cv_coeff = [
        # no constant, rows denote values for 1%, 5%, and 10% level, respectively
        -2.56574 -02.2358 -03.627    0.000
        -1.94100 -00.2686 -03.365   31.223
        -1.61682  00.2656 -02.714   25.364
        # constant
        -3.43035 -06.5393 -16.786 -079.433
        -2.86154 -02.8903 -04.234 -040.040
        -2.56677 -01.5384 -02.809    0.000
        # constant and trend
        -3.95877 -09.0531 -28.428 -134.155
        -3.41049 -04.3904 -09.036 -045.374
        -3.12705 -02.5856 -03.925 -022.380
        # constant, trend, and trend squared
        -4.37113 -11.5882 -35.819 -334.047
        -3.83239 -05.9057 -12.490 -118.284
        -3.55326 -03.6596 -05.293 -063.559
    ]

    # Equation A.1 in MacKinnon (2010)
    if deterministic == :none
        adf_crit_val = cv_coeff[1:3, 1] + cv_coeff[1:3, 2] ./ nobs +
                          cv_coeff[1:3, 3] ./ (nobs^2) + cv_coeff[1:3, 4] ./ (nobs^3)
    elseif deterministic == :constant
        adf_crit_val = cv_coeff[4:6, 1] + cv_coeff[4:6, 2] ./ nobs +
                          cv_coeff[4:6, 3] ./ (nobs^2) + cv_coeff[4:6, 4] ./ (nobs^3)
    elseif deterministic == :trend
        adf_crit_val = cv_coeff[7:9, 1] + cv_coeff[7:9, 2] ./ nobs +
                          cv_coeff[7:9, 3] ./ (nobs^2) + cv_coeff[7:9, 4] ./ (nobs^3)
    elseif deterministic == :squared_trend
        adf_crit_val = cv_coeff[10:12, 1] + cv_coeff[10:12, 2] ./ nobs +
                          cv_coeff[10:12, 3] ./ (nobs^2) + cv_coeff[10:12, 4] ./ (nobs^3)
    else
        throw(ArgumentError("deterministic = $(deterministic) is invalid"))
    end

end

function adf_pv_aux(adf_stat::Float64, deterministic::Symbol)
    # helper function for p-value computation
    #
    # based on James G. MacKinnon, "Approximate Asymptotic Distribution Functions for
    # Unit-Root and Cointegration Tests", Journal of Business & Economic Statistics,
    # Vol. 12, No. 2 (Apr., 1994), pp. 167-176, http://www.jstor.org/stable/1391481
    #
    # A slightly more accurate approach could be based on an algorithm described in
    # James G. MacKinnon, "Numerical Distribution Functions for Unit Root and Cointegration
    # Tests", Journal of Applied Econometrics, Vol. 11, No. 6 (Nov.-Dec., 1996),
    # pp. 601-618, http://www.jstor.org/stable/2285154
    # and made available under GPL as Fortran routines (used e.g. by R's urca package) at
    # http://qed.econ.queensu.ca/pub/faculty/mackinnon/numdist/

    pv_coeff_smallp = [
        0.6344 1.2378 0.032496 -1.04 -19.04
        2.1659 1.4412 0.038269 -1.61 -18.83
        3.2512 1.6047 0.049588 -2.89 -16.18
        4.0003 1.6580 0.048288 -3.21 -17.17
    ]

    pv_coeff_largep = [
        0.4797 0.93557 -0.06999  0.033066 -1.04 Inf
        1.7339 0.93202 -0.12745 -0.010368 -1.61 2.74
        2.5261 0.61654 -0.37956 -0.060285 -2.89 0.70
        3.0778 0.49529 -0.41477 -0.059359 -3.21 0.54
    ]

    if deterministic == :none
        tab_row = 1
    elseif deterministic == :constant
        tab_row = 2
    elseif deterministic == :trend
        tab_row = 3
    elseif deterministic == :squared_trend
        tab_row = 4
    else
        throw(ArgumentError("deterministic = $(deterministic) is invalid"))
    end

    if adf_stat < pv_coeff_smallp[tab_row, 5]
        aux_var = -Inf
    elseif adf_stat > pv_coeff_largep[tab_row, 6]
        aux_var = Inf
    else
        if adf_stat < pv_coeff_smallp[tab_row, 4]
            aux_var = pv_coeff_smallp[tab_row, 1] + pv_coeff_smallp[tab_row, 2] * adf_stat +
                        pv_coeff_smallp[tab_row, 3] * (adf_stat^2)
        else
            aux_var = pv_coeff_largep[tab_row, 1] + pv_coeff_largep[tab_row, 2] * adf_stat +
                        pv_coeff_largep[tab_row, 3] * (adf_stat^2) +
                        pv_coeff_largep[tab_row, 4] * (adf_stat^3)
        end

    end
end

testname(::ADFTest) = "Augmented Dickey-Fuller unit root test"
population_param_of_interest(x::ADFTest) =
    ("coefficient on lagged non-differenced variable", 0, x.coef)

function show_params(io::IO, x::ADFTest, ident)
    println(io, ident, "sample size in regression:          ", x.n)
    println(io, ident, "number of lags:                     ", x.lag)
    println(io, ident, "ADF statistic:                      ", x.stat)
    print(io, ident, "Critical values at 1%, 5%, and 10%: ")
    show(io, x.cv')
    println(io)
end

pvalue(x::ADFTest) = HypothesisTests.pvalue(Normal(0, 1), adf_pv_aux(x.stat, x.deterministic); tail=:left)


text:
using Test
using Printf
using MathOptInterface

const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

import GAMS

const OPTIMIZER_CONSTRUCTOR = MOI.OptimizerWithAttributes(GAMS.Optimizer, MOI.Silent() => true)
const OPTIMIZER = MOI.instantiate(OPTIMIZER_CONSTRUCTOR)
OPTIMIZER.gamswork = GAMS.GAMSWorkspace()

const CACHING_OPTIMIZER = MOIU.CachingOptimizer(MOIU.Model{Float64}(), OPTIMIZER);
const OPTIMIZER_SPLITINTERVAL = MOIB.Constraint.SplitInterval{Float64}(
    MOIU.CachingOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()), OPTIMIZER))
const BRIDGED = MOI.instantiate(OPTIMIZER_CONSTRUCTOR, with_bridge_type = Float64)

const TEST_SOLVERS = (
   "alphaecp",
   "antigone",
   "baron",
   "bdmlp",
   "bonmin",
   "bonminh",
   "cbc",
   "conopt",
   "conopt3",
   "conopt4",
   "couenne",
   "cplex",
   "dicopt",
   "glomiqo",
   "gurobi",
   "ipopt",
   "ipopth",
   "knitro",
   "lgo",
   "lindo",
   "lindoglobal",
   "localsolver",
   "localsolver70",
   "minos",
   "mosek",
   "msnlp",
   "path",
   "pathc",
   "pathnlp",
   "quadminos",
   "sbb",
   "shot",
   "scip",
   "snopt",
   "soplex",
   "xa",
   "xpress"
)

atol = 1e-5
rtol = 1e-5
const CONFIG = MOIT.TestConfig(atol=atol, rtol=rtol)
const CONFIG_LOCAL = MOIT.TestConfig(atol=atol, rtol=rtol, optimal_status = MOI.LOCALLY_SOLVED)
const CONFIG_LOCAL_NODUAL = MOIT.TestConfig(atol=atol, rtol=rtol, duals = false, optimal_status = MOI.LOCALLY_SOLVED)
const CONFIG_NODUAL = MOIT.TestConfig(atol=atol, rtol=rtol, duals = false)

println("Unit Tests" * "_"^19)
testname = @sprintf("%-15s", "default")
@testset "$testname" begin
    @test MOI.get(OPTIMIZER, MOI.SolverName()) == "GAMS"
    @test MOIU.supports_default_copy_to(OPTIMIZER, false)
    @test !MOIU.supports_default_copy_to(OPTIMIZER, true)
    exclude = [
        "delete_variable",                          # deleting not supported
        "delete_variables",                         # deleting not supported
        "solve_affine_deletion_edge_cases",         # deleting not supported
        "solve_affine_lessthan",                    # constraint names for variable bounds not supported
        "solve_objbound_edge_cases",                # constraint names for variable bounds not supported
        "getconstraint",                            # constraint names for variable bounds not supported
        "getvariable",                              # constraint names for variable bounds not supported
        "solve_with_upperbound",                    # constraint names for variable bounds not supported
        "solve_blank_obj",                          # constraint names for variable bounds not supported
        "solve_affine_equalto",                     # constraint names for variable bounds not supported
        "solve_single_variable_dual_min",           # constraint names for variable bounds not supported
        "solve_integer_edge_cases",                 # constraint names for variable bounds not supported
        "solve_constant_obj",                       # constraint names for variable bounds not supported
        "solve_singlevariable_obj",                 # constraint names for variable bounds not supported
        "solve_with_lowerbound",                    # constraint names for variable bounds not supported
        "solve_zero_one_with_bounds_1",             # constraint names for variable bounds not supported
        "solve_zero_one_with_bounds_2",             # constraint names for variable bounds not supported
        "solve_zero_one_with_bounds_3",             # constraint names for variable bounds not supported
        "get_objective_function",                   # function getters not supported
        "update_dimension_nonnegative_variables",   # function getters not supported
        "delete_nonnegative_variables",             # function getters not supported
        "solve_result_index",                       # DualObjectiveValue not supported
        "delete_soc_variables",                     # second order cone not supported
        "solve_affine_interval",                    # get constraint index not supported
        "solve_affine_greaterthan",                 # get constraint index not supported
        "solve_qp_edge_cases",                      # conopt finds only local optimal solution
        "solve_qcp_edge_cases",                     # conopt finds only local optimal solution
        "solve_qp_zero_offdiag",                    # conopt finds only local optimal solution
        "solve_farkas_equalto_upper",               # definition of duals for infeasible models differs
        "solve_farkas_equalto_lower",               # definition of duals for infeasible models differs
        "solve_farkas_variable_lessthan",           # definition of duals for infeasible models differs
        "solve_farkas_variable_lessthan_max",       # definition of duals for infeasible models differs
        "solve_farkas_interval_upper",              # definition of duals for infeasible models differs
        "solve_farkas_lessthan",                    # definition of duals for infeasible models differs
    ]
    config = MOIT.TestConfig(atol=1e-3, rtol=1e-3)
    MOIT.unittest(BRIDGED, config, exclude)
end
println()

println("Continuous Linear" * "_"^12)
for solver in TEST_SOLVERS
    if ! GAMS.check_solver(OPTIMIZER.gamswork, solver, GAMS.MODEL_TYPE_LP)
        continue
    end
    testname_s = @sprintf("%-15s", "$solver")
    MOI.set(OPTIMIZER_SPLITINTERVAL, GAMS.Solver(), solver)
    exclude = [
        "linear7",  # VectorAffineFunction not supported
        "linear15", # VectorAffineFunction not supported
        "linear10", # DualObjectiveValue not supported
        "linear2",  # DualObjectiveValue not supported
        "linear1",  # DualObjectiveValue not supported
        "linear14", # DualObjectiveValue not supported
        "linear12", # dual behavior in infeasible case doesn't match test
        "linear8a", # dual behavior in infeasible case doesn't match test
    ]
    @testset "$testname_s" begin
        if solver in ("xa",)
            exclude_solver = [
                "linear8c", # wrong primal variables
            ]
            MOIT.contlineartest(OPTIMIZER_SPLITINTERVAL, CONFIG, union(exclude, exclude_solver))

        elseif solver in ("knitro", "lgo")
            exclude_solver = [
                "linear8b", # unable to detect unboundedness
                "linear8c", # unable to detect unboundedness
            ]
            MOIT.contlineartest(OPTIMIZER_SPLITINTERVAL, CONFIG, union(exclude, exclude_solver))

        elseif solver in ("conopt4",)
            exclude_solver = [
                "partial_start", # wrong objective (bug)
                "linear9",       # wrong objective (bug)
            ]
            MOIT.contlineartest(OPTIMIZER_SPLITINTERVAL, CONFIG, union(exclude, exclude_solver))

        elseif solver in ("ipopt", "ipopth")
            exclude_solver = [
                "linear13", # wrong dual
                "linear8c", # does not terminated unbounded
            ]
            MOIT.contlineartest(OPTIMIZER_SPLITINTERVAL, CONFIG, union(exclude, exclude_solver))

        elseif solver in ("localsolver", "localsolver70")
            # needs extra license

        else
            MOIT.contlineartest(OPTIMIZER_SPLITINTERVAL, CONFIG, exclude)
        end
    end
end
println()

MOI.empty!(OPTIMIZER)

println("Integer Linear" * "_"^15)
MOI.set(OPTIMIZER, GAMS.ModelType(), "MIP")
for solver in TEST_SOLVERS
    if ! GAMS.check_solver(OPTIMIZER.gamswork, solver, GAMS.MODEL_TYPE_MIP)
        continue
    end
    testname_s = @sprintf("%-15s", "$solver")
    MOI.set(OPTIMIZER, GAMS.Solver(), solver)
    exclude = [
        "indicator1",   # indicator constraints not supported
        "indicator2",   # indicator constraints not supported
        "indicator3",   # indicator constraints not supported
        "indicator4",   # indicator constraints not supported
    ]
    @testset "$testname_s" begin
        if solver in ("baron", "mosek")
            exclude_solver = [
                "semiconttest", # semicont constraints not supported
                "semiinttest",  # semiint constraints not supported
                "int2",         # sos constraints not supported
            ]
            MOIT.intlineartest(OPTIMIZER_SPLITINTERVAL, CONFIG, union(exclude, exclude_solver))

        elseif solver in ("lindo", "lindoglobal")
            exclude_solver = [
                "semiconttest", # semicont constraints not supported
                "semiinttest",  # semiint constraints not supported
                "int2",         # finds different optimal solution
            ]
            MOIT.intlineartest(OPTIMIZER_SPLITINTERVAL, CONFIG, union(exclude, exclude_solver))

        elseif solver in ("bdmlp", "cbc", "cplex", "gurobi")
            exclude_solver = [
                "semiconttest", # does not provide objective bound
                "semiinttest",  # does not provide objective bound
            ]
            MOIT.intlineartest(OPTIMIZER_SPLITINTERVAL, CONFIG, union(exclude, exclude_solver))

        elseif solver in ("xa",)
            exclude_solver = [
                "semiconttest", # does not provide objective bound
                "semiinttest",  # does not provide objective bound
                "int2",         # sos constraints not supported
            ]
            MOIT.intlineartest(OPTIMIZER_SPLITINTERVAL, CONFIG, union(exclude, exclude_solver))

        elseif solver in ("scip",)
            exclude_solver = [
                "int1",         # finds suboptimal solution
                "semiconttest", # does not provide objective bound
                "semiinttest",  # does not provide objective bound
            ]
            MOIT.intlineartest(OPTIMIZER_SPLITINTERVAL, CONFIG, union(exclude, exclude_solver))

        elseif solver in ("xpress",)
            exclude_solver = [
                "int3",         # finds suboptimal solution
                "semiconttest", # does not provide objective bound
                "semiinttest",  # does not provide objective bound
            ]
            MOIT.intlineartest(OPTIMIZER_SPLITINTERVAL, CONFIG, union(exclude, exclude_solver))

        elseif solver in ("localsolver", "localsolver70")
            # needs extra license

        else
            MOIT.intlineartest(OPTIMIZER_SPLITINTERVAL, CONFIG, exclude)
        end
    end
end
println()

MOI.empty!(OPTIMIZER)

println("Continuous Quadratic" * "_"^9)
MOI.set(OPTIMIZER, GAMS.ModelType(), "QCP")
for solver in TEST_SOLVERS
    if ! GAMS.check_solver(OPTIMIZER.gamswork, solver, GAMS.MODEL_TYPE_QCP)
        continue
    end
    testname_s = @sprintf("%-15s", "$solver")
    MOI.set(OPTIMIZER, GAMS.Solver(), solver)
    exclude = [
        "qcp1",   # VectorAffineFunction not supported
    ]
    @testset "$testname_s" begin
        if solver in ("antigone", "glomiqo")
            exclude_solver = [
                "ncqcp1", # finds only local optimal solution
                "socp1",  # finds only local optimal solution
            ]
            MOIT.contquadratictest(CACHING_OPTIMIZER, CONFIG_NODUAL, union(exclude, exclude_solver))

        elseif solver in ("conopt", "conopt3", "dicopt", "minos", "snopt", "pathnlp")
            exclude_solver = [
                "socp1",  # fails to find feasible solution
                "ncqcp2", # fails to find feasible solution
                "qp3",    # manages to find global optimal solution
            ]
            MOIT.contquadratictest(CACHING_OPTIMIZER, CONFIG_LOCAL_NODUAL, union(exclude, exclude_solver))

        elseif solver in ("baron",)
            exclude_solver = [
                "qp1",    # finds only local optimal solution
                "qp3",    # finds only local optimal solution
                "qp2",    # fails to prove optimality
                "qcp2",   # fails to prove optimality
                "qcp3",   # fails to prove optimality
                "qcp4",   # fails to prove optimality
                "qcp5",   # fails to prove optimality
                "socp1",  # fails to prove optimality
                "ncqcp1", # fails to prove optimality
            ]
            MOIT.contquadratictest(CACHING_OPTIMIZER, CONFIG_NODUAL, union(exclude, exclude_solver))

        elseif solver in ("couenne",)
            exclude_solver = [
                "qp1",    # finds different solution
                "qp2",    # takes too much time
                "qcp3",   # fails to prove optimality
                "qp3",    # fails to prove optimality
            ]
            MOIT.contquadratictest(CACHING_OPTIMIZER, CONFIG_NODUAL, union(exclude, exclude_solver))

        elseif solver in ("cplex",)
            exclude_solver = [
                "socp1",  # not positive semi-definite
                "ncqcp1", # not positive semi-definite
                "ncqcp2", # not positive semi-definite
            ]
            cplex_config = MOIT.TestConfig(atol=1e-3, rtol=1e-3, duals = false)
            MOIT.contquadratictest(CACHING_OPTIMIZER, cplex_config, union(exclude, exclude_solver))

        elseif solver in ("conopt4",)
            exclude_solver = [
                "qp3",    # wrong objective (bug)
                "ncqcp1", # wrong objective (bug)
                "ncqcp2", # manages to find global optimal solution
            ]
            MOIT.contquadratictest(CACHING_OPTIMIZER, CONFIG_LOCAL_NODUAL, union(exclude, exclude_solver))

        elseif solver in ("gurobi",)
            exclude_solver = [
                "ncqcp1", # not positive semi-definite
            ]
            MOIT.contquadratictest(CACHING_OPTIMIZER, CONFIG_NODUAL, union(exclude, exclude_solver))

        elseif solver in ("lindo", "lindoglobal")
            exclude_solver = [
                "ncqcp1", # not positive semi-definite
                "ncqcp2", # only local optimal
            ]
            MOIT.contquadratictest(CACHING_OPTIMIZER, CONFIG_NODUAL, union(exclude, exclude_solver))

        elseif solver in ("ipopt", "ipopth", "knitro", "lgo", "msnlp")
            exclude_solver = [
                "qp3",    # manages to find global optimal solution
            ]
            MOIT.contquadratictest(CACHING_OPTIMIZER, CONFIG_LOCAL_NODUAL, union(exclude, exclude_solver))

        elseif solver in ("mosek", "xpress")
            exclude_solver = [
                "ncqcp1", # not supported
                "ncqcp2", # not supported
            ]
            MOIT.contquadratictest(CACHING_OPTIMIZER, CONFIG_NODUAL, union(exclude, exclude_solver))

        elseif solver in ("scip",)
            MOIT.contquadratictest(CACHING_OPTIMIZER, CONFIG_NODUAL, exclude)

        elseif solver in ("localsolver", "localsolver70")
            # needs extra license

        else
            MOIT.contquadratictest(CACHING_OPTIMIZER, CONFIG, exclude)
        end
    end
end
println()

MOI.empty!(OPTIMIZER)

println("Continuous Nonlinear" * "_"^9)
MOI.set(OPTIMIZER, GAMS.ModelType(), "NLP")
for solver in TEST_SOLVERS
    if ! GAMS.check_solver(OPTIMIZER.gamswork, solver, GAMS.MODEL_TYPE_NLP)
        continue
    end
    testname_s = @sprintf("%-15s", "$solver")
    MOI.set(OPTIMIZER, GAMS.Solver(), solver)
    @testset "$testname_s" begin
        if solver in ("antigone", "baron", "couenne")
            exclude = [
                "hs071",            # finds only local optimal solution / fails to prove optimality
                "hs071_no_hessian", # finds only local optimal solution / fails to prove optimality
            ]
            MOIT.nlptest(OPTIMIZER, CONFIG, exclude)

        elseif solver in ("scip",)
            exclude = [
                "hs071",            # fails to prove optimality
                "hs071_no_hessian", # fails to prove optimality
                "nlp_objective_and_moi_objective", # finds different solution
            ]
            MOIT.nlptest(OPTIMIZER, CONFIG, exclude)

        elseif solver in ("mosek",)
            # not supported

        elseif solver in ("lindo", "lindoglobal")
            exclude = [
                "hs071",            # finds global optimal solution
                "hs071_no_hessian", # finds global optimal solution
            ]
            MOIT.nlptest(OPTIMIZER, CONFIG_LOCAL, exclude)

        elseif solver in ("xpress",)
            exclude = [
                "hs071",            # finds different point
                "hs071_no_hessian", # finds different point
            ]
            MOIT.nlptest(OPTIMIZER, CONFIG_LOCAL, exclude)

        elseif solver in ("localsolver", "localsolver70")
            # needs extra license

        else
            MOIT.nlptest(OPTIMIZER, CONFIG_LOCAL)
        end
    end
end

return


text:
using LinearAlgebra

function deigen(x)
    copyto!(x, Symmetric(x, :U))
    xd = Elemental.DistMatrix(eltype(x))
    Elemental.zeros!(xd, size(x)...)
    copy_local_symm!(xd, x)
    @debug @assert x ≈ Array(xd)
    w, Z = Elemental.eigHermitian(Elemental.UPPER, xd)
    vals = vec(Array(w))
    vecs = fill!(similar(x), 0)
    copy_local_symm!(vecs, Z)
    vecs = permutedims(vecs)
    @debug @assert vecs ≈ Array(Z)
    return vals, vecs
end

function copy_local_symm!(xd::Elemental.DistMatrix, x)
    is = part(1:size(x, 1))
    for i in is, j in 1:size(x, 2)
        Elemental.queueUpdate(xd, i, j, x[i, j])
    end
    Elemental.processQueues(xd)
    return xd
end

function copy_local_symm!(x, xd::Elemental.DistMatrix)
    is = part(1:size(x, 1))
    xl = view(x, :, is)
    for i in is, j in 1:size(x, 2)
        Elemental.queuePull(xd, i, j)
    end
    Elemental.processPullQueue(xd, unsafe_wrap(Array, pointer(xl), size(xl)))
    counts = MPI.Allgather(Cint(length(xl)), worldcomm())
    MPI.Allgatherv!(vec(x), counts, worldcomm())
    return x
end


text:
# These constants were manually copied from the file mpi.h in the Microsoft
# MPI SDK v7

const MPI_BYTE = Int32(0x4c00010d)
const MPI_WCHAR = Int32(0x4c00020e)
const MPI_INT8_T = Int32(0x4c000133)
const MPI_UINT8_T = Int32(0x4c000137)
const MPI_INT16_T = Int32(0x4c000234)
const MPI_UINT16_T = Int32(0x4c000238)
const MPI_INT32_T = Int32(0x4c000435)
const MPI_UINT32_T = Int32(0x4c000439)
const MPI_INT64_T = Int32(0x4c000836)
const MPI_UINT64_T = Int32(0x4c00083a)
const MPI_REAL4 = Int32(0x4c000427)
const MPI_REAL8 = Int32(0x4c000829)
const MPI_COMPLEX8 = Int32(0x4c000828)
const MPI_COMPLEX16 = Int32(0x4c00102a)
const MPI_INTEGER1 = Int32(0x4c00012d)
const MPI_INTEGER2 = Int32(0x4c00022f)
const MPI_INTEGER4 = Int32(0x4c000430)
const MPI_INTEGER8 = Int32(0x4c000831)
const MPI_COMM_NULL = Int32(67108864)
const MPI_COMM_SELF = Int32(1140850689)
const MPI_COMM_WORLD = Int32(1140850688)
const MPI_OP_NULL = Int32(402653184)
const MPI_BAND = Int32(1476395014)
const MPI_BOR = Int32(1476395016)
const MPI_BXOR = Int32(1476395018)
const MPI_LAND = Int32(1476395013)
const MPI_LOR = Int32(1476395015)
const MPI_LXOR = Int32(1476395017)
const MPI_MAX = Int32(1476395009)
const MPI_MAXLOC = Int32(1476395020)
const MPI_MIN = Int32(1476395010)
const MPI_MINLOC = Int32(1476395019)
const MPI_PROD = Int32(1476395012)
const MPI_REPLACE = Int32(1476395021)
const MPI_SUM = Int32(1476395011)
const MPI_REQUEST_NULL = Int32(738197504)
const MPI_STATUS_SIZE = Int32(5)
const MPI_ERROR = Int32(5)
const MPI_SOURCE = Int32(3)
const MPI_TAG = Int32(4)
const MPI_ANY_SOURCE = Int32(-2)
const MPI_ANY_TAG = Int32(-1)
const MPI_TAG_UB = Int32(1681915906)
const MPI_UNDEFINED = Int32(-32766)
const HAVE_MPI_COMM_C2F = false

const MPI_INIT = (:MPI_INIT, "msmpi.dll")
const MPI_COMM_RANK = (:MPI_COMM_RANK, "msmpi.dll")
const MPI_COMM_SIZE = (:MPI_COMM_SIZE, "msmpi.dll")
const MPI_BARRIER = (:MPI_BARRIER, "msmpi.dll")
const MPI_FINALIZE = (:MPI_FINALIZE, "msmpi.dll")
const MPI_BCAST = (:MPI_BCAST, "msmpi.dll")
const MPI_REDUCE = (:MPI_REDUCE, "msmpi.dll")
const MPI_IRECV = (:MPI_IRECV, "msmpi.dll")
const MPI_RECV = (:MPI_RECV, "msmpi.dll")
const MPI_ISEND = (:MPI_ISEND, "msmpi.dll")
const MPI_WAITALL = (:MPI_WAITALL, "msmpi.dll")
const MPI_ALLGATHER = (:MPI_ALLGATHER, "msmpi.dll")
const MPI_ALLGATHERV = (:MPI_ALLGATHERV, "msmpi.dll")
const MPI_ALLTOALL = (:MPI_ALLTOALL, "msmpi.dll")
const MPI_ALLTOALLV = (:MPI_ALLTOALLV, "msmpi.dll")
const MPI_INITIALIZED = (:MPI_INITIALIZED, "msmpi.dll")
const MPI_FINALIZED = (:MPI_FINALIZED, "msmpi.dll")
const MPI_OP_CREATE = (:MPI_OP_CREATE, "msmpi.dll")
const MPI_OP_FREE = (:MPI_OP_FREE, "msmpi.dll")
const MPI_SCATTER = (:MPI_SCATTER, "msmpi.dll")
const MPI_SCATTERV = (:MPI_SCATTERV, "msmpi.dll")
const MPI_SEND = (:MPI_SEND, "msmpi.dll")
const MPI_SCAN = (:MPI_SCAN, "msmpi.dll")
const MPI_EXSCAN = (:MPI_EXSCAN, "msmpi.dll")
const MPI_GATHER = (:MPI_GATHER, "msmpi.dll")
const MPI_GATHERV = (:MPI_GATHERV, "msmpi.dll")
const MPI_COMM_DUP = (:MPI_COMM_DUP, "msmpi.dll")
const MPI_IPROBE = (:MPI_IPROBE, "msmpi.dll")
const MPI_PROBE = (:MPI_PROBE, "msmpi.dll")
const MPI_COMM_FREE = (:MPI_COMM_FREE, "msmpi.dll")
const MPI_GET_COUNT = (:MPI_GET_COUNT, "msmpi.dll")
const MPI_TEST = (:MPI_TEST, "msmpi.dll")
const MPI_TESTSOME = (:MPI_TESTSOME, "msmpi.dll")
const MPI_TESTANY = (:MPI_TESTANY, "msmpi.dll")
const MPI_TESTALL = (:MPI_TESTALL, "msmpi.dll")
const MPI_TYPE_CREATE_STRUCT = (:MPI_TYPE_CREATE_STRUCT, "msmpi.dll")
const MPI_TYPE_COMMIT = (:MPI_TYPE_COMMIT, "msmpi.dll")
const MPI_WAIT = (:MPI_WAIT, "msmpi.dll")
const MPI_WAITSOME = (:MPI_WAITSOME, "msmpi.dll")
const MPI_WAITANY = (:MPI_WAITANY, "msmpi.dll")
const MPI_CANCEL = (:MPI_CANCEL, "msmpi.dll")
const MPI_ALLREDUCE = (:MPI_ALLREDUCE, "msmpi.dll")

bitstype 32 CComm


text:
module aux

function logMeanExp( x )
  xmax = maximum(x)
  return xmax + log.(mean(exp.(x - xmax)))
end

function proximitySequence( x )
  m = length(x)
  seq = collect(1:m)
  for i = 1:m-1
    j = i + indmin(abs.(x[seq[i+1:m]] - x[seq[i]]))
    seq[[i+1,j]] = seq[[j,i+1]]
  end
  return seq[2:end]
end

function SumOfDeviationsPerBlock( y::Matrix{T}, ym::Matrix{T}, b::Integer ) where T<:AbstractFloat
  (m,n) = size(y)
  Δy = y .- ym
  B = Matrix{T}(m-b+1,n)
  B[1,:] = sum(Δy[1:b,:],1)
  for j = 1:m-b
    B[j+1,:] = B[j,:] + Δy[j+b,:] - Δy[j,:]
  end
  return B
end

function compute!( u0, P, π, f, u )
  g = (f + log.(π))'
  for i = 1:length(P)
    a = g .- u[i]
    max = maximum(a,2)
    numer = exp.(a .- max)
    denom = sum(numer,2)
    P[i] = numer ./ denom
    u0[i] = -(max + log.(denom))
  end
end

function overlapSampling( u::Vector{Matrix{Float64}} )
  m = length(u)
  f = zeros(m)
  seq = aux.proximitySequence( [mean(u[i][:,i]) for i=1:m] )
  i = 1
  for j in seq
    f[j] = f[i] + aux.logMeanExp(0.5(u[j][:,j] - u[j][:,i])) -
                  aux.logMeanExp(0.5(u[i][:,i] - u[i][:,j]))
    i = j
  end
  return f
end

function info( msg::String, val )
  const msg_color = "\033[1;36m"
  const val_color = "\033[0;36m"
  const no_color  = "\033[0m"
  print( STDOUT, msg_color, msg )
  if (isa(val,AbstractArray))
    println( STDOUT, val_color )
    Base.showarray( STDOUT, val, false; header=false )
  elseif (val != nothing)
    print( STDOUT, val_color, val )
  end
  println( STDOUT, no_color )
end

str(a) = join( string.(a), ", " )

end


text:
@testset "variable" begin
    @test x == x
    @test x ≠ y
end

@testset "free" begin
    @test a == a
    @test @term(a) == @term(a)
    @test a ≠ b
    @test f(a, a) == f(a, a)
    @test f(a, b) == f(a, b)
    @test f(a, b) ≠ f(b, a)
    @test f(x, b) == f(x, b)
    @test f(x, b) ≠ f(b, x)
end

@testset "c" begin
    @test p(a, b) == p(a, b)
    @test p(a, b) == p(b, a)
    @test p(a, b) ≠ p(a, c)
    @test p(p(a, b), c) ≠ p(a, p(b, c))
    @test p(a, p(b, c)) == p(p(c, b), a)
    @test p(p(a, b), q(a, b)) == p(q(a, b), p(a, b))
    @test p(p(x, a), q(x, b)) ≠ p(p(b, x), q(x, a))
    @test p(q(x, a), q(x, b)) == p(q(b, x), q(x, a))
    @test p(f(a, b), f(a, c)) == p(f(a, c), f(a, b))
end

@testset "ac" begin
    @test s(a, b) == s(a, b)
    @test s(a, b) == s(b, a)
    @test s(a, b, c) == s(c, a, b)
    @test s(a, x, b) == s(b, a, x)
    @test s(a, s(b, c)) == s(a, b, c)
    @test s(a, t(b, c)) ≠ s(a, b, c)
    @test s(a, s(b, c, a), s(c)) == s(a, a, b, c, c)
end


text:
struct WalkerIdealParam <: EoSParam
    Mw::SingleParam{Float64}
    Nrot::SingleParam{Int}
    theta1::SingleParam{Float64}
    theta2::SingleParam{Float64}
    theta3::SingleParam{Float64}
    theta4::SingleParam{Float64}
    deg1::SingleParam{Int}
    deg2::SingleParam{Int}
    deg3::SingleParam{Int}
    deg4::SingleParam{Int}
end

abstract type WalkerIdealModel <: IdealModel end
@newmodelgc WalkerIdeal WalkerIdealModel WalkerIdealParam

export WalkerIdeal
function WalkerIdeal(components::Array{String,1}; userlocations::Array{String,1}=String[], verbose=false)
    groups = GroupParam(components,["ideal/WalkerIdeal_Groups.csv"], verbose=verbose)
    params = getparams(groups, ["ideal/WalkerIdeal.csv"]; userlocations=userlocations, verbose=verbose)
    Mw = params["Mw"]
    Nrot = params["Nrot"]
    theta1 = params["theta_1"]
    theta2 = params["theta_2"]
    theta3 = params["theta_3"]
    theta4 = params["theta_4"]
    deg1 = params["deg_1"]
    deg2 = params["deg_2"]
    deg3 = params["deg_3"]
    deg4 = params["deg_4"]
    packagedparams = WalkerIdealParam(Mw, Nrot, theta1, theta2, theta3, theta4, deg1, deg2, deg3, deg4)
    references = ["10.1021/acs.jced.0c00723"]
    model = WalkerIdeal(packagedparams, groups, BasicIdeal, references=references, verbose=verbose)
    return model
end

function a_ideal(model::WalkerIdealModel,V,T,z)
    Mw = model.params.Mw.values
    Nrot = model.params.Nrot.values
    θ1 = model.params.theta1.values
    θ2 = model.params.theta2.values
    θ3 = model.params.theta3.values
    θ4 = model.params.theta4.values
    g1 = model.params.deg1.values
    g2 = model.params.deg2.values
    g3 = model.params.deg3.values
    g4 = model.params.deg4.values
    θ_vib = [θ1, θ2, θ3, θ4]
    g_vib = [g1, g2, g3, g4]
    n = model.groups.n_flattenedgroups
    res = zero(V+T+first(z))
    Σz = sum(z)
    x = z/Σz
    @inbounds for i in @comps
        ni = n[i]
        Mwi = sum(ni[k]*Mw[k] for k in @groups(i))
        Nroti = sum(ni[k]*Nrot[k] for k in @groups(i))/sum(ni[k] for k in @groups(i))
        Λ   = h/√(k_B*T*Mwi/N_A)
        res +=x[i]*(log(z[i]*N_A/V*Λ^3)-Nroti/2*log(T)+sum(ni[k]*sum(g_vib[v][k]*(θ_vib[v][k]/2/T+log(1-exp(-θ_vib[v][k]/T))) for v in 1:4) for k in @groups(i)))
    end
    return res - 1.
end


text:
using Printf
using LinearAlgebra
using Random

include("../common/PWGrid_v01.jl")
include("../common/wrappers_fft.jl")
include("../common/ortho_gram_schmidt.jl")
include("../common/gen_lattice_pwscf.jl")

include("gen_dr.jl")
include("init_pot_harm_3d.jl")
include("op_K.jl")
include("op_Vpot.jl")
include("op_H.jl")

include("calc_rho.jl")
include("calc_Etot.jl")
include("Kprec.jl")
include("diag_lobpcg.jl")

function test_main( ns1::Int,ns2::Int,ns3::Int )
    #
    Ns = [ns1,ns2,ns3]
    LatVecs = gen_lattice_sc(6.0)
    #
    pw = PWGrid( Ns, LatVecs )
    #
    Npoints = pw.Npoints
    Ω = pw.Ω
    r = pw.r
    G = pw.G
    G2 = pw.G2
    #
    # Generate array of distances
    #
    center = sum(LatVecs,dims=2)/2
    dr = gen_dr( r, center )
    #
    # Setup potential
    #
    Vpot = init_pot_harm_3d( pw, dr )
    println("sum(Vpot)*Ω/Npoints = ", sum(Vpot)*Ω/Npoints);
    #
    Nstates = 4
    Random.seed!(2222)
    psi = randn(ComplexF64,Npoints,Nstates)
    psi = ortho_gram_schmidt(psi)
    #
    evals, evecs = diag_lobpcg( pw, Vpot, psi, verbose=true, tol_avg=1e-10 )

    for ist = 1:Nstates
        @printf("State %3d, Energy = %18.10f\n", ist, real(evals[ist]))
    end
end

@time test_main( 30,30,30 )