text:
module Test_dplyap

using LinearAlgebra
using MatrixEquations
using Test

@testset "Testing positive discrete Lyapunov equation solvers" begin

n = 30
m = 5
p = 3
Ty = Float64

@testset "Positive discrete Lyapunov equations" begin

reltol = eps(float(100))
a = -.1+0.2im; b = 2im; @time u = plyapd(a,b)
@test abs(a*u*u'*a'-u*u'+b*b') < reltol

reltol = eps(float(100f0))
a = -.1f0+-.1f0*im; b = 2f0im; @time u = plyapd(a,b)
@test abs(a*u*u'*a'-u*u'+b*b')  < reltol

for Ty in (Float64, Float32)

ar = rand(Ty,n,n)
ar = ar/(one(Ty) + norm(ar))
ac = rand(Ty,n,n)+im*rand(Ty,n,n)
ac = ac/(one(Ty) + norm(ac))
br = rand(Ty,n,m)
brw = rand(Ty,n,n+m)
bc = br+im*rand(Ty,n,m)
bcw = brw+im*rand(Ty,n,n+m)
cr = rand(Ty,p,n)
crt = rand(Ty,n+p,n)
cc = cr+im*rand(Ty,p,n)
cct = crt+im*rand(Ty,n+p,n)
Ty == Float64 ? reltol = eps(float(100)) : reltol = eps(100*n*one(Ty))

@time u = plyapd(ar,br);
x = u*u'; @test norm(ar*x*ar'-x+br*br')/norm(x)/max(1.,norm(ar)^2) < reltol

@time u = plyapd(ar,brw);
x = u*u'; @test norm(ar*x*ar'-x+brw*brw')/norm(x)/max(1.,norm(ar)^2) < reltol

@time u = plyapd(ar,I,br);
x = u*u'; @test norm(ar*x*ar'-x+br*br')/norm(x)/max(1.,norm(ar)^2) < reltol

@time u = plyapd(ar',cr');
x = u'*u; @test norm(ar'*x*ar-x+cr'*cr)/norm(x)/max(1.,norm(ar)^2) < reltol

@time u = plyapd(ar',crt');
x = u'*u; @test norm(ar'*x*ar-x+crt'*crt)/norm(x)/max(1.,norm(ar)^2) < reltol

@time u = plyapd(ac,bc);
x = u*u'; @test norm(ac*x*ac'-x+bc*bc')/norm(x)/max(1.,norm(ac)^2) < reltol

@time u = plyapd(ac,bcw);
x = u*u'; @test norm(ac*x*ac'-x+bcw*bcw')/norm(x)/max(1.,norm(ac)^2) < reltol

@time u = plyapd(ac',cc');
x = u'*u; @test norm(ac'*x*ac-x+cc'*cc)/norm(x)/max(1.,norm(ac)^2) < reltol

@time u = plyapd(ac',cct');
x = u'*u; @test norm(ac'*x*ac-x+cct'*cct)/norm(x)/max(1.,norm(ac)^2) < reltol

@time u = plyapd(ac',cr');
x = u'*u; @test norm(ac'*x*ac-x+cr'*cr)/norm(x)/max(1.,norm(ac)^2) < reltol

@time u = plyapd(ar,bc);
x = u*u'; @test norm(ar*x*ar'-x+bc*bc')/norm(x)/max(1.,norm(ar)^2) < reltol

@time u = plyapd(ar',cc');
x = u'*u; @test norm(ar'*x*ar-x+cc'*cc)/norm(x)/max(1.,norm(ar)^2) < reltol

end
end

@testset "Positive discrete generalized Lyapunov equations" begin

reltol = eps(float(100))
a = -1+im; ee = 4; b = 2im; @time u = plyapd(a,ee,b)
@test abs(a*u*u'*a'-ee*u*u'*ee'+b*b') < reltol

reltol = eps(float(100f0))
a = -1f0+1f0*im; ee = 4f0; b = 2f0im; @time u = plyapd(a,ee,b)
@test abs(a*u*u'*a'-ee*u*u'*ee'+b*b')  < reltol

for Ty in (Float64, Float32)

ar = rand(Ty,n,n)
ar = ar/(one(Ty) + norm(ar))
er = rand(Ty,n,n)
ar = er*ar
ac = rand(Ty,n,n)+im*rand(Ty,n,n)
ac = ac/(one(Ty) + norm(ac))
ec = er+im*rand(Ty,n,n)
ac = ec*ac
br = rand(Ty,n,m)
brw = rand(Ty,n,n+m)
bc = br+im*rand(Ty,n,m)
bcw = brw+im*rand(Ty,n,n+m)
cr = rand(Ty,p,n)
crt = rand(Ty,n+p,n)
cc = cr+im*rand(Ty,p,n)
cct = crt+im*rand(Ty,n+p,n)
Ty == Float64 ? reltol = eps(float(100)) : reltol = eps(100*n*one(Ty))

@time u = plyapd(ar',er',cr');
x = u'*u; @test norm(ar'*x*ar-er'*x*er+cr'*cr)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time u = plyapd(ar',er',crt');
x = u'*u; @test norm(ar'*x*ar-er'*x*er+crt'*crt)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time u = plyapd(ar,er,br);
x = u*u'; @test norm(ar*x*ar'-er*x*er'+br*br')/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time u = plyapd(ar,er,brw);
x = u*u'; @test norm(ar*x*ar'-er*x*er'+brw*brw')/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time u = plyapd(ac,ec,bc);
x = u*u'; @test norm(ac*x*ac'-ec*x*ec'+bc*bc')/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

@time u = plyapd(ac,ec,bcw);
x = u*u'; @test norm(ac*x*ac'-ec*x*ec'+bcw*bcw')/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

@time u = plyapd(ac',ec',cc');
x = u'*u; @test norm(ac'*x*ac-ec'*x*ec+cc'*cc)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

@time u = plyapd(ac',ec',cct');
x = u'*u; @test norm(ac'*x*ac-ec'*x*ec+cct'*cct)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

@time u = plyapd(ar,er,bc);
x = u*u'; @test norm(ar*x*ar'-er*x*er'+bc*bc')/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time u = plyapd(ar',er',cc');
x = u'*u; @test norm(ar'*x*ar-er'*x*er+cc'*cc)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time u = plyapd(ac,ec,br);
x = u*u'; @test norm(ac*x*ac'-ec*x*ec'+br*br')/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time u = plyapd(ac',ec',cr');
x = u'*u; @test norm(ac'*x*ac-ec'*x*ec+cr'*cr)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol
end
end

@testset "Positive discrete 2x2 Lyapunov equations" begin

A = [-1.1 1.; -1. -1.]/10
E = [1. 1.; 0. 1.]
R = [1. 1.; 0. 1.]
reltol = eps(float(100))

U, scale, β, α = plyap2(A, R, adj = true, disc = true)
X = U'*U; @test norm(A'*X*A-X+scale^2*R'*R)/max(1,norm(X))/norm(A) < reltol &&
                norm(β*U-U*A)/max(1,norm(U))/norm(A) < reltol &&
                norm(α*U/scale - scale*R)/max(1,norm(scale*R)) < reltol

U, scale, β, α = plyap2(A, R, adj = false, disc = true)
X = U*U'; @test norm(A*X*A'-X+scale^2*R*R')/max(1,norm(X))/norm(A) < reltol &&
                norm(U*β-A*U)/max(1,norm(U))/norm(A) < reltol &&
                norm(U*α/scale - scale*R)/max(1,norm(scale*R)) < reltol

U, scale, β, α = pglyap2(A, E, R, adj = true, disc = true)
X = U'*U; @test norm(A'*X*A-E'*X*E+scale^2*R'*R)/max(1,norm(X))/norm(A) < reltol &&
                norm(β*U*E-U*A)/max(1,norm(U))/norm(A) < reltol &&
                norm(α*U*E/scale - scale*R)/max(1,norm(scale*R)) < reltol

U, scale, β, α = pglyap2(A, E, R, adj = false, disc = true)
X = U*U'; @test norm(A*X*A'-E*X*E'+scale^2*R*R')/max(1,norm(X))/norm(A) < reltol &&
                norm(A*U-E*U*β)/max(1,norm(U))/norm(A) < reltol &&
                norm(E*U*α/scale - scale*R)/max(1,norm(scale*R)) < reltol

A = [-1.1 1.; -1. -1.]/10; A = convert(Matrix{Float32},A)
E = [1. 1.; 0. 1.]; E = convert(Matrix{Float32},E)
R = [1. 1.; 0. 1.]; R = convert(Matrix{Float32},R)
reltol = eps(100f0)

U, scale, β, α = plyap2(A, R, adj = true, disc = true)
X = U'*U; @test norm(A'*X*A-X+scale^2*R'*R)/max(1,norm(X))/norm(A) < reltol &&
                norm(β*U-U*A)/max(1,norm(U))/norm(A) < reltol &&
                norm(α*U/scale - scale*R)/max(1,norm(scale*R)) < reltol

U, scale, β, α = plyap2(A, R, adj = false, disc = true)
X = U*U'; @test norm(A*X*A'-X+scale^2*R*R')/max(1,norm(X))/norm(A) < reltol &&
                norm(U*β-A*U)/max(1,norm(U))/norm(A) < reltol &&
                norm(U*α/scale - scale*R)/max(1,norm(scale*R)) < reltol

U, scale, β, α = pglyap2(A, E, R, adj = true, disc = true)
X = U'*U; @test norm(A'*X*A-E'*X*E+scale^2*R'*R)/max(1,norm(X))/norm(A) < reltol &&
                norm(β*U*E-U*A)/max(1,norm(U))/norm(A) < reltol &&
                norm(α*U*E/scale - scale*R)/max(1,norm(scale*R)) < reltol

U, scale, β, α = pglyap2(A, E, R, adj = false, disc = true)
X = U*U'; @test norm(A*X*A'-E*X*E'+scale^2*R*R')/max(1,norm(X))/norm(A) < reltol &&
                norm(A*U-E*U*β)/max(1,norm(U))/norm(A) < reltol &&
                norm(E*U*α/scale - scale*R)/max(1,norm(scale*R)) < reltol

end

@testset "Positive discrete Lyapunov equations - Schur form" begin

for Ty in (Float64, Float32)

ar = rand(Ty,n,n)
ar = ar/(one(Ty) + norm(ar))
as,  = schur(ar)
ac = rand(Ty,n,n)+im*rand(Ty,n,n)
ac = ac/(one(Ty) + norm(ac))
acs,  = schur(ac)
br = rand(Ty,n,m)
bc = br+im*rand(Ty,n,m)
cr = rand(Ty,p,n)
cc = cr+im*rand(Ty,p,n)
Ty == Float64 ? reltol = eps(float(100)) : reltol = eps(100*n*one(Ty))

@time u = plyaps(as,br,disc = true);
x = u*u'; @test norm(as*x*as'-x+br*br')/norm(x)/norm(as) < reltol

@time u = plyaps(as,ar,disc = true);
x = u*u'; @test norm(as*x*as'-x+ar*ar')/norm(x)/norm(as) < reltol

@time u = plyaps(as',cr',disc = true);
x = u'*u; @test norm(as'*x*as-x+cr'*cr)/norm(x)/norm(as) < reltol

@time u = plyaps(as',ar',disc = true);
x = u'*u; @test norm(as'*x*as-x+ar'*ar)/norm(x)/norm(as) < reltol

@time u = plyaps(acs,bc,disc = true);
x = u*u'; @test norm(acs*x*acs'-x+bc*bc')/norm(x)/norm(acs) < reltol

@time u = plyaps(acs',cc',disc = true);
x = u'*u; @test norm(acs'*x*acs-x+cc'*cc)/norm(x)/norm(as) < reltol

R = UpperTriangular(rand(Ty,n,n))
U = copy(R)
@time plyapds!(as, U, adj = false)
X = U*U'; @test norm(as*X*as'-X+R*R')/max(1,norm(X))/norm(as) < reltol

R = UpperTriangular(rand(Ty,n,n))
U = copy(R)
@time plyapds!(as,U,adj = true)
X = U'*U; @test norm(as'*X*as-X+R'*R)/max(1,norm(X))/norm(as) < reltol

R = UpperTriangular(rand(Complex{Ty},n,n))
U = copy(R)
@time plyapds!(acs,U,adj = false)
X = U*U'; @test norm(acs*X*acs'-X+R*R')/max(1,norm(X))/norm(acs) < reltol

R = UpperTriangular(complex(rand(Ty,n,n)))
U = copy(R)
@time plyapds!(acs,U,adj = true)
X = U'*U; @test norm(acs'*X*acs-X+R'*R)/max(1,norm(X))/norm(acs) < reltol

er = rand(Ty,n,n)
ar = er*ar
as, es = schur(ar,er)
ec = er+im*rand(Ty,n,n)
ac = ec*ac
acs, ecs = schur(ac,ec)

@time u = plyaps(as,es,br,disc = true);
x = u*u'; @test norm(as*x*as'-es*x*es'+br*br')/norm(x)/norm(as) < reltol

@time u = plyaps(as,es,ar,disc = true);
x = u*u'; @test norm(as*x*as'-es*x*es'+ar*ar')/norm(x)/norm(as) < reltol

@time u = plyaps(as',es',cr',disc = true);
x = u'*u; @test norm(as'*x*as-es'*x*es+cr'*cr)/norm(x)/norm(as) < reltol

@time u = plyaps(as',es',ar',disc = true);
x = u'*u; @test norm(as'*x*as-es'*x*es+ar'*ar)/norm(x)/norm(as) < reltol

@time u = plyaps(acs,ecs,bc,disc = true);
x = u*u'; @test norm(acs*x*acs'-ecs*x*ecs'+bc*bc')/norm(x)/norm(as) < reltol

@time u = plyaps(acs',ecs',cc',disc = true);
x = u'*u; @test norm(acs'*x*acs-ecs'*x*ecs+cc'*cc)/norm(x)/norm(as) < reltol

R = UpperTriangular(rand(Ty,n,n))
U = copy(R)
@time plyapds!(as,es,U,adj = false)
X = U*U'; @test norm(as*X*as'-es*X*es'+R*R')/max(1,norm(X))/norm(as) < reltol

R = UpperTriangular(rand(Ty,n,n))
U = copy(R)
@time plyapds!(as,es,U,adj = true)
X = U'*U; @test norm(as'*X*as-es'*X*es+R'*R)/max(1,norm(X))/norm(as) < reltol

R = UpperTriangular(rand(Complex{Ty},n,n)) # error
U = copy(R)
@time plyapds!(acs,ecs,U,adj = false)
X = U*U'; @test norm(acs*X*acs'-ecs*X*ecs'+R*R')/max(1,norm(X))/norm(acs) < reltol

R = UpperTriangular(rand(Complex{Ty},n,n)) # error
U = copy(R)
@time plyapds!(acs,ecs,U,adj = true)
X = U'*U; @test norm(acs'*X*acs-ecs'*X*ecs+R'*R)/max(1,norm(X))/norm(acs) < reltol
end
end
end

end


text:
"""
        Helper function for Smac reading of binary images
"""

"""
    struct Smac1ImageInfo

Stores the information in a Smac binary image header.
"""
struct Smac1ImageInfo

    snap::Int32                 # number of input snapshot
    z::Float32                  # redshift of snapshot
    m_vir::Float32              # virial mass of halo
    r_vir::Float32              # virial radius of halo
    xcm::Float32                # x coordinate of image center
    ycm::Float32                # y coordinate of image center
    zcm::Float32                # z coordinate of image center
    z_slice_kpc::Float32        # depth of the image in kpc
    boxsize_kpc::Float32        # xy-size of the image in kpc
    boxsize_pix::Float32        # xy-size of the image in pixels
    pixsize_kpc::Float32        # size of one pixel in kpc
    xlim::Array{Float64,1}      # x limits of image
    ylim::Array{Float64,1}      # y limits of image
    zlim::Array{Float64,1}      # z limits of image
    units::String               # unitstring of image

    function Smac1ImageInfo(snap::Int32, z::Float32, m_vir::Float32, r_vir::Float32,
                       xcm::Float32, ycm::Float32, zcm::Float32,
                       z_slice_kpc::Float32,
                       boxsize_kpc::Float32, boxsize_pix::Float32, pixsize_kpc::Float32,
                       units::String)

        xlim = [xcm - boxsize_kpc/2.0, xcm + boxsize_kpc/2.0]
        ylim = [ycm - boxsize_kpc/2.0, ycm + boxsize_kpc/2.0]
        zlim = [zcm - z_slice_kpc/2.0, zcm + z_slice_kpc/2.0]

        new(snap, z, m_vir, r_vir,
            xcm, ycm, zcm,
            z_slice_kpc,
            boxsize_kpc, boxsize_pix, pixsize_kpc,
            xlim, ylim, zlim, units)
    end

end

"""
    read_smac1_binary_image(filename::String)

Reads a binary image file from Smac and returns a Matrix with the pixel values.
"""
function read_smac1_binary_image(filename::String)
    f = open(filename)
    first = read(f, Int32)
    n_pixels = Int64(sqrt.(first/4.0))
    a = read!(f, Array{Float32,2}(undef,n_pixels,n_pixels))
    last = read(f, Int32)
    close(f)
    if first == last
        return collect(transpose(a))
    else
        error("Read error: Incorrect image format!")
    end
end

"""
    read_smac1_binary_info(filename::String)

Returns the image info in a `Smac1ImageInfo` struct.
"""
function read_smac1_binary_info(filename::String)

    f = open(filename)

    # skip image
    image = read(f, Int32)
        seek(f, position(f) + image + 8)

    # read info block
    snap = read(f, Int32)
        seek(f, position(f) + 8)
    z = read(f, Float32)
        seek(f, position(f) + 8)
    m_vir = read(f, Float32)
        seek(f, position(f) + 8)
    r_vir = read(f, Float32)
        seek(f, position(f) + 8)
    xcm = read(f, Float32)
        seek(f, position(f) + 8)
    ycm = read(f, Float32)
        seek(f, position(f) + 8)
    zcm = read(f, Float32)
        seek(f, position(f) + 8)
    z_slice_kpc = read(f, Float32)
        seek(f, position(f) + 8)
    boxsize_kpc = read(f, Float32)
        seek(f, position(f) + 8)
    boxsize_pix = read(f, Float32)
        seek(f, position(f) + 8)
    pixsize_kpc = read(f, Float32)
        seek(f, position(f) + 4)
    # read length of unit string
    unitstring_length = read(f, Int32)
    # read unti string
    unit = String(Char.((read!(f, Array{Int8,1}(undef,unitstring_length)))))
    close(f)

    return Smac1ImageInfo(snap, z, m_vir, r_vir,
                    xcm, ycm, zcm,
                    z_slice_kpc,
                    boxsize_kpc, boxsize_pix, pixsize_kpc,
                    unit)
end

"""
    write_smac1_par([...])

Writes a Smac parameter file. Not all relevant parameters implemented yet!
"""
function write_smac1_par(path, use_keys, out_dir, snap_base, snap_file, image_prefix,
                         halo_id, snap_start, snap_end, file_format, kernel,
                         xy_size, z_size, out_effect, out_subeffect,
                         cr_bins, cr_pmin, cr_pmax, cr_subsamples,
                         projection, x0, y0, z0)

    open(path * "smac1.par", "w") do f
        write(f,
        """#=============================================================
        # Input file
        #=============================================================

        SIM_FORMAT = 4

        #**** Input file
        USE_KEYS = $use_keys
        OUTPUT_DIR = $out_dir
        SNAP_BASE = $snap_base

        SNAP_FILE = $snap_file

        #**** Prefix of the output files:
        PREFIX_OUT = $image_prefix

        #**** Halo ID
        HALO_ID = $halo_id

        SNAP_START = $snap_start
        SNAP_END   = $snap_end

        INITIAL_Z = 0.0
        IMAGE_Z   = 0.0

        FILE_FORMAT = $file_format

        FILE_HEADER = 1

        KERNEL_TYPE = $kernel

        PART_DISTR = 1

        IMG_XY_UNITS = 2
        IMG_XY_SIZE  = $xy_size

        IMG_Z_UNITS = 2
        IMG_Z_SIZE  = $z_size
        FILL_CARROT = 0.

        IMG_SIZE = 1024

        SMOOTH_UNITS = 0
        SMOOTH_FWHM  = 5.

        NSIDE = 8

        MIN_DIST =   5000.
        MAX_DIST = 420000.

        OUTPUT_MAP = $out_effect
        OUTPUT_SUB = $out_subeffect

        METAL_SPECIES = 0

        TEMP_UNIT = 1

        TEMP_CUT = 0.5

        XRAY_E0 = 0.5
        XRAY_E1 = 2.

        X_TABLES =

        FRQ_ZS = 300.

        #**** Input for Radio Maps
        #XCRP_TABLE_PATH gives the change of XCRP over radius.
        ECRP_MIN = 1e9
        FRQ_Pnu = 1.4e9
        XCRP = 0.01
        XCRP_TABLE_PATH = ~/
        B_TABLE_PATH = ~/
        KERNEL_TABLE_PATH =
        GAM_nu = 1.25
        #XCRP_TABLE_FILE_FMT = (1F8.6,2X,1F14.6)

        #energies in GeV
        E_gam = 1e1
        E_gam_min = 1e-2
        E_gam_max = 30e3

        # Number of momentum bins for CR model
        CR_nbins = $cr_bins
        CR_pmin = $cr_pmin
        CR_pmax = $cr_pmax
        CR_subsamples = $cr_subsamples
        CR_DSlope = 1.0e-6

        GIVE_MORE_INFO = 0

        PROJECT       = $projection

        CENTER_MOTION = 0

        REMOVE_LOCAL_GROUP_VEL = 0

        CENTER = 1

        CENTER_X = $x0
        CENTER_Y = $y0
        CENTER_Z = $z0

        PERIODIC = 1

        LIGHTCONE = 0
        X_ORIGIN = 0.
        Y_ORIGIN = 0.
        Z_ORIGIN = 0.
        OPEN_ANGLE = 1.0

        MAIN_PROG = /afs/mpa/project/hydrosims/Hutt/g676/csf/Post/main_prog.a_gas.gv

        HUBBLE = 0.7
        OMEGA  = 0.3
        LAMBDA = 0.7

        """
        )
    end
end

ext: jl

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
function plot_states(ts, xs, nx; figure_name = "states", clear_figure = true, save_figure = true, fig_name = "states.png", color = "blue", stl = "-", label = "")
    fig_states = PP.figure(figure_name)
    if clear_figure; fig_states.clear(); end

    fig_states.add_subplot(231, ylabel = LS.L"x[m]")
    fig_states.add_subplot(234, ylabel = LS.L"y[m]", xlabel = LS.L"t[s]")
    fig_states.add_subplot(232, ylabel = LS.L"v_x[m/s]")
    fig_states.add_subplot(235, ylabel = LS.L"v_y[m/s]", xlabel = LS.L"t[s]")
    fig_states.add_subplot(233, ylabel = LS.L"\theta[^\circ]")
    fig_states.add_subplot(236, ylabel = LS.L"\dot{\theta}[^\circ/s]", xlabel = LS.L"t[s]")

    axs = fig_states.get_axes()
    for (i,ax) in enumerate(axs)
        ax.plot(ts, index_states(xs, nx, i), color = color, stl, label = label)
    end

    fig_states.tight_layout()
    fig_states.set_size_inches(11, 7, forward=true)

    if save_figure
        if !isdir(joinpath(topdirectory, "data", "plots", TODAY)); mkpath(joinpath(topdirectory, "data", "plots", TODAY)); end
        fig_states.savefig(joinpath(topdirectory, "data", "plots", TODAY, fig_name))
    end
end

max_stars_repo_name: byuflowlab/MultifidelityTrajectories

max_issues_repo_name: byuflowlab/MultifidelityTrajectories

max_forks_repo_name: byuflowlab/MultifidelityTrajectories


text:
using Test

import Base.Order: Forward, Reverse

# every argument can flip the integer order by passing the right value. Here,
# we enumerate a few of these combinations and check that all these flips
# compound so that in total we either have an increasing or decreasing sort.
for (s1, rev) in enumerate([true, false])
    for (s2, lt) in enumerate([>, <, (a, b) -> a - b > 0, (a, b) -> a - b < 0])
        for (s3, by) in enumerate([-, +])
            for (s4, order) in enumerate([Reverse, Forward])
                if iseven(s1 + s2 + s3 + s4)
                    target = [1, 2, 3]
                else
                    target = [3, 2, 1]
                end
                @test target == sort([2, 3, 1], rev=rev, lt=lt, by=by, order=order)
            end
        end
    end
end

@test [1 => 3, 2 => 5, 3 => 1] ==
            sort([1 => 3, 2 => 5, 3 => 1]) ==
            sort([1 => 3, 2 => 5, 3 => 1], by=first) ==
            sort([1 => 3, 2 => 5, 3 => 1], rev=true, order=Reverse) ==
            sort([1 => 3, 2 => 5, 3 => 1], lt= >, order=Reverse)

@test [3 => 1, 1 => 3, 2 => 5] ==
            sort([1 => 3, 2 => 5, 3 => 1], by=last) ==
            sort([1 => 3, 2 => 5, 3 => 1], by=last, rev=true, order=Reverse) ==
            sort([1 => 3, 2 => 5, 3 => 1], by=last, lt= >, order=Reverse)

struct SomeOtherOrder <: Base.Order.Ordering end

@test_throws ErrorException sort([1, 2, 3], lt=(a, b) -> a - b < 0, order=SomeOtherOrder())

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
#==#
#Ratios.BisectSort()
@propagate_inbounds function div(n=10,m=6)   #inline's maximum is slower its counterparty by factor of 2.6928
   k=m
    @inbounds for i ∈ 1:k
    n /= m
    end
end
#=
maybe the point of inline is not the speed at all, but
only the stability of your code? #research
=#
#inline is not the best  for speed
  @inline function div2(n=10,m=6)    # the faster # depreciate
   k=m
    @inbounds for i ∈ 1:k
    n /= m
    end
end

@propagate_inbounds function mul(n=10,m=6)
    k = m
   @inbounds for i ∈ 1:k
        n *= m
    end
end
#-----------------------
t1 = @benchmark div()
t2 = @benchmark div2()

maxt1 = maximum(t1)
maxt2   = maximum(t2)

maxt1.time/mint1.time #
mint1.time/maxt1.time

maxratio = maxt1.time / maxt2.time #max time of div1  is larger  than max time of div2 by factor of 2.6928
minratio = mint1.time/mint2.time  # minimum is equal for div1 & div2
#Idea: min-max matrix
#vcat([maxt1/maxt2,maxt1/mint2])

#meanDiv2ns = median(t2)
#modeDiv2ns = mode(t2)
#stdDiv2ns = std(t2)

mint1.time = minimum(t1)

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

avg = 0
tmul = @benchmark mul()#10^6

maxi = 100
count=0
tmul = @benchmark mul()#

#---------------------------------------------------
@inbounds for i ∈ 1:maxi #is 10_000
    tmul = @benchmark mul()#
    tmul
    maxm = maximum(tmul.times)
    count += maxm
end
count #count = 55975
avg = count/maxi #population MostLikelihoodEstimae =  E(x)=avg =55.975
sampleMax=maximum(tmul)
#tmul.median
sampletime = sampleMax.time #143.3 # is it higher or lower of the 'Normal'?
sampletime/avg #2.56 Infer: sample is higher by factor of 2.56 # standard deviant got to be high as well
dispersion = sampletime- avg# 87.32499
dispersion/sampletime/avg #idea 0.01 (pvalue significant)

mint2 = minimum(t2)

mint1.time/mint2.time
mint1.time ≈ mint2.time
mint1.time == mint2.time

#=
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):  1.400 ns … 143.300 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.500 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.637 ns ±   2.125 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ██ ▃  ▂ ▁                                                   ▂
  ██▁█▁██▁█▁▇▁▇▄▁▅▁▄▁▄▃▁▁▁▁▁▃▁▁▁▁▄▄▁▁▁▃▁▅▃▁▄▁▁▁▃▃▁▁▁▃▁▄▃▁▄▁▁▃ █
  1.4 ns       Histogram: log(frequency) by time       4.8 ns <

julia> dump(tmul)

  params: BenchmarkTools.Parameters
.....
     seconds: Float64 5.0
    time_tolerance: Float64 0.05
    memory_tolerance: Float64 0.01
....
julia>     gcsample: Bool falseBenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):  1.400 ns … 143.300 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.500 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.637 ns ±   2.125 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ██ ▃  ▂ ▁                                                   ▂
  ██▁█▁██▁█▁▇▁▇▄▁▅▁▄▁▄▃▁▁▁▁▁▃▁▁▁▁▄▄▁▁▁▃▁▅▃▁▄▁▁▁▃▃▁▁▁▃▁▄▃▁▄▁▁▃ █
  1.4 ns       Histogram: log(frequency) by time       4.8 ns <
vial info:
    seconds: Float64 5.0

    time_tolerance: Float64 0.05
    memory_tolerance: Float64 0.01

    -----------
Lesson learned:population (of the poulation) : returned time could be dubious, depending on the distribution of variance (which is nasty indee d)
#=stil, giving emphasis on just the mean, could mean turning a blind eye on possible meaningful anomalies
most data is ranged differently according to its interval (interval analysis ?)
    a fast fourier  would enhance analysis but the best approach would be a "weighted kernel approach"
that too may not possible even explain the existence of such anomalies, in general-  not only in this one
mint2 = minimum(t2)

mint1.time/mint2.time
mint1.time ≈ mint2.time
mint1.time == mint2.time

------------------

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
"""
    Kernel(perturbation_function::Function,
    pdf_function::Function,
    calculate_kernel_parameters::Function)

Create a parameter perturbation kernel. Required inputs are 3 functions.
First is the `perturbation_function` which should take 2 parameters,
the parameter to be perturbed and any kernel specific parameter
(for example the standard deviation of a normal distribution if this is the kernel of choice).
Second function is the `pdf_function`, that requires 4 inputs:
    1) the newparticle
    2) the old particle
    3) kernel specific parameters and
    4) an index i.
The third function is `calculate_kernel_parameters` which given an array of
particles should calculate the kernel specific parameters for the next population.
Should you wish to keep the same parameters throughout you can just write a function that returns a number(s).
"""
mutable struct Kernel
    perturbation_function::Function
    pdf_function::Function
    calculate_kernel_parameters::Function
    kernel_parameters

    Kernel(perturbation_function::Function,
    pdf_function::Function,
    calculate_kernel_parameters::Function) =
    new(perturbation_function,
    pdf_function,
    calculate_kernel_parameters,
    [])
end

perturbation_function(prevparameter, scale) = rand(Uniform(prevparameter - scale, prevparameter + scale))

pdf_function(newparticle, prevparticle, scales, i) = pdf(Uniform(prevparticle.params[i] - scales[i], prevparticle.params[i] + scales[i]), newparticle.params[i])
function calculate_kernel_parameters(particles)
  #calculate the range of parameter values (ie the scale) to use for the
  #perturbation kernel
  parameters = hcat(map(x -> x.params, particles)...)'
  scales = ((maximum(parameters, dims = 1) -
                  minimum(parameters, dims = 1)) ./ 2.0)[:]
  return scales
end

uniformkernel = Kernel(
    perturbation_function,
    pdf_function,
    calculate_kernel_parameters
)

perturbation_function_gauss(currentparameter, stdev) = rand(Normal(currentparameter, stdev))
pdf_function_gauss(newparticle, prevparticle, stdev, i) = pdf(Normal(prevparticle.params[i], stdev[i]), newparticle.params[i])
function calculate_kernel_parameters_gauss(particles)
  #calculate the standard deviation of previous population parameters,
  # returns an array for the standard deviation of each parameter
  std(hcat(map(x -> x.params, particles)...), dims = 2)
end

gaussiankernel = Kernel(
    perturbation_function_gauss,
    pdf_function_gauss,
    calculate_kernel_parameters_gauss
)

function modelselection_kernel(ABCsetup, particles)
    #calculate the parameters for the perturbation kernel for each model

    #get the indeces of particles that are a member of each model
    modelindex = trues(ABCsetup.nparticles, ABCsetup.nmodels)
    for i in 1:ABCsetup.nmodels
        modelindex[:, i] = map(x -> x.model, particles) .== i
    end

    #calculate the frequency of each model
    modelfreq = sum(modelindex, dims = 1)

    #calculate the kernel parameters for each model
    for i in 1:ABCsetup.nmodels
      if modelfreq[i] == 0
        ABCsetup.Models[i].kernel.kernel_parameters = [0.0]
      elseif modelfreq[i] == 1
        ABCsetup.Models[i].kernel.kernel_parameters = ABCsetup.Models[i].kernel.kernel_parameters
      else
        ABCsetup.Models[i].kernel.kernel_parameters =
        ABCsetup.Models[i].kernel.calculate_kernel_parameters(particles[modelindex[:, i]])
      end
    end
    return ABCsetup
end


text:
# mobility.jl

"""
----------------------------------------------------------------------
The Moblity of the Polaron.
----------------------------------------------------------------------
"""

"""
polaron_mobility_ac(Ω::Float64, β::Float64, α::Float64, v::Float64, w::Float64)

    Calculate the ac moblity μ(Ω) of the polaron at finite temperature (equation (1) in Hellwarth 1999) for a given frequency Ω. β is the thermodynamic beta. v and w are the variational polaron parameters that minimise the free energy, for the supplied α Frohlich coupling.
"""
function polaron_mobility_ac(Ω, β, α, v, w)
    Ω / imag(χ(Ω, β, α, v, w))
end

"""
polaron_mobility_ac(Ω::Float64, α::Float64, v::Float64, w::Float64)

    Calculate the ac moblity μ(Ω) of the polaron at zero-temperatures (equation (46) in FHIP 1962) for a given frequency Ω. v and w are the variational polaron parameters that minimise the free energy, for the supplied α Frohlich coupling.
"""
function polaron_mobility_ac(Ω, α, v, w)
    Ω / imag(χ(Ω, α, v, w))
end

"""
polaron_mobility_dc(β::Float64, α::Float64, v::Float64, w::Float64)

    Calculate the dc moblity μ(β) of the polaron at finite temperatures (equation (46) in FHIP 1962) at zero frequenc. v and w are the variational polaron parameters that minimise the free energy, for the supplied α Frohlich coupling.
"""
function polaron_mobility_dc(β, α, v, w)
    1 / imag(χ_dc(β, α, v, w))
end


text:
module ArithmeticKoans

function sum(x, y)
    return x + y
end

function difference(x, y)
    return x - y
end

function multiplication(x, y)
    return x * y
end

function division(x, y)
    return x / y
end

function integer_division(x, y)
    return div(x, y)
end

function power(x, y)
    return x^y
end

function modulo(x, y)
    return x % y
end

end


text:
module MasslistFunctions
	using HDF5
	using DelimitedFiles

	export IntegrationBorders, masslistPos, createMassList, loadMasslist, createCompound, massFromComposition, massFromCompositionArray, massFromCompositionArrayList, isotopesFromComposition, isotopesFromCompositionArray, sumFormulaStringFromCompositionArray, sumFormulaStringListFromCompositionArrayList, filterMassListByContribution1, filterMassListByContribution2, findClosestMassIndex, inCompositions, findInCompositions

	struct IntegrationBorders
	    centerMass::Array{Float64}
	    lowMass::Array{Float64}
	    highMass::Array{Float64}

	    function IntegrationBorders(masslist::Array{Float64,1}, resolution)
		borders = new(similar(masslist), similar(masslist),similar(masslist))
		if length(masslist) < 2
		    # return just one
		    borders.lowMass[1] = masslist[1]-masslist[1]/resolution
		    borders.highMass[1] = masslist[1]+masslist[1]/resolution
		    borders.centerMass[1] = masslist[1]
		    return borders
		end
		#set center masses
		#1st and last don't have two neighbors
		borders.lowMass[1] = masslist[1]-(masslist[1]/resolution)
		borders.highMass[end] = masslist[end]+(masslist[end]/resolution)
		borders.centerMass[1] = masslist[1]
		borders.centerMass[end] = masslist[end]

		for i=1:length(masslist)-1 # working on gaps between masses, not on masses
		    borders.centerMass[i] = masslist[i]
		    if (2*masslist[i]/resolution) < (masslist[i+1]-masslist[i]) # distance to next mass is bigger than 2*resolution --> use resolution border
		        borders.highMass[i] = masslist[i]+(masslist[i]/resolution)
		        borders.lowMass[i+1] = masslist[i+1]-(masslist[i]/resolution)
		    else # use half the distance as border
		        center = (masslist[i] + masslist[i+1])/2
		        borders.highMass[i] = center
		        borders.lowMass[i+1] = center
		    end
		end
		return borders
	    end
	    function IntegrationBorders(masslist::Array{Float64,1})
		return IntegrationBorders(masslist, 3000.0)
	    end
	end

	massC=12
	massC13=13.00335
	massH=1.00783
	massHplus=1.007276
	massN=14.00307
	massO=15.99492
	massO18=17.99916
	massS=31.97207
	nC13=0
	nO18=0

	abundanceC13 = 0.010816
	abundanceO18 = 0.002005

	function masslistPos(x)
	  if (x>0)
	    return x
	  else
	    return 0
	  end
	end

	function inCompositions(composition, compositionList)
	  for i=1:size(compositionList,2)
	    if composition == compositionList[:,i]
	      return true
	    end
	  end
	  return false
	end

	masslistElements = ["C", "C(13)", "H", "H+", "N", "O", "O(18)", "S"]
	masslistElementMasses = [massC, massC13, massH, massHplus, massN, massO, massO18, massS]
	function createMassList(; C=0:1:0, O=0:0, N=0:0, S=0:0, nHplus=1, allowRadicals=false)
	  masses = [37.02]
	  masslistCompositions = [[0 0 4 1 0 2 0 0]]
	  for nC in C
	    for nN in N
	      for nS in S
		for nO in O[O.<=nC+1]
		  if allowRadicals
		    possibleH = MasslistFunctions.masslistPos(round(nC/2)*2-4):1:2*nC+2+nN*3
		  else
		    possibleH = MasslistFunctions.masslistPos(round(nC/2)*2-4):2:2*nC+2+nN*3
		  end
		  if (isodd(nN))
		    possibleH = possibleH .+ 1
		  end
		  for nH in possibleH
		    append!(masses,nS*massS+nC*massC+nN*massN+nO*massO+nH*massH+nHplus*massHplus)
		    push!(masslistCompositions, [nC nC13 nH nHplus nN nO nO18 nS])
		  end
		end
	      end
	    end
	  end
	  sortIndices = sortperm(masses)
	  return masses[sortIndices],masslistElements,masslistElementMasses, masslistCompositions[sortIndices]
	end

	function createMassList2(; C=0:1:1, O=0:1:1, H=0:2:10, N=0:1, S=0:1, nHplus=1)
	  masses = [37.02]
	  masslistCompositions = [[0 0 4 1 0 2 0 0]]
	  for nC in C
	    for nN in N
	      for nS in S
		for nO in O[O.<=nC+1]
		  for nH in H
		    append!(masses,nS*massS+nC*massC+nN*massN+nO*massO+nH*massH+nHplus*massHplus)
		    push!(masslistCompositions, [nC nC13 nH nHplus nN nO nO18 nS])
		  end
		end
	      end
	    end
	  end
	  sortIndices = sortperm(masses)
	  return masses[sortIndices],masslistElements,masslistElementMasses, masslistCompositions[sortIndices]
	end

	function loadMasslist(filename)
	    fileIsValidResultfile = (endswith(filename, ".hdf5") || endswith(filename, ".h5"))
	    if fileIsValidResultfile
		println("Trying to load masslist from result hdf5 file ($filename)!")
		masses =  HDF5.h5read(filename, "MassList")
		masslistElementsLoaded = HDF5.h5read(filename,"ElementNames")
		compRaw = HDF5.h5read(filename, "ElementalCompositions")
		masslistCompositions = []
		for i=1:size(compRaw,2)
		    push!(masslistCompositions,compRaw[:,i])
		end
		return masses,masslistElementsLoaded, masslistElementMasses, masslistCompositions
	    else
	      println("Trying to load masslist from csv file ($filename)!")
	      list = readdlm(filename, '\t', skipstart=7)
	      masses = Array{Float64,1}()
	      masslistCompositions = []
	      print("Unidentified Peaks found in masslist $filename: ")
	      for i=1:size(list,1)
		if sum(list[i,1:8]) > 0
		  push!(masses,massFromCompositionArray(list[i,1:8]))
		  push!(masslistCompositions, list[i,1:8])
		else
		    if list[i,9] > 0
		    print("$(list[i,9]) ")
		    push!(masses,list[i,9])
		    push!(masslistCompositions, list[i,1:8])
		    end
		end
	      end
	      sortIndices = sortperm(masses)
	      return masses[sortIndices],masslistElements, masslistElementMasses, masslistCompositions[sortIndices]
	  end
	end

	function createCompound(; C=0, C13=0, H=0, Hplus=1, N=0, O=0, O18=0, S=0)
	  mass = C*massC + O*massO + N*massN +H*massH + S*massS + Hplus*massHplus + C13*massC13 + O18*massO18
	  composition = [C, C13, H, Hplus, N, O, O18, S]
	  return mass, masslistElements, masslistElementMasses, composition
	end

	function massFromComposition(; C=0, C13=0, H=0, Hplus=1, N=0, O=0, O18=0, S=0)
	  mass = C*massC + O*massO + N*massN +H*massH + S*massS + Hplus*massHplus + C13*massC13 + O18*massO18
	  #composition = [C C13 H Hplus N O O18 S]
	  return mass
	end

	function isotopesFromComposition(; C=0, C13=0, H=0, Hplus=1, N=0, O=0, O18=0, S=0)
	  masses = Array{Float64,1}()
	  compositions = []
	  abundances = Array{Float64,1}()
	  for carbIsotopeNr = 0:minimum([C,2])
	    mass = (C-carbIsotopeNr)*massC + O*massO + N*massN +H*massH + S*massS + Hplus*massHplus + (C13+carbIsotopeNr)*massC13 + O18*massO18
	    abundance = (abundanceC13^carbIsotopeNr)*binomial(Int64(C),Int64(carbIsotopeNr))
	    composition = [C-carbIsotopeNr C13+carbIsotopeNr H Hplus N O O18 S]
	    push!(masses,mass)
	    push!(compositions, composition)
	    push!(abundances, abundance)
	  end
	  if (length(masses) == 1 && masses[1] == 0)
	    println("Strange composition found, maybe unidentified compound: $(compositions[1])")
	  end
	  return masses, masslistElements, compositions, abundances
	end

	function isotopesFromCompositionArray(composition)
	  #[C C13 H Hplus N O O18 S]
	  return isotopesFromComposition(C=composition[1],
	  C13=composition[2],
	  H=composition[3],
	  Hplus=composition[4],
	  N=composition[5],
	  O=composition[6],
	  O18=composition[7],
	  S=composition[8])
	end

	function massFromCompositionArray(composition)
	  #mass = composition[1]*massC + composition[2]*massC13 +composition[3]*massH + composition[4]*massHplus
	  #+ composition[5]*massN + composition[6]*massO + composition[7]*massO18 + composition[8]*massS
	  mass = sum(masslistElementMasses .* composition)
	  #composition = [C C13 H Hplus N O O18 S]
	  return mass
	end

	function massFromCompositionArrayList(compositions)
	  ret = Array{Float32,1}()
	  for i=1:size(compositions,2)
	    push!(ret,massFromCompositionArray(compositions[:,i]))
	  end
	  return ret
	end

	function sumFormulaStringFromCompositionArray(composition)
	  if (composition[1]>1)
	    CString = "C$(composition[1])"
	  elseif composition[1] == 1
	    CString = "C"
	  else
	    CString = ""
	  end

	  if (composition[3]>1)
	    HString = "H$(composition[3])"
	  elseif composition[3] == 1
	    HString = "H"
	  else
	    HString = ""
	  end

	  if (composition[6]>1)
	    OString = "O$(composition[6])"
	  elseif composition[6] == 1
	    OString = "O"
	  else
	    OString = ""
	  end

	  if (composition[5]>1)
	    NString = "N$(composition[5])"
	  elseif composition[5] == 1
	    NString = "N"
	  else
	    NString = ""
	  end

	  if (composition[8]>1)
	    SString = "S$(composition[8])"
	  elseif composition[8] == 1
	    SString = "S"
	  else
	    SString = ""
	  end

	  return "$CString$HString$OString$NString$SString"
	end

	function sumFormulaStringListFromCompositionArrayList(compositions; showMass = false)
	  ret = Array{String,1}()
	  for i=1:size(compositions,2)
	    if showMass
	      push!(ret,"$(round(massFromCompositionArray(compositions[:,i]),4)) - $(sumFormulaStringFromCompositionArray(compositions[:,i]))")
	    else
	      push!(ret,sumFormulaStringFromCompositionArray(compositions[:,i]))
	    end
	  end
	  return ret
	end

	function filterMassListByContribution1(masses, countrates, relThreshold)
	  selected = Array{Bool}(undef,length(masses))
	  for i=1:ceil(maximum(masses))
	    sel = (masses.>i-0.3) & (masses.<i+0.7)
	    quant = maximum(countrates[sel])*threshold
	    selected[sel] = countrates[sel].>quant
	  end
	  return selected
	end

	function filterMassListByContribution2(masses, countrates, resolution, relThreshold)
	  selected = Array{Bool}(undef,length(masses))
	  fill!(selected,true)

	  for i=2:length(masses)-1
	    if abs(masses[i]-masses[i+1]) < masses[i]/4*resolution # if neighbor is within resolution
	      th =  relThreshold * masses[i]/(resolution*abs(masses[i]-masses[i+1]))
	      print("\nMass $(masses[i]) close to neighbor, th=$th  ")
	      if countrates[i] < (countrates[i+1] * th) # deselect if contribution is lower than relThreshold*neighbor
		print("REMOVED")
		selected[i] = false
	      end
	    end
	    if abs(masses[i]-masses[i-1]) < masses[i]/4*resolution # if neighbor is within resolution
	      th = relThreshold * masses[i]/(resolution*abs(masses[i]-masses[i-1]))
	      print("\nMass $(masses[i]) close to neighbor, th=$th  ")
	      if countrates[i] < (countrates[i-1] * th) # deselect if contribution is lower than relThreshold*neighbor
		print("REMOVED")
		selected[i] = false
	      end
	    end
	    if countrates[i] < 0
	      selected[i] = false
	    end
	  end
	  return selected
	end

	function OScFromComposition(composition)
	  return (2*composition[6]-composition[3])/composition[1] # OSc = 2*O/C - (H/C)
	end
	function OScFromCompositionArray(compositions)
	  ret = Array{Float64,1}()
	  for i=1:size(compositions,2)
	      push!(ret,OScFromComposition(compositions[:,i]))
	  end
	end

	function findClosestMassIndex(mass, masses)
	    closestIndex::Int64 = 1
	    for i=1:length(masses)
		if abs(masses[i]-mass)<abs(masses[closestIndex]-mass)
		    closestIndex = i
		end
	    end
	    return closestIndex
	end

	function findInCompositions(compositions, composition)
	    for i=1:size(compositions,2)
		if(compositions[:,i] == composition)
		    return i
		end
	    end
	    return 0
	end

end


text:
function avgfilter1d_constprop_setup(n, w = 3)
    xs = rand(SVector{w,UInt8}, n)
    ys = zero(xs)
    return (; ys, xs)
end

@inline function avg_at(xs0, i, N)
    @static if USE_ALIASSCOPE
        xs = Const(xs0)
    else
        xs = xs0
    end
    s = zero(eltype(xs))
    for k in 1:N
        s = s .+ @inbounds xs[i + k - 1]
    end
    return s .÷ N
end

@inline function avgfilter1d_constprop_seq!(ys, xs, N = 4)
    @assert axes(ys) == axes(xs)
    for i0 in 0:length(xs) - N
        i = firstindex(xs) + i0
        @inbounds ys[i] = avg_at(xs, i, N)
    end
    return ys
end

@inline function avgfilter1d_constprop_threads!(ys, xs, N = 4)
    @assert axes(ys) == axes(xs)
    Threads.@threads for i0 in 0:length(xs) - N
        i = firstindex(xs) + i0
        @inbounds ys[i] = avg_at(xs, i, N)
    end
    return ys
end

# Using `@inline` to nudge that the 3-arg method sees the constant 4 set by the
# "outer" 2-arg function.
@inline function avgfilter1d_constprop_tapir_dac!(ys, xs, N = 4)
    @assert axes(ys) == axes(xs)
    Tapir.@par dac for i0 in 0:length(xs) - N
        i = firstindex(xs) + i0
        @inbounds ys[i] = avg_at(xs, i, N)
        @grainsize 131072
    end
    return ys
end

@inline function avgfilter1d_constprop_tapir_seq!(ys, xs, N = 4)
    @assert axes(ys) == axes(xs)
    Tapir.@par seq for i0 in 0:length(xs) - N
        i = firstindex(xs) + i0
        @inbounds ys[i] = avg_at(xs, i, N)
        @grainsize 131072
    end
    return ys
end


text:
@testset "Slow Manifold System" begin
  function slow_manifold(du, u, p, t)
    du[1] = p[1]*u[1]
    du[2] = p[2]*(u[2]-u[1]^2)
  end

  u0 = [3.0; -2.0]
  tspan = (0.0, 5.0)
  p = [-0.8; -0.7]

  problem = ODEProblem(slow_manifold, u0, tspan, p)
  solution = solve(problem, Tsit5(), saveat = 0.01)

  prob = ContinuousDataDrivenProblem(solution)

  @variables u[1:2]
  Ψ = Basis([u; u[1]^2], u)
  for alg in [DMDPINV(); DMDSVD(); TOTALDMD()]
    res = solve(prob, Ψ, alg, digits = 1)
    b = result(res)
    m = metrics(res)
    @test isapprox(eigvals(b), [2*p[1]; p[1]; p[2]], atol = 1e-1)
    @test all(m[:L₂] .< 1e-10)

    _prob = ODEProblem((args...)->b(args...), u0, tspan, parameters(res))
    _sol = solve(_prob, Tsit5(), saveat = solution.t)
    @test norm(solution - _sol)/size(solution, 2) < 1e-1
  end
end

@testset "Slow Manifold Discrete System" begin
  function slow_manifold(du, u, p, t)
    du[1] = p[1]*u[1]
    du[2] = p[2]*u[2]+(p[1]^2 - p[2])*u[1]^2
  end

  u0 = [3.0; -2.0]
  tspan = (0.0, 10.0)
  p = [-0.8; -0.7]

  problem = DiscreteProblem(slow_manifold, u0, tspan, p)
  solution = solve(problem, FunctionMap())

  ddprob = DiscreteDataDrivenProblem(solution)

  @variables u[1:2]

  Ψ = Basis([u[1]; u[1]^2; u[2]-u[1]^2], u)

  for alg in [DMDPINV(); DMDSVD(); TOTALDMD()]
    res = solve(ddprob, Ψ, alg, digits = 1)
    b = result(res)
    m = metrics(res)
    @test isapprox(eigvals(b), [p[1]; p[2]; p[1]^2], atol = 1e-1)
    @test all(m[:L₂] .< 3e-1)
  end
end


text:
# adapted from https://github.com/JuliaPlots/StatsPlots.jl/blob/master/src/ecdf.jl
# to avoid dependency and to make it a series plot instead of a user plot
# not exported to avoid clashing with StatsPlots.ecdfplot

"""
    ecdfplot(x; kwargs...)
    ecdfplot!(x; kwargs...)

Plot the empirical cumulative distribution function (ECDF) of the data `x`.

# Positional Arguments
- `x`: Values to plot

# Keyword Arguments
- `weights=nothing`: Vector of weights to be used for a weighted ECDF
- `kwargs`: Additional attributes understood by Plots.jl
"""
ecdfplot(x; kwargs...) = RecipesBase.plot(x; seriestype=:ecdfplot, kwargs...)
ecdfplot!(x; kwargs...) = RecipesBase.plot!(x; seriestype=:ecdfplot, kwargs...)

RecipesBase.@recipe function f(::Type{Val{:ecdfplot}}, x, y, z; weights=nothing)
    ecdf = weights === nothing ? StatsBase.ecdf(y) : StatsBase.ecdf(y; weights=weights)
    xnew, ynew = _ecdf_points(ecdf)
    if isrotated(plotattributes)
        xnew, ynew = ynew, xnew
    end
    seriestype := :steppost
    x := xnew
    y := ynew
    return ()
end

function _ecdf_points(ecdf::StatsBase.ECDF)
    x = [ecdf.sorted_values[1]; ecdf.sorted_values]
    if hasproperty(ecdf, :weights) && !isempty(ecdf.weights)
        y = cumsum(ecdf.weights)
        y ./= y[end]
        prepend!(y, [zero(eltype(y))])
    else
        T = float(eltype(ecdf.sorted_values))
        y = range(zero(T), one(T); length=length(x))
    end
    return x, y
end


text:
import FLOWFarm; const ff = FLOWFarm

# set initial turbine x and y locations
diam = 126.4
data = readdlm("inputfiles/layout_amalia.txt",  ' ', skipstart=1)
turbine_x = data[:, 1]
nturbines = length(turbine_x)
turbine_y = data[:, 2]
turbine_z = zeros(nturbines)

turbine_x = turbine_x .- turbine_x[1]
turbine_y = turbine_y .- turbine_y[1]

# calculate the number of turbines
nturbines = length(turbine_x)

# set turbine base heights
turbine_z = zeros(nturbines) .+ 0.0

# set turbine yaw values
turbine_yaw = zeros(nturbines)

# set turbine design parameters
rotor_diameter = zeros(nturbines) .+ diam # m
hub_height = zeros(nturbines) .+ 90.0   # m
cut_in_speed = zeros(nturbines) .+3.  # m/s
cut_out_speed = zeros(nturbines) .+25.  # m/s
rated_speed = zeros(nturbines) .+11.4  # m/s
rated_power = zeros(nturbines) .+5.0E6  # W
generator_efficiency = zeros(nturbines) .+ 0.944

# rotor swept area sample points (normalized by rotor radius)
rotor_points_y = [0.0]
rotor_points_z = [0.0]

# set flow parameters
data = readdlm("inputfiles/windrose_amalia_72dirs.txt",  ' ', skipstart=1)
winddirections = data[:, 1].*pi/180.0
windspeeds = data[:,2]
windprobabilities = data[:, 3]
nstates = length(windspeeds)

air_density = 1.1716  # kg/m^3
ambient_ti = 0.077
shearexponent = 0.15
ambient_tis = zeros(nstates) .+ ambient_ti
measurementheight = zeros(nstates) .+ hub_height[1]

# load power curve
cpctdata = readdlm("inputfiles/NREL5MWCPCT.txt",  ' ', skipstart=1)
velpoints = cpctdata[:,1]
cppoints = cpctdata[:,2]

# initialize power model
power_model = ff.PowerModelCpPoints(velpoints, cppoints)
power_models = Vector{typeof(power_model)}(undef, nturbines)
for i = 1:nturbines
    power_models[i] = power_model
end

# load thrust curve
ctpoints = cpctdata[:,3]

# initialize thurst model
ct_model = ff.ThrustModelCtPoints(velpoints, ctpoints)
ct_models = Vector{typeof(ct_model)}(undef, nturbines)
for i = 1:nturbines
    ct_models[i] = ct_model
end

# initialize wind shear model
wind_shear_model = ff.PowerLawWindShear(shearexponent)

# get sorted indecies
sorted_turbine_index = sortperm(turbine_x)

# initialize the wind resource definition
windresource = ff.DiscretizedWindResource(winddirections, windspeeds, windprobabilities, measurementheight, air_density, ambient_tis, wind_shear_model)

# set up wake and related models
wakedeficitmodel = ff.JensenTopHat()

wakedeflectionmodel = ff.JiminezYawDeflection()
wakecombinationmodel = ff.SumOfSquaresFreestreamSuperposition()
localtimodel = ff.LocalTIModelNoLocalTI()

# initialize model set
model_set = ff.WindFarmModelSet(wakedeficitmodel, wakedeflectionmodel, wakecombinationmodel, localtimodel)


text:
# Autogenerated wrapper script for fzf_jll for armv7l-linux-gnueabihf
export fzf

JLLWrappers.@generate_wrapper_header("fzf")
JLLWrappers.@declare_executable_product(fzf)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        fzf,
        "bin/fzf",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "mindomain.jl VariableSelection" begin
    @testset "MinDomainVariableSelection{TakeObjective=true}" begin
        trailer = SeaPearl.Trailer()
        cpmodel = SeaPearl.CPModel(trailer)
        x1 = SeaPearl.IntVar(1, 4, "x1", trailer)
        x2 = SeaPearl.IntVar(1, 5, "x2", trailer)
        y = SeaPearl.IntVar(1, 6, "y", trailer)

        SeaPearl.addVariable!(cpmodel, x1)
        SeaPearl.addVariable!(cpmodel, x2)
        SeaPearl.addVariable!(cpmodel, y)
        SeaPearl.addObjective!(cpmodel,y)

        variableselection = SeaPearl.MinDomainVariableSelection{true}()

        @test variableselection(cpmodel) == x1

        SeaPearl.removeAbove!(x2.domain, 3)
        @test variableselection(cpmodel) == x2

        SeaPearl.removeAbove!(y.domain, 2)
        @test variableselection(cpmodel) == y

        # With branchable branchable variables
        trailer = SeaPearl.Trailer()
        cpmodel = SeaPearl.CPModel(trailer)
        x1 = SeaPearl.IntVar(1, 4, "x1", trailer)
        x2 = SeaPearl.IntVar(1, 5, "x2", trailer)
        y = SeaPearl.IntVar(1, 6, "y", trailer)

        SeaPearl.addVariable!(cpmodel, x1)
        SeaPearl.addVariable!(cpmodel, x2; branchable=false)
        SeaPearl.addVariable!(cpmodel, y)
        SeaPearl.addObjective!(cpmodel,y)

        variableselection = SeaPearl.MinDomainVariableSelection{true}()

        @test variableselection(cpmodel) == x1

        SeaPearl.removeAbove!(x2.domain, 3)
        @test variableselection(cpmodel) == x1

        SeaPearl.removeAbove!(y.domain, 2)
        @test variableselection(cpmodel) == y
    end
    @testset "MinDomainVariableSelection{TakeObjective=false}" begin
        trailer = SeaPearl.Trailer()
        cpmodel = SeaPearl.CPModel(trailer)
        x1 = SeaPearl.IntVar(1, 4, "x1", trailer)
        x2 = SeaPearl.IntVar(1, 5, "x2", trailer)
        y = SeaPearl.IntVar(1, 6, "y", trailer)

        SeaPearl.addVariable!(cpmodel, x1)
        SeaPearl.addVariable!(cpmodel, x2)
        SeaPearl.addVariable!(cpmodel, y)
        SeaPearl.addObjective!(cpmodel,y)

        variableselection = SeaPearl.MinDomainVariableSelection{false}()

        @test variableselection(cpmodel) == x1

        SeaPearl.removeAbove!(x2.domain, 3)
        @test variableselection(cpmodel) == x2

        SeaPearl.removeAbove!(y.domain, 2)
        @test variableselection(cpmodel) == x2

        # With branchable variables
        trailer = SeaPearl.Trailer()
        cpmodel = SeaPearl.CPModel(trailer)
        x1 = SeaPearl.IntVar(1, 4, "x1", trailer)
        x2 = SeaPearl.IntVar(1, 5, "x2", trailer)
        y = SeaPearl.IntVar(1, 6, "y", trailer)

        SeaPearl.addVariable!(cpmodel, x1; branchable=false)
        SeaPearl.addVariable!(cpmodel, x2)
        SeaPearl.addVariable!(cpmodel, y)
        SeaPearl.addObjective!(cpmodel,y)

        variableselection = SeaPearl.MinDomainVariableSelection{false}()

        @test variableselection(cpmodel) == x2

        SeaPearl.removeAbove!(x2.domain, 3)
        @test variableselection(cpmodel) == x2

        SeaPearl.removeAbove!(y.domain, 2)
        @test variableselection(cpmodel) == x2
    end
end


text:
using Test

m = page_model()
include("../src/components/N2Ocycle.jl")

add_comp!(m, n2ocycle)

set_param!(m, :n2ocycle, :e_globalN2Oemissions, readpagedata(m,"test/validationdata/e_globalN2Oemissions.csv"))
set_param!(m, :n2ocycle, :y_year, [2009.,2010.,2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.]) #real value
set_param!(m, :n2ocycle, :y_year_0, 2008.) #real value
set_param!(m, :n2ocycle, :rtl_g_landtemperature, readpagedata(m,"test/validationdata/rtl_g_landtemperature.csv"))

##running Model
run(m)

conc=m[:n2ocycle,  :c_N2Oconcentration]
conc_compare=readpagedata(m,"test/validationdata/c_n2oconcentration.csv")

@test conc ≈ conc_compare rtol=1e-4


text:
using UnivariateFunctions
using StochasticIntegrals
using DataFrames
using Dates
using LinearAlgebra
using Distributions: Normal, MersenneTwister, quantile, Exponential
using Statistics: var, mean, cov
using Sobol
using Random
const tol = 10*eps()
twister = MersenneTwister(22)

USD_mean_reversion  = 0.15
USD_vol = 0.01
GBP_mean_reversion = 0.1
GBP_vol = 0.001
GBP_FX_Vol = 0.2

today    = Date(2016,12,1)
tommorow = Date(2016,12,2)

Barclays_vol = PE_Function(0.1, 0.0, 0.0, 0)
CS_vol       = PE_Function(0.2, 0.0, 0.0, 0)

brownian_corr_matrix = Hermitian([1.0 0.75;
                                  0.75 1.0])
assets = [ :CS, :BARC]
BARC_ito    = ItoIntegral(:BARC, Barclays_vol)
CS_ito    = ItoIntegral(:CS, CS_vol)
ito_integrals = Dict([:CS, :BARC] .=> [CS_ito, BARC_ito])

ito_set_ = ItoSet(brownian_corr_matrix, assets, ito_integrals)
covar = ForwardCovariance(ito_set_, years_from_global_base(today), years_from_global_base(tommorow))

T = Float64
conditioning_draws = Dict{Symbol,T}(:CS => 0.02)
conditional_mean, conditional_covar = generate_conditioned_distribution(covar, conditioning_draws)
conditional_mean[1] > 100*eps()

# Testing ito processes
stock_processes = Dict(assets .=> map(a -> ItoProcess(0.0, 0.0, PE_Function(0.00, 0.0, 0.0, 0), ito_integrals[a]), assets))
update_rates = Dict(assets .=> Exponential.(    rand(twister, length(assets))))

ts_syncronous  = make_ito_process_syncronous_time_series(stock_processes, covar, mean(a-> update_rates[a].θ, assets),1000; ito_twister = twister)
nrow(ts_syncronous) == 2000
ts_asyncronous = make_ito_process_non_syncronous_time_series(stock_processes, covar, update_rates, 1000; timing_twister = twister, ito_twister = twister)
nrow(ts_asyncronous) == 1000


text:
"""
    module ManageDatabase

Module for reading `coefs` from files and making database as needed for calculations
"""
module ManageDatabase

#using FileIO
#using JLD2

using ..CrystalMod:crystal
using ..CalcTB:coefs

using ..ThreeBodyTB:DATSDIR1
using ..ThreeBodyTB:DATSDIR2
using ..CalcTB:read_coefs

datdir1 = DATSDIR1 #primary directory
datdir2 = DATSDIR2 #backup directory

database_list = Set()

database_cached = Dict()
database_cached["SCF"] = true
database_cached["scf"] = true

"""
    function prepare_database(c::crystal)

Get ready database of precalculated `coefs` for `crystal`
"""
function prepare_database(c::crystal; directory = missing)
#    println("prepare c ", c.types)
    prepare_database(c.types, directory = directory)
end

"""
    function clear_database()

Remove loaded databasea
"""
function clear_database()

    if length(database_list) > 0
        for a in database_list
            delete!(database_list, a)
        end
    end

    for a in keys(database_cached)
        delete!(database_cached, a)
    end
    database_cached["SCF"] = true
    database_cached["scf"] = true

end

"""
    function prepare_database(at_list)
"""
function prepare_database(at_list; directory = missing)

    println("prepare atoms ", at_list)
    at_list = Symbol.(at_list)
    #    println("database_list ", database_list)
    s = Set(at_list)
    for s1 in s
        for s2 in s
            for s3 in s
                if ! ( Set([s1,s2,s3]) in database_list)
                    add_to_database(Set([s1,s2,s3]), directory = directory)
                end
            end
        end
    end
end

"""
    function add_to_database(s::Set)

Load elements or twobody terms from precalcuated `coefs` from files.
"""
function add_to_database(s::Set; directory = missing)#

    println("add_to_database $s")

    dirlist = []
    if !ismissing(directory)
        push!(dirlist, directory)
    end
    push!(dirlist, datdir1)
    push!(dirlist, datdir2)

    at_arr = collect(s)
    if length(s) == 1
        a1 = Symbol(at_arr[1])
        if !haskey(database_cached , (a1, a1))
            loaded = false
            for d in dirlist
                f =  "$d/els/coef.el.2bdy.$a1.xml.gz"

                if isfile(f) || isfile(f*".gz")
                    try
                        dat = read_coefs(f)
                        database_cached[(a1, a1)] = dat
                        println("added to cache ", (a1, a1))
                        loaded=true
                    catch
                        println("WARNING - error loading $f")
                    end
                end
                if loaded == true
                    break
                end

            end
            if loaded==false
                println("WARNING, FAILED TO LOAD 2bdy ", a1)
            end

        end

        if !haskey(database_cached , (a1, a1, a1))
#            f = "$defaultdatdir/v0.1_dat_3body_scf_pbesol_el.$a1.xml"
#            f =  "$datdir1/coef.el.3bdy.$a1.xml.gz"
#            f2 = "$datdir2/coef.el.3bdy.$a1.xml.gz"

            loaded = false
            for d in dirlist
                f = "$d/els/coef.el.3bdy.$a1.xml.gz"

                if isfile(f) || isfile(f*".gz")
                    try
                        #                    jldopen(f)
                        dat = read_coefs(f)
                        database_cached[(a1, a1, a1)] = dat
                        println("added to cache ", (a1, a1, a1))
                        loaded == true
                    catch
                        println("WARNING - error loading $f")
                    end
                    loaded=true
                end
                if loaded == true
                    break
                end

            end
            if loaded == false
                println("WARNING, FAILED LOADING 3bdy ", a1)
            end

#                elseif isfile(f2) || isfile(f2*".gz")
#                try
#                    jldopen(f)
#                    dat = read_coefs(f2)
#                    database_cached[(a1, a1, a1)] = dat
#                    println("added to cache ", (a1, a1, a1))
#                catch
#                    println("WARNING - error loading $f2")
#                end
#            else
#                println("WARNING, no file for 3bdy database $a1")
#            end
        end

    elseif length(s) == 2

        a1 = Symbol(at_arr[1])
        a2 = Symbol(at_arr[2])
        if !haskey(database_cached , (a1, a2))

#            fab = "$defaultdatdir/v0.1_dat_2body_scf_pbesol_binary.$a1.$a2.xml"
#            fba = "$defaultdatdir/v0.1_dat_2body_scf_pbesol_binary.$a2.$a1.xml"

#            fab =  "$datdir1/coef.el.2bdy.$a1.$a2.xml.gz"
#            fba =  "$datdir1/coef.el.2bdy.$a2.$a1.xml.gz"

#            fab2 = "$datdir2/coef.el.2bdy.$a1.$a2.xml.gz"
#            fba2 = "$datdir2/coef.el.2bdy.$a2.$a1.xml.gz"

            loaded = false
            for d in dirlist

                fab =  "$d/binary/coef.el.2bdy.$a1.$a2.xml.gz"
                fba =  "$d/binary/coef.el.2bdy.$a2.$a1.xml.gz"

                f=missing
                if isfile(fab) || isfile(fab*".gz")
                    f = fab
                elseif isfile(fba) || isfile(fba*".gz")
                    f = fba
                end
#            elseif isfile(fab2) || isfile(fab2*".gz")
#                f = fab2
#            elseif isfile(fba2) || isfile(fba2*".gz")
#                f = fba2
#            else
#                f = missing
#                println("WARNING - binary file missing ")
#            end

                if !ismissing(f)
                    try
                        dat = read_coefs(f)

                        database_cached[(a1, a2)] = dat
                        database_cached[(a2, a1)] = dat
                        println("added to cache ", (a1, a2), " twobody ")
                        loaded = true
                    catch
                        println("WARNING - error loading binary file $f")
                    end
                end
                if loaded == true
                    break
                end
            end
            if loaded == false
                println("WARNING, FAILED LOAD 2BDY ", a1, " ", a2)
            end

#################

#            fab = "$defaultdatdir/v0.1_dat_3body_scf_pbesol_binary.$a1.$a2.xml"
#            fba = "$defaultdatdir/v0.1_dat_3body_scf_pbesol_binary.$a2.$a1.xml"

#            fab =  "$datdir1/coef.el.3bdy.$a1.$a2.xml.gz"
#            fba =  "$datdir1/coef.el.3bdy.$a2.$a1.xml.gz"

#            fab2 = "$datdir2/coef.el.3bdy.$a1.$a2.xml.gz"
#            fba2 = "$datdir2/coef.el.3bdy.$a2.$a1.xml.gz"

            loaded = false
            for d in dirlist
                fab =  "$d/binary/coef.el.3bdy.$a1.$a2.xml.gz"
                fba =  "$d/binary/coef.el.3bdy.$a2.$a1.xml.gz"

                f=missing
                if isfile(fab) || isfile(fab*".gz")
                    f = fab
                elseif isfile(fba) || isfile(fba*".gz")
                    f = fba
                end
#                elseif isfile(fab2) || isfile(fab2*".gz")
#                    f = fab2
#                elseif isfile(fba2) || isfile(fba2*".gz")
#                f = fba2

#            else
#                f = missing
 #               println("WARNING - binary file missing ")
 #           end

                if !ismissing(f)
                    try
#                    jldopen(f)
                        dat = read_coefs(f)

                        database_cached[(a1,a1,a2)] = dat
                        database_cached[(a1,a2,a1)] = dat
                        database_cached[(a2,a1,a1)] = dat

                        database_cached[(a1,a2,a2)] = dat
                        database_cached[(a2,a1,a2)] = dat
                        database_cached[(a2,a2,a1)] = dat
                        println("added to cache ", (a1, a2), " threebody ")
                        loaded=true
                    catch
                        println(" MISSING FILE - WARNING loading binary 3body $f ")
                    end
                end
                if loaded == true
                    break
                end

            end
            if loaded==false
                println("WARNING, FAILED TO LOAD 3bdy ", a1, " " , a2)
            end

        end
    elseif length(s) == 3

        println("warning, ternary not currently supported")
        return

        a1 = Symbol(at_arr[1])
        a2 = Symbol(at_arr[2])
        a3 = Symbol(at_arr[3])

        ats  = sort([a1,a2,a3])

        a1 = ats[1]
        a2 = ats[2]
        a3 = ats[3]

        if !haskey(database_cached , Tuple(ats))

            loaded = false
            for d in dirlist

                f =  "$d/tern_$a1/coef.el.3bdy.$a1.$a2.$a3.xml.gz"
                if  isfile(f) || isfile(f*".gz")
                    try
                        dat = read_coefs(f)
                        database_cached[(a1,a2,a3)] = dat
                        database_cached[(a1,a3,a2)] = dat
                        database_cached[(a2,a1,a3)] = dat
                        database_cached[(a2,a3,a1)] = dat
                        database_cached[(a3,a1,a2)] = dat
                        database_cached[(a3,a2,a1)] = dat
                        println("added to cache ", (a1, a2, a3), " threebody ")
                        loaded=true
                    catch
                        println("failed to load $f")
                    end
                end
                if loaded == true
                    break
                end
            end
            if loaded == false
                println("WARNING, FAILED TO LOAD ", [a1, a2, a3])
            end
        end

    else
        println("WARNING - not setup for quaternary+ ", s)
    end

    push!(database_list, s)

end

end #end module


text:
export CGPChromo, node_genes, get_positions

# vanilla CGP

struct CGPChromo <: Chromosome
    genes::Array{Float64}
    nodes::Array{CGPNode}
    outputs::Array{Int64}
    nin::Int64
    nout::Int64
end

function CGPChromo(genes::Array{Float64}, nin::Int64, nout::Int64)::CGPChromo
    num_nodes = Int64(ceil((length(genes)-nin-nout)/4))
    rgenes = reshape(genes[(nin+nout+1):end], (4, num_nodes))'
    connections = Array{Int64}(undef, 2, nin+num_nodes)
    connections[:, 1:nin] = zeros(2, nin)
    positions = collect(1:(nin+num_nodes))/(1.0*(nin+num_nodes))
    fc = deepcopy(hcat(zeros(2, nin), [rgenes[:, 1]'; rgenes[:, 2]']))
    e = (1.0 .- positions) .* Config.recurrency .+ positions
    for j in 1:size(fc)[1]
        fc[j, :] .*= e
    end
    if Config.recurrency == 0
        for i in (nin+1):length(positions)
            for j in 1:size(fc)[1]
                fc[j, i] = max(fc[j,i], positions[i-1])
            end
        end
    end
    connections = snap(fc, positions)
    functions = Array{Function}(undef, nin+num_nodes)
    functions[1:nin] .= Config.f_input
    functions[(nin+1):end] = map(i->Config.index_in(Config.functions, i), rgenes[:, 3])
    outputs = Int64.(ceil.(genes[nin .+ (1:nout)]*(nin+num_nodes)))
    active = find_active(nin, outputs, connections)
    params = [zeros(nin); 2.0*rgenes[:, 4] .- 1.0]
    nodes = Array{CGPNode}(undef, nin+num_nodes)
    for i in 1:(nin+num_nodes)
        nodes[i] = CGPNode(connections[:, i], functions[i], active[i], params[i])
    end
    CGPChromo(genes, nodes, outputs, nin, nout)
end

function CGPChromo(nin::Int64, nout::Int64)::CGPChromo
    n_nodes = Config.static_node_size
    if Config.bloat()
        n_nodes = Config.starting_nodes
    end
    CGPChromo(rand(nin+nout+4*n_nodes), nin, nout)
end

function CGPChromo(c::CGPChromo)::CGPChromo
    gene_mutate(c)
end

function node_genes(c::CGPChromo)
    4
end

function get_positions(c::CGPChromo)
    collect(1:length(c.nodes))/(1.0*length(c.nodes))
end

function add_subtree(c::CGPChromo)
    add_nodes(c)
end

function delete_subtree(c::CGPChromo)
    delete_nodes(c)
end


text:
module JuSpglibTests

using Test

using JuSpglib

@testset "Get symmetry of a simple lattices" begin
    cell = Cell([0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0], transpose([0 0 0; 0.25 0.25 0.25]), [:Si, :Si])
    symops = get_symmetry(cell)
    @test length(symops) == 48
    symbol = get_international(cell)
    @test symbol == "Fd-3m"
    symbol = get_schoenflies(cell)
    @test symbol == "Oh^7"
end

@testset "Zincblende" begin
    cell = Cell([0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0], transpose([0 0 0; 0.25 0.25 0.25]), [:Zn, :S])
    symops = get_symmetry(cell)
    @test length(symops) == 24
    symbol = get_international(cell)
    @test symbol == "F-43m"
    symbol = get_schoenflies(cell)
    @test symbol == "Td^2"
end

end


text:
function buildlodf(branches::Array{T}, nodes::Array{Bus}, dist_slack::Array{Float64}=[0.1] ) where {T<:Branch}
    linecount = length(branches)
    ptdf , a = PowerSystems.buildptdf(branches,nodes,dist_slack)
    ptdf = ptdf.data; a = a.data;
    H =  gemm('N','N',ptdf,a)
    ptdf_denominator = H;
    for iline = 1:linecount
        if (1.0 - ptdf_denominator[iline,iline] ) < 1.0E-06
            ptdf_denominator[iline,iline] = 0.0;
        end
    end
    (Dem, dipiv, dinfo) = getrf!(Matrix{Float64}(LinearAlgebra.I,linecount,linecount) - Array(LinearAlgebra.Diagonal(ptdf_denominator)))
    LODF = gemm('N','N',H,getri!(Dem,dipiv))
    LODF = LODF - Array(LinearAlgebra.Diagonal(LODF)) - Matrix{Float64}(LinearAlgebra.I,linecount,linecount)
    return LODF
end


text:
using Printf

using BenchmarkTools
using CUDA
using DataFrames
using PrettyTables

using Oceananigans
using Oceananigans.Architectures
using JULES

using BenchmarkTools: prettytime, prettymemory

Archs = [CPU]
@hascuda Archs = [CPU, GPU]

Ns = [32, 192]
Tvars = [Energy, Entropy]
Gases = [DryEarth, DryEarth3]

suite = BenchmarkGroup()

sync_step!(model) = time_step!(model, 1)
sync_step!(model::CompressibleModel{GPU}) = CUDA.@sync time_step!(model, 1)

for Arch in Archs, N in Ns, Gas in Gases, Tvar in Tvars
    @info "Running static atmosphere benchmark [$Arch, N=$N, $Tvar, $Gas]..."

    grid = RegularCartesianGrid(size=(N, N, N), halo=(2, 2, 2), extent=(1, 1, 1))
    model = CompressibleModel(architecture=Arch(), grid=grid, thermodynamic_variable=Tvar(), gases=Gas())

    sync_step!(model) # warmup

    b = @benchmark sync_step!($model) samples=10
    display(b)

    key = (Arch, N, Gas, Tvar)
    suite[key] = b
end

function benchmarks_to_dataframe(suite)
    df = DataFrame(architecture=[], size=[], gases=[],
                   thermodynamic_variable=[], min=[], median=[],
                   mean=[], max=[], memory=[], allocs=[])

    for Arch in Archs, N in Ns, Gas in Gases, Tvar in Tvars
        b = suite[Arch, N, Gas, Tvar]

        d = Dict(
            "architecture" => Arch,
            "size" => "$(N)³",
            "gases" => Gas,
            "thermodynamic_variable" => Tvar,
            "min" => minimum(b.times) |> prettytime,
            "median" => median(b.times) |> prettytime,
            "mean" => mean(b.times) |> prettytime,
            "max" => maximum(b.times) |> prettytime,
            "memory" => prettymemory(b.memory),
            "allocs" => b.allocs
        )

        push!(df, d)
    end

    return df
end

header = ["Arch" "Size" "Gases" "ThermoVar" "min" "median" "mean" "max" "memory" "allocs"]

df = benchmarks_to_dataframe(suite)
pretty_table(df, header, title="Static atmosphere benchmarks", nosubheader=true)

function gpu_speedups(suite)
    df = DataFrame(size=[], gases=[], thermodynamic_variable=[], speedup=[])

    for N in Ns, Gas in Gases, Tvar in Tvars
        b_cpu = suite[CPU, N, Gas, Tvar] |> median
        b_gpu = suite[GPU, N, Gas, Tvar] |> median
        b_ratio = ratio(b_cpu, b_gpu)

        d = Dict(
            "size" => "$(N)³",
            "gases" => Gas,
            "thermodynamic_variable" => Tvar,
            "speedup" => @sprintf("%.3fx", b_ratio.time)
        )

        push!(df, d)
    end

    return df
end

header = ["Size" "Gases" "ThermoVar" "speedup"]

df = gpu_speedups(suite)
pretty_table(df, header, title="Static atmosphere speedups", nosubheader=true)


text:
import JuMP as J
using Ipopt

# Dados de entrada

Sy=350.0;
r1=0.4;
E=2000.0;
P=200.0;
B=1.0;

f(r2,H) = pi*(r2^2)*(H^2+B^2)^(1/2)

g1(r2,H) = P*B/(pi*H*(r1^2)*Sy) - 1;
g2(r2,H) = P*((B^2+H^2)^(1/2))/(H*pi*r2^2*Sy)-1;
g3(r2,H) = 4*P*(B^2+H^2)^(3/2)/(H*pi^3*r2^4*E)-1;

model = J.Model(Ipopt.Optimizer)

J.set_optimizer_attribute(model, "print_level",0)

#registrando o número de variáveis

J.register(model, :g1, 2, g1; autodiff = true)
J.register(model, :g2, 2, g2; autodiff = true)
J.register(model, :g3, 2, g3; autodiff = true)
J.register(model, :f, 2, f; autodiff = true)

J.@variable(model, 0.4 <= r2 <= 0.8);
J.@variable(model, 1.0 <= H <= 6.0);

J.@NLobjective(model, Min, f(r2,H));

J.@NLconstraint(model, g1(r2,H) <= 0);
J.@NLconstraint(model, g2(r2,H) <= 0);
J.@NLconstraint(model, g3(r2,H) <= 0);

println(model)
J.optimize!(model)

println("Termination status: ", J.termination_status(model))
println("Primal status: ", J.primal_status(model))
println("f*: ", J.objective_value(model))
println("r2* : ", J.value(r2), ", H* : ", J.value(H))

a = 4;


text:
import MapImages: MapImage
import CSV

geolist78 = Sario.load_geolist_from_h5(
    "/Users/scott/Documents/Learning/masters-thesis/zoom_images78/stack_max700_alpha100.h5",
    "stack/1",
);
dfmt = dateformat"y-m-d H:M:S.s"
stack78 = MapImage(
    "/Users/scott/Documents/Learning/masters-thesis/zoom_images78/stack_max700_alpha100.h5",
    "stack/1",
);
eqs15 = CSV.read("TXAR_Q12events_2015.csv", dateformat = dfmt);
eqs17 = CSV.read("TXAR_Q12events_2017.csv", dateformat = dfmt);
eqs17texnet = CSV.read("texnet_events_2017.csv", dateformat = dfmt);

@time include("/Users/scott/repos/InsarTimeseries.jl/scripts/plotting.jl");

animate_imgs_vs_pts(
    stack78,
    eqs17texnet,
    geolist78;
    latcol = :Latitude,
    loncol = :Longitude,
    sizecol = :Magnitude,
    datecol = :date,
    alpha = 0.2,
    vm = 6,
    c = "r",
    outname = "animation_texnet78.gif",
)
animate_imgs_vs_pts(
    stack78,
    [eqs15; eqs17],
    geolist78;
    latcol = :Latitude,
    loncol = :Longitude,
    sizecol = :Magnitude,
    datecol = :date,
    alpha = 0.2,
    vm = 6,
    c = "r",
    outname = "animation_TXAR78.gif",
)


text:
export
    qpdf,
    qpdf!

"""
    qpdf([T=Float64], d::GaussianStateBHD, θ::Real, x::Real)
    qpdf([T=Float64], d::GaussianStateBHD, θs::AbstractRange, xs::AbstractRange)
    qpdf([T=Float64], d::QuantumStateBHD, θ::Real, x::Real)
    qpdf([T=Float64], d::QuantumStateBHD, θs::AbstractRange, xs::AbstractRange)
    qpdf([T=Float64], ρ::AbstractArray, θ::Real, x::Real)
    qpdf([T=Float64], ρ::AbstractArray, θs::AbstractRange, xs::AbstractRange)

Quadrature prabability in intensity-to-measurement-phase quadrature coordinate.

## Arguments

* `T`: Data type, default as Float64.
* `state`: State can be `GaussianStateBHD` distribution, `QuantumStateBHD` distribution or density matrix.
* `θ`: Measurement phase, can be `Real` or `AbstractRange`.
* `x`: Intensity in quadrature coordinate, can be `Real` or `AbstractRange`.

``p(\\rho, \\theta, x) = tr(\\hat{\\Pi}(\\theta, x) \\rho)``
"""
qpdf(state, θ, x) = qpdf(Float64, state, θ, x)

function qpdf(T::Type{<:Real}, d::GaussianStateBHD, θ::Real, x::Real)
    μ = QuantumStateDistributions.mean(d, θ)
    σ = QuantumStateDistributions.std(d, θ)

    return T(pdf(Normal(μ, σ), x))
end

function qpdf(T::Type{<:Real}, d::GaussianStateBHD, θs::AbstractRange, xs::AbstractRange)
    m, n = length(θs), length(xs)

    μs = Vector{T}(undef, m)
    σs = Vector{T}(undef, m)
    gaussians = Vector{Normal}(undef, m)
    𝐩 = Matrix{T}(undef, m, n)

    return qpdf!(μs, σs, gaussians, 𝐩, d, θs, xs)
end

function qpdf!(μs, σs, gaussians, 𝐩, d::GaussianStateBHD, θs::AbstractRange, xs::AbstractRange)
    μs .= QuantumStateDistributions.mean(d, θs)
    σs .= QuantumStateDistributions.std(d, θs)

    gaussians .= Normal.(μs, σs)

    for i in 1:length(θs)
        for (j, x) in enumerate(xs)
		    𝐩[i, j] = pdf(gaussians[i], x)
        end
    end

    return 𝐩
end

qpdf(T::Type{<:Real}, d::QuantumStateBHD, θ, x) = qpdf(T, d.ρ, θ, x)

function qpdf(T::Type{<:Real}, ρ::AbstractMatrix, θ::Real, x::Real)
    dim = size(ρ, 1)
    𝛑̂_res = Matrix{Complex{T}}(undef, dim, dim)

    return qpdf!(𝛑̂_res, ρ, θ, x)
end

function qpdf!(𝛑̂_res::AbstractMatrix, ρ::AbstractMatrix, θ::Real, x::Real)
    dim = size(ρ, 1)

    return real_tr_mul(𝛑̂!(𝛑̂_res, θ, x, dim=dim), ρ)
end

function qpdf(T::Type{<:Real}, ρ::AbstractMatrix, θs::AbstractRange, xs::AbstractRange)
    dim = size(ρ, 1)
    𝛑̂_res_vec = [Matrix{Complex{T}}(undef, dim, dim) for _ in 1:Threads.nthreads()]
    𝐩 = Matrix{T}(undef, length(θs), length(xs))

    return qpdf!(𝛑̂_res_vec, 𝐩, ρ, θs, xs)
end

function qpdf!(
    𝛑̂_res_vec::AbstractVector{Matrix{Complex{T}}}, 𝐩::Matrix{T},
    ρ::AbstractMatrix, θs::AbstractRange, xs::AbstractRange
) where {T}
    @sync for (j, x) in enumerate(xs)
        for (i, θ) in enumerate(θs)
            Threads.@spawn 𝐩[i, j] = qpdf!(𝛑̂_res_vec[Threads.threadid()], ρ, θ, x)
        end
    end

    return 𝐩
end

##### for arb. state in intensity-to-measurement-phase quadrature coordinate #####

"""
    ψₙ(n::Integer, θ::Real, x::Real)

Eigenstate of BHD measurement operator.

``\\psi_n(\\theta, x) = \\langle n | \\theta, x \\rangle``
"""
function ψₙ(n::Integer, θ::Real, x::Real)
    # |θ, x⟩ = ∑ₙ |n⟩ ⟨n|θ, x⟩ = ∑ₙ ψₙ(θ, x) |n⟩
    # ⟨n|θ, x⟩ = ψₙ(θ, x) = exp(im n θ) (2/π)^(1/4) exp(-x^2) Hₙ(√2 x)/√(2^n n!)

    return (2/π)^(1/4) * exp(im*n*θ - x^2) * hermiteh(n, sqrt(2)x) / sqrt(2^n * factorial(n))
end

"""
    𝛑̂(θ::Real, x::Real; dim::Integer)

BHD measurement operator.

``\\hat{\\Pi}_{m, n}(\\theta, x) = \\langle m | \\hat{\\Pi}(\\theta, x) | n \\rangle = \\langle m | \\theta, x \\rangle \\langle \\theta, x | n \\rangle``
"""
𝛑̂(θ::Real, x::Real; dim) = 𝛑̂(ComplexF64, θ, x, dim=dim)

function 𝛑̂(T::Type{<:Complex}, θ::Real, x::Real; dim)
    result = Matrix{T}(undef, dim, dim)

    return 𝛑̂!(result, θ, x, dim=dim)
end

function 𝛑̂!(result::AbstractMatrix{<:Complex}, θ::Real, x::Real; dim)
    view(result, :, 1) .= ψₙ.(big.(0:dim-1), θ, x)
    result .= view(result, :, 1) * view(result, :, 1)'

    return result
end

# #########
# # utils #
# #########

real_tr_mul(𝐚, 𝐛) = sum(real(𝐚[i, :]' * 𝐛[:, i]) for i in 1:size(𝐚, 1))


text:
InputFolder = './Images/NucleiCytoo/';
OutputFolder = './Results/Images/NucleiCytoo/';
Rescale = 1;
Fill = 1;

@iA = '*C00*.tif';		% Image filter

@fxg_sSBFRegMax [iA] > [M];
params.UnderSamp = 2;
params.Rmin = 7;
params.Rmax = 13;
params.Rstep = 1;
params.Band = 1;
params.lambdaMag = 1;
params.LocMaxThr = 4;
params.BckLvl = 20;
params.Sym = 1;
/endf

/show iA > M;
/keep M > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
@testset "group" begin
    @test group(identity, 11:20) == Dict(Pair.(11:20, (x->[x]).(11:20)))
    @test group(iseven, 1:10) == Dict(true => [2,4,6,8,10], false => [1,3,5,7,9])

    @test group(iseven, x -> x*2, 1:10) == Dict(true => [4,8,12,16,20], false => [2,6,10,14,18])

    @test group((x,y) -> iseven(x+y), (x,y) -> x, 1:10, [1,3,4,2,5,6,4,2,3,9]) == Dict(true => [1,4,5,6,8,9], false => [2,3,7,10])
end

@testset "groupinds" begin
    @test groupinds(identity, 11:20) == Dict(Pair.(11:20, (x->[x]).(1:10)))
    @test groupinds(iseven, 11:20) == Dict(true => [2,4,6,8,10], false => [1,3,5,7,9])
end

@testset "groupview" begin
    @test groupview(identity, 11:20)::Groups == group(identity, 11:20)::Dict
    @test groupview(iseven, 11:20)::Groups == group(iseven, 11:20)::Dict
end

@testset "groupreduce" begin
    @test groupreduce(identity, +, 1:10) == Dict(Pair.(1:10, 1:10))
    @test groupreduce(iseven, +, 1:10) == Dict(true => 30, false => 25)

    @test groupreduce(iseven, x -> x*2, +, 1:10) == Dict(true => 60, false => 50)

    @test groupreduce(iseven, x -> x*2, +, 1:10; init=10) == Dict(true => 70, false => 60)

    @test groupreduce((x,y) -> iseven(x+y), (x,y) -> x+y, +, 1:10, 1:10; init=10) == Dict(true => 120)
    @test groupreduce((x,y) -> iseven(x+y), (x,y) -> x+y, +, 1:10, [1,3,4,2,5,6,4,2,3,9]; init=10) == Dict(true => 62, false => 52)

    @test groupcount(iseven, 1:10) == Dict(true => 5, false => 5)
    @test groupsum(iseven, 1:10) == Dict(true => 30, false => 25)
    @test groupprod(iseven, 1:10) == Dict(true => 2*4*6*8*10, false => 1*3*5*7*9)
end


text:
include("adaptive_parameters.jl")
include("center_of_mass.jl")

mutable struct ECA <: AbstractParameters
    η_max::Float64
    K::Int
    N::Int
    N_init::Int
    p_exploit::Float64
    p_bin::Float64
    p_cr::Array{Float64}
    ε::Float64
    adaptive::Bool
    resize_population::Bool
end

"""
    ECA(;
        η_max = 2.0,
        K = 7,
        N = 0,
        N_init = N,
        p_exploit = 0.95,
        p_bin = 0.02,
        p_cr = Float64[],
        adaptive = false,
        resize_population = false,
        information = Information(),
        options = Options()
    )

Parameters for the metaheuristic ECA: step-size `η_max`,`K` is number of vectors to
generate the center of mass, `N` is the population size.

# Example

```jldoctest
julia> f(x) = sum(x.^2)
f (generic function with 1 method)

julia> optimize(f, [-1 -1 -1; 1 1 1.0], ECA())
+=========== RESULT ==========+
  iteration: 1429
    minimum: 3.3152400000000004e-223
  minimizer: [4.213750597785841e-113, 5.290977430907081e-112, 2.231685329262638e-112]
    f calls: 29989
 total time: 0.1672 s
+============================+

julia> optimize(f, [-1 -1 -1; 1 1 1.0], ECA(N = 10, η_max = 1.0, K = 3))
+=========== RESULT ==========+
  iteration: 3000
    minimum: 0.000571319
  minimizer: [-0.00017150889316537758, -0.007955828028420616, 0.022538733289139145]
    f calls: 30000
 total time: 0.1334 s
+============================+
```

"""
function ECA(;
    η_max::Float64 = 2.0,
    K::Int = 7,
    N::Int = 0,
    N_init::Int = N,
    p_exploit::Float64 = 0.95,
    p_bin::Float64 = 0.02,
    p_cr::Array{Float64} = Float64[],
    ε::Float64 = 0.0,
    adaptive::Bool = false,
    resize_population::Bool = false,
    information = Information(),
    options = Options(),
)

    N_init = N

    parameters = ECA(
        η_max,
        K,
        N,
        N_init,
        p_exploit,
        p_bin,
        p_cr,
        ε,
        adaptive,
        resize_population,
    )
    Algorithm(
        parameters,
        information = information,
        options = options,
    )

end

function update_state!(
    status,
    parameters::ECA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
)
    K = parameters.K
    I = randperm(parameters.N)
    D = size(problem.bounds, 2)
    is_multiobjective = typeof(status.population[1]) == xFgh_indiv

    parameters.adaptive && (Mcr_fail = zeros(D))

    stop = false
    a = problem.bounds[1, :]
    b = problem.bounds[2, :]

    if is_multiobjective
        shuffle!(status.population)
    end

    # For each elements in Population
    for i = 1:parameters.N
        p = status.f_calls / options.f_calls_limit

        # current
        x = status.population[i].x

        # generate U masses
        U = getU(status.population, parameters.K, I, i, parameters.N)

        # generate center of mass
        c, u_worst, u_best = center(U)

        # stepsize
        η = parameters.η_max * rand()

        # u: worst element in U
        u = U[u_worst].x

        # current-to-center/bin
        if p < parameters.p_exploit
            # u: worst element in U
            u = U[u_worst].x

            # current-to-center/bin
            y = x .+ η .* (c .- u)
        elseif parameters.p_exploit < 0
            y =
                x .+ (1 - p^5) * η * (c .- u) .+
                (p^5) * η * (status.best_sol.x .- c)
        else
            # current-to-best/bin
            y = x .+ η .* (status.best_sol.x .- c)
        end

        # binary crossover
        y, M_current = crossover(U[u_best].x, y, parameters.p_cr)

        evo_boundary_repairer!(y, c, problem.bounds)

        sol = create_solution(y, problem)

        # replace worst element
        if is_better(sol, status.population[i])
            wi = argworst(status.population)
            status.population[wi] = sol
            if is_better(sol, status.best_sol)
                status.best_sol = sol
            end
        else
            if is_multiobjective && !is_better(status.population[i], sol)
                push!(status.population, sol)
            end
            parameters.adaptive && (Mcr_fail += M_current)
        end

        stop_criteria!(status, parameters, problem, information, options)
        status.stop && break
    end

    status.f_calls = problem.f_calls

    if status.stop
        return
    end

    if parameters.adaptive
        parameters.p_cr =
            adaptCrossover(parameters.p_cr, Mcr_fail / parameters.N)
    end

    if is_multiobjective && length(status.population) > parameters.N
        status.population = truncate_population!(status.population, parameters.N, (w,z) -> is_better(w, z))
    end

    p = status.f_calls / options.f_calls_limit

    if parameters.resize_population
        K = parameters.K

        # new size
        parameters.N = 2K .+ round(Int, (1 - p) * (parameters.N_init .- 2K))

        if parameters.N < 2K
            parameters.N = 2K
        end

        status.population = resizePop!(
            status.population,
            parameters.N,
            K
        )
    end
end

function initialize!(
        status,
        parameters::ECA,
        problem::AbstractProblem,
        information::Information,
        options::Options,
        args...;
        kargs...
    )
    D = size(problem.bounds, 2)

    if parameters.N <= parameters.K
        parameters.N = parameters.K * D
    end

    if options.f_calls_limit == 0
        options.f_calls_limit = 10000D
        options.debug &&
            @warn( "f_calls_limit increased to $(options.f_calls_limit)")
    end

    if options.iterations == 0
        options.iterations = div(options.f_calls_limit, parameters.N) + 1
    end

    N_init = parameters.N

    if parameters.adaptive
        parameters.p_cr = rand(D)
    else
        parameters.p_cr = parameters.p_bin .* ones(D)
    end

    # initialize!(problem, nothing, parameters, status, information, options)
    st = gen_initial_state(problem,parameters,information,options, status)
    st

end

function final_stage!(
    status,
    parameters::ECA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
    )
    status.final_time = time()

    # compute Pareto front if it is a multiobjective problem
    if typeof(status.population[1].f) <: Array
        options.debug && @info "Computing Pareto front..."
        status.best_sol = get_pareto_front(status.population, is_better)
    end
end


text:
using Oceananigans: AbstractModel
using Oceananigans.Grids
using Oceananigans.Utils: tupleit

struct KernelComputedField{X, Y, Z, A, S, D, G, T, K, B, F, P} <: AbstractDataField{X, Y, Z, A, G, T, 3}
                     data :: D
             architecture :: A
                     grid :: G
                   kernel :: K
      boundary_conditions :: B
    computed_dependencies :: F
               parameters :: P
                   status :: S

    function KernelComputedField{X, Y, Z}(kernel::K, arch::A, grid::G;
                                          boundary_conditions::B = ComputedFieldBoundaryConditions(grid, (X, Y, Z)),
                                           computed_dependencies = (),
                                                   parameters::P = nothing,
                                                         data::D = new_data(arch, grid, (X, Y, Z)),
                                                recompute_safely = true) where {X, Y, Z, A, D, B, G, K, P}

        computed_dependencies = tupleit(computed_dependencies)

        # Use FieldStatus if we want to avoid always recomputing
        status = recompute_safely ? nothing : FieldStatus(0.0)

        S = typeof(status)
        F = typeof(computed_dependencies)
        T = eltype(grid)

        return new{X, Y, Z, A, S, D, G, T, K, B, F, P}(
            data, arch, grid, kernel, boundary_conditions, computed_dependencies, parameters, status)
    end
end

"""
    KernelComputedField(X, Y, Z, kernel, model;
                        boundary_conditions = ComputedFieldBoundaryConditions(grid, (X, Y, Z)),
                        computed_dependencies = (),
                        parameters = nothing,
                        data = nothing,
                        recompute_safely = true)

Builds a `KernelComputedField` at `X, Y, Z` computed with `kernel` and `model.architecture` and `model.grid`, with `boundary_conditions`.

`computed_dependencies` are an iterable of `AbstractField`s or other objects on which `compute!` is called prior to launching `kernel`.

`data` is a three-dimensional `OffsetArray` of scratch space where the kernel computation is stored.

If `data=nothing` (the default) then additional memory will be allocated to store the `data` of `KernelComputedField`.

If `isnothing(parameters)`, `kernel` is launched with the function signature

`kernel(data, grid, computed_dependencies...)`

Otherwise, `kernel` is launched with the function signature

`kernel(data, grid, computed_dependencies..., parameters)`

`recompute_safely` (default: `true`) determines whether the `KernelComputedField` is "recomputed" if embedded in the expression
tree of another operation.
    - If `recompute_safely=true`, the `KernelComputedField` is always recomputed.
    - If `recompute_safely=false`, the `KernelComputedField` will not be recomputed if its status is up-to-date.

Example
=======

```julia
using KernelAbstractions: @index, @kernel
using Oceananigans.Fields: AveragedField, KernelComputedField, compute!
using Oceananigans.Grids: Center, Face

@inline ψ′²(i, j, k, grid, ψ, Ψ) = @inbounds (ψ[i, j, k] - Ψ[i, j, k])^2
@inline ψ′²(i, j, k, grid, ψ, Ψ::Number) = @inbounds (ψ[i, j, k] - Ψ)^2

@kernel function compute_variance!(var, grid, ϕ, Φ)
    i, j, k = @index(Global, NTuple)

    @inbounds var[i, j, k] = ψ′²(i, j, k, grid, ϕ, Φ)
end

u, v, w = model.velocities

U = AveragedField(u, dims=(1, 2))
V = AveragedField(v, dims=(1, 2))

u′² = KernelComputedField(Face, Center, Center, compute_variance!, model; computed_dependencies=(u, U))
v′² = KernelComputedField(Center, Face, Center, compute_variance!, model; computed_dependencies=(v, V))
w′² = KernelComputedField(Center, Center, Face, compute_variance!, model; computed_dependencies=(w, 0))

compute!(u′²)
compute!(v′²)
compute!(w′²)
```
"""
KernelComputedField(X, Y, Z, kernel, model::AbstractModel; kwargs...) =
    KernelComputedField{X, Y, Z}(kernel, model.architecture, model.grid; kwargs...)

KernelComputedField(X, Y, Z, kernel, arch::AbstractArchitecture, grid::AbstractGrid; kwargs...) =
    KernelComputedField{X, Y, Z}(kernel, arch, grid; kwargs...)

function compute!(kcf::KernelComputedField{X, Y, Z}) where {X, Y, Z}

    for dependency in kcf.computed_dependencies
        compute!(dependency)
    end

    arch = architecture(kcf)

    args = isnothing(kcf.parameters) ?
        tuple(kcf.data, kcf.grid, kcf.computed_dependencies...) :
        tuple(kcf.data, kcf.grid, kcf.computed_dependencies..., kcf.parameters)

    event = launch!(arch, kcf.grid, :xyz, kcf.kernel, args...;
                    location=(X, Y, Z), include_right_boundaries=true)

    wait(device(arch), event)

    fill_halo_regions!(kcf, arch)

    return nothing
end

compute!(field::KernelComputedField{X, Y, Z, <:FieldStatus}, time) where {X, Y, Z} =
    conditional_compute!(field, time)

Adapt.adapt_structure(to, kcf::KernelComputedField) = Adapt.adapt(to, kcf.data)


text:
function test_sqg_problemtype(dev, T)
  prob = SurfaceQG.Problem(dev; T=T)

  A = ArrayType(dev)

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

function test_sqg_stochasticforcedproblemconstructor(dev::Device=CPU())

  function calcF!(Fqh, sol, t, clock, vars, params, grid)
    Fqh .= Ffh
    return nothing
  end

  prob = SurfaceQG.Problem(dev; calcF=calcF!, stochastic=true)

  return typeof(prob.vars.prevsol) == typeof(prob.sol)
end

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

Tests the advection term in the SurfaceQG module by timestepping a
test problem with timestep dt and timestepper identified by the string stepper.
The test problem is derived by picking a solution bf (with associated
streamfunction ψf) for which the advection term J(ψf, bf) is non-zero. Next, a
forcing Ff is derived according to Ff = ∂bf/∂t + J(ψf, bf) - ν∇²bf. One solution
to the buoyancy equation forced by this Ff is then bf. (This solution may not
be realized, at least at long times, if it is unstable.)
"""
function test_sqg_advection(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=1e-2, nν=1)
  n, L  = 128, 2π
  ν, nν = 1e-2, 1
  tf = 0.5  # SQG piles up energy at small scales so running for t ⪆ 0.5 brings instability
  nt = round(Int, tf/dt)

  grid = TwoDGrid(dev, n, L)
  x, y = gridpoints(grid)

  ψf = @.             sin(2x)*cos(2y) +            2sin(x)*cos(3y)
  bf = @. - sqrt(8) * sin(2x)*cos(2y) - sqrt(10) * 2sin(x)*cos(3y)

  Ff = @. (
    - ν*(16*sqrt(2) * sin(2x)*cos(2y) + 10*sqrt(10) * 2sin(x)*cos(3y) )
    - 4*sqrt(2)*(sqrt(5)-2) * (cos(x)*cos(3y)*sin(2x)*sin(2y)
      - 3cos(2x)*cos(2y)*sin(x)*sin(3y))
  )

  Ffh = rfft(Ff)

  # Forcing
  function calcF!(Fh, sol, t, clock, vars, params, grid)
    Fh .= Ffh
    return nothing
  end

  prob = SurfaceQG.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, dt=dt, stepper=stepper, calcF=calcF!, stochastic=false)

  SurfaceQG.set_b!(prob, bf)

  stepforward!(prob, nt)

  SurfaceQG.updatevars!(prob)

  isapprox(prob.vars.b, bf, rtol=rtol_surfaceqg)
end

function test_sqg_kineticenergy_buoyancyvariance(dev::Device=CPU())
  nx, Lx  = 128, 2π
  ny, Ly  = 128, 3π

  grid = TwoDGrid(dev, nx, Lx, ny, Ly)
  x, y = gridpoints(grid)

  k₀, l₀ = 2π/grid.Lx, 2π/grid.Ly # fundamental wavenumbers
  ψ₀ = @. sin(2k₀*x)*cos(2l₀*y) + 2sin(k₀*x)*cos(3l₀*y)
  b₀ = @. - sqrt(8) * sin(2k₀*x)*cos(2l₀*y) - sqrt(10) * 2sin(k₀*x)*cos(3l₀*y)

  kinetic_energy_calc = 6
  buoyancy_variance_calc = 12

  prob = SurfaceQG.Problem(dev; nx=nx, Lx=Lx, ny=ny, Ly=Ly, stepper="ForwardEuler")

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

  SurfaceQG.set_b!(prob, b₀)
  SurfaceQG.updatevars!(prob)

  kinetic_energy_b₀ = SurfaceQG.kinetic_energy(prob)
  buoyancy_variance_b₀ = SurfaceQG.buoyancy_variance(prob)

  params = SurfaceQG.Params(p.ν, p.nν)

  return (isapprox(kinetic_energy_b₀, kinetic_energy_calc, rtol=rtol_surfaceqg) &&
   isapprox(buoyancy_variance_b₀, buoyancy_variance_calc, rtol=rtol_surfaceqg))
end

function test_sqg_paramsconstructor(dev::Device=CPU())
  n, L = 128, 2π
  ν, nν = 1e-3, 4

  prob = SurfaceQG.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, stepper="ForwardEuler")

  params = SurfaceQG.Params(ν, nν)

  return prob.params == params
end

function test_sqg_noforcing(dev::Device=CPU())
  n, L = 16, 2π

  prob_unforced = SurfaceQG.Problem(dev; nx=n, Lx=L, stepper="ForwardEuler")

  SurfaceQG.addforcing!(prob_unforced.timestepper.N, prob_unforced.sol, prob_unforced.clock.t, prob_unforced.clock, prob_unforced.vars, prob_unforced.params, prob_unforced.grid)

  function calcF!(Fh, sol, t, clock, vars, params, grid)
    Fh .= 2 * ArrayType(dev)(ones(size(sol)))
    return nothing
  end

  prob_forced = SurfaceQG.Problem(dev; nx=n, Lx=L, stepper="ForwardEuler", calcF=calcF!)

  SurfaceQG.addforcing!(prob_forced.timestepper.N, prob_forced.sol, prob_forced.clock.t, prob_forced.clock, prob_forced.vars, prob_forced.params, prob_forced.grid)

  return prob_unforced.timestepper.N == Complex.(ArrayType(dev)(zeros(size(prob_unforced.sol)))) && prob_forced.timestepper.N == Complex.(2*ArrayType(dev)(ones(size(prob_unforced.sol))))
end

function test_sqg_deterministicforcing_buoyancy_variance_budget(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, tf=10.0)
  n, L  = 256, 2π
  ν, nν = 1e-7, 2
  dt, tf = 0.005, 10
  nt = round(Int, tf/dt)

  grid  = TwoDGrid(dev, n, L)
  x, y = gridpoints(grid)

  # buoyancy forcing = 0.01cos(4x)cos(5y)cos(2t)
  f = @. 0.01cos(4x) * cos(5y)
  fh = rfft(f)
  function calcF!(Fh, sol, t, clock, vars, params, grid)
    @. Fh = fh * cos(2t)
    return nothing
  end

  prob = SurfaceQG.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, dt=dt, stepper="RK4",
    calcF=calcF!, stochastic=false)

  SurfaceQG.set_b!(prob, 0*x)

  B = Diagnostic(SurfaceQG.buoyancy_variance,    prob, nsteps=nt)
  D = Diagnostic(SurfaceQG.buoyancy_dissipation, prob, nsteps=nt)
  W = Diagnostic(SurfaceQG.buoyancy_work,        prob, nsteps=nt)
  diags = [B, D, W]

  stepforward!(prob, diags, nt)

  SurfaceQG.updatevars!(prob)

  dBdt_numerical = (B[3:B.i] - B[1:B.i-2]) / (2 * prob.clock.dt)
  dBdt_computed  = W[2:B.i-1] - D[2:B.i-1]

  return isapprox(dBdt_numerical, dBdt_computed, rtol = 1e-4)
end

function test_sqg_stochasticforcing_buoyancy_variance_budget(dev::Device=CPU(); n=256, L=2π, dt=0.005, ν=1e-7, nν=2, tf=2.0)
  nt = round(Int, tf/dt)

  # Forcing parameters
  kf, dkf = 12.0, 2.0
  εᵇ = 0.01

  grid = TwoDGrid(dev, n, L)
  x, y = gridpoints(grid)

  Kr = ArrayType(dev)([CUDA.@allowscalar grid.kr[i] for i=1:grid.nkr, j=1:grid.nl])

  forcing_spectrum = ArrayType(dev)(zero(grid.Krsq))
  @. forcing_spectrum = exp(-(sqrt(grid.Krsq) - kf)^2 / (2 * dkf^2))
  @. forcing_spectrum = ifelse(grid.Krsq < 2^2, 0, forcing_spectrum)
  @. forcing_spectrum = ifelse(grid.Krsq > 20^2, 0, forcing_spectrum)
  @. forcing_spectrum = ifelse(Kr < 2π/L, 0, forcing_spectrum)
  εᵇ0 = parsevalsum(forcing_spectrum, grid) / (grid.Lx * grid.Ly)
  forcing_spectrum .= εᵇ / εᵇ0 * forcing_spectrum

  Random.seed!(1234)

  function calcF!(Fh, sol, t, clock, vars, params, grid)
    eta = ArrayType(dev)(exp.(2π * im * rand(Float64, size(sol))) / sqrt(clock.dt))
    CUDA.@allowscalar eta[1, 1] = 0.0
    @. Fh = eta * sqrt(forcing_spectrum)
    nothing
  end

  prob = SurfaceQG.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, dt=dt, stepper="RK4",
    calcF=calcF!, stochastic=true)

  SurfaceQG.set_b!(prob, 0*x)

  B = Diagnostic(SurfaceQG.buoyancy_variance,    prob, nsteps=nt)
  D = Diagnostic(SurfaceQG.buoyancy_dissipation, prob, nsteps=nt)
  W = Diagnostic(SurfaceQG.buoyancy_work,        prob, nsteps=nt)
  diags = [B, D, W]

  stepforward!(prob, diags, nt)

  SurfaceQG.updatevars!(prob)

  dBdt_numerical = (B[2:B.i] - B[1:B.i-1]) / prob.clock.dt

  dBdt_computed = W[2:B.i] - D[1:B.i-1]

  return isapprox(dBdt_numerical, dBdt_computed, rtol = 1e-4)
end


text:
# To create all plots execute this file from within this folder
# "julia create_plots.jl"

include("plot_baseline_empirical.jl")
include("plot_policy_comparison.jl")
include("plot_policy_comparison_dynamics.jl")
include("plot_variance_single_runs.jl")


text:
macro check_pValues(pValues)
    :(if any(($pValues .< 0.0) | ($pValues .> 1.0))
      throw(DomainError())
    end)
end

macro step_up_adjustment(multiplier, sortedPValues)
  # applies the step up adjustment as e.g. described here:
  # http://www.unc.edu/courses/2007spring/biol/145/001/docs/lectures/Nov12.html
  # but a bit more efficiently as the array of pValues is ordered
  :(
    for i in 0:(n - 1)
      if i == 0
        $sortedPValues[n - i] *= $multiplier
      else
        $sortedPValues[n - i] = min($sortedPValues[n - i + 1], $sortedPValues[n - i] * $multiplier)
      end
    end
  )
end

macro step_down_adjustment(multiplier, sortedPValues)
  # applies the step down adjustment as e.g. described here:
  # http://en.wikipedia.org/w/index.php?title=Holm%E2%80%93Bonferroni_method&oldid=615389617#Adjusted_P-value
  # but a bit more efficiently as the array of pValues is ordered
  :(
    for i in 1:n
      if i == 1
        $sortedPValues[i] *= $multiplier
      else
        $sortedPValues[i] = max(sortedPValues[i - 1], sortedPValues[i] * $multiplier)
      end
    end
  )
end

function order(vector)
  indexes = collect(1:length(vector))
  sortedIndexes = sort(indexes, by = (i) -> vector[i])
  sortedIndexes
end