text:
@info("loading libraries...")
using NBodyIPs, StaticArrays, BenchmarkTools, JuLIP, Test, LinearAlgebra

using NBodyIPs.Polys
using NBodyIPs.EnvIPs
using NBodyIPs.EnvIPs: EnvPoly

# TRANSFORM = PolyTransform(3, 2.9)
# DEGREES = [14, 12, 10]
# ENVDEGS = [3, 2, 1]
# r0 = rnn(:W)
# RCUT = [3.1*r0, 2.5*r0, 2.1*r0]
# RCUT = [2.6*r0, 2.2*r0, 1.8*r0]

## fast tests
TRANSFORM = PolyTransform(3, 2.9)
DEGREES = [14, 12, 10]
ENVDEGS = [3, 2, 1]
r0 = rnn(:W)
RCUT = [2.6*r0, 2.2*r0, 1.8*r0]

function random_ip(DT)
   DD = [ DT(TRANSFORM, CosCut(rcut-1.0, rcut))  for rcut in RCUT ]
   rn = 1.45 * rnn(:W)
   Vn = ("1.0/(1.0 + exp(1.0 / ($rn - r + 1e-2)))", rn)
   BB = [ envpolys(n+1, DD[n], DEGREES[n], Vn, ENVDEGS[n]) for n = 1:3 ]
   B = vcat(BB...)
   c = rand(length(B))
   IP = NBodyIP(B, c)
   return IP
end

##
@info("Generate random EnvIP...")
IP = random_ip(BondLengthDesc)
IPf = fast(IP)
IPcl = random_ip(ClusterBLDesc)
IPclf = fast(IPcl)

@info("Generate faster EnvIP => EnvPoly ...")
V2 = IPf.components[1:4]
V3 = IPf.components[5:7]
V4 = IPf.components[8:9]
IPff = NBodyIP( [ EnvPoly([v.Vr for v in V2], V2[1].Vn, V2[1].str_Vn),
                  EnvPoly([v.Vr for v in V3], V3[1].Vn, V3[1].str_Vn),
                  EnvPoly([v.Vr for v in V4], V4[1].Vn, V4[1].str_Vn) ] )

IPff2 = NBodyIPs.Experimental.faster_envpot(IP)

V2cl = IPclf.components[1:4]
V3cl = IPclf.components[5:7]
V4cl = IPclf.components[8:9]
IPclff = NBodyIP( [ EnvPoly([v.Vr for v in V2cl], V2cl[1].Vn, V2cl[1].str_Vn),
                    EnvPoly([v.Vr for v in V3cl], V3cl[1].Vn, V3cl[1].str_Vn),
                    EnvPoly([v.Vr for v in V4cl], V4cl[1].Vn, V4cl[1].str_Vn) ] )

@info("Test Correctness...")
at = rattle!(bulk(:W, cubic=true) * 3, 0.02)
@show energy(IP, at)
@show energy(IP, at) - energy(IPf, at)
@show energy(IPff, at) - energy(IPf, at)
@show energy(IPff2, at) - energy(IPf, at)
@show energy(IPcl, at)
@show energy(IPcl, at) - energy(IPclf, at)
# @show energy(IPclff, at) - energy(IPclf, at)

@info("Test evaluation time...")
at = rattle!(bulk(:W, cubic=true) * 8, 0.02)
print("nlist : ")
neighbourlist(at, 3.0); @time begin
      neighbourlist(at, RCUT[1]);
      neighbourlist(at, RCUT[2]);
      neighbourlist(at, RCUT[3])
end
GC.gc()
print("IP    : ")
energy(IP, at); @time energy(IP, at); GC.gc()
print("IPf   : ")
energy(IPf, at); @time energy(IPf, at);GC.gc()
print("IPff  : ")
energy(IPff, at); @time energy(IPff, at);GC.gc()
print("IPff2 : ")
energy(IPff2, at); @time energy(IPff, at);GC.gc()
print("IPcl  : ")
energy(IPcl, at); @time energy(IPcl, at);GC.gc()
print("IPclf : ")
energy(IPclf, at); @time energy(IPclf, at);GC.gc()
print("IPclff: ")
energy(IPclff, at); @time energy(IPclff, at);GC.gc()

@info("For comparison some simpler potentials:")
IP0 = NBodyIP( [V2[1], V3[1], V4[1]] )
IP1 = NBodyIP( [V2[1].Vr, V3[1].Vr, V4[1].Vr] )
print("IP0: EnvIP deg = 0 : ")
energy(IP0, at); @time energy(IP0, at);
print("IP1:  4B IP        : ")
energy(IP1, at); @time energy(IP1, at);

@info("Test Error in Forces")
@show maximum(norm.(forces(IP, at) - forces(IPf, at)))
@show maximum(norm.(forces(IPff, at) - forces(IPf, at)))
@show maximum(norm.(forces(IPff2, at) - forces(IPf, at)))

@info("Test evaluation time for forces...")
at = rattle!(bulk(:W, cubic=true) * 8, 0.02)
print("IP  : ")
forces(IP, at); @time forces(IP, at);
GC.gc()
print("IPf  : ")
forces(IPf, at); @time forces(IPf, at);
GC.gc()
print("IPff  : ")
forces(IPff, at); @time forces(IPff, at);
print("IPff2 : ")
forces(IPff2, at); @time forces(IPff, at);
GC.gc()
print("IP0  : ")
forces(IP0, at); @time forces(IP0, at);
GC.gc()
print("IP1  : ")
forces(IP1, at); @time forces(IP1, at);

# BACKUP PROFILING CODE
# at = rattle!(bulk(:W, cubic=true) * 7, 0.02)
# Profile.clear(); @profile forces(IPff, at);
# Profile.print()
#
# at = rattle!(bulk(:W, cubic=true) * 3, 0.02)
# V4 = IPff.components[end]
# @code_warntype forces(V4, at);

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
using Pkg

# Collection of sources required to build LCIOWrapBuilder
sources = [
    DirectorySource("LCIO_Julia_Wrapper")
]

julia_version = v"1.5.3"

name = "LCIO_Julia_Wrapper"
version = get(ENV, "TRAVIS_TAG", "")
if version == ""
	version = v"0.13.1-alpha1"
else
	version = VersionNumber(version)
end

# Bash recipe for building across all platforms
script = raw"""
	mkdir build && cd build
	cmake -DJulia_PREFIX=${prefix} -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DCMAKE_FIND_ROOT_PATH=${prefix} ..
	VERBOSE=ON cmake --build . --config Release --target install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = Platform[
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(cxxstring_abi=:cxx11)),
]

# The products that we will ensure are always built
products = [
    LibraryProduct("liblciowrap", :lciowrap)
]

# Dependencies that must be installed before this package can be built
dependencies = [
	Dependency(PackageSpec(name="libcxxwrap_julia_jll",version=v"0.8.5")),
	Dependency(PackageSpec(name="LCIO_jll", version=v"2.15.4")),
	BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version))
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"7", julia_compat = "^$(julia_version.major).$(julia_version.minor)")

ext: jl


text:
using Tissue
import Tissue as T

const NUM_ITERATIONS = 5
const TARGET_PERIOD = 0.5

mutable struct TimestampSourceCalculator
    last::Int
    max_gen::Int
    TimestampSourceCalculator() = new(1, NUM_ITERATIONS)
end

function T.process(c::TimestampSourceCalculator)
    if c.last > c.max_gen
        return nothing
    end

    c.last += 1

    return time()
end

mutable struct PeriodCalculator
    last_ts::Float64

    PeriodCalculator() = new(-1)
end

function T.process(calc::PeriodCalculator, new_ts::Float64)
    if calc.last_ts == -1
        calc.last_ts = new_ts
        return new_ts
    end

    global latest_period
    latest_period = new_ts - calc.last_ts

    calc.last_ts = new_ts

    new_ts
end

struct SleeperCalculator
    "Sleep time in seconds"
    sleep_time::Float64
end

function T.process(calc::SleeperCalculator, in_num)
    sleep(calc.sleep_time)

    in_num
end

@graph SleeperGraph begin
    @calculator source = TimestampSourceCalculator()
    @calculator period = PeriodCalculator()
    @calculator sleeper = SleeperCalculator(TARGET_PERIOD)

    @bindstreams period (new_ts = source)
    @bindstreams sleeper (in_num = source)
end

graph = SleeperGraph()

T.start(graph)
T.wait_until_done(graph)

@test isapprox(latest_period, TARGET_PERIOD; rtol=0.1)


text:
using Test, Distributed

@testset "Tests" begin
    tests = [
        # Cleaned
        "PlotsX/PlotsX.jl",
        "Datasets/Datasets.jl",
        "Scripting/Scripting.jl",
        # To clean
        "Neural/Neural.jl",
        "MonteCarlo/MonteCarlo.jl",
        "utility.jl",
        # "special.jl",
        "transformations.jl",
        # Deprecated
        # "distributions/distributions.jl",
        # "neural/neural.jl",
    ]

    pmap(tests) do t
        @eval module $(Symbol("Test", t))
            include($t)
        end
        return
    end
end


text:
using TopOpt, LinearAlgebra, StatsFuns
using TopOpt.TrussTopOptProblems
using TopOptMakie
using VTKDataIO
using PyCall
using VTKDataTypes
using DashVtk
using Dash
using DashHtmlComponents

#dash = pyimport("dash")
#dash_html = pyimport("dash_html_components")
dash_vtk = pyimport("dash_vtk")


text:
@deprecate in_edges in_neighbors
@deprecate out_edges out_neighbors

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
"""
    image(cmd0::String="", arg1=[]; kwargs...)

Place images or EPS files on maps.

Full option list at [`psimage`](http://gmt.soest.hawaii.edu/doc/latest/psimage.html)

Parameters
----------

- $(GMT.opt_B)
- **D** : **ref_point** : -- Str --

    Sets reference point on the map for the image using one of four coordinate systems.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/psimage.html#d)
- **F** : **box** : -- Str or [] --

    Without further options, draws a rectangular border around the image using MAP_FRAME_PEN.
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/psimage.html#f)
- **I** : **invert_1bit** : -- Number or Str --

    Invert 1-bit image before plotting.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/psimage.html#i)
- $(GMT.opt_J)
- $(GMT.opt_Jz)
- **M** : **monochrome** : -- Bool or [] --

    Convert color image to monochrome grayshades using the (television) YIQ-transformation.
    [`-M`](http://gmt.soest.hawaii.edu/doc/latest/psimage.html#m)
- $(GMT.opt_R)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function image(cmd0::String="", arg1=[]; K=false, O=false, first=true, kwargs...)

	length(kwargs) == 0 && return monolitic("psimage", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)
	output, opt_T, fname_ext = fname_out(d)		# OUTPUT may have been an extension only

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12c/12c")
	cmd, = parse_JZ(cmd, d)
	cmd  = parse_UVXY(cmd, d)
	cmd, = parse_p(cmd, d)
	cmd, = parse_t(cmd, d)
	cmd  = parse_params(cmd, d)

	cmd = add_opt(cmd, 'D', d, [:D :ref_point])
	cmd = add_opt(cmd, 'F', d, [:F :box])
	cmd = add_opt(cmd, 'G', d, [:G :bit_color])
	cmd = add_opt(cmd, 'I', d, [:I :invert_1bit])
	cmd = add_opt(cmd, 'M', d, [:M :monochrome])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)		# Find how data was transmitted

	cmd = finish_PS(d, cmd, output, K, O)
	return finish_PS_module(d, cmd, "", output, fname_ext, opt_T, K, "psimage", arg1)
end

# ---------------------------------------------------------------------------------------------------
image!(cmd0::String="", arg1=[]; kw...) = image(cmd0, arg1; K=true, O=true,  first=false, kw...)

# ---------------------------------------------------------------------------------------------------
psimage  = image			# Alias
psimage! = image!			# Alias


text:
function evaluate(modl::KrigingModel, PointToExamine::Array{T,1}; k = correlation_vector_of_a_point(PointToExamine,  modl.X,  modl.cov_func, modl.hyper)) where T<:Real
    yhat =  modl.mu + (k') *  modl.invK * modl.y_m_mu
    return yhat
end

function predicted_error(modl::KrigingModel, PointToExamine::Array{T,1}; k = correlation_vector_of_a_point(PointToExamine,  modl.X,  modl.cov_func, modl.hyper)) where T<:Real
    # Equation 9 of Jones Schonlau, Welch
    dim = size(modl.invK)[1]
    Kinvk = modl.invK * k
    bracketed_term = max(eps(), 1 + ((1 - sum(Kinvk))^2) / sum(modl.invK * ones(dim)) - ((k') * Kinvk))
    return sqrt(modl.sigma2 * bracketed_term)
end

function (modl::KrigingModel)(PointToExamine::Array{T,1}) where T<:Real
    k = correlation_vector_of_a_point(PointToExamine,  modl.X,  modl.cov_func, modl.hyper)
    value = evaluate(modl, PointToExamine; k = k)
    error = predicted_error(modl, PointToExamine; k = k)
    return (value = value, standard_error = error)
end


text:
module Relax

using ..SCF:scf_energy
using ..Force_Stress:get_energy_force_stress_fft
using ..Force_Stress:safe_mode_energy
using ..CrystalMod:get_grid
using ..CrystalMod:crystal
using ..CrystalMod:print_with_force_stress
using ..CrystalMod:makecrys
using ..CrystalMod:write_axsf
using ..CalcTB:calc_tb_fast
using ..Atomdata:atom_radius

using LinearAlgebra
using ..Force_Stress:inv_reshape_vec
using ..Force_Stress:reshape_vec

using Optim
using LineSearches
using ..ManageDatabase:prepare_database
using ..ManageDatabase:database_cached
using ..MyOptim:conjgrad

export relax_structure

"""
    function relax_structure(crys::crystal, database; smearing = 0.01, grid = missing, mode="vc-relax", nsteps=50, update_grid=true, conv_thr=2e-4)

Relax structure. Primary user function is relax_structure in TightlyBound.jl, which calls this one.
"""
function relax_structure(crys::crystal, database; smearing = 0.01, grid = missing, mode="vc-relax", nsteps=50, update_grid=true, conv_thr = 1e-2, filename="t.axsf")

    if update_grid==false
        grid = get_grid(crys)
    end

    eden = missing #starts off missing

    #do this ahead of first iteration, to get memory in correct place
    tbc = calc_tb_fast(deepcopy(crys), database)
    energy_tot, efermi, e_den, dq, VECTS, VALS, error_flag, tbcx  = scf_energy(tbc, smearing=smearing, grid=grid, e_den0=eden, conv_thr=1e-7)

    if error_flag
        println("warning error computing scf in relax_structure, zeroth iteration")
    end

    eden = deepcopy(tbc.eden)

    A0 = deepcopy(crys.A)
    strain = zeros(3,3)

    x0 = inv_reshape_vec(crys.coords, strain, crys.nat, strain_mode=true)
    crys_working = deepcopy(crys)

    fcall = 0
    firstiter = true

    nat = crys.nat

    function fix_strain(s)
        for i = 1:3
            for j = 1:3
                s[i,j] = min(s[i,j], 0.60)
                s[i,j] = max(s[i,j], -0.60)
            end
        end
        return 0.5*(s'+s)
    end

    energy_global = -99.0

    CRYSTAL = []
    FORCES = []

    function fn(x)
#        println("CALL FN", x)
        coords, strain = reshape_vec(x, nat, strain_mode=true)

        strain=fix_strain(strain)
        coords = coords .% 1.0

        A = A0 * (I(3) + strain)

        crys_working.coords = coords
        if  mode == "vc-relax"
            crys_working.A = A
        end

        if crys_working == tbc.crys
            #            println("SKIP")
            #we already have energy from calling grad, we don't need to call again.
            return energy_global, true
        end

        #        println("FN crys")
#        println(crys_working)
#        println("before short ")

        tooshort, energy_short = safe_mode_energy(crys_working, database)

        if tooshort
            return energy_short, false
        end

        if crys_working != tbc.crys
#                println("yes calc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
#                println(crys_working)
            tbc = calc_tb_fast(deepcopy(crys_working), database, verbose=false)
        else
#                println("nocalc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
#                println(crys_working)
#                println(tbc.crys)

        end
#            println(crys_working)
        energy_tot, efermi, e_den, dq, VECTS, VALS, error_flag, tbcx  = scf_energy(tbc, smearing=smearing, grid=grid, e_den0=eden, verbose=false, conv_thr=1e-6)
        eden = deepcopy(tbcx.eden)

#        println("fn $energy_tot fnffnfnfnffnfnfnfnfnfnfnfnfnfnfffffff")

#        println("ENERGY $energy_tot $energy_global")

        energy_global=energy_tot
#        println("FN end")

        return energy_tot, !error_flag

    end

    f_cart_global = []
    stress_global = []

#    function grad(x)
    function grad(storage, x)

#        println("CALL GRAD", x)

#        println("typeof x ", typeof(x))
        fcall += 1

        coords, strain = reshape_vec(x, nat, strain_mode=true)

        strain=fix_strain(strain)
        coords = coords .% 1.0

        A = A0 * (I(3) + strain)

        crys_working.coords = coords

        if  mode == "vc-relax"
            crys_working.A = A
        end

#        println("GRAD crys")
#        println(crys_working)

#        println("crys_working $fcall")
#        println(crys_working)

        tooshort, energy_short = safe_mode_energy(crys_working, database)

#        println("too short ", tooshort)

        tbc = calc_tb_fast(deepcopy(crys_working), database, verbose=false)

        #        if crys_working != tbc.crys && !tooshort
        if !tooshort
#            println("yes calc forces -----------------------------------------------------------------------------------------")
#            println(crys_working)

            energy_tot, efermi, e_den, dq, VECTS, VALS, error_flag, tbcx  = scf_energy(tbc, smearing=smearing, grid=grid, e_den0=eden, verbose=false, conv_thr=1e-6)

            eden = deepcopy(tbcx.eden)

            #        else
#            println("no calc forces ------------------------------------------------------------------------------------------")
#            println(crys_working)
#            println(tbc.crys)
#            println("__")
        elseif tooshort
            energy_tot=energy_short
        end

        energy_global=energy_tot

        energy_tmp,  f_cart, stress =  get_energy_force_stress_fft(tbc, database; do_scf = false, smearing = smearing, grid = grid, vv=[VECTS, VALS, efermi] )

        #energy_tmp,  f_cart, stress =  get_energy_force_stress_fft(tbc, database; do_scf = true, smearing = smearing, grid = grid )

        push!(CRYSTAL, deepcopy(tbc.crys))
        push!(FORCES, f_cart)

        f_cart_global = f_cart
        stress_global = stress

        fsum = sum((f_cart).^2)^0.5
        ssum= sum((stress).^2)^0.5

#        println("f_cart")
#        println(f_cart)

        println()
        println("FCALL $fcall en:  $energy_tot (Ryd)  fsum:  $fsum  ssum:  $ssum    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
        println()

        print_with_force_stress(crys_working, f_cart_global, stress_global)
        println()
#        println("grad stress: ", stress)

#        stress_units = cell_force( crys_working.A, stress)
        #        stress_units = (stress * abs(det(crys_working.A)) *  inv(crys_working.A))'
        #stress_units = (stress * abs(det(crys_working.A)) )'

#        println("grad stress units: ", -stress_units)

        #        f_crys = f_cart * inv(crys_working.A)

        f_crys = f_cart * crys_working.A'

        g = inv_reshape_vec(-f_crys, -stress* abs(det(crys_working.A))  , crys_working.nat, strain_mode=true)

#        println("g ", g)

        if mode != "vc-relax" #don't need stress
            g[3*nat+1:end] = zeros(6)
        end

        storage[:] = g

        neaten(storage)

        return energy_global, storage
#        println("ret storage $fcall")
#        println(storage)
    end

    println("starting vec ")
    println(x0)
    println()

#    opts = Optim.Options(g_tol = 7e-4,f_tol = 7e-4, x_tol = 7e-4,
#                         iterations = nsteps,
#                         store_trace = true,
#                         show_trace = false)

    opts = Optim.Options(g_tol = conv_thr,f_tol = conv_thr, x_tol = conv_thr,
                             iterations = nsteps,
                             store_trace = true,
                             show_trace = false)

    #res = optimize(fn,grad, x0, ConjugateGradient(), opts)

    #=
    #preconditioner, guess for inv Hess is based on the  metric. see qe bfgs_module.f90
    function init(x)
        num = 3*nat + 6

#        factor=0.1
        factor = 1.0

        P = zeros(eltype(x), num, num)

        A = crys.A
        g = A'*A
        ginv = inv(g)

        for a = 1:nat
            aa = (a-1)*3
            for i = 1:3
                for j = 1:3
                    P[aa+i ,aa+j] = g[i,j] * factor
                end
            end
        end

        vol = abs(det(A))

#        if mode == "vc-relax"
        for b = 1:2
            bb = 3 * nat + (b - 1)*3
            for i = 1:3
                for j = 1:3
                    P[bb+i,bb+j] = 0.04 * vol * ginv[i,j] * factor
                end
            end
        end
#    end

        return inv(P)
    end
=#

    res = missing
    if false

        #        res = optimize(fn,grad, x0, BFGS(  initial_invH = init, linesearch=LineSearches.MoreThuente() ), opts)
        #res = optimize(fn,grad, x0, BFGS(initial_invH = init ), opts)
        #        res = optimize(fn,grad, x0, BFGS(initial_invH = init, linesearch=LineSearches.MoreThuente() ), opts)

        #    try
        #res = optimize(fn,grad, x0, BFGS( linesearch=LineSearches.MoreThuente() ), opts)

        #    res = optimize(fn,grad, x0, ConjugateGradient( linesearch=LineSearches.MoreThuente() ) , opts)

        res = optimize(fn,grad, x0, ConjugateGradient( linesearch=LineSearches.BackTracking( maxstep=0.05  ) ) , opts)

        #    res = optimize(fn,grad, x0, ConjugateGradient( linesearch=LineSearches.BackTracking( maxstep=0.05) ) , opts)
        #    res = optimize(fn,grad, x0, ConjugateGradient( linesearch=LineSearches.MoreThuente(alphamax=0.05) ) , opts)

        minvec = Optim.minimizer(res)

        coords, strain = reshape_vec(minvec, nat, strain_mode=true)
        coords= coords .% 1.0

        A = A0 *( I(3) + strain)

        crys = makecrys(A, coords, crys.types, units="Bohr")

        tbc = calc_tb_fast(deepcopy(crys), database)
        energy_tot, efermi, e_den, dq, VECTS, VALS, error_flag, tbcx  = scf_energy(tbc, smearing=smearing, grid=grid, e_den0=eden, conv_thr=1e-7)

        eden = deepcopy(tbc.eden)

        A0 = deepcopy(crys.A)
        strain = zeros(3,3)

        x0 = inv_reshape_vec(crys.coords, strain, crys.nat, strain_mode=true)
        crys_working = deepcopy(crys)

        fcall = 1000
        firstiter = true

        println("RELAX Round 2")
        #    res = optimize(fn,grad, x0, ConjugateGradient( linesearch=LineSearches.MoreThuente(alphamax=0.05) ) , opts)
        res = optimize(fn,grad, x0, ConjugateGradient(  ) , opts)

        nat = crys.nat

        #    res = optimize(fn,grad, x0, ConjugateGradient(  ) , opts)

        #res = optimize(fn,grad, x0, LBFGS(m=8 ) , opts)

        #    catch e
        #        if e isa InterruptException
        #            println("user interrupt")
        #            return tbc.crys
        #        else
        #            println("unknown error")
        #            return tbc.crys
        #        end
        #    end

        #res = optimize(fn,grad, x0, BFGS(  initial_invH = init, linesearch=LineSearches.HagerZhang() ), opts)

        energy = res.minimum

        # bad : res = optimize(fn,grad, x0, ConjugateGradient(linesearch=LineSearches.BackTracking(order=3) ), opts)
        #linesearch=LineSearches.BackTracking(order=2),
        #LineSearches.MoreThuente()
        #LineSearches.HagerZhang()
        #linesearch=LineSearches.BackTracking(order=2)
        #linesearch=LineSearches.StrongWolfe()

        println()
        println("res")
        println(res)

        minvec = Optim.minimizer(res)

    else

        minvec, energy, grad = conjgrad(fn, grad, x0; maxstep=5.0, niters=nsteps, conv_thr = conv_thr)

    end

#    println("minvec")
#    println(minvec)

    coords, strain = reshape_vec(minvec, nat, strain_mode=true)
    coords = coords .% 1.0

    strain = fix_strain(strain)

    A = A0 *( I(3) + strain)

    cfinal = makecrys(A, coords, crys.types, units="Bohr")

    if !ismissing(filename)
        write_axsf(CRYSTAL, filename=filename, FORCES=FORCES)
    end

    return cfinal, tbc, energy, f_cart_global, stress_global

#    return res

end

"""
    function neaten(storage)

Detect/Enforce symmetries in forces/stress, tight tolerance. Deal with minor numerical issues causing symmetry breaking
"""
function neaten(storage)
#    println("n before ", storage)
    n = length(storage)
    for i in 1:n-6
        if abs(storage[i]) < 1e-8
            storage[i] = 0.0
        end
    end

    for i in n-2:n
        if abs(storage[i]) < 1e-7
            storage[i] = 0.0
        end
    end

    for i in 1:n
        for j in i+1:n
            if abs(storage[j] - storage[i]) < 1e-8
                storage[j] = storage[i]
            elseif abs(storage[j] + storage[i]) < 1e-8
                storage[j] = -storage[i]
            elseif abs(storage[j] + 2.0*storage[i]) < 1e-8
                storage[j] = -2.0*storage[i]
            elseif abs(storage[j] + 0.5*storage[i]) < 1e-8
                storage[j] = -0.5*storage[i]
            elseif abs(storage[j] - 2.0*storage[i]) < 1e-8
                storage[j] = 2.0*storage[i]
            elseif abs(storage[j] - 0.5*storage[i]) < 1e-8
                storage[j] = 0.5*storage[i]
            end
        end
    end
#    println("n after ", storage)
end

#=
function finite_diff_grad(fn, x0, step)
    x = deepcopy(x0)
    g = zeros(size(x))
    for i in 1:length(x)
        x[:] = x0[:]
        x[i] += step
        a = fn(x)

        x[:] = x0[:]
        x[i] -= step
        b = fn(x)

        g[i] = (a-b)/(2*step)
    end

    return g

end
=#

"
    function make_random_crystal(types)

Make a random crystal with types t
"
function make_random_crystal(types)

    nat = length(types)

    prepare_database(types)

    total_vol = 0.0
    for t in types
        total_vol += 4 * pi / 3 * (atom_radius[t]/100 / 0.529177)^3 * 1.0
    end
    println("total_vol $total_vol")

    for i = 1:30

        A = zeros(3,3)
        c = rand(nat, 3)
        c[1,:] = [0.0 0.0 0.0]

        A[1,1] = 1.0

        l2 = 1.0 + (rand(1)[1] - 0.5)
        th2 =  pi*(1.0 + 0.5*(rand(1)[1] - 0.5))

        l3 = 1.0 + (rand(1)[1] - 0.5)
        th3 =  pi*(1.0 + 0.5*(rand(1)[1] - 0.5))

        phi = pi*(0.25*(rand(1)[1] - 0.5))

        A[2,:] = [l2 * sin(th2)  l2 * cos(th2) 0.0]
        A[3,:] = [sin(phi)  l3 * sin(th3) *cos(phi)  l3 * cos(th3)* cos(phi)]

#        println("A before ", det(A), "  ", total_vol^(1/3) * 1.5 )
        A = A / (abs(det(A)))^(1/3) * total_vol^(1/3) * 1.2
#        println("A after ", det(A))

        A = A * det(A) / abs(det(A))

        crys = makecrys(A, c, types, units="Bohr")

        within_fit = calc_tb_fast(crys, database_cached, check_only=true)
        if within_fit
            return crys * 1.05
        end

        total_vol = total_vol * 1.05

    end
    println("failed to generate crystal")
    return missing

end

end #end module


text:
using .util_convexopt

function _recon2d_tv_primaldual!(u, A, b0, niter, w_tv, sigmas, tau)
    At = sparse(A')
    H, W = size(u)

    b = vec(b0)

    ubar = deepcopy(u)
    u_prev = similar(u)

    p1 = zeros(size(b))
    p2 = zeros(H, W, 2)

    p_adjoint = zeros(H, W)

    du = similar(p2)
    divp2 = similar(p2)

    data1 = similar(b)
    u_temp = similar(u)

    for it=1:niter
        u_prev .= u

        # dual update
        data1 .= mul!(data1, A, vec(ubar)) .- b
        # p1 = proj_dual_l1(p1_ascent, w_data) # l1 norm
        p1 .= (p1 .+ sigmas[1] * data1) ./ (sigmas[1] + 1.0) # l2 norm
        mul!(view(p_adjoint, :), At, p1)
        # p_adjoint .= reshape(At * p1, H, W)

        grad!(du, ubar)
        p2 .+= sigmas[2] .* du
        util_convexopt.proj_dual_iso!(p2, u_temp, w_tv)
        # util_convexopt.proj_dual_l1!(p2, w_tv) #anisotropic TV

        p_adjoint .-= div!(divp2, p2) # p1_adjoint + p2_adjoint

        # primal update
        u .= max.(u .- tau .* p_adjoint, 0.0) # positivity constraint

        # acceleration
        ubar .= 2 .* u .- u_prev

        # compute primal energy (optional)
        if it % 20 == 0
            energy = sum(data1.^2) / length(data1) + sum(abs.(du)) / length(du)
            println("$it, approx. primal energy: $energy")
        end
    end
    return u
end

"""
    recon2d_tv_primaldual!(u::Array{T, 2}, A, b, niter::Int, w_tv::T, c=10.0)

Reconstruct a 2d image by TV-L2 model using Primal Dual optimization method

# Args
u : Initial guess of image
A : Forward opeartor
b : Projection data (2 dimension or 1 dimension)
niter: number of iterations
w_tv: weight for TV term
c : See 61 page in 2016_Chambolle,Pock_An_introduction_to_continuous_optimization_for_imagingActa_Numerica
"""
function recon2d_tv_primaldual!(u::Array{T, 2}, A, b, niter::Int, w_tv::T, c=10.0) where {T <: AbstractFloat}
    @time op_A_norm = util_convexopt.compute_opnorm(A)
    println("@ opnorm of forward projection operator: $op_A_norm")
    ops_norm = [op_A_norm, sqrt(8)]

    sigmas = zeros(length(ops_norm))
    for i=1:length(ops_norm)
        sigmas[i] = 1.0 / (ops_norm[i] * c)
    end

    tau = c / sum(ops_norm)
    println("@ step sizes sigmas: ", sigmas, ", tau: $tau")

    return _recon2d_tv_primaldual!(u, A, b, niter, w_tv, sigmas, tau)
end

function _recon2d_stack_tv_primaldual!(u::Array{T, 3}, A, b0::Array{T, 3}, niter::Int, w_tv::T, sigmas, tau::T) where {T<:AbstractFloat}
    H, W, nslice = size(u)

    At = sparse(A') # this significatnly improves the performance

    ubar = deepcopy(u)
    u_prev = similar(u)

    b_axWxH = permutedims(b0, [3, 2, 1]) # PermutedDimsArray(A, (3,1,2));
    b = reshape(b_axWxH, :, nslice)

    halfslice = Int(floor(nslice/2))

    p1 = zeros(size(b))
    p2 = zeros(H, W, nslice, 3)

    p_adjoint = zeros(H, W, nslice)
    du = similar(p2)
    divu = similar(p2)
    temp_residual = zeros(size(A, 1), nslice)

    u_temp = similar(u)

    for it=1:niter
        u_prev .= u

        # dual update
        Threads.@threads  for slice=1:nslice
            ubar_slice = view(ubar, :, :, slice)
            p_adjoint_slice = vec(view(p_adjoint, :, :, slice))

            # for l2 norm
            @views temp_residual[:, slice] .= mul!(temp_residual[:, slice], A, vec(ubar_slice)) .- b[:, slice]
            @views p1[:, slice] .= (p1[:, slice] .+ sigmas[1] .* temp_residual[:, slice]) ./ (sigmas[1] + 1.0)

            @views mul!(p_adjoint_slice, At, p1[:, slice])
        end

        # p2: 3d gradient and divergence
        grad!(du, ubar)
        p2 .+= sigmas[2] .* du
        util_convexopt.proj_dual_iso!(p2, u_temp, w_tv)
        # util_convexopt.proj_dual_l1!(p2, w_tv) # anisotropic TV

        # if p_adjoint .++ -view(), memory allocation increases. WHY??
        p_adjoint .-= div!(divu, p2)

        # primal update
        u .= max.( u .- tau .* p_adjoint, 0.0) # positivity constraint

        # acceleration
        ubar .= 2.0 .* u .- u_prev

        # compute primal energy (optional)
        if it % 10 == 0
            # energy = sum(data1.^2) / length(data1) + sum(abs.(data2)) / length(data2)
            println("iter: $it, residua: $(sum(temp_residual .^ 2 ) / length(temp_residual))")
        end
    end
    residual = sum(temp_residual .^ 2 ) / length(temp_residual)
    return residual
end

"""
    recon2d_tv_primaldual!(u::Array{T, 2}, A, b::Array{T, 2}, niter::Int, w_tv::T, c=1.0)

Reconstruct a 2d image by TV-L2 model using Primal Dual optimization method

# Args
- u : Initial guess of 3d image
- A : Forward opeartor
- b : Projection data [nangles x detheight x detwidth]
- niter: number of iterations
- w_tv: weight for TV term
- c : See 61 page in 2016_Chambolle,Pock_An_introduction_to_continuous_optimization_for_imagingActa_Numerica
"""
function recon2d_stack_tv_primaldual!(u::Array{T, 3}, A, b::Array{T, 3}, niter::Int, w_tv::T, c=10.0) where {T <: AbstractFloat}
    @time op_A_norm = util_convexopt.compute_opnorm(A)
    println("@ opnorm of forward projection operator: $op_A_norm")
    ops_norm = [op_A_norm, sqrt(8)]
    println("TODO! in 3d, sqrt(8) would be wrong")

    sigmas = zeros(length(ops_norm))
    for i=1:length(ops_norm)
        sigmas[i] = 1.0 / (ops_norm[i] * c)
    end

    tau = c / sum(ops_norm)
    println("@ step sizes sigmas: ", sigmas, ", tau: $tau")

    residual = _recon2d_stack_tv_primaldual!(u, A, b, niter, w_tv, sigmas, tau)
end

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
blank = ""
nonblank = "!"

println("The length of blank is ", length(blank))
println("That blank is empty is ", isempty(blank))
println("That blank is not empty is ", !isempty(blank))

println()
println("The length of nonblank is ", length(nonblank))
println("That nonblank is empty is ", isempty(nonblank))
println("That nonblank is not empty is ", !isempty(nonblank))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module MLDatasets

include("io/CoNLL.jl")

include("CIFAR10.jl")
include("CIFAR100.jl")
include("MNIST/MNIST.jl")
include("PTBLM.jl")
include("UD_English.jl")

end


text:
using Distributions
using VariantsNMF
using Test

#### simulate
V = get_one_simulated_V()
V = round.(V .* 100)

@testset "nmf_MU" begin
    #### params for rank selection
    rs_params = RSParams(
        K_min     = 1,
        K_max     = 10,
        n_iter    = 50,
        pert_meth = :multinomial,
        seed      = 123
    )

    #### nmf algorithm
    nmf_global_params = NMFParams(
        init          = :random,
        dist          = Uniform(0, 1),
        max_iter      = 1_000,
        stopping_crit = :conn,
        stopping_iter = 10,
        verbose       = false,
    )

    nmf_local_params = MUParams(
        β            = 1,
        div          = :β,
        scale_W_iter = false,
        scale_W_last = false,
        alg          = :mu
    )

    nmf = NMF(
        solver        = nmf_MU,
        global_params = nmf_global_params,
        local_params  = nmf_local_params
    )

    #### run rank selection by stability procedure
    @time rs_results = rank_by_stability(V, rs_params, nmf)

    @test "rank" in names(rs_results.df_metrics)
    @test "stab_avg" in names(rs_results.df_metrics)
    @test "stab_std" in names(rs_results.df_metrics)
    @test "fid_avg" in names(rs_results.df_metrics)
    @test "fid_std" in names(rs_results.df_metrics)
end


text:
import Term: Measure
import Term: remove_markup

@testset "\e[34mMeasure - str" begin
    @test Measure("a"^10).w == 10

    strings = (
        "asadasda"^2,
        "~~±±||__ASdaSDvxcvxc\nasdfsiaudfhsdiufhndskjv",
        "dasdadsa\nsdfsd"^7,
        "asdasdsa|dfvcxvashfusdn\nfidsuhfsdf",
    )
    for string in strings
        m = Measure(string)
        @test m.w == lw(string)
        @test m.h == nlines(string)
    end

    strings = (
        "[red]is my color"^2,
        "~~±±||__[red on_green]aasdas[bold]asdsad[/bold]asdas[/red on_green]xc\nasdfsiaudfhsdiufhndskjv",
        "dasda[#ffffff]dsa[#ffffff]\nsdfsd"^7,
        "asdasdsa|dfvcxvashf[bold]usdn\nfid[/bold]suhfsdf",
    )
    for string in strings
        m = Measure(string)
        @test m.w == lw(remove_markup(string))
        @test m.h == nlines(remove_markup(string))
    end
end


text:
module OpenQASM_Map

IBM_qelib1.inc="""
// Quantum Experience (QE) Standard Header
// file: qelib1.inc

// --- QE Hardware primitives ---

// 3-parameter 2-pulse single qubit gate
gate u3(theta,phi,lambda) q { U(theta,phi,lambda) q; }
// 2-parameter 1-pulse single qubit gate
gate u2(phi,lambda) q { U(pi/2,phi,lambda) q; }
// 1-parameter 0-pulse single qubit gate
gate u1(lambda) q { U(0,0,lambda) q; }
// controlled-NOT
gate cx c,t { CX c,t; }
// idle gate (identity)
gate id a { U(0,0,0) a; }

// --- QE Standard Gates ---

// Pauli gate: bit-flip
// Pauli gate: bit and phase flip
gate y a { u3(pi,pi/2,pi/2) a; }
// Pauli gate: phase flip
gate z a { u1(pi) a; }
// Clifford gate: Hadamard
gate h a { u2(0,pi) a; }
// Clifford gate: sqrt(Z) phase gate
gate s a { u1(pi/2) a; }
// Clifford gate: conjugate of sqrt(Z)
gate sdg a { u1(-pi/2) a; }
// C3 gate: sqrt(S) phase gate
gate t a { u1(pi/4) a; }
// C3 gate: conjugate of sqrt(S)
gate tdg a { u1(-pi/4) a; }

// --- Standard rotations ---
// Rotation around X-axis
gate rx(theta) a { u3(theta,-pi/2,pi/2) a; }
// rotation around Y-axis
gate ry(theta) a { u3(theta,0,0) a; }
// rotation around Z axis
gate rz(phi) a { u1(phi) a; }

// --- QE Standard User-Defined Gates  ---

// controlled-Phase
gate cz a,b { h b; cx a,b; h b; }
// controlled-Y
gate cy a,b { sdg b; cx a,b; s b; }
// controlled-H
gate ch a,b {
h b; sdg b;
cx a,b;
h b; t b;
cx a,b;
t b; h b; s b; x b; s a;
}
// C3 gate: Toffoli
gate ccx a,b,c
{
  h c;
  cx b,c; tdg c;
  cx a,c; t c;
  cx b,c; tdg c;
  cx a,c; t b; t c; h c;
  cx a,b; t a; tdg b;
  cx a,b;
}
// controlled rz rotation
gate crz(lambda) a,b
{
  u1(lambda/2) b;
  cx a,b;
  u1(-lambda/2) b;
  cx a,b;
}
// controlled phase rotation
gate cu1(lambda) a,b
{
  u1(lambda/2) a;
  cx a,b;
  u1(-lambda/2) b;
  cx a,b;
  u1(lambda/2) b;
}
// controlled-U
gate cu3(theta,phi,lambda) c, t
{
  // implements controlled-U(theta,phi,lambda) with  target t and control c
  u1((lambda-phi)/2) t;
  cx c,t;
  u3(-theta/2,0,-(phi+lambda)/2) t;
  cx c,t;
  u3(theta/2,phi,0) t;
}

"""

function x(;args...)
    return "x $(args[:target]);"
end
function y(;args...)
    return "y $(args[:target]);"
end
function z(;args...)
    return "z $(args[:target]);"
end
function h(;args...)
    return "h $(args[:target]);"
end
function c_x(;args...)
    return "cx $(args[:ctrl]), $(args[:target]);"
end
function c_y(;args...)
    return "cy $(args[:ctrl]), $(args[:target]);"
end
function c_z(;args...)
    return "cz $(args[:ctrl]), $(args[:target]);"
end
function r_x(;args...)
    return "rx($(args[:param])) $(args[:target]);"
end
function r_y(;args...)
    return "ry($(args[:param])) $(args[:target]);"
end
function r_z(;args...)
    return "rz($(args[:param])) $(args[:target]);"
end
function r_phase(;args...)
    return "u1($(args[:param])) $(args[:target]);"
end

function c_r_z(;args...)
    return "crz($(args[:param])) $(args[:ctrl]),$(args[:target]);"
end
function c_r_phase(;args...)
    return "cu1($(args[:param])) $(args[:ctrl]),$(args[:target]);"
end

Dict(   :x=>x,
        :y=>y,
        :z=>z,
        :h=>h,
        :c_x=>c_x,
        :c_y=>c_y,
        :c_z=>c_z,
        :r_x=>r_x,
        :r_y=>r_y,
        :r_z=>r_z,
        :c_r_z=>c_r_z,
        :c_r_phase=>c_r_phase
)

end


text:
function rand_weighted_index(cumulative_weights::Array{Float64, 1})
	val::Float64 = rand()
	idx::Int64 = 1
	while cumulative_weights[idx] < val
		idx += 1
	end
	return idx
end


text:
function plotCabidaOptima(setCeldasActivas, pisoActivoCeldas, pos_x, pos_y, largo_x, largo_y, R, r, dca, fig_, ax_, ax_mat_)

    global fig, ax, ax_mat

    fig=fig_; ax=ax_; ax_mat=ax_mat_;
    for j in setCeldasActivas

        p1_vec_ = (R[r].mat * ([pos_x[j]'; pos_y[j]'] - R[r].cr) + R[r].cr)';
        p2_vec_ = (R[r].mat * ([pos_x[j]' + largo_x[j,r]'; pos_y[j]'] - R[r].cr) + R[r].cr)';
        p3_vec_ = (R[r].mat * ([pos_x[j]' + largo_x[j,r]'; pos_y[j]' + largo_y[j,r]'] - R[r].cr) + R[r].cr)';
        p4_vec_ = (R[r].mat * ([pos_x[j]'; pos_y[j]' + largo_y[j,r]'] - R[r].cr) + R[r].cr)';

        V_tapa = [p1_vec_; p2_vec_; p3_vec_; p4_vec_]
        for k = 1:round(sum(pisoActivoCeldas[j,:]))
            V_edif = [[p1_vec_ dca.ALTURAPISO * (k - 1)];
                  [p2_vec_ dca.ALTURAPISO * (k - 1)];
                  [p3_vec_ dca.ALTURAPISO * (k - 1)];
                  [p4_vec_ dca.ALTURAPISO * (k - 1)];
                  [p1_vec_ dca.ALTURAPISO * k];
                  [p2_vec_ dca.ALTURAPISO * k];
                  [p3_vec_ dca.ALTURAPISO * k];
                  [p4_vec_ dca.ALTURAPISO * k] ]
            fig, ax, ax_mat = polyShape.plotPolyshape3d_v2(PolyShape([V_edif],1), dca.ALTURAPISO * (k - 1), fig, ax, ax_mat, "teal", 1)
            fig, ax, ax_mat = polyShape.plotPolyshape3d_v2(PolyShape([V_tapa],1), dca.ALTURAPISO * k, fig, ax, ax_mat, "teal", 1)
        end

    end

    return fig, ax, ax_mat

end


text:
# Custom composite types.

### 1-D input kernels for modeling a scalar-valued state function of
#   time.

abstract type ODEKernelType end

struct GaussianKernelODEType{T} <: ODEKernelType
    ϵ_sq::T
    s::T # amplitude scale.
end

### elemntary Kernels.

abstract type StationaryKernelType end

# first number is dimension, second is differentiability.
struct Spline12KernelType{T} <: StationaryKernelType
    a::T
end

struct Spline32KernelType{T} <: StationaryKernelType
    a::T
end

struct Spline34KernelType{T} <: StationaryKernelType
    a::T
end

struct RationalQuadraticKernelType{T} <: StationaryKernelType
    a::T
end

struct TunableRationalQuadraticKernelType{T} <: StationaryKernelType
    a::T
    w::T # variance.
end

struct ModulatedSqExpKernelType{T} <: StationaryKernelType
    ϵ_sq::T # decay.
    ν::T # frequency.
end

struct ModulatedMultivariateSqExpKernelType{T}
    ϵ_sq::T # decay.
    ν::Vector{T} # frequency.
end

struct AdaptiveModulatedSqExpKernelType{T} <: StationaryKernelType
    ϵ_sq::T # decay.
    ν::Vector{T} # frequency.
    ϕ::Function # warp map for Sq Exp.
    ψ::Function # warp map for cosine.
end

### Non-canonical kernels
struct KRWarpKernelType{KTv,KTx}
    θ_v::KTv
    θ_x::KTx
end

struct GaussianKernel1DType{T} <: StationaryKernelType
    ϵ_sq::T
end

struct DoubleProductKernelType{KT}
    θ_base::KT
    M::Int
end

abstract type BrownianBridgeKernelType end

struct BrownianBridge10{T} <: BrownianBridgeKernelType
    a::T
end

struct BrownianBridge20{T} <: BrownianBridgeKernelType
    a::T
end

struct BrownianBridge20CompactDomain{T} <: BrownianBridgeKernelType
    a::Vector{T} # lower bound of domain.
    b::Vector{T} # upper bound of domain.
end

# struct BrownianBridgeβ0 <: BrownianBridgeKernelType
#     β::Int
# end

struct BrownianBridge1ϵ{T} <: BrownianBridgeKernelType
    ϵ::T
end

struct BrownianBridge2ϵ{T} <: BrownianBridgeKernelType
    ϵ::T
end

struct BrownianBridgeSemiInfDomain{BT <: BrownianBridgeKernelType}
    θ_base::BT
end

struct BrownianBridgeCompactDomain{BT <: BrownianBridgeKernelType, T}
    θ_base::BT
    a::T # common lower bound for all input dimensions.
    b::T # common upper bound for all input dimensions.
end

### adaptive kernel.

struct AdaptiveKernelDPPType{KT}
    canonical_params::KT # Canonical kernel parameters.
    warpfunc::Function
end

struct AdaptiveKernelType{KT}
    canonical_params::KT # Canonical kernel parameters.
    warpfunc::Function
end

struct FastAdaptiveKernelType{KT, T, D}
    canonical_kernel::KT # Canonical kernel parameters.
    warpfuncs::Vector{Function}
    w_X::Array{T,D} # pre-computed warp map evaluations at training positions.
    s::Vector{T} # weights for each warp function.
end

struct AdaptiveKernelMultiWarpType{KT, T}
    canonical_params::KT # Canonical kernel parameters.
    warpfuncs::Vector{Function} # length M.

    # length M. Positive weights for the m-th warpfunc
    #   to the diagonal entries of the kernel matrix.
    a::Vector{T}
end

struct AdaptiveKernelMultiWarpDPPType{KT, T}
    canonical_params::KT # Canonical kernel parameters.
    warpfuncs::Vector{Function} # length M.

    # length M. Positive weights for the m-th warpfunc
    #   to the diagonal entries of the kernel matrix.
    a::Vector{T}
    self_gain::T
end

# sum of kernels. Weights should be positive.
"""
out = θ.base_gain * evalkernel(p, q, θ.base_kernel_parameters) + δ(p,q) * additive_function(p).
δ(p,q) is kronecker delta.

Type fields:
base_kernel_parameters::KT
additive_function::Function # takes input on 𝓧 to ℝ_{++}.
base_gain::T
zero_tol::T
"""
struct AdditiveVarianceKernelType{KT, T}
    base_kernel_parameters::KT
    additive_function::Function # takes input on 𝓧 to ℝ_{++}.
    base_gain::T
    zero_tol::T
end

struct AdaptiveKernelΨType{KT}
    canonical_params::KT # Canonical kernel parameters.
    ψ::Function
end

struct BrownianBridge20Adaptive
    β::Int # set to 2 for now.
    warpfunc::Function # must be between 0 and 1
end

struct SOSKernelType{KT}
    base_params::KT
end

struct ElementaryKRKernelType{KT,T}
    θ_a::AdaptiveKernelType{KT}
    θ_canonical::RationalQuadraticKernelType{T}
end

struct KRKernelType{KT,T}
    θ_EKR::Vector{ElementaryKRKernelType{KT,T}}
end

### for fitting conditional pdf, univaraite.
# p(A | B), B could be multivariate, A is univariate.
struct CPDF1DKernelType{T, KT}
    A::RationalQuadraticKernelType{T}
    B::AdaptiveKernelType{KT}
end

# struct CPDF1DKernelType{KT_B,KT_A}
#     B::KT_B
#     A::KT_A
# end

### Squared kernel, from multinomial theorem.

struct SquaredFunctionsKernelType{KT}
    base_params::KT
end

### approximation kernel.
# struct ApproxKernelType{Kernel_Type1, Kernel_Type2, T, X_Type}
#     c::Vector{Vector{T}} # solution of the problem.
#     X::Vector{Vector{X_Type}} # sampling locations.
#     θ1::Kernel_Type1 # outer kernel
#     θ2::Kernel_Type2 # inner kernel
#     σ²::T
# end

### warpmaps.

### For RKHS regularization problems.

# classic RKHS
struct RKHSProblemType{Kernel_Type,T,X_Type}
    c::Vector{T} # solution of the problem.
    X::Vector{X_Type} # sampling locations.
    θ::Kernel_Type
    σ²::T
end

struct AdaptiveRKHSProblemType{Kernel_Type,T,X_Type}
    c::Vector{T} # solution of the problem.
    X::Vector{X_Type} # sampling locations.
    θ::Kernel_Type
    σ²::Vector{T}
end

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
struct PhyloTree
    node::Taxon
    children::Vector{PhyloTree}
end

AbstractTrees.children(ptree::PhyloTree) = ptree.children
AbstractTrees.nodetype(::PhyloTree) = PhyloTree
AbstractTrees.printnode(io::IO, ptree::PhyloTree) = print(io, "$(ptree.node)")
Base.show(io::IO,ptree::PhyloTree) = print_tree(io, ptree)

function topolgoy(taxa::Vector{Taxon}; intermediate=false)
    root = lca(taxa)
    lineages = map(Lineage, taxa)
    all_taxon = union(lineages...)
    branches = Dict{Taxon, Taxon}()
    for node in all_taxon
        parent_node = parent(node)
        if parent_node !== nothing
            branches[node] = parent_node
        end
    end
    return _Phylotree(root, branches; intermediate=intermediate)
end

function _Phylotree(root::Taxon,branches::Dict{Taxon,Taxon}; intermediate=false)
    children_node = findall(isequal(root),branches)
    if (!intermediate) &&  (length(children_node) == 1)
        root = children_node[1]
        return _Phylotree(root, branches; intermediate=intermediate)
    end
    children = map(x -> _Phylotree(x,branches;intermediate=intermediate),children_node)
    return PhyloTree(root, children)
end


text:
using SingularSpectrumAnalysis, Test
L = 20
K = 10
N = K*L;
t = 1:N;
T = 20;
y = sin.(2pi/T*t);            # Add seasons
y .+= (0.5sin.(2pi/T*4*t)).^2 # Add seasons
y .+= LinRange(0,1,N)         # Add trend
e = 0.1randn(N);
yn = y+e;                     # Add noise

yt, ys = analyze(yn, L) # trend and seasons
using ControlSystemIdentification
pd  = PredictionData(yt,ys, trend_order=1, ar_order=2)
yth = trend(pd)
ysh = seasons(pd)
@test yth isa Vector{Float64}
@test ysh isa Vector{<:Vector}
@test length(yth) == N
@test length(ysh[1]) == N-2
@test length(ysh) >= 1
plot(pd)

pd  = pred(pd,2) # pd now contains extended fields with the predictions in the end
plot(pd)
yth = trend(pd)
ysh = seasons(pd)

@test yth isa Vector{Float64}
@test ysh isa Vector{<:Vector}
@test length(yth) == N+2
@test length(ysh[1]) == N

@test trend(pd,2) == trend(pd)[2]
@test seasons(pd,2) == getindex.(seasons(pd),2)


text:
#  _____ _____ ____ _____ ____        ____        ____  ___   ___  ____
# |_   _| ____/ ___|_   _/ ___|      |  _ \      / ___|/ _ \ / _ \|  _ \
#   | | |  _| \___ \ | | \___ \ _____| |_) |____| |  _| | | | | | | | | |
#   | | | |___ ___) || |  ___) |_____|  _ <_____| |_| | |_| | |_| | |_| |
#   |_| |_____|____/ |_| |____/      |_| \_\     \____|\___/ \___/|____/

# NOTE TO READERS: These tests validate against the json generation
# functions. For tests dealing directly with DynamoDB, check out the
# integration tests.

## GET_ITEM

# helper to make get_item_query_dict match get_item's iterface
get_item_dict(table :: DynamoTable, key, range=nothing;
              consistant_read=true, only_returning=nothing :: Union{Void, Array{DynamoReference}}) =
    DynamoDB.get_item_query_dict(table, key, range, consistant_read, only_returning)

@test get_item_dict(foo_basic, "asdf") ==
    Dict("TableName" => "foo_basic", "ConsistentRead" => true,
         "Key" => Dict("a" => Dict("S" => "asdf")))

@test get_item_dict(foo_range, "asdf", 3) ==
    Dict("TableName" => "foo_range", "ConsistentRead" => true,
         "Key" => Dict("a" => Dict("S" => "asdf"), "b" => Dict("N" => "3")))

@test get_item_dict(foo_basic, "asdf"; consistant_read=false) ==
    Dict("TableName" => "foo_basic", "ConsistentRead" => false,
         "Key" => Dict("a" => Dict("S" => "asdf")))

@test get_item_dict(foo_basic, "asdf"; consistant_read=false, only_returning=[attr("a")]) ==
    Dict("TableName" => "foo_basic", "ConsistentRead" => false, "ProjectionExpression" => "#1",
         "ExpressionAttributeNames" => Dict("#1" => "a"), "Key" => Dict("a" => Dict("S" => "asdf")))

## BATCH_GET_ITEM

@test DynamoDB.batch_get_item_dict([DynamoDB.batch_get_item_part(foo_basic, 1, 2)]) ==
    Dict("RequestItems" => Dict("foo_basic" => Dict("ConsistentRead" => true,
                                                    "Keys" => [Dict("a" => Dict("N" => "1")),
                                                               Dict("a" => Dict("N" => "2"))])))

@test DynamoDB.batch_get_item_dict([DynamoDB.batch_get_item_part(foo_range, (1, 2), (3, 4))]) ==
    Dict("RequestItems" => Dict("foo_range" =>
                                Dict("ConsistentRead" => true,
                                     "Keys" => [Dict("a" => Dict("N" => "1"), "b" => Dict("N" => "2")),
                                                Dict("a" => Dict("N" => "3"), "b" => Dict("N" => "4"))])))

@test DynamoDB.batch_get_item_dict([DynamoDB.batch_get_item_part(foo_basic, 1, 2),
                                    DynamoDB.batch_get_item_part(foo_range, (1, 2), (3, 4))]) ==
    Dict("RequestItems"=>Dict("foo_basic"=>Dict("Keys"=>[Dict("a"=>Dict("N"=>"1")),
                                                         Dict("a"=>Dict("N"=>"2"))],
                                                "ConsistentRead"=>true),
                              "foo_range"=>Dict("Keys"=>[Dict("a"=>Dict("N"=>"1"),"b"=>Dict("N"=>"2")),
                                                         Dict("a"=>Dict("N"=>"3"),"b"=>Dict("N"=>"4"))],
                                                "ConsistentRead"=>true)))

## PUT ITEM

@test DynamoDB.put_item_dict(foo_basic, Foo(1, 2)) ==
    Dict("TableName" => "foo_basic",
         "Item" => Dict("a" => Dict("N" => "1"),
                        "b" => Dict("N" => "2")))

@test DynamoDB.put_item_dict(foo_basic, Foo(1, 2); conditions=attr("b") < 2) ==
    Dict("TableName" => "foo_basic",
         "Item" => Dict("a" => Dict("N" => "1"),
                        "b" => Dict("N" => "2")),
         "ConditionExpression" => "(#1) < (:2)",
         "ExpressionAttributeNames" => Dict("#1" => "b"),
         "ExpressionAttributeValues" => Dict(":2" => Dict("N" => "2")))

@test DynamoDB.put_item_dict(foo_basic, Foo(1, 2); conditions=attr("b") < 2, return_old=true) ==
    Dict("TableName" => "foo_basic",
         "Item" => Dict("a" => Dict("N" => "1"),
                        "b" => Dict("N" => "2")),
         "ReturnValues" => "ALL_OLD",
         "ConditionExpression" => "(#1) < (:2)",
         "ExpressionAttributeNames" => Dict("#1" => "b"),
         "ExpressionAttributeValues" => Dict(":2" => Dict("N" => "2")))

## BATCH WRITE ITEM

@test DynamoDB.batch_write_item_dict([DynamoDB.batch_delete_part(foo_basic, 1, 2)]) ==
    [Dict("foo_basic" => [Dict("DeleteRequest" =>
                               Dict("Key" => Dict("a" => Dict("N" => "1")))),
                          Dict("DeleteRequest" =>
                               Dict("Key" => Dict("a" => Dict("N" => "2"))))])]

@test DynamoDB.batch_write_item_dict([DynamoDB.batch_put_part(foo_basic, Foo(1, 2), Foo(3, 4))]) ==
    [Dict("foo_basic" => [Dict("PutRequest" =>
                               Dict("Item" => Dict("a" => Dict("N" => "1"),
                                                   "b" => Dict("N" => "2")))),
                          Dict("PutRequest" =>
                               Dict("Item" => Dict("a" => Dict("N" => "3"),
                                                   "b" => Dict("N" => "4"))))])]

@test DynamoDB.batch_write_item_dict([DynamoDB.batch_delete_part(foo_basic, 1, 2),
                                      DynamoDB.batch_put_part(foo_basic, Foo(1, 2), Foo(3, 4))]) ==
    [Dict("foo_basic" => [Dict("DeleteRequest" =>
                               Dict("Key" => Dict("a" => Dict("N" => "1")))),
                          Dict("DeleteRequest" =>
                               Dict("Key" => Dict("a" => Dict("N" => "2")))),
                          Dict("PutRequest" =>
                               Dict("Item" => Dict("a" => Dict("N" => "1"),
                                                   "b" => Dict("N" => "2")))),
                          Dict("PutRequest" =>
                               Dict("Item" => Dict("a" => Dict("N" => "3"),
                                                   "b" => Dict("N" => "4"))))])]

## UPDATE ITEM

@test DynamoDB.update_item_dict(foo_basic, 1, nothing, [assign(attr("a"), 22)]) ==
    Dict("TableName" => "foo_basic",
         "Key" => Dict("a" => Dict("N" => "1")),
         "ReturnValues" => "NONE",
         "UpdateExpression" => "SET #1 = :2",
         "ExpressionAttributeNames" => Dict("#1" => "a"),
         "ExpressionAttributeValues" => Dict(":2" => Dict("N" => "22")))

@test DynamoDB.update_item_dict(foo_basic, 1, nothing, [assign(attr("a"), 22)]; conditions=attr("b") <= 2) ==
    Dict("TableName" => "foo_basic",
         "Key" => Dict("a" => Dict("N" => "1")),
         "ReturnValues" => "NONE",
         "UpdateExpression" => "SET #1 = :2",
         "ConditionExpression" => "(#3) <= (:4)",
         "ExpressionAttributeNames" => Dict("#1" => "a", "#3" => "b"),
         "ExpressionAttributeValues" => Dict(":2" => Dict("N" => "22"), ":4" => Dict("N" => "2")))

## DELETE ITEM

@test DynamoDB.delete_item_dict(foo_basic, 1) ==
    Dict("TableName" => "foo_basic",
         "Key" => Dict("a" => Dict("N" => "1")))

@test DynamoDB.delete_item_dict(foo_basic, 1; conditions=attr("b") > 2) ==
    Dict("TableName" => "foo_basic",
         "Key" => Dict("a" => Dict("N" => "1")),
         "ConditionExpression" => "(#1) > (:2)",
         "ExpressionAttributeNames" => Dict("#1" => "b"),
         "ExpressionAttributeValues" => Dict(":2" => Dict("N" => "2")))

## QUERY

@test DynamoDB.query_dict(foo_range, 77, attr("b") > 17) ==
    Dict("TableName" => "foo_range",
         "KeyConditionExpression" => "((#1) = (:2)) AND ((#3) > (:4))",
         "ExpressionAttributeNames" => Dict("#1" => "a", "#3" => "b"),
         "ExpressionAttributeValues" => Dict(":2" => Dict("N" => "77"),
                                             ":4" => Dict("N" => "17")))

@test DynamoDB.query_dict(foo_range, 77, attr("b") > 17;
                 filter=attr("c") != "cat", scan_index_forward=false, limit=100) ==
    Dict("TableName" => "foo_range",
         "KeyConditionExpression" => "((#1) = (:2)) AND ((#3) > (:4))",
         "FilterExpression" => "(#5) <> (:6)",
         "ScanIndexForward" => false,
         "Limit" => 100,
         "ExpressionAttributeNames" => Dict("#1" => "a", "#3" => "b", "#5" => "c"),
         "ExpressionAttributeValues" => Dict(":2" => Dict("N" => "77"),
                                             ":4" => Dict("N" => "17"),
                                             ":6" => Dict("S" => "cat")))

@test DynamoDB.query_dict(foo_range, 77, attr("b") > 17;
                 projection = [attr("one"), attr("two")], consistant_read=false, index_name="MyIndex",
                 filter=attr("c") != "cat", scan_index_forward=false, limit=100) ==
    Dict("TableName" => "foo_range",
         "KeyConditionExpression" => "((#1) = (:2)) AND ((#3) > (:4))",
         "FilterExpression" => "(#7) <> (:8)",
         "ProjectionExpression" => "#5, #6",
         "ScanIndexForward" => false, "ConsistentRead" => false,
         "IndexName" => "MyIndex", "Limit" => 100,
         "ExpressionAttributeNames" => Dict("#1" => "a", "#3" => "b", "#5" => "one", "#6" => "two", "#7" => "c"),
         "ExpressionAttributeValues" => Dict(":2" => Dict("N" => "77"),
                                             ":4" => Dict("N" => "17"),
                                             ":8" => Dict("S" => "cat")))

## SCAN

@test DynamoDB.scan_dict(foo_range, attr("c") > 17) ==
    Dict("TableName" => "foo_range",
         "FilterExpression" => "(#1) > (:2)",
         "ExpressionAttributeNames" => Dict("#1" => "c"),
         "ExpressionAttributeValues" => Dict(":2" => Dict("N" => "17")))

@test DynamoDB.scan_dict(foo_range, attr("c") == 17;
                projection=[attr("atty1"), attr("atty2")], consistant_read=true, scan_index_forward=false,
                limit=31, segment=3, total_segments=17) ==
    Dict("TableName" => "foo_range",
         "ConsistentRead" => true,
         "ScanIndexForward" => false,
         "Segment" => 3,
         "TotalSegments" => 17,
         "Limit" => 31,
         "FilterExpression" => "(#1) = (:2)",
         "ProjectionExpression" => "#3, #4",
         "ExpressionAttributeNames" => Dict("#1" => "c", "#3" => "atty1", "#4" => "atty2"),
         "ExpressionAttributeValues" => Dict(":2" => Dict("N" => "17")))


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.
# This file includes parts from Combinatorics.jl in modified form.

# Version of: UTC 2019-05-26 10:30:37
# 89df4600-7f90-11e9-3bf6-c39fb02ccd92

# Do not edit this file, it is generated from the modules and will be overwritten!
# Edit the modules in the src directory and build this file with BuildSequences.jl!

tstdir = realpath(joinpath(dirname(@__FILE__)))
srcdir = joinpath(dirname(tstdir), "src")
tstdir ∉ LOAD_PATH && push!(LOAD_PATH, tstdir)
srcdir ∉ LOAD_PATH && push!(LOAD_PATH, srcdir)
module perftests
using IntegerSequences, Dates, InteractiveUtils
InteractiveUtils.versioninfo()
start = Dates.now()
# +++ Abundant.jl +++
GC.gc()
println("\nTEST: L005101(5000)")
@time L005101(5000)
# +++ AndreNumbers.jl +++
GC.gc()
println("\nTEST: (for m in 1:20, n in 0:100 André(m, n) end)")
@time (for m in 1:20, n in 0:100 André(m, n) end)
# +++ BellNumbers.jl +++
GC.gc()
BellNumberList(5)
println("\nTEST: BellNumberList(1000)")
@time BellNumberList(1000)
println("\nTEST: BellTriangle(100, [1 for _ in 1:100])")
@time BellTriangle(100, [1 for _ in 1:100])
# +++ BernoulliNumbers.jl +++
GC.gc()
println("\nTEST: (for m in 1:10, n in 0:499 BernoulliInt(m, n) end)")
@time (for m in 1:10, n in 0:499 BernoulliInt(m, n) end)
println("\nTEST: (for m in 1:10 BernoulliIntList(m, 500) end)")
@time (for m in 1:10 BernoulliIntList(m, 500) end)
GC.gc()
println("\nTEST: (for n in 0:1000 Bernoulli(n) end)")
@time (for n in 0:1000 Bernoulli(n) end)
println("\nTEST: BernoulliList(1000)")
@time BernoulliList(1000)
println("\nTEST: (for n in 0:1000 V027641(n) end)")
@time (for n in 0:1000 V027641(n) end)
GC.gc()
println("\nTEST: (for n in 0:10000 V195441(n) end)")
@time (for n in 0:10000 V195441(n) end)
# +++ BinaryInteger.jl +++
println("\nTEST: [V295513(k) for k in 0:100000]")
@time [V295513(k) for k in 0:100000]
# +++ BinaryQF.jl +++
GC.gc()
println("\nTEST: L034017(10000)")
@time L034017(10000)
# +++ CantorMachines.jl +++
# +++ CarmichaelNumbers.jl +++
println("\nTEST: L = L002997(30)")
@time L = L002997(30)
# +++ ClausenNumbers.jl +++
GC.gc()
println("\nTEST: ClausenNumberList(10000)")
@time ClausenNumberList(10000)
println("\nTEST: (for n in 0:10000 ClausenNumber(n) end)")
@time (for n in 0:10000 ClausenNumber(n) end)
# +++ CombinationsIterator.jl +++
# +++ Compositions.jl +++
GC.gc()
println("\nTEST: L097805(500)")
@time L097805(500)
# +++ Counts.jl +++
println("\nTEST: [A257993(n) for n in 1:10000]")
@time [A257993(n) for n in 1:10000]
GC.gc()
println("\nTEST: PrimePiList(10000)")
@time PrimePiList(10000)
# +++ CyclotomicBinaryForms.jl +++
println("\nTEST: L299214(100)")
@time L299214(100)
# +++ DedekindEta.jl +++
println("\nTEST: PartitionNumberList(10000)")
@time PartitionNumberList(10000)
println("\nTEST: L000731(10000)")
@time L000731(10000)
println("\nTEST: RamanujanTauList(10000)")
@time RamanujanTauList(10000)
# +++ DelehamDelta.jl +++
GC.gc()
println("\nTEST: T225478(100)")
@time T225478(100)
println("\nTEST: T055883(100)")
@time T055883(100)
# +++ EulerTransforms.jl +++
println("\nTEST: (for n in 0:150 V107895(n) end)")
@time (for n in 0:150 V107895(n) end)
# +++ Fibonacci.jl +++
GC.gc()
println("\nTEST: (for n in 1:1000 V000045(n) end)")
@time (for n in 1:1000 V000045(n) end)
println("\nTEST: (for fib in I000045(10000) end)")
@time (for fib in I000045(10000) end)
println("\nTEST: L000045(10000)")
@time L000045(10000)
# +++ FigurativeNumbers.jl +++
println("\nTEST: (for n in 1:1000 V000326(n) end)")
@time (for n in 1:1000 V000326(n) end)
# +++ GaussFactorials.jl +++
# +++ GeneralizedBinomial.jl +++
println("\nTEST: (for n in 0:10000 Binomial(2n, n) end)")
@time (for n in 0:10000 Binomial(2n, n) end)
println("\nTEST: (for n in -100:100, k in -100:100 Binomial(n, k) end)")
@time (for n in -100:100, k in -100:100 Binomial(n, k) end)
println("\nTEST: (for k in -10000:10000 Binomial(-5, k) end)")
@time (for k in -10000:10000 Binomial(-5, k) end)
# +++ HighlyAbundant.jl +++
GC.gc()
println("\nTEST: (for n in 1:100 V034885(n) end)")
@time (for n in 1:100 V034885(n) end)
# +++ Hyper1F1.jl +++
GC.gc()
println("\nTEST: (for n in 0:150 V000262(n) end)")
@time (for n in 0:150 V000262(n) end)
# +++ JacobiTheta.jl +++
println("\nTEST: L000143(100000)")
@time L000143(100000)
println("\nTEST: L035016(100000)")
@time L035016(100000)
# +++ Kolakoski.jl +++
println("\nTEST: I000002(100000)")
@time I000002(100000)
println("\nTEST: KolakoskiList(100000)")
@time KolakoskiList(100000)
# +++ NarayanaCows.jl +++
println("\nTEST: (for cow in NarayanasCows(10000) end)")
@time (for cow in NarayanasCows(10000) end)
# +++ NumberTheory.jl +++
println("\nTEST: [Divisors(n) for n in 1:10000]")
@time [Divisors(n) for n in 1:10000]
println("\nTEST: [Radical(n)  for n in 1:10000]")
@time [Radical(n)  for n in 1:10000]
# +++ OEISUtils.jl +++
# +++ OrthoPolynomials.jl +++
GC.gc()
println("\nTEST: T111062(500)")
@time T111062(500)
GC.gc()
println("\nTEST: T066325(500)")
@time T066325(500)
GC.gc()
println("\nTEST: T053120(500)")
@time T053120(500)
# +++ Partitions.jl +++
# +++ PrimesIterator.jl +++
println("\nTEST: primepi = PrimePi(1000000)")
@time primepi = PrimePi(1000000)
# +++ PrimeSwingFactorial.jl +++
GC.gc()
println("\nTEST: Base.factorial(BigInt(1000000))")
@time Base.factorial(BigInt(1000000))
println("\nTEST: PSfactorial(1000000)")
@time PSfactorial(1000000)
# +++ Products.jl +++
GC.gc()
println("\nTEST: (for n in 1:10000 F!(n) end)")
@time (for n in 1:10000 F!(n) end)
GC.gc()
println("\nTEST: (for n in 1:1000 V000407(n) end)")
@time (for n in 1:1000 V000407(n) end)
GC.gc()
println("\nTEST: (for n in 1:200, k in 1:200 RisingFactorial(n, k) end)")
@time (for n in 1:200, k in 1:200 RisingFactorial(n, k) end)
# +++ QueensProblems.jl +++
println("\nTEST: L319284(15)")
@time L319284(15)
# +++ RationalTrees.jl +++
GC.gc()
println("\nTEST: EuclidTree(100)")
@time EuclidTree(100)
GC.gc()
println("\nTEST: [SchinzelSierpinskiEncoding(l) for l in EuclidTree(100)]")
@time [SchinzelSierpinskiEncoding(l) for l in EuclidTree(100)]
# +++ RecordSearch.jl +++
# +++ RiordanSquares.jl +++
# +++ SelfConvolutive.jl +++
GC.gc()
println("\nTEST: L005411(500)")
@time L005411(500)
# +++ SeqUtils.jl +++
# +++ SeriesExpansion.jl +++
println("\nTEST: L068875(1000)")
@time L068875(1000)
# +++ StirlingLahNumbers.jl +++
GC.gc()
println("\nTEST: (for row in I132393(500) end)")
@time (for row in I132393(500) end)
GC.gc()
println("\nTEST: (for row in I048993(500) end)")
@time (for row in I048993(500) end)
GC.gc()
println("\nTEST: (for row in I271703(500) end)")
@time (for row in I271703(500) end)
# +++ SwingFactorial.jl +++
GC.gc()
n = 1000000
println("\nTEST: Base.factorial(BigInt(n))")
@time Base.factorial(BigInt(n))
println("\nTEST: Sfactorial(n)")
@time Sfactorial(n)
println("n = $n")
# +++ Triangles.jl +++
# +++ UlamNumbers.jl +++
GC.gc()
for n in 6:10
println("\nTEST: UlamList(2^n)")
@time UlamList(2^n)
end
# +++ ZumkellerNumbers.jl +++
GC.gc()
println("\nTEST: (for n in 1:2000 isZumkeller(n) end)")
@time (for n in 1:2000 isZumkeller(n) end)
println("\nTEST: L083207(500)")
@time L083207(500)
stop = Dates.now()
tdiff = stop - start
println("\nJulia version: " * string(VERSION) )
println(start)
println("Total test time: ", tdiff)
end # module


text:
function (store::CSVExporter{<:AbstractSchema})(path::Path, li::LabelledIntervals)
    # Convert labelled intervals to a dataframe.
    labelled_intervals = get_labelled_intervals(li)
    df = DataFrame(label = String[], start = Int64[], stop = Int64[], first = Float64[], step = Float64[], last  = Float64[])
    for (key, labelled_interval) in pairs(labelled_intervals)
        nested_interval = labelled_interval.nested_interval
        start = get_start(nested_interval)
        stop = get_stop(nested_interval)
        interval = get_interval(nested_interval)
        push!(df, [key start stop first(interval) step(interval) last(interval)])
    end
    Base.display(df)
    store(path, df)
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/domain_reduction.jl
# Contains subroutines used for domain reduction.
#############################################################################

"""
$(FUNCTIONNAME)

Tightens the bounds of the `_current_node` using the current global upper bound
and the duality information obtained from the relaxation.
"""
function variable_dbbt!(n::NodeBB, mult_lo::Vector{Float64}, mult_hi::Vector{Float64},
                        LBD::Float64, UBD::Float64, nx::Int64)

    delta = UBD - LBD
    lvbs = n.lower_variable_bounds
    uvbs = n.upper_variable_bounds
    if LBD <= UBD
        for i = 1:nx
            ml = @inbounds mult_lo[i]
            if ml > 0.0
                cut = @inbounds lvbs[i] + delta/ml
                if cut < @inbounds uvbs[i]
                    @inbounds uvbs[i] = cut
                end
            else
                mh = @inbounds mult_hi[i]
                if mh > 0.0
                    cut = @inbounds uvbs[i] - delta/mh
                    if cut > @inbounds lvbs[i]
                        @inbounds lvbs[i] = cut
                    end
                end
            end
         end
    end

    return nothing
end

"""
$(FUNCTIONNAME)

Excludes OBBT on variable indices that are tight for the solution of the relaxation.
"""
function trivial_filtering!(m::Optimizer, n::NodeBB)

    obbt_tolerance = m._parameters.obbt_tolerance
    m._preprocess_termination_status = MOI.get(m.relaxed_optimizer, MOI.TerminationStatus())
    m._preprocess_result_status = MOI.get(m.relaxed_optimizer, MOI.PrimalStatus())
    valid_flag, feasible_flag = is_globally_optimal(m._preprocess_termination_status,
                                                    m._preprocess_result_status)

    if valid_flag && feasible_flag
        for j = 1:length(m._obbt_working_lower_index)
            if @inbounds m._obbt_working_lower_index[j]
                vi = @inbounds m._relaxed_variable_index[j]
                diff = MOI.get(m.relaxed_optimizer, MOI.VariablePrimal(), vi)
                diff -= @inbounds n.lower_variable_bounds[j]
                if abs(diff) <= obbt_tolerance
                    @inbounds m._obbt_working_lower_index[j] = false
                end
            end
        end
        for j = 1:length(m._obbt_working_upper_index)
            if @inbounds m._obbt_working_upper_index[j]
                vi = @inbounds m._relaxed_variable_index[j]
                diff = -MOI.get(m.relaxed_optimizer, MOI.VariablePrimal(), vi)
                diff += @inbounds n.upper_variable_bounds[j]
                if abs(diff) <= obbt_tolerance
                    @inbounds m._obbt_working_upper_index[j] = false
                end
            end
        end
    end

    return
end

"""
$(FUNCTIONNAME)

Utility function used to set vector of booleans z to x & ~y. Avoids the
generation of conversion of the BitArray created by broadcasting logical operators.
"""
function bool_indx_diff(z::Vector{Bool},x::Vector{Bool}, y::Vector{Bool})
    for i = 1:length(z)
        @inbounds z[i] = (x[i] & ~y[i])
    end
    return
end

"""
$(FUNCTIONNAME)

Excludes OBBT on variable indices after a search in a filtering direction.
"""
function aggressive_filtering!(m::Optimizer, n::NodeBB)

    # Initial filtering vector (negative one direction per remark in Gleixner2017)
    variable_number = m._working_problem._variable_count
    v = -ones(variable_number)

    # Copy prior index set (ignores linear and binary terms)
    obbt_variable_count = m._obbt_variable_count
    copyto!(m._old_low_index, m._obbt_working_lower_index)
    copyto!(m._old_upp_index, m._obbt_working_upper_index)
    copyto!(m._new_low_index, m._obbt_working_lower_index)
    copyto!(m._new_upp_index, m._obbt_working_upper_index)

    # Exclude unbounded directions
    for i = 1:obbt_variable_count
        if @inbounds m._new_low_index[i] && @inbounds n.lower_variable_bounds[i] === -Inf
            @inbounds m._new_low_index[i] = false
        end
        if @inbounds m._new_low_index[i] && @inbounds n.upper_variable_bounds[i] === Inf
            @inbounds m._new_low_index[i] = false
        end
    end

    # Begin the main algorithm
    for k = 1:m._parameters.obbt_aggressive_max_iteration

        # Set index differences and vector for filtering direction
        bool_indx_diff(m._lower_indx_diff, m._old_low_index, m._new_low_index)
        bool_indx_diff(m._upper_indx_diff, m._old_upp_index, m._new_upp_index)

        for i = 1:obbt_variable_count
            vi = @inbounds v[i]
            if @inbounds m._lower_indx_diff[i] && vi < 0.0
                @inbounds v[i] = 0.0
            end
            if @inbounds m._upper_indx_diff[i] && vi > 0.0
                @inbounds v[i] = 0.0
            end
        end

        # Termination Condition
        ((~any(m._new_low_index) & ~any(m._new_upp_index)) || (iszero(v))) && break
        if k >= 2
            if (count(m._lower_indx_diff) + count(m._upper_indx_diff)) < m._parameters.obbt_aggressive_min_dimension
                break
            end
        end

        # Set objective in OBBT problem to filtering vector
        MOI.set(m.relaxed_optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        saf = SAF(SAT.(v, m._relaxed_variable_index), 0.0)
        MOI.set(m.relaxed_optimizer, MOI.ObjectiveFunction{SAF}(), saf)

        # Optimizes the problem and if successful filter additional bounds
        MOI.optimize!(m.relaxed_optimizer)

        m._preprocess_termination_status = MOI.get(m.relaxed_optimizer, MOI.TerminationStatus())
        m._preprocess_result_status = MOI.get(m.relaxed_optimizer, MOI.PrimalStatus())
        valid_flag, feasible_flag = is_globally_optimal(m._preprocess_termination_status,
                                                        m._preprocess_result_status)

        if valid_flag && feasible_flag
            variable_primal = MOI.get(m.relaxed_optimizer, MOI.VariablePrimal(), m._relaxed_variable_index)
            copyto!(m._new_low_index, m._old_low_index)
            copyto!(m._new_upp_index, m._old_upp_index)
            for i = 1:obbt_variable_count
                vp_value =  @inbounds variable_primal[i]
                if @inbounds m._old_low_index[i] && vp_value == @inbounds n.lower_variable_bounds[i]
                    @inbounds m._new_low_index[i] = false
                end
                if @inbounds m._old_upp_index[i] && vp_value == @inbounds n.upper_variable_bounds[i]
                    @inbounds m._new_upp_index[i] = false
                end
            end
        else
            return false
        end
    end
    copyto!(m._obbt_working_lower_index, m._new_low_index)
    copyto!(m._obbt_working_upper_index, m._new_upp_index)
    return true
end

"""
$(FUNCTIONNAME)
"""
function set_node_flag!(m::Optimizer)
    for constr in m._working_problem._nonlinear_constr
        set_node_flag!(constr)
    end

    if m._working_problem._objective_type === NONLINEAR
        set_node_flag!(m._working_problem._objective_nl)
    end

    return nothing
end

"""
$(FUNCTIONNAME)
"""
function set_reference_point!(m::Optimizer)

    evaluator = m._working_problem._relaxed_evaluator
    evaluator_x = evaluator.x
    current_xref = m._current_xref

    new_reference_point = false
    for node_i = 1:m._branch_variable_count
        solution_i = m._branch_to_sol_map[node_i]

        node_x = current_xref[node_i]
        solution_x = evaluator_x[solution_i]

        if node_x !== solution_x
            evaluator_x[solution_i] = node_x
            new_reference_point = true
        end
    end

    if new_reference_point

        for constr in m._working_problem._nonlinear_constr
            constr.has_value = false
        end

        if m._working_problem._objective_type === NONLINEAR
            m._working_problem._objective_nl.has_value = false
        end
    end
    fill!(evaluator.subexpressions_eval, false)

    return nothing
end

"""
$(FUNCTIONNAME)

Performs OBBT with filtering and greedy ordering as detailed in:
Gleixner, A.M., Berthold, T., Müller, B. et al. J Glob Optim (2017) 67: 731.
https://doi.org/10.1007/s10898-016-0450-4
"""
function obbt!(m::Optimizer)

    feasibility = true

    n = m._current_node
    branch_to_sol_map = m._branch_to_sol_map
    relaxed_optimizer = m.relaxed_optimizer

    # set node and reference point if necessary then solve initial problem to
    # feasibility. This is repeated `obbt_repetitions` number of times in the
    # following fashion. Relax the problem, populate affine constraints, run
    # obbt which contracts variable bounds, delete affine constraints...
    # update variable bounds and repeat. TODO: Keep track of which variables
    # participate in which functions and only delete a constraint if a
    # variable participating in a nonlinear term changes it bounds.
    if m._obbt_performed_flag
        reset_relaxation!(m)
        set_first_relax_point!(m)
    end

    update_relaxed_problem_box!(m)
    if m._nonlinear_evaluator_created
        set_node!(m._working_problem._relaxed_evaluator, n)
        set_node_flag!(m)
        set_reference_point!(m)
    end
    relax_constraints!(m, 1)
    relax_objective!(m, 1)
    MOI.set(relaxed_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(relaxed_optimizer)

    # Sets indices to attempt OBBT on
    obbt_variable_count = m._obbt_variable_count
    fill!(m._obbt_working_lower_index, true)
    fill!(m._obbt_working_upper_index, true)

    # Filters out any indicies with active bounds on variables
    # determined by solving the feasibility problem
    trivial_filtering!(m, n)

    # Applies an aggressive filter to eliminate indices that
    # cannot be tightened by obbt
    if m._parameters.obbt_aggressive_on
        feasibility = aggressive_filtering!(m, n)
    end

    # extracts info from relaxed problem (redundant if aggressive obbt is called)
    m._preprocess_termination_status = MOI.get(relaxed_optimizer, MOI.TerminationStatus())
    m._preprocess_result_status = MOI.get(relaxed_optimizer, MOI.PrimalStatus())
    valid_flag, feasible_flag = is_globally_optimal(m._preprocess_termination_status,
                                                    m._preprocess_result_status)

    if valid_flag && feasible_flag
        xLP = MOI.get(relaxed_optimizer, MOI.VariablePrimal(), m._relaxed_variable_index)
    else
        return false
    end

    # continue tightening bounds by optimization until all indices have been checked
    # or the node is empty and the problem is thus proven infeasible
    while (any(m._obbt_working_lower_index) || any(m._obbt_working_upper_index)) && !isempty(n)

        # Get lower value
        lower_indx = -1;     upper_indx = -1
        lower_value = Inf;   upper_value = Inf

        # min of xLP - yL on active
        if any(m._obbt_working_lower_index)
            for i = 1:obbt_variable_count
                if @inbounds m._obbt_working_lower_index[i]
                    sol_indx = branch_to_sol_map[i]
                    temp_value = @inbounds xLP[sol_indx] - n.lower_variable_bounds[i]
                    # Need less than or equal to handle unbounded cases
                    if temp_value <= lower_value
                        lower_value = temp_value
                        lower_indx = i
                    end
                end
            end
        end

        # min of yU - xLP on active
        if any(m._obbt_working_upper_index)
            for i = 1:obbt_variable_count
                if @inbounds m._obbt_working_upper_index[i]
                    sol_indx = branch_to_sol_map[i]
                    temp_value = @inbounds n.upper_variable_bounds[i] - xLP[sol_indx]
                    if temp_value <= upper_value
                        upper_value = temp_value
                        upper_indx = i
                    end
                end
            end
        end

        # default to upper bound if no lower bound is found, use maximum distance otherwise
        if lower_value <= upper_value && lower_indx > 0

            @inbounds m._obbt_working_lower_index[lower_indx] = false
            var = SV(m._relaxed_variable_index[lower_indx])

            MOI.set(relaxed_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
            MOI.set(relaxed_optimizer, MOI.ObjectiveFunction{SV}(), var)

            MOI.optimize!(m.relaxed_optimizer)
            m._preprocess_termination_status = MOI.get(relaxed_optimizer, MOI.TerminationStatus())
            m._preprocess_result_status = MOI.get(relaxed_optimizer, MOI.PrimalStatus())
            valid_flag, feasible_flag = is_globally_optimal(m._preprocess_termination_status,
                                                            m._preprocess_result_status)

            if valid_flag && feasible_flag
                xLP .= MOI.get(relaxed_optimizer, MOI.VariablePrimal(), m._relaxed_variable_index)

                node_index = branch_to_sol_map[lower_indx]
                updated_value = xLP[node_index]
                previous_value = n.lower_variable_bounds[lower_indx]

                # if bound is improved update node and corresponding constraint update
                # the node bounds and the single variable bound in the relaxation
                # we assume branching does not occur on fixed variables and interval
                # constraints are internally bridged by EAGO. So the only L <= x
                # constraint in the model is a GreaterThan.
                if updated_value > previous_value && (updated_value - previous_value) > 1E-6
                    ci_list = MOI.get(relaxed_optimizer, MOI.ListOfConstraintIndices{SAF,LT}())
                    sv_geq_ci = m._node_to_sv_geq_ci[lower_indx]
                    MOI.set(relaxed_optimizer, MOI.ConstraintSet(), sv_geq_ci, GT(updated_value))
                    @inbounds n.lower_variable_bounds[lower_indx] = updated_value
                end

                if isempty(n)
                    feasibility = false
                    break
                end

            elseif valid_flag && !feasible_flag
                feasibility = false
                break

            else
                break
            end

        elseif upper_indx > 0

            m._obbt_working_upper_index[upper_indx] = false
            var = SV(m._relaxed_variable_index[upper_indx])
            MOI.set(relaxed_optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)
            MOI.set(relaxed_optimizer, MOI.ObjectiveFunction{SV}(), var)
            MOI.optimize!(relaxed_optimizer)
            m._preprocess_termination_status = MOI.get(m.relaxed_optimizer, MOI.TerminationStatus())
            m._preprocess_result_status = MOI.get(m.relaxed_optimizer, MOI.PrimalStatus())
            valid_flag, feasible_flag = is_globally_optimal(m._preprocess_termination_status,
                                                            m._preprocess_result_status)

            if valid_flag && feasible_flag
                xLP .= MOI.get(relaxed_optimizer, MOI.VariablePrimal(), m._relaxed_variable_index)
                node_index = branch_to_sol_map[upper_indx]
                updated_value = xLP[node_index]
                previous_value = n.upper_variable_bounds[upper_indx]

                # if bound is improved update node and corresponding constraint update
                # the node bounds and the single variable bound in the relaxation
                # we assume branching does not occur on fixed variables and interval
                # constraints are internally bridged by EAGO. So the only U => x
                # constraint in the model is a LessThan.
                if updated_value < previous_value && (previous_value - updated_value) > 1E-6
                    sv_leq_ci = m._node_to_sv_leq_ci[upper_indx]
                    MOI.set(relaxed_optimizer, MOI.ConstraintSet(), sv_leq_ci, LT(updated_value))
                    @inbounds n.upper_variable_bounds[upper_indx] = updated_value
                end

                if isempty(n)
                    feasibility = false
                    break
                end

            elseif valid_flag && !feasible_flag
                feasibility = false
                break

            else
                break
            end

        else
            break
        end
        trivial_filtering!(m, n)
    end

    return feasibility
end

"""
$(FUNCTIONNAME)
"""
function load_fbbt_buffer!(m::Optimizer)

    n = m._current_node
    sol_to_branch = m._sol_to_branch_map
    lower_variable_bounds = n.lower_variable_bounds
    upper_variable_bounds = n.upper_variable_bounds

    for i = 1:m._working_problem._variable_count
        if @inbounds m._branch_variables[i]
            indx = @inbounds sol_to_branch[i]
            @inbounds m._lower_fbbt_buffer[i] = lower_variable_bounds[indx]
            @inbounds m._upper_fbbt_buffer[i] = upper_variable_bounds[indx]

        else
            @inbounds m._lower_fbbt_buffer[i] = m._working_problem._variable_info[i].lower_bound
            @inbounds m._upper_fbbt_buffer[i] = m._working_problem._variable_info[i].upper_bound

        end
    end

    return nothing
end

"""
$(FUNCTIONNAME)
"""
function unpack_fbbt_buffer!(m::Optimizer)

    n = m._current_node
    sol_to_branch = m._sol_to_branch_map
    lower_variable_bounds = n.lower_variable_bounds
    upper_variable_bounds = n.upper_variable_bounds

    for i = 1:m._working_problem._variable_count
        if m._branch_variables[i]
            indx = sol_to_branch[i]
            if m._lower_fbbt_buffer[i] > lower_variable_bounds[indx]
                lower_variable_bounds[indx] = m._lower_fbbt_buffer[i]
            end
            if upper_variable_bounds[indx] > m._upper_fbbt_buffer[i]
                upper_variable_bounds[indx] = m._upper_fbbt_buffer[i]
            end
        end
    end

    return nothing
end

"""
$(FUNCTIONNAME)

Performs feasibility-based bound tightening on a back-end constraint and returns `true` if it is feasible or
`false` if it is infeasible.
"""
function fbbt! end

function fbbt!(m::Optimizer, f::AffineFunctionIneq)

    # compute full sum
    lower_bounds = m._lower_fbbt_buffer
    upper_bounds = m._upper_fbbt_buffer

    terms = f.terms
    temp_sum = -f.constant

    for k = 1:f.len

        aik, indx_k = @inbounds terms[k]
        if aik !== 0.0
            aik_xL = aik*(@inbounds lower_bounds[indx_k])
            aik_xU = aik*(@inbounds upper_bounds[indx_k])
            temp_sum -= min(aik_xL, aik_xU)
        end
    end

    # subtract extra term, check to see if implied bound is better, if so update the node and
    # the working sum if the node is now empty then break

    for k = 1:f.len

        aik, indx_k = @inbounds terms[k]
        if aik !== 0.0

            xL = @inbounds lower_bounds[indx_k]
            xU = @inbounds upper_bounds[indx_k]

            aik_xL = aik*xL
            aik_xU = aik*xU

            temp_sum += min(aik_xL, aik_xU)
            xh = temp_sum/aik

            if aik > 0.0
                (xh < xL) && return false
                if xh > xL
                    @inbounds upper_bounds[indx_k] = xh
                end

            elseif aik < 0.0
                (xh > xU) && return false
                if xh < xU
                    @inbounds lower_bounds[indx_k] = xh
                end

            else
                temp_sum -= min(aik_xL, aik_xU)
                continue

            end

            aik_xL = aik*(@inbounds lower_bounds[indx_k])
            aik_xU = aik*(@inbounds upper_bounds[indx_k])
            temp_sum -= min(aik_xL, aik_xU)

        end
    end

    return true
end

function fbbt!(m::Optimizer, f::AffineFunctionEq)

    # compute full sum
    lower_bounds = m._lower_fbbt_buffer
    upper_bounds = m._upper_fbbt_buffer

    terms = f.terms
    temp_sum_leq = -f.constant
    temp_sum_geq = -f.constant

    for k = 1:f.len
        aik, indx_k = @inbounds terms[k]

        if aik !== 0.0
            aik_xL = aik*(@inbounds lower_bounds[indx_k])
            aik_xU = aik*(@inbounds upper_bounds[indx_k])
            temp_sum_leq -= min(aik_xL, aik_xU)
            temp_sum_geq -= max(aik_xL, aik_xU)

        end
    end

    # subtract extra term, check to see if implied bound is better, if so update the node and
    # the working sum if the node is now empty then break
    for k = 1:f.len

        aik, indx_k = @inbounds terms[k]
        if aik !== 0.0

            xL = @inbounds lower_bounds[indx_k]
            xU = @inbounds upper_bounds[indx_k]

            aik_xL = aik*xL
            aik_xU = aik*xU

            temp_sum_leq += min(aik_xL, aik_xU)
            temp_sum_geq += max(aik_xL, aik_xU)

            xh_leq = temp_sum_leq/aik
            xh_geq = temp_sum_geq/aik

            if aik > 0.0
                (xh_leq < xL) && return false
                if xh_leq > xL
                    @inbounds upper_bounds[indx_k] = xh_leq
                end
                (xh_geq > xU) && return false
                if (xh_geq < xU)
                    @inbounds lower_bounds[indx_k] = xh_geq
                end

            elseif aik < 0.0
                (xh_leq > xU) && return false
                if xh_leq < xU
                    @inbounds lower_bounds[indx_k] = xh_leq
                end
                (xh_geq < xL) && return false
                if (xh_geq > xL)
                    @inbounds upper_bounds[indx_k] = xh_geq
                end

            else
                temp_sum_leq -= min(aik_xL, aik_xU)
                temp_sum_geq -= max(aik_xL, aik_xU)
                continue

            end
            aik_xL = aik*(@inbounds lower_bounds[indx_k])
            aik_xU = aik*(@inbounds upper_bounds[indx_k])

            temp_sum_leq -= min(aik_xL, aik_xU)
            temp_sum_geq -= max(aik_xL, aik_xU)
        end
    end

    return true
end

cp_condition(m::Optimizer) = false

"""
Performs bound tightening based on forward/reverse interval and/or McCormick passes. This routine
resets the current node with new interval bounds.
"""
function set_constraint_propagation_fbbt!(m::Optimizer)
    feasible_flag = true

    if m._nonlinear_evaluator_created
        evaluator = m._working_problem._relaxed_evaluator
        set_node!(m._working_problem._relaxed_evaluator, m._current_node)
        set_node_flag!(m)
        set_reference_point!(m)

        m._working_problem._relaxed_evaluator.is_first_eval = m._new_eval_constraint
        for constr in m._working_problem._nonlinear_constr
            if feasible_flag
                forward_pass!(evaluator, constr)
                feasible_flag &= reverse_pass!(evaluator, constr)
                evaluator.interval_intersect = true
            end
        end

        m._working_problem._relaxed_evaluator.is_first_eval = m._new_eval_constraint
        for constr in m._working_problem._nonlinear_constr
            if feasible_flag
                forward_pass!(evaluator, constr)
            end
        end

        evaluator.is_post = m._parameters.subgrad_tighten

        m._working_problem._relaxed_evaluator.is_first_eval = m._new_eval_objective
        if feasible_flag && (m._working_problem._objective_type === NONLINEAR)
            obj_nonlinear = m._working_problem._objective_nl
            set_node_flag!(obj_nonlinear)
            forward_pass!(evaluator, obj_nonlinear)
            feasible_flag &= reverse_pass!(evaluator, obj_nonlinear)
            evaluator.interval_intersect = true
        end

        if feasible_flag && (m._working_problem._objective_type === NONLINEAR)
            obj_nonlinear = m._working_problem._objective_nl
            set_node_flag!(obj_nonlinear)
            forward_pass!(evaluator, obj_nonlinear)
        end

        m._new_eval_constraint = false
        m._new_eval_objective = false

        retrieve_x!(m._current_xref, evaluator)
        m._current_node = retrieve_node(evaluator)
    end

    return feasible_flag
end


text:
module RobustAndOptimalControl

using ControlSystems: issiso
using LinearAlgebra, Statistics
using RecipesBase
using ControlSystems
import ControlSystems: ss, ssdata, ninputs, noutputs, nstates, isdiscrete, iscontinuous, to_matrix, timeevol, _string_mat_with_headers, PartionedStateSpace, common_timeevol

using ComponentArrays

using MonteCarloMeasurements, Optim, UnPack
import Distributions: Uniform
import IntervalArithmetic
import IntervalArithmetic: Interval

import MatrixPencils, MatrixEquations, DescriptorSystems
using DescriptorSystems: dss

using ChainRulesCore

export show_construction, vec2sys
include("utils.jl")

export dss, hinfnorm2, h2norm, hankelnorm, nugap, νgap, baltrunc2
include("descriptor.jl")

export ExtendedStateSpace, system_mapping, performance_mapping, noise_mapping, ssdata_e, partition, ss
include("ExtendedStateSpace.jl")

export δ, δr, δc, δss, nominal, UncertainSS, uss, blocksort
include("uncertainty_interface.jl")

export makeweight, neglected_delay, gain_and_delay_uncertainty, neglected_lag, fit_complex_perturbations
include("weights.jl")

export hinfsynthesize, hinfassumptions, hinfpartition, hinfsignals, bilinearc2d, bilineard2c, fudge_inv, hinfgrad
include("hinfinity_design.jl")

include("plotting.jl")

export frequency_weighted_reduction, hsvd
include("reduction.jl")

export h2synthesize
include("h2_design.jl")

export LQGProblem, sensitivity, input_sensitivity, output_sensitivity, comp_sensitivity, input_comp_sensitivity, output_comp_sensitivity, controller, ff_controller, extended_controller, closedloop, static_gain_compensation, G_PS, G_CS, gangoffour, loopgain, stabilityrobustness, returndifference
include("lqg.jl")

export NamedStateSpace, named_ss, expand_symbol, measure, connect, sumblock, splitter
include("named_systems2.jl")

export find_lft, δ
include("find_lft.jl")

export add_disturbance, add_low_frequency_disturbance, add_measurement_disturbance, add_resonant_disturbance, add_output_differentiator, add_output_integrator, add_input_integrator, add_input_differentiator
include("model_augmentation.jl")

export glover_mcfarlane, hanus
include("glover_mcfarlane.jl")

export diskmargin, Diskmargin, Disk, sim_diskmargin, loop_diskmargin, structured_singular_value, broken_feedback, robstab, loop_scaling, loop_scale
include("diskmargin.jl")
include("mimo_diskmargin.jl")

end


text:
module Advection

using Interpolations: GriddedInterpolation

export advect_rk4

"""
    advect_rk4(dt, pos, wrap, u_prev, v_prev, u_inter, v_inter, u_next, v_next)

RK4 advection kernel

Advect a particle located at position (x, y) with a timestep dt. Velocity fields
are functions returning the velocity evaluated (or interpolated) to a given location.
The velocity fields are specified at the current timestep and at the following timestep.
"""
function advect_rk4(dt, pos, wrap,
                    u_prev::GriddedInterpolation, v_prev::GriddedInterpolation,
                    u_next::GriddedInterpolation, v_next::GriddedInterpolation)
    u_inter(pos) = (u_prev.(pos[:,1], pos[:,2]) + u_next.(pos[:,1], pos[:,2])) / 2
    v_inter(pos) = (v_prev.(pos[:,1], pos[:,2]) + v_next.(pos[:,1], pos[:,2])) / 2

    pos = permutedims(pos, (2, 1))

    vel1 = hcat(u_prev.(pos[:,1], pos[:,2]),
                v_prev.(pos[:,1], pos[:,2]))
    pos1 = wrap(pos .+ 0.5 * vel1 * dt)

    vel2 = hcat(u_inter(pos1), v_inter(pos1))
    pos2 = wrap(pos .+ 0.5 * vel2 * dt)

    vel3 = hcat(u_inter(pos2), v_inter(pos2))
    pos3 = wrap(pos .+ vel3 * dt)

    vel4 = hcat(u_next.(pos3[:,1], pos3[:,2]),
                v_next.(pos3[:,1], pos3[:,2]))

    permutedims(wrap(pos .+ (vel1 .+ 2*vel2 .+ 3*vel3 .+ vel4) / 6*dt), (2, 1))
end

end # module


text:
@testset "CLBLAS.L1" begin
    @testset "CLBLAS gemm" begin
        for T in testtypes
            alpha = T(2.0)
            beta = T(0.0)
            hA = rand(T, 64, 32)
            hB = rand(T, 32, 64)
            hC = rand(T, 64, 64)
            A = CLArray(q, hA)
            B = CLArray(q, hB)
            C = CLArray(q, hC)

            cl.wait(gemm!('N', 'N', alpha, A, B, beta, C; queue=q))
            gemm!('N', 'N', alpha, hA, hB, beta, hC)

            @test cl.to_host(C) ≈ hC
        end
    end
end


text:
using Test

@testset verbose = true "TDD" begin
    include("graph_test.jl");
    include("interval_test.jl")
end;