text:
using CSyntax
using Dates
using GLFW
using ModernGL

const DEFAULT_R = 0.8
const FRAGMENT_SHADER_BW = "fragment_shader_bw.glsl"
const FRAGMENT_SHADER_COLORFUL = "fragment_shader_colorful.glsl"
const VERTEX_SHADER = "vertex_shader.glsl"

function glsetup()
    # set up OpenGL context version
    # it seems OpenGL 4.1 is the highest version supported by MacOS.
    @static if Sys.isapple()
        VERSION_MAJOR = 4
        VERSION_MINOR = 1
    end

    @static if Sys.isapple()
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, VERSION_MAJOR)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, VERSION_MINOR)
        GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
        GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)
    else
        GLFW.DefaultWindowHints()
    end
end

function get_radius(args)
    if length(args) == 0
        return DEFAULT_R
    end
    try
        r = parse(Float64, args[1])
        return if r > 0 r else DEFAULT_R end
    catch
        return DEFAULT_R
    end
end

function get_fragment_shader(args)
    if length(args) >= 2 && args[2] == "bw"
        return FRAGMENT_SHADER_BW
    end
    return FRAGMENT_SHADER_COLORFUL
end

function create_shader_program(fs_source, vs_source)
    vs = glCreateShader(GL_VERTEX_SHADER)
    glShaderSource(vs, 1, Ptr{GLchar}[pointer(vs_source)], C_NULL)
    glCompileShader(vs)
    fs = glCreateShader(GL_FRAGMENT_SHADER)
    glShaderSource(fs, 1, Ptr{GLchar}[pointer(fs_source)], C_NULL)
    glCompileShader(fs)
    program = glCreateProgram()
    glAttachShader(program, vs)
    glAttachShader(program, fs)
    glLinkProgram(program)
    return program
end

mutable struct Controls
    right::Bool
    up::Bool
    left::Bool
    down::Bool
    zoomin::Bool
    zoomout::Bool
    paused::Bool
    backward::Bool
    forward::Bool
    max_iter_increase::Bool
    max_iter_decrease::Bool
    print::Bool
    quit::Bool
end

Controls() = Controls(false, false, false, false, false, false, false, false,
    false, false, false, false, false)

function make_key_listener(controls)
    function key_listener(window, key, a, action, b)
        if action == GLFW.PRESS
            if key == GLFW.KEY_I
                controls.zoomin = true
            elseif key == GLFW.KEY_O
                controls.zoomout = true
            elseif key == GLFW.KEY_D
                controls.right = true
            elseif key == GLFW.KEY_W
                controls.up = true
            elseif key == GLFW.KEY_A
                controls.left = true
            elseif key == GLFW.KEY_S
                controls.down = true
            elseif key == GLFW.KEY_SPACE
                controls.paused = !controls.paused
            elseif key == GLFW.KEY_LEFT
                controls.backward = true
            elseif key == GLFW.KEY_RIGHT
                controls.forward = true
            elseif key == GLFW.KEY_UP
                controls.max_iter_increase = true
            elseif key == GLFW.KEY_DOWN
                controls.max_iter_decrease = true
            elseif key == GLFW.KEY_P
                controls.print = true
            elseif key == GLFW.KEY_Q
                controls.quit = true
            end
        elseif action == GLFW.RELEASE
            if key == GLFW.KEY_I
                controls.zoomin = false
            elseif key == GLFW.KEY_O
                controls.zoomout = false
            elseif key == GLFW.KEY_D
                controls.right = false
            elseif key == GLFW.KEY_W
                controls.up = false
            elseif key == GLFW.KEY_A
                controls.left = false
            elseif key == GLFW.KEY_S
                controls.down = false
            elseif key == GLFW.KEY_LEFT
                controls.backward = false
            elseif key == GLFW.KEY_RIGHT
                controls.forward = false
            elseif key == GLFW.KEY_UP
                controls.max_iter_increase = false
            elseif key == GLFW.KEY_DOWN
                controls.max_iter_decrease = false
            end
        end
    end
end

function main()

    glsetup()

    window = GLFW.CreateWindow(720, 720, "Julia Set Explorer")
    @assert window != C_NULL "could not open window with GLFW3."
    GLFW.MakeContextCurrent(window)

    glEnable(GL_DEPTH_TEST)
    glDepthFunc(GL_LESS)

    program = create_shader_program(read(get_fragment_shader(ARGS), String),
        read(VERTEX_SHADER, String))

    screen_rect = GLfloat[
        -1.0, -1.0, 0.0,
        1.0, -1.0, 0.0,
        1.0, 1.0, 0.0,
        -1.0, -1.0, 0.0,
        1.0, 1.0, 0.0,
        -1.0, 1.0, 0.0
    ]

    vbo = GLuint(0)
    @c glGenBuffers(1, &vbo)
    glBindBuffer(GL_ARRAY_BUFFER, vbo)
    glBufferData(GL_ARRAY_BUFFER, sizeof(screen_rect), screen_rect,
        GL_STATIC_DRAW)

    vao = GLuint(0)
    @c glGenVertexArrays(1, &vao)
    glBindVertexArray(vao)
    glBindBuffer(GL_ARRAY_BUFFER, vbo)
    glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, C_NULL)
    glEnableVertexAttribArray(0)

    glClearColor(0.0, 0.0, 0.0, 1.0)

    # angular frequency
    w = 0.1 * 2 * pi
    r = get_radius(ARGS)

    controls = Controls()
    GLFW.SetKeyCallback(window, make_key_listener(controls))
    t0 = now()
    tm = t0
    c = r * exp(0im)
    center = 0 + 0im
    centerdot = 2.0
    scale = 0.5
    scalebase = 3.0
    max_iter = 200.0
    max_iter_dot = 20.0
    while !GLFW.WindowShouldClose(window)

        if controls.quit
            break
        end

        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

        glUseProgram(program)

        tn = now()
        t = (tn - t0).value / 1000.0
        dt = (tn - tm).value / 1000.0
        if !controls.paused && !controls.backward && !controls.forward
            c = c * exp(im * w * dt)
        end
        if controls.forward
            c = c * exp(2 * im * w * dt)
        end
        if controls.backward
            c = c * exp(-2 * im * w * dt)
        end
        if controls.zoomin
            scale += scale * log(scalebase) * dt
        end
        if controls.zoomout
            scale -= scale * log(scalebase) * dt
        end
        if controls.right
            center += centerdot * dt / scale
        end
        if controls.up
            center += im * centerdot * dt / scale
        end
        if controls.left
            center -= centerdot * dt / scale
        end
        if controls.down
            center -= im * centerdot * dt / scale
        end
        if controls.max_iter_increase
            max_iter += max_iter_dot * dt
        end
        if controls.max_iter_decrease
            max_iter -= max_iter_dot * dt
        end
        if controls.print
            println("c = $c")
            println("Showing the section between " *
                "$(round(center - (1 + 1im)/scale, digits=5)) and " *
                "$(round(center + (1 + 1im)/scale, digits=5))")
            controls.print = false
        end
        tm = tn

        loc = glGetUniformLocation(program, "max_iter")
        glUniform1i(loc, floor(max_iter))
        loc = glGetUniformLocation(program, "c")
        glUniform2f(loc, real(c), imag(c))
        loc = glGetUniformLocation(program, "center")
        glUniform2f(loc, real(center), imag(center))
        loc = glGetUniformLocation(program, "scale")
        glUniform1f(loc, scale)

        glBindVertexArray(vao)
        glDrawArrays(GL_TRIANGLES, 0, 6)

        GLFW.PollEvents()
        GLFW.SwapBuffers(window)
    end

    GLFW.DestroyWindow(window)

end

main()

ext: jl

max_stars_repo_name: antunmagdic/computer-graphics-lab

max_issues_repo_name: antunmagdic/computer-graphics-lab

max_forks_repo_name: antunmagdic/computer-graphics-lab


text:
mutable struct APWLOVars
    # energy step used for numerical calculation of energy derivatives
    deapwlo::Float64
    # maximum allowable APW order
    maxapword::Int64 # parameter=4
    # APW order
    apword::Array{Int64,2} #(0:maxlapw,maxspecies)
    # maximum of apword over all angular momenta and species
    apwordmax::Int64
    # total number of APW coefficients (l, m and order) for each species
    lmoapw::Vector{Int64} # (maxspecies)
    # polynomial order used for APW radial derivatives
    npapw::Int64
    # APW initial linearisation energies
    apwe0::Array{Float64,3} #(maxapword,0:maxlapw,maxspecies)
    # APW linearisation energies
    apwe::Array{Float64,3} # (:,:,:)
    # APW derivative order
    apwdm::Array{Int64,3} # (maxapword,0:maxlapw,maxspecies)
    # apwve is .true. if the linearisation energies are allowed to vary
    apwve::Array{Bool,3} # (maxapword,0:maxlapw,maxspecies)
    # APW radial functions
    apwfr::Array{Float64,5}
    # derivate of radial functions at the muffin-tin surface
    apwdfr::Array{Float64,3}
    # maximum number of local-orbitals
    maxlorb::Int64 # parameter=200
    # maximum allowable local-orbital order
    maxlorbord::Int64 # parameter=5
    # number of local-orbitals
    nlorb::Vector{Int64} #(maxspecies)
    # maximum nlorb over all species
    nlomax::Int64
    # total number of local-orbitals
    nlotot::Int64
    # local-orbital order
    lorbord::Array{Int64,2} #(maxlorb,maxspecies)
    # maximum lorbord over all species
    lorbordmax::Int64
    # polynomial order used for local-orbital radial derivatives
    nplorb::Int64
    # local-orbital angular momentum
    lorbl::Array{Float64,2} #(maxlorb,maxspecies)
    # maximum lorbl over all species
    lolmax::Int64
    # (lolmax+1)^2
    lolmmax::Int64
    # local-orbital initial energies
    lorbe0::Array{Float64,3} # (maxlorbord,maxlorb,maxspecies)
    # local-orbital energies
    lorbe::Array{Float64,3} #(:,:,:)
    # local-orbital derivative order
    lorbdm::Array{Int64,3} #(maxlorbord,maxlorb,maxspecies)
    # lorbve is .true. if the linearisation energies are allowed to vary
    lorbve::Array{Bool,3} # (maxlorbord,maxlorb,maxspecies)
    # local-orbital radial functions
    lofr::Array{Float64,4} # (:,:,:,:)
    # band energy search tolerance
    epsband::Float64
    # maximum allowed change in energy during band energy search; enforced only if
    # default energy is less than zero
    demaxbnd::Float64
    # minimum default linearisation energy over all APWs and local-orbitals
    e0min::Float64
    # if autolinengy is .true. then the fixed linearisation energies are set to the
    # Fermi energy minus dlefe
    autolinengy::Bool
    # difference between linearisation and Fermi energies when autolinengy is .true.
    dlefe::Float64
    # lorbcnd is .true. if conduction state local-orbitals should be added
    lorbcnd::Bool
    # conduction state local-orbital order
    lorbordc::Bool
    # excess order of the APW and local-orbital functions
    nxoapwlo::Int64
    # excess local orbitals
    nxlo::Int64
end

# effective size for 0:maxlapw -> 0:lmaxapw
# Using indices starting from 1
function APWLOVars(Nspecies::Int64, maxlapw::Int64)

    deapwlo = 0
    maxapword = 4
    apword = zeros(Int64,maxlapw+1,Nspecies)
    apwordmax = 0
    lmoapw = zeros(Int64,Nspecies)
    npapw = 0

    apwe0 = zeros(Float64, maxapword, maxlapw+1, Nspecies)
    apwe = zeros(Float64, 1, 1, 1)
    apwdm = zeros(Int64,maxapword, maxlapw+1, Nspecies)

    apwve = zeros(Bool, maxapword, maxlapw+1, Nspecies)
    apwfr = zeros(Float64,1,1,1,1,1)
    apwdfr = zeros(Float64,1,1,1)

    maxlorb = 200
    maxlorbord = 5
    nlorb = zeros(Int64,Nspecies)
    nlomax = 0
    nlotot = 0
    lorbord = zeros(Int64,maxlorb,Nspecies)
    lorbordmax = 0
    nplorb = 0

    lorbl = zeros(Float64,maxlorb,Nspecies)
    lolmax = 0
    lolmmax = 0
    lorbe0 = zeros(Float64,maxlorbord,maxlorb,Nspecies)
    lorbe = zeros(Float64,1,1,1)
    lorbdm = zeros(Int64,maxlorbord,maxlorb,Nspecies)
    lorbve = zeros(Bool,maxlorbord,maxlorb,Nspecies)
    lofr = zeros(Float64,1,1,1,1)
    epsband = 0.0
    demaxbnd = 0.0
    e0min = 0.0
    autolinengy = false
    dlefe = 0.0
    lorbcnd = false
    lorbordc = 0
    nxoapwlo = 0
    nxlo = 0

    return APWLOVars(
        deapwlo, maxapword, apword, apwordmax,
        lmoapw, npapw, apwe0, apwe, apwdm, apwve, apwfr, apwdfr,
        maxlorb, maxlorbord, nlorb, nlomax, nlotot, lorbord,
        lorbordmax, nplorb, lorbl, lolmax,
        lolmmax, lorbe0, lorbe, lorbdm, lorbve,
        lofr, epsband, demaxbnd, e0min, autolinengy,
        dlefe, lorbcnd, lorbordc, nxoapwlo, nxlo
    )

end

#apwlo_vars = APWLOVars(2,50)

# indexing for array allocated with lmaxo
# allocate array((-lmaxo-1:lmaxo+2))
# from lmaxo idx to the usual 1-based indexing
function lmaxo2idx(idx, lmaxo)
    return idx + lmaxo + 2
end


text:
struct DistributionOperator <: AbstractOperator
    domain
    target
    pindices
end

domain(op::DistributionOperator) = op.domain
target(op::DistributionOperator) = op.target

function apply(op::DistributionOperator, (val_full, val_distribute))
    res = similar(val_full)
    @avx for i in eachindex(val_full)
        res[i] = val_full[i] * val_distribute[op.pindices[i]]
    end
    return res
end

PowerDistributor(pdom::PowerDomain) = begin
    DistributionOperator((codomain(pdom), pdom), codomain(pdom), pdom._pindices)
end


text:
using FASTX
using Random
using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--proteins", "-p"
        required = true
        arg_type = String
        help = "Path to the FASTA/tsv file containing the protein sequences."
    "--peptides", "-e"
        required = true
        arg_type = String
        help = "Path to the FASTA/tsv file containing the peptide sequences."
    "--output", "-o"
        required = true
        arg_type = String
        help = "Path to the tsv file listing the interactions to score for PIPR."
end

args = parse_args(ARGS, s)

PROTEINS = args["proteins"]
PEPTIDES = args["peptides"]
OUTPUT = args["output"]
PROTEIN_FASTA = occursin("fasta", PROTEINS) ? true : false
PEPTIDES_FASTA = occursin("fasta", PEPTIDES) ? true : false

if PROTEIN_FASTA
    protein_ids = [FASTA.identifier(record) for record in FASTA.Reader(open(PROTEINS))]
else
    protein_ids = [split(line, "\t")[1] for line in readlines(PROTEINS)]
end

if PEPTIDES_FASTA
    peptide_ids = [FASTA.identifier(record) for record in FASTA.Reader(open(PEPTIDES))]
else
    peptide_ids = [split(line, "\t")[1] for line in readlines(PEPTIDES)]
end

open(OUTPUT, "w") do io
    write(io, "p1\tp2\tlabel\n")
	for peptide in peptide_ids
            for protein in protein_ids
		write(io, "$peptide\t$protein\t$(rand([0,1]))\n")
            end
	end
end

max_stars_repo_name: GreenCUBIC/PeptideScreening

max_issues_repo_name: GreenCUBIC/PeptideScreening

max_forks_repo_name: GreenCUBIC/PeptideScreening


text:
# Testing precomposition by nonlinear operator

b = randn(10)
g = SqrNormL2(3.0)
G = AffineAdd(AbstractOperators.Sigmoid((10,), 1.0),b,false)
f = StructuredOptimization.PrecomposeNonlinear(g, G)

x = randn(10)

grad_f_x, f_x = gradient(f, x)

@test size(grad_f_x) == size(x)
@test abs(f_x - 3.0/2 * norm(1.0 ./ (1.0 .+ exp.(-x)) - b)^2) <= 1e-10
expx = exp.(x)
expmx = 1.0./expx
grad_f_x_ref = 3.0 * ( expx ./ (1 .+ expx).^2 ) .* (1.0 ./ (1.0 .+ expmx) - b)
@test norm(grad_f_x - grad_f_x_ref) <= 1e-10

## with compose
#with vectors
l,m1,m2,n1,n2 = 2,3,4,5,6
x = ArrayPartition(randn(m1,m2),randn(n1,n2))
A = MatrixOp(randn(l,m1),m2)
B = MatrixOp(randn(m2,n1),n2)
r = randn(l,n2)

b = randn(l,n2)
G = AffineAdd(Ax_mul_Bx(
              HCAT(A,Zeros(codomainType(B), size(B,2), size(A,1) )),
              HCAT(Zeros(codomainType(A), size(A,2), size(B,1) ),B)
                       ),
              b,false)

g = SqrNormL2(3.0)
f = StructuredOptimization.PrecomposeNonlinear(g, G)

x = ArrayPartition(randn(m1,m2),randn(n1,n2))

grad_f_x, f_x = gradient(f, x)

r = G*x
grad_f_x2, f_x2 = gradient(g, r)
grad_f_x2 = jacobian(G,x)'*grad_f_x2

@test norm(f_x-f_x2) < 1e-8
@test norm(grad_f_x2.-grad_f_x2) < 1e-8


text:
struct EvalDfs end
# TODO doc sync, not real DFS

function fire(eve::EvalDfs, d::Dag, src_nid, v)
    # simple approach: DFS style triggering
    # TODO doc async behaviour
    # TODO able to fire multiple concurrently
    targets = Base.get(d.links, src_nid, [])
    # println("-- trigger $src_nid -> $targets")
    @sync for (dst_nid, f) in targets
        # TODO ugly hack to register external callback
        if dst_nid == 0
            @async f(v)
            continue
        end

        # TODO store dst directly in targets?
        dst = d.nodes[dst_nid]

        @async begin
            tv = f(v)
            if !isnull(tv)
                fire(eve, d, dst.nid, get(tv))
            end
        end
    end
end


text:
using TaskHeritableStorage  # Load the module (overrides Task() constructor)

using Test

@testset "TaskHeritableStorage" begin
    include("TaskHeritableStorage.jl")

    # Examples
    include("ThreadSafeAccumulators.jl")
    include("Tracing.jl")
    # NOTE: MultiTest.jl currently modifies Test.@testset behavior, so keep it last:
    include("MultiTest.jl")
end


text:
using TomoForward
using XfromProjections

img = zeros(100, 100)
img[30:51, 40:81] .= 1

nangles = 90
detcount = Int(floor(size(img,1)*1.4))
proj_geom = ProjGeom(1.0, detcount, LinRange(0,pi,nangles+1)[1:nangles])

# test line projection model
A = fp_op_parallel2d_line(proj_geom, size(img, 1), size(img, 2))
@time p = A * vec(img);
p = reshape(Array(p), (:, detcount));
q = filter_proj(p)

@time bp = A' * vec(p)
fbp = A' * vec(q) .* (pi / nangles)

fbp_img = reshape(fbp, size(img))
bp_img = reshape(bp, size(img))

# test strip projection model
A_strip = fp_op_parallel2d_strip(proj_geom, size(img, 1), size(img, 2))
@time p = A_strip * vec(img);
p = reshape(Array(p), (:, detcount));
q = filter_proj(p)

p_ = vec(p)
q_ = vec(q)

bp_strip = A' * p_
@time fbp_strip = A' * q_ .* (pi / nangles)

fbp_img_strip = reshape(fbp_strip, size(img))
bp_img_strip = reshape(bp_strip, size(img))

using PyPlot

ax00 = plt.subplot2grid((2,3), (0,0))
ax01 = plt.subplot2grid((2,3), (0,1))
ax02 = plt.subplot2grid((2,3), (0,2))
ax10 = plt.subplot2grid((2,3), (1,0))
ax11 = plt.subplot2grid((2,3), (1,1))
ax12 = plt.subplot2grid((2,3), (1,2))
ax00.imshow(p); ax00.set_title("sinogram")
ax01.imshow(bp_img); ax01.set_title("bp")
ax02.imshow(bp_img_strip); ax02.set_title("bp_strip")
ax10.imshow(img); ax10.set_title("ground truth")
ax11.imshow(fbp_img); ax11.set_title("fbp")
ax12.imshow(fbp_img_strip); ax12.set_title("fbp_strip")
plt.show()

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
using SpecialFunctions, QuadGK, ForwardDiff

function fluence_spatial_freq(z, s, μa, μsp, l)

    #get numerical instabillity at large s
    D = @. 1/(3(μsp + μa))
    α = @. sqrt((D*s^2 + μa)/D)
    A = 0.493
    zb = 2*A*D[1]
    z0 = 1/(μsp[1] + μa[1])
    tolr = 1e-50

    ## equation 11
    if z < z0

        # have to account for exp blowing up to inf
        if exp(α[1]*z0 - α[1]*z) > 1e8
            ϕ = 0.0
        else
            ϕ = D[1]*α[1]*(1 + exp(-2*α[1]*(l - z))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l - z)))
            ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
            ϕ = ϕ*(1 + exp(-2*α[1]*(zb + z0)))/(2*D[1]*α[1])
            ϕ = ϕ*exp(α[1]*z0 - α[1]*z)
            ϕ = ϕ - ((exp(α[1]*(z0 - z))*(1 - exp(-2*α[1]*(z0 - z))))/(2*D[1]*α[1]))
        end

    ## equation 12
    elseif z > z0 && z < l
        ϕ = D[1]*α[1]*(1 + exp(-2*α[1]*(l - z))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l - z)))
        ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
        ϕ = ϕ*(1 + exp(-2*α[1]*(zb + z0)))/(2*D[1]*α[1])
        ϕ = ϕ*exp(α[1]*z0 - α[1]*z)

    #equation 13
    elseif z >= l

        ϕ = exp(α[1]*(z0 - l) + α[2]*(l - z))
        ϕ = ϕ*(1 - exp(-2*α[1]*(zb + z0)))
        ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) + D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
    end

    return ϕ
end

# compute equation 14 integrand
function compute_integral(ρ, z, s, μa, μsp, l)
    return fluence_spatial_freq(z, s, μa, μsp, l)*s*besselj0(ρ*s)
end

# numerically integrate equation 14
function fluence_steadystate(ρ, z, μa, μsp, l)
    return quadgk((s) -> compute_integral(ρ, z, s, μa, μsp, l), 0, Inf, rtol=1e-8)[1]/2π
 end

#compare semi-infinite fluence

function ss_compare(ρ, z, μa, μsp)

    D = 1/(3(μsp + μa))
    A = 0.493
    z0 = 1/(μsp + μa)
    zb = 2*A*D
    μe = sqrt(3*μa*(μa + μsp))

    ϕ = exp(-μe*sqrt((z - z0)^2 + ρ^2))/sqrt((z - z0)^2 + ρ^2)
    ϕ = ϕ - exp(-μe*sqrt((z + z0 + 2*zb)^2 + ρ^2))/sqrt((z + z0 + 2*zb)^2 + ρ^2)
    ϕ = ϕ/(4*π*D)
end

# Eqn 16 to compute reflectance (use auto diff)
function refl_steadystate(ρ, μa, μsp, l)

    f(z) = fluence_steadystate(ρ, z, μa, μsp, l)
    g(z) = ForwardDiff.derivative(z -> f(z), z)

    return 0.118*fluence_steadystate(ρ, 0, μa, μsp, l) + 0.306/(3*(μa[1] + μsp[1]))*g(0)
end

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
using Base.Test
import ArchGDAL; const AG = ArchGDAL

@testset "Test methods for dataset" begin
    AG.registerdrivers() do
        AG.read("data/utmsmall.tif") do dataset
            AG.createcopy(dataset, "tmp/utmcopy.tif","GTiff") do copydataset
                @test AG.ngcp(copydataset) == 0
                @test AG.noverview(AG.getband(copydataset,1)) == 0
                AG.buildoverviews!(copydataset, Cint[2,4,8])
                @test AG.noverview(AG.getband(copydataset,1)) == 3
                AG.copywholeraster(dataset, copydataset,
                                   progressfunc=GDAL.C.GDALTermProgress)
            end
            AG.copyfiles("GTiff", "tmp/utmcopy2.tif", "tmp/utmcopy.tif")
            AG.update("tmp/utmcopy2.tif") do copydataset
                AG.copywholeraster(dataset, copydataset, ["COMPRESS=LZW"])
            end
        end
    end
    rm("tmp/utmcopy.tif")
    rm("tmp/utmcopy2.tif")
end

# untested: AG.deletelayer!(copydataset, 0)


text:
# TODO: Add ci()
# TODO: Add one-sided hypothesis tests

immutable BinomialTest <: HypothesisTest
    test_name::ASCIIString
    # Number of Successes
    statistic::Int
    # Test parameter: p
    parameter::Float64
    p_value::Float64
    confidence_interval::Vector{Float64}
    estimate::Float64
    null_value::Float64
    alternative_hypothesis::ASCIIString
    method::ASCIIString
end

function p_value(::Type{BinomialTest},
                 x::Integer,
                 n::Integer,
                 p::Real)
    d = Binomial(n, p)
    expectation = mean(d)

    if x == expectation
        return 1.0
    end

    if (x != 0 && p == 0.0) || (x != n && p == 1.0)
        return 0.0
    end

    mirror = expectation + (expectation - x) - 1.0

    pv = pmf(d, x)

    if x < expectation
        # TODO: Perform a binary search here
        if pmf(d, ifloor(mirror)) > pv
            lower_bound = ifloor(mirror)
        else
            lower_bound = iceil(expectation)
        end
        for i in lower_bound:n
            if pmf(d, i) <= pv
                return cdf(d, x) + ccdf(d, i - 1)
            end
        end
        return cdf(d, x)
    else
        # TODO: Perform a binary search here
        if pmf(d, iceil(mirror)) > pv
            upper_bound = iceil(mirror)
        else
            upper_bound = ifloor(expectation)
        end
        for i in upper_bound:-1:0
            if pmf(d, i) <= pv
                return cdf(d, i) + ccdf(d, x - 1)
            end
        end
        return ccdf(d, x - 1)
    end
end

function p_value(::Type{BinomialTest},
                 s::Real,
                 n::Real,
                 p::Real)
    p_value(BinomialTest, int(s), int(n), p)
end

function p_value{T <: Real}(::Type{BinomialTest},
                            x::Vector{T},
                            p::Real)
    p_value(BinomialTest, int(sum(x)), length(x), p)
end

function Base.run(::Type{BinomialTest},
                  s::Int,
                  n::Int;
                  p::Float64 = 0.5)
    return BinomialTest("Binomial Test",
                        s,
                        p,
                        p_value(BinomialTest, s, n, p),
                        [NaN, NaN],
                        s / n,
                        p,
                        "Two-sided",
                        "Exact binomial test")
end

function Base.run(::Type{BinomialTest},
                  s::Real,
                  n::Real;
                  p::Real = 0.5)
    Base.run(BinomialTest, int(s), int(n), p = float64(p))
end

function Base.run{T <: Real}(::Type{BinomialTest},
                             x::Vector{T};
                             p::Real = 0.5)
    Base.run(BinomialTest, int(sum(x)), length(x), p = float64(p))
end


text:
@testset "mle" begin
    # gen data
    ρ = SqueezedState(0.8, π/4, Matrix, dim=100)
    d = GaussianStateBHD(ρ)
    n = 8192
    data = rand(d, n)

    # mle
    t = 50
    ppit = PositivePreservingIterator(data, t, dim=35)
    run!(ppit)
end


text:
# run this file using "include" after adding Plots to the default Julia environment
using GlobalEnergyGIS, Plots, Plots.PlotMeasures

function plotlines(models, startday; variables=defaultvariables, era_year=2018)
    normdemand, normdemand_predicted, regionlist = getcvdata(models; variables=defaultvariables, era_year=era_year)
    regionlist[5] = "Bosnia & Herz."
    plot(normdemand, linecolor = "#0000FFFF", size=(1800,1350), lw=4, layout = (4,11), label="", xticks=[0,4000,8000],
    		title=reshape(regionlist, (1,44)), ylims=(0.5,1.7), titlefont=(13,:Calibri), tickfont=(11,:Calibri), bottom_margin = 10px)
    display(plot!(normdemand_predicted, linecolor = "#FF00007F", lw=4, label=""))
    # alternative colors "#00FC" and "#D008"
    plot(normdemand[startday*24 .+ (1:24*7), :], linecolor = "#0000FFFF", size=(1800,1350), layout = (4,11), label="",
    		title=reshape(regionlist, (1,44)), ylims=(0.5,1.7), titlefont=(13,:Calibri), tickfont=(11,:Calibri), bottom_margin = 10px)
    plot!(normdemand_predicted[startday*24 .+ (1:24*7), :], linecolor = "#FF00007F", label="")
end

plotly()

models = crossvalidate(variables=defaultvariables, nrounds=1000, max_depth=7, eta=0.005, subsample=0.05, metrics=["mae"])
plotlines(models, 123)

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module Modia3D.Composition

Structuring of objects moving in 3D. Most important constructors
(dof are the `degrees-of-freedom`):

| Function                                    | dof  | Description                                  |
|:--------------------------------------------|:----:|:---------------------------------------------|
| [`@assembly`](@ref)` Name(..) begin .. end` |   -  | Return a `Name` constructor for Object3Ds    |
| [`Object3D`](@ref)`([data];..)`             |   0  | Return a reference Object3D                  |
| [`Object3D`](@ref)`(parent [, data];..)`    |  0,6 | Return Object3D fixed/moving w.r.t. `parent` |
| [`Modia3D.Revolute`](@ref)`(obj1,obj2;..)`  |   1  | Return a revolute joint                      |
| [`Modia3D.Prismatic`](@ref)`(obj1,obj2;..)` |   1  | Return a prismatic joint                     |

The optional `data` associated with an `Object3D`
can be one of the following:

| `data`                            | Description                                          |
|:----------------------------------|:-----------------------------------------------------|
| `::`[`Modia3D.Solid`](@ref)       | Solids with geometry, mass, visual/contact material  |
| `<:Modia3D.AbstractVisualElement` | Visual elements ([`Modia3D.Graphics`](@ref))         |

# Main developers
Andrea Neumayr and Martin Otter, [DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
"""
module Composition

export @assembly

export initialize, initAnalysis!, performAnalysis!, closeAnalysis!, visualize!, visualizeAssembly!, visualizeWorld!
export updatePosition!, update!, driveJoint!
export Object3D, AssemblyInternal, initAssemblyInternal!, initAssemblyComponent!, Scene
export Part
export SimulationModel, Model
export BodyWithTwoFrames, ContactBox, ContactBox2, ContactBox3
export printObject3DAndAllChildren, writeObject3DAndAllChildrenOnJsonFile
export set_r!, set_q!, set_r_abs!
export Fixed
export Revolute, setAngle!, connect, addRevoluteTorqueObject, computeTorque
export Prismatic, setDistance!
export FreeMotion
export getModelResidues!, model_f1!, model_f2!
export distanceAndAngles, distance, planarRotationAngle

export Flange, RevoluteFlange

export SignalToFlangeAngle, SignalToFlangeTorque
export computeSignal

export AdaptorForceElementToFlange, AdaptorForceElementToPFlange

export setVariable_phi!, setVariable_r_rel!, setVariable_R_rel!

export SceneOptions, DLR_Visualization, ContactDetectionMPR
export NoGravityField, UniformGravityField, PointGravityField, gravityAcceleration
export G, EarthMass, EarthRadius

export DLR_Visualization_renderer

export ContactDetectionMPR_handler
export initializeContactDetection!, selectContactPairsNoEvent!, selectContactPairsWithEvent!, getDistances!, setComputationFlag, closeContactDetection!

export responseCalculation, print_ModelVariables
export contactStart, contactEnd, regularize, resultantCoefficientOfRestitution, resultantDampingCoefficient
export ElasticContactPairResponseMaterial, elasticContactPairCoefficients, normalRelativeVelocityAtContact

using StaticArrays
using DataStructures

import Base
import ModiaMath
import Modia3D
import Modia3D.Basics
import Modia3D.Solids
import Modia3D.Graphics
import Modia3D.Signals
import JSON
import Printf

const NOTHING = Nothing

const MaterialOrNothing                = Union{Graphics.Material,NOTHING}
const MassPropertiesOrNothing          = Union{Solids.MassProperties,NOTHING}
const AbstractContactMaterialOrNothing = Union{Modia3D.AbstractContactMaterial,NOTHING}

include("flange.jl")
include("object3D.jl")
include(joinpath("responseCalculation", "elasticCollisionResponse.jl"))
include(joinpath("joints", "FixedJoint.jl"))
include(joinpath("joints", "FreeMotion.jl"))
include(joinpath("superObjects.jl"))
include("scene.jl")
include("assembly.jl")
include("sensors.jl")
include("joints.jl")
include(joinpath("joints", "Fixed.jl"))
include(joinpath("joints", "Revolute.jl"))
include(joinpath("joints", "Prismatic.jl"))
include("adaptors.jl")
include("assemblies.jl")
include("assignObjects.jl")
include("massComputation.jl")
include("handler.jl")
include("responseCalculation.jl")
include("dynamics.jl")

end


text:
using Distributed
@everywhere using FileTrees, DataFrames, CSV, Test

@testset "taxi tests" begin
    dir = joinpath(@__DIR__, "..", "page", "taxi-data")
    global taxi_dir = FileTree(dir)
    @test string(name(taxi_dir)) == dir
    @test string(path(taxi_dir["2019"])) == joinpath(dir, "2019")

    @test string(path(taxi_dir["2019"])) == joinpath(dir, "2019")
    @test name(FileTrees.rename(taxi_dir, "taxi-data")) == "taxi-data"
end

@testset "loading" begin
    global dfs = FileTrees.load(taxi_dir) do file
        DataFrame(CSV.File(path(file)))
    end

    global lazy_dfs = FileTrees.load(taxi_dir, lazy=true) do file
        DataFrame(CSV.File(path(file)))
    end

    @testset "count files" begin
        @test reducevalues(+, mapvalues(x -> (@test x isa DataFrame; 1), dfs)) == 8
        @test reducevalues(+, mapvalues(x -> (@test x isa DataFrame; 1), dfs, lazy=true)) isa FileTrees.Thunk
        @test exec(reducevalues(+, mapvalues(x -> (@test x isa DataFrame; 1), dfs, lazy=true))) == 8

        @test reducevalues(+, mapvalues(x -> (@test x isa DataFrame; 1), lazy_dfs)) isa FileTrees.Thunk
        @test exec(reducevalues(+, mapvalues(x -> 1, lazy_dfs))) == 8
        # lazy=false
        @test reducevalues(+, mapvalues(x -> x isa FileTrees.Thunk,
                                        lazy_dfs, lazy=false)) == 8
    end

end

@testset "mv" begin
    t2 = mv(dfs, r"^([^/]*)/([^/]*)/yellow.csv$", s"yellow/\1/\2.csv")

    # move */*/green.csv to green/*/*.csv
    t2 = mv(t2, r"^([^/]*)/([^/]*)/green.csv$", s"green/\1/\2.csv")

    # All files must get moved
    @test isempty(setdiff(name.(children(t2)), ["yellow", "green"]))

    # this should throw
    @test_throws ErrorException mv(dfs, r"^([^/]*)/([^/]*)/yellow.csv$", s"yellow.csv")

    yellow = mv(dfs, r"^([^/]*)/([^/]*)/yellow.csv$", s"yellow.csv", combine=vcat)["yellow.csv"]
    @test get(yellow) isa DataFrame

    # correct?
    @test Set(DataFrames.Tables.rowtable(reducevalues(vcat, dfs[glob"*/*/yellow.csv"]))) ==
    Set(DataFrames.Tables.rowtable(get(yellow)))
end

@testset "metadata" begin
    using FileTrees, DataFrames, CSV

    # Lazy-load the data
    data = FileTrees.load(taxi_dir;lazy=true) do file
        DataFrame(CSV.File(path(file)))
    end

    metadata = mapvalues(data) do df
        (unique(df.RatecodeID)...,)
    end |> exec

    only_5 = filter(x->FileTrees.hasvalue(x) && (5 in get(x)), metadata, dirs=false)

    df_only_5 = data[only_5]

    reference = data[glob"2020/01/*.csv"]

    @test isempty(diff(df_only_5, reference))

    @test isequal(exec(reducevalues(tuple, reference)),
                  exec(reducevalues(tuple, df_only_5)))
end


text:
const MAX_SHOW_COMPONENTS = 10
const MAX_SHOW_FORECASTS = 10
const MAX_SHOW_FORECAST_INITIAL_TIMES = 1

function Base.summary(components::Components)
    return "$(typeof(components)): $(get_num_components(components))"
end

function Base.show(io::IO, components::Components)
    i = 1
    for component in iterate_components(components)
        if i <= MAX_SHOW_COMPONENTS
            show(io, component)
            println(io)
        end
        i += 1
    end

    if i > MAX_SHOW_COMPONENTS
        num = i - MAX_SHOW_COMPONENTS
        println(io, "\n***Omitted $num components***\n")
    end
end

function Base.show(io::IO, ::MIME"text/plain", components::Components)
    df = create_components_df(components)
    num_components = get_num_components(components)
    println(io, "Components")
    println(io, "==========")
    println(io, "Num components: $num_components")
    if num_components > 0
        println(io)
        show(io, df)
    end
end

function Base.show(io::IO, ::MIME"text/html", components::Components)
    df = create_components_df(components)
    num_components = get_num_components(components)
    println(io, "
Components
")
    println(io, "
Num components: $num_components
")
    if num_components > 0
        withenv("LINES" => 100, "COLUMNS" => 200) do
            show(io, MIME"text/html"(), df)
        end
    end
end

function Base.summary(container::TimeSeriesContainer)
    return "$(typeof(container)): $(length(container))"
end

function Base.show(io::IO, ::MIME"text/plain", container::TimeSeriesContainer)
    println(io, summary(container))
    for key in keys(container.data)
        println(io, "$(key.time_series_type): name=$(key.name)")
    end
end

function Base.summary(time_series::TimeSeriesData)
    return "$(typeof(time_series)) time_series ($(length(time_series)))"
end

function Base.summary(time_series::TimeSeriesMetadata)
    return "$(typeof(time_series)) time_series"
end

function Base.show(io::IO, data::SystemData)
    show(io, data.components)
    println(io, "\n")
    show(io, data.time_series_params)
end

function Base.show(io::IO, ::MIME"text/plain", data::SystemData)
    component_count, ts_count, forecast_count = get_time_series_counts(data)
    show(io, MIME"text/plain"(), data.components)
    println(io, "\n")

    println(io, "TimeSeriesContainer")
    println(io, "===================")
    println(io, "Components with time series data: $component_count")
    println(io, "Total StaticTimeSeries: $ts_count")
    println(io, "Total Forecasts: $forecast_count")
    if component_count == 0
        return
    end

    res = get_time_series_resolution(data)
    res = res <= Dates.Minute(1) ? Dates.Second(res) : Dates.Minute(res)
    println(io, "Resolution: $res")
    if forecast_count > 0
        initial_times = get_forecast_initial_times(data)
        println(io, "First initial time: $(first(initial_times))")
        println(io, "Last initial time: $(last(initial_times))")
        println(io, "Horizon: $(get_forecast_horizon(data))")
        println(io, "Interval: $(Dates.Minute(get_forecast_interval(data)))")
        println(io, "Forecast window count: $(get_forecast_window_count(data))")
    end
end

function Base.show(io::IO, ::MIME"text/html", data::SystemData)
    show(io, MIME"text/html"(), data.components)
    println(io, "\n")

    component_count, ts_count, forecast_count = get_time_series_counts(data)
    println(io, "
TimeSeriesContainer
")
    println(io, "
Components with time series data: $component_count
")
    println(io, "
Total StaticTimeSeries: $ts_count
")
    println(io, "
Total Forecasts: $forecast_count
")

    res = get_time_series_resolution(data)
    res = res <= Dates.Minute(1) ? Dates.Second(res) : Dates.Minute(res)
    println(io, "
Resolution: $(res)
")
    if forecast_count > 0
        initial_times = get_forecast_initial_times(data)
        window_count = get_forecast_window_count(data)
        println(io, "
First initial time: $(first(initial_times))
")
        println(io, "
Last initial time: $(last(initial_times))
")
        println(io, "
Horizon: $(get_forecast_horizon(data))
")
        println(io, "
Interval: $(Dates.Minute(get_forecast_interval(data)))
")
        println(io, "
Forecast window count: $(window_count)
")
    end
end

function Base.summary(ist::InfrastructureSystemsComponent)
    # All InfrastructureSystemsComponent subtypes are supposed to implement get_name.
    # Some don't.  They need to override this function.
    return "$(get_name(ist)) ($(typeof(ist)))"
end

function Base.show(io::IO, ::MIME"text/plain", system_units::SystemUnitsSettings)
    print(io, summary(system_units), ":")
    for name in fieldnames(typeof(system_units))
        val = getfield(system_units, name)
        print(io, "\n      ", name, ": ", val)
    end
end

function Base.show(io::IO, ::MIME"text/plain", ist::InfrastructureSystemsComponent)
    print(io, summary(ist), ":")
    for (name, field_type) in zip(fieldnames(typeof(ist)), fieldtypes(typeof(ist)))
        obj = getfield(ist, name)
        if obj isa InfrastructureSystemsInternal
            continue
        elseif obj isa TimeSeriesContainer || obj isa InfrastructureSystemsType
            val = summary(getfield(ist, name))
        elseif obj isa Vector{<:InfrastructureSystemsComponent}
            val = summary(getfield(ist, name))
        else
            val = getfield(ist, name)
        end
        # Not allowed to print `nothing`
        if isnothing(val)
            val = "nothing"
        end
        print(io, "\n   ", name, ": ", val)
    end
end

function Base.show(io::IO, ist::InfrastructureSystemsComponent)
    print(io, string(nameof(typeof(ist))), "(")
    is_first = true
    for (name, field_type) in zip(fieldnames(typeof(ist)), fieldtypes(typeof(ist)))
        if field_type <: TimeSeriesContainer || field_type <: InfrastructureSystemsInternal
            continue
        else
            val = getfield(ist, name)
        end
        if is_first
            is_first = false
        else
            print(io, ", ")
        end
        print(io, val)
    end
    print(io, ")")
end

function Base.show(io::IO, ::MIME"text/plain", it::FlattenIteratorWrapper)
    println(io, "$(eltype(it)) Counts: ")
    for (ctype, count) in _get_type_counts(it)
        println(io, "$ctype: $count")
    end
end

function _get_type_counts(it::FlattenIteratorWrapper)
    data = SortedDict()
    for component in it
        ctype = string(typeof(component))
        if !haskey(data, ctype)
            data[ctype] = 1
        else
            data[ctype] += 1
        end
    end

    return data
end

function create_components_df(components::Components)
    counts = Dict{String, Int}()
    rows = []

    for (subtype, values) in components.data
        type_str = strip_module_name(string(subtype))
        counts[type_str] = length(values)
        parents = [strip_module_name(string(x)) for x in supertypes(subtype)]
        row = (
            ConcreteType = type_str,
            SuperTypes = join(parents, " <: "),
            Count = length(values),
        )
        push!(rows, row)
    end

    isempty(rows) && return DataFrames.DataFrame()

    sort!(rows, by = x -> x.ConcreteType)

    return DataFrames.DataFrame(rows)
end

## This function takes in a time period or date period and returns a compound period
function convert_compound_period(period::Union{Dates.TimePeriod, Dates.DatePeriod})
    period = time_period_conversion(period)

    milli_weeks = period - (period % Dates.Millisecond(604800000))
    weeks = convert(Dates.Week, milli_weeks)
    period -= milli_weeks

    milli_days = period - (period % Dates.Millisecond(86400000))
    days = convert(Dates.Day, milli_days)
    period -= milli_days

    milli_hours = period - (period % Dates.Millisecond(3600000))
    hours = convert(Dates.Hour, milli_hours)
    period -= milli_hours

    milli_minutes = period - (period % Dates.Millisecond(60000))
    minutes = convert(Dates.Minute, milli_minutes)
    period -= milli_minutes

    seconds = period - (period % Dates.Millisecond(1000)) # finding the seconds
    seconds = convert(Dates.Second, seconds)
    remainder = period % Dates.Millisecond(1000) #finding the remainding milliseconds
    total = weeks + days + hours + minutes + seconds + remainder
    return total
end


text:
include("CubicSolvers.jl")
include("HirotaSolvers.jl")
include("SSSolvers.jl")

"""
    solve!(sim::Sim)

Solves the object `sim::Sim`.

See also: [`Sim`](@ref)
"""
function solve!(sim::Sim)
    # Set initial condition
    sim.ψ[:, 1] = sim.ψ₀
    # Check for pruning and calculate indices
    ind_p = zeros(Int64, 1)
    if sim.β > 0
        ind_p = [i for i in 2:(sim.box.Nₜ÷2+1) if (i-1)%sim.box.n_periods != 0]
        ind_p = sort([ind_p; sim.box.Nₜ.-ind_p.+2])
        @debug "Computed pruning indices $ind_p"
    end

    ops = Operators(sim)

    @info "Starting evolution"
    # Define these tuples as global consts
    if sim.T̂ ∈ (T1A!, T2A!, T4A_TJ!, T6A_TJ!, T8A_TJ!, T4A_SF!, T4A_SF!, T6A_SF!, T8A_SF!, T4A_CMP!, T6A_CMP!, T8A_CMP!, T6A_Ss14!, T6A_Ys7!, T6A_KLs9!, T8A_Ss15!, T8A_KLs17!, T1A_H!, T2A_H!, T1A_SS!)
        soln_loop_A(sim, ops, ind_p)
    elseif sim.T̂ ∈ (T1B!, T2B!, T4B_TJ!, T6B_TJ!, T8B_TJ!, T4B_SF!, T6B_SF!, T8B_SF!, T4B_CMP!, T6B_CMP!, T8B_CMP!, T6B_Ss14!, T6B_Ys7!, T6B_KLs9!, T8B_Ss15!, T8B_KLs17!)
        soln_loop_B(sim, ops, ind_p)
    else
        throw(ArgumentError("Unknown algorithm $(sim.T̂)"))
    end

    @info "Computation Done!"
end #solve

"""
    soln_loop_A(sim, ops, ind_p))

Propagate in `x` for integrators of type `A` (i.e. starting with nonlinear step, thus in real space).

See also: [`Sim`](@ref), [`solve!`](@ref)
"""
function soln_loop_A(sim, ops, ind_p)
    for i = 1:sim.box.Nₓ-1
        @views sim.T̂(sim.ψ[:, i+1], sim.ψ[:, i], sim.box.dx, ops)
        # Pruning
        if sim.β > 0
            ops.F̂*view(sim.ψ,:,i+1)
            for j in ind_p
                sim.ψ[j, i+1] *= exp(-sim.β*abs(sim.ψ[j, i+1]))
            end
            ops.F̃̂*view(sim.ψ,:,i+1)
        end # if
    end # for
    @info "Computing Spectrum"
    sim.ψ̃ .= fftshift(fft(sim.ψ, 1), 1)./sim.box.Nₜ
end

"""
    soln_loop_B(sim, ops, ind_p))

Propagate in `x` for integrators of type `B` (i.e. starting with dispersion step, thus in frequency space).

See also: [`Sim`](@ref), [`solve!`](@ref)
"""
function soln_loop_B(sim, ops, ind_p)
    F̃̂ = plan_ifft(@view sim.ψ[:, 1])
    F̂ = plan_fft(@view sim.ψ[:, 1])
    sim.ψ̃[:, 1] .= F̂*view(sim.ψ, :, 1)
    for i = 1:sim.box.Nₓ-1
        @views sim.T̂(sim.ψ̃[:, i+1], sim.ψ̃[:, i], sim.box.dx, ops)
        # Pruning
        if sim.β > 0
            for j in ind_p
                sim.ψ̃[j, i+1] *= exp(-sim.β*abs(sim.ψ̃[j, i+1]))
            end
        end
        # Compute the actual ψ
        #sim.ψ[:, i+1] .= F̃̂*view(sim.ψ̃, :, i+1)
    end
    # Shift the FT
    @info "Computing Spectrum"
    sim.ψ .= ifft(sim.ψ̃, 1)
    sim.ψ̃ .= fftshift(sim.ψ̃, 1)./sim.box.Nₜ
end


text:
function BestHit(leaves::Dict{Taxon,Int})
    return last(findmax(leaves))
end

function freeLCA(leaves::Dict{Taxon,BlastResult}, minimal::Float64,ranks::Vector{Symbol},precision::Dict{Symbol, Float64})
    besthitscore = first(findmax([last(l).bitscore for l in leaves]))
    filter!(x -> last(x).bitscore > besthitscore*minimal, leaves)
    taxa = collect(keys(leaves))
    current_lca = topolgoy(taxa)
    corrected_taxon = cut_by_precision(current_lca, ranks, precision, leaves)
    return corrected_taxon
end

function weightedLCA(leaves::Dict{Taxon,BlastResult}, minimal::Float64, cutoff::Float64, ranks::Vector{Symbol},precision::Dict{Symbol, Float64})
    @assert cutoff > 0.5 && cutoff < 1
    besthitscore = first(findmax([last(l).bitscore for l in leaves]))
    filtered_leaves = filter(x -> last(x).bitscore > besthitscore*minimal, leaves)

    taxa = collect(keys(filtered_leaves))
    tree = topolgoy(taxa)

    total_bitscore = sum([last(l).bitscore for l in filtered_leaves])
    threshold_bitscore = cutoff * total_bitscore

    next = Stack{PhyloTree}()
    push!(next,tree)
    current_lca = tree
    while true
        current_node = pop!(next)
        sub_bitscore = sum([leaves[leave.node].bitscore for leave in Leaves(current_node)])
        if sub_bitscore >= threshold_bitscore
            current_lca = current_node
            children_tree =  current_node.children
            for ctree in children_tree; push!(next,ctree); end
        end
        if isempty(next)
            break
        end
    end
    corrected_taxon = cut_by_precision(current_lca, ranks, precision, leaves)
    return corrected_taxon
end

function cut_by_precision(current_lca::PhyloTree, ranks::Vector{Symbol}, precision::Dict{Symbol, Float64}, leaves::Dict{Taxon,BlastResult})
    max_sub_pident = maximum([leaves[leave.node].pident for leave in Leaves(current_lca)])
    lineage = Lineage(current_lca.node)
    for taxon in lineage
        r = rank(taxon)
        if ! (r in ranks) || !haskey(precision, r)
            continue
        end
        if precision[r] > max_sub_pident
            return taxon
        end
    end
    return lineage[end]
end


text:
function DecimateSym!(var::PlmVar, alg::PlmAlg, dec::DecVar{1})

    Jvec,pseudolike = MinimizePLSym(alg, var)
    Nc2 = binomial(var.N,2)
    blockdecimate = dec.blockdecimate
    LL = blockdecimate ? Nc2 : Nc2 * var.q2
    iteration = 1
    while true
        SetAllCouplingToZero!(Jvec,  var, dec)
        nzeros = sum(!,dec.dmask) # number of couplings to be set to zero
        fraczeros = nzeros / LL
        println("# set to zero ", nzeros ," over ", LL ," equal to ", 100.0*fraczeros," % of the couplings. Converging ...")
        pseudolike, Jvec = MinimizePLSym(alg, var, dec, Jvec, x0=Jvec)
        println("# Converged!")
#        push!(results,(iteration,1.0-fraczeros, fmin,fminloc))
        fraczeros >= 1.0 && break
        fraczeros >= dec.fracmax && break
        iteration += 1
    end
    return Jvec, pseudolike
end

function SetAllCouplingToZero!(vecJ::Array{Float64,1}, var::PlmVar, dec::DecVar{1})

    N = var.N
    q  = var.q
    q2 = var.q2
    blockdecimate = dec.blockdecimate

    if blockdecimate
        LL = binomial(N,2)
    else
        LL = length(vecJ)
        LL == (binomial(N,2)*q2 + N*q) || error("LL = $LL should be ", binomial(N,2)*q2 + N*q)
    end

    vecval = Tuple{Int,Float64}[]
    sizehint!(vecval,LL)
    if blockdecimate   # compute the FN of the block
        for i=1:LL
            if dec.dmask[i]
                fnvecJ = 0.0
                for j in ((i-1)*q2 + 1):i*q2
                    fnvecJ += vecJ[j]^2
                end
                push!(vecval,(i,fnvecJ))
            end
        end
    else
        for i=1:LL-N*q
            dec.dmask[i] && push!(vecval,(i, vecJ[i]^2))
        end
    end

    lvec = length(vecval)
    idxperm = sortperm(Float64[ vecval[i][2] for i=1:lvec])

    ndec = blockdecimate ? round(Int, dec.fracdec * (LL-N*q)) : round(Int, dec.fracdec*LL)

    for i=1:ndec
        site = vecval[idxperm[i]][1]
        dec.dmask[site] = false
    end
end

function MinimizePLSym(alg::PlmAlg, var::PlmVar, dec::DecVar{1}, Jvec::Array{Float64,1}; x0 = zeros(Float64,length(Jvec)))

    N  = var.N
    q  = var.q
    q2 = var.q2
    Nc2 = binomial(N,2)

    opt = Opt(alg.method, length(x0))
    ftol_abs!(opt, alg.epsconv)
    lb,ub = ComputeConstraintsSym!(dec, var.q, var.N, x0)
    lower_bounds!(opt, lb)
    upper_bounds!(opt, ub)

    maxeval!(opt, alg.maxit)
    min_objective!(opt, (x,g)->optimfunwrapper(x,g,var))
    elapstime = @elapsed  (minf, minx, ret) = optimize(opt, x0)

    alg.verbose && @printf("pl = %.4f\t time = %.4f\n",  minf, elapstime)
    return minf,minx
end

function ComputeConstraintsSym!(decvar::DecVar{1}, q::Int, N::Int, x0::Array{Float64,1} )

    q2 = q*q
    dmask = decvar.dmask
    blockdecimate = decvar.blockdecimate

    Nc2 = binomial(N,2)
    lb = -Inf * ones(Float64, Nc2*q2 + N*q)
    ub =  Inf * ones(Float64, Nc2*q2 + N*q)
    LL = length(x0)

    if blockdecimate
        for i=1:Nc2
            if dmask[i] == false
                for j in ((i-1)*q2 + 1):i*q2
                    lb[ j ] = -1e-12
                    ub[ j ] = 1e-12
                    x0[ j ] = 0.0
                end
            else
                for j in ((i-1)*q2 + 1):i*q2
                    x0[ j ] += (1.0 - 2.0*rand())*1e-3
                end
            end
        end
    else
        for i=1:LL-N*q
            if dmask[i] == false
                lb[ i ] = -1e-6
                ub[ i ] =  1e-6
                x0[ i ] =  0.0
            else
                x0[ i ] += (1.0 - 2.0*rand())*1e-3
            end
        end
    end
    return lb,ub
end

max_stars_repo_name: YoshitakaMo/PlmDCA

max_issues_repo_name: YoshitakaMo/PlmDCA

max_forks_repo_name: YoshitakaMo/PlmDCA


text:
module ForcedFlowFixedSlip

using Printf

using Oceananigans
using Oceananigans.BoundaryConditions
using Oceananigans.OutputWriters
using Oceananigans.Fields

# Functions that define the forced flow problem

 ξ(t) = 1 + sin(t^2)
ξ′(t) = 2t * cos(t^2)

 f(x, t) =   cos(x - ξ(t))
fₓ(x, t) = - sin(x - ξ(t))

 g(y) = y^3 - y^2
g′(y) = 3y^2 - 2y

H₁(y) = y^4 / 4 - y^3 / 3 - 3y^2 + 2y
H₂(y) = 3y^4 - 4y^3 + 2y^2
H₃(y) = y^4 / 4 - y^3 / 3 - 6y^2 + 4y

Fᵘ(x, y, t) = ξ′(t) * fₓ(x, t) * H₁(y) + f(x, t) * fₓ(x, t) * H₂(y) - f(x, t) * H₃(y)

Fᵛ(x, y, t) = 3y^5 - 5y^4 + 2y^3

u(x, y, t) = f(x, t) * g′(y)
v(x, y, t) = - fₓ(x, t) * g(y)

function print_progress(simulation)
    model = simulation.model
    i, t = model.clock.iteration, model.clock.time
    progress = 100 * (i / simulation.stop_iteration)
    @info @sprintf("[%05.2f%%] iteration: %d, time: %.5e", progress, i, t)
    return nothing
end

const DATA_DIR = joinpath(@__DIR__, "..", "data")

function setup_xy_simulation(; Nx, Δt, stop_iteration, architecture=CPU(), dir=DATA_DIR)

    u_forcing(x, y, z, t) = Fᵘ(x, y, t)
    v_forcing(x, y, z, t) = Fᵛ(x, y, t)

    grid = RegularRectilinearGrid(size=(Nx, Nx, 1), x=(0, 2π), y=(0, 1), z=(0, 1),
                                topology=(Periodic, Bounded, Bounded))

    # "Fixed slip" boundary conditions (eg, no-slip on south wall, finite slip on north wall)."
    u_bcs = UVelocityBoundaryConditions(grid, north = ValueBoundaryCondition((x, y, t) -> f(x, t)),
                                              south = ValueBoundaryCondition(0))

    model = IncompressibleModel(       architecture = CPU(),
                                               grid = grid,
                                           coriolis = nothing,
                                           buoyancy = nothing,
                                            tracers = nothing,
                                            closure = IsotropicDiffusivity(ν=1),
                                boundary_conditions = (u=u_bcs,),
                                            forcing = (u=u_forcing, v=v_forcing))

    set!(model, u = (x, y, z) -> u(x, y, 0),
                v = (x, y, z) -> v(x, y, 0))

    simulation = Simulation(model, Δt=Δt, stop_iteration=stop_iteration, progress=print_progress, iteration_interval=20)

    outputs = Dict()
    pressure_output = model -> parent(model.pressures.pHY′) .+ parent(model.pressures.pNHS)
    outputs = merge((p=pressure_output,), model.velocities)

    simulation.output_writers[:fields] = JLD2OutputWriter(model, outputs,
                                                          dir = dir, force = true,
                                                          prefix = @sprintf("forced_fixed_slip_xy_Nx%d_Δt%.1e", Nx, Δt),
                                                          field_slicer = nothing,
                                                          schedule = TimeInterval(stop_iteration * Δt / 2))

    return simulation
end

function setup_and_run_xy(args...; kwargs...)
    simulation = setup_xy_simulation(args...; kwargs...)
    @info "Running forced flow fixed slip simulation with Nx = Ny = $(kwargs[:Nx]), Δt = $(kwargs[:Δt])"
    @time run!(simulation)
    return simulation
end

function setup_xz_simulation(; Nx, Δt, stop_iteration, architecture=CPU(), dir=DATA_DIR)

    u_forcing(x, y, z, t) = Fᵘ(x, z, t)
    w_forcing(x, y, z, t) = Fᵛ(x, z, t)

    grid = RegularRectilinearGrid(size=(Nx, 1, Nx), x=(0, 2π), y=(0, 1), z=(0, 1),
                                topology=(Periodic, Bounded, Bounded))

    # "Fixed slip" boundary conditions (eg, no-slip on bottom and finite slip on top)."
    u_bcs = UVelocityBoundaryConditions(grid,    top = ValueBoundaryCondition((x, z, t) -> f(x, t)),
                                              bottom = ValueBoundaryCondition(0))

    model = IncompressibleModel(       architecture = CPU(),
                                               grid = grid,
                                           coriolis = nothing,
                                           buoyancy = nothing,
                                            tracers = nothing,
                                            closure = IsotropicDiffusivity(ν=1),
                                boundary_conditions = (u=u_bcs,),
                                            forcing = (u=u_forcing, w=w_forcing))

    set!(model, u = (x, y, z) -> u(x, z, 0),
                w = (x, y, z) -> v(x, z, 0))

    simulation = Simulation(model, Δt=Δt, stop_iteration=stop_iteration, progress=print_progress, iteration_interval=20)

    outputs = Dict()
    pressure_output = model -> parent(model.pressures.pHY′) .+ parent(model.pressures.pNHS)
    outputs = merge((p=pressure_output,), model.velocities)

    simulation.output_writers[:fields] = JLD2OutputWriter(model, outputs,
                                                          dir = dir, force = true,
                                                          prefix = @sprintf("forced_fixed_slip_xz_Nx%d_Δt%.1e", Nx, Δt),
                                                          field_slicer = nothing,
                                                          schedule = TimeInterval(stop_iteration * Δt / 2))

    return simulation
end

function setup_and_run_xz(args...; kwargs...)
    simulation = setup_xz_simulation(args...; kwargs...)
    @info "Running forced flow fixed slip simulation with Nx = Nz = $(kwargs[:Nx]), Δt = $(kwargs[:Δt])"
    @time run!(simulation)
    return simulation
end

end # module


text:
### This code was taken and modified from Optim.jl.
# Optim.jl is licensed under the MIT License:

# > Copyright (c) 2012: John Myles White, Tim Holy, and other contributors.
# > Copyright (c) 2016: Patrick Kofod Mogensen, John Myles White, Tim Holy,
# >                     and other contributors.
# > Copyright (c) 2017: Patrick Kofod Mogensen, Asbjørn Nilsen Riseth,
# >                     John Myles White, Tim Holy, and other contributors.
# >
# > Permission is hereby granted, free of charge, to any person
# > obtaining a copy of this software and associated documentation files
# > (the "Software"), to deal in the Software without restriction,
# > including without limitation the rights to use, copy, modify, merge,
# > publish, distribute, sublicense, and/or sell copies of the Software,
# > and to permit persons to whom the Software is furnished to do so,
# > subject to the following conditions:
# >
# > The above copyright notice and this permission notice shall be
# > included in all copies or substantial portions of the Software.
# >
# > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# > EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# > MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# > NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
# > BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# > ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# > CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# > SOFTWARE.

# mutable struct
abstract type Manifold end
struct ManifoldObjective{M <: Manifold, P, T, O <:AbstractDifferentiableObject{P,T}} <: AbstractDifferentiableObject{P,T}
    inner_obj::O
end
function ManifoldObjective(::M, inner_obj::O) where {M <: Manifold, P, T, O <: AbstractDifferentiableObject{P,T}}
    ManifoldObjective{M,P,T,O}(inner_obj)
end
# @inline function Base.getproperty(MO::ManifoldObjective{M}, s::Symbol) where M
#     if s == :inner_obj
#         return MO.inner_obj
#     elseif s == :manifold
#         return M
#     end
#     throw("ManifoldObjective has no field $s.")
# end

"""Flat Euclidean space {R,C}^N, with projections equal to the identity."""
struct Flat <: Manifold end
ManifoldObjective(::Flat, inner_obj::AbstractDifferentiableObject) = inner_obj
@inline retract(::Flat, x) = x
@inline retract!(::Flat,x) = x
@inline project_tangent(::Flat, g, x) = g
@inline project_tangent!(::Flat, g, x) = g
"""Spherical manifold {|x| = 1}."""
struct Sphere <: Manifold end
@inline retract!(::Sphere, x) = normalize!(x)
@inline project_tangent!(::Sphere, g, x) = (g .-= real(dot(x,g)).*x)
# @inline gradient(m::ManifoldObjective) = gradient(m.inner_obj)
function value!(obj::ManifoldObjective{M}, x) where M
    xin = retract(M, x)
    value!(obj.inner_obj, xin)
end
function value(obj::ManifoldObjective{M}) where M
    value(obj.inner_obj)
end
function gradient(obj::ManifoldObjective{M}) where M
    gradient(obj.inner_obj)
end
function gradient(obj::ManifoldObjective{M}, i::Int) where M
    gradient(obj.inner_obj,i)
end
function gradient!(obj::ManifoldObjective{M},x) where M
    xin = retract(M, x)
    gradient!(obj.inner_obj,xin)
    project_tangent!(M,gradient(obj.inner_obj),xin)
    return gradient(obj.inner_obj)
end
function value_gradient!(obj::ManifoldObjective{M}, x) where M
    xin = retract(M, x)
    value_gradient!(obj.inner_obj,xin)
    project_tangent!(M,gradient(obj.inner_obj),xin)
    return value(obj.inner_obj)
end

# abstract type AbstractOptimizer{M<:Manifold} end
# abstract type AbstractOptimizerState end
# abstract type FirstOrderOptimizer{M} <: AbstractOptimizer{M} end

# struct BFGS{M, IL, L, H<:Function} <: FirstOrderOptimizer{M}
#     alphaguess!::IL# = LineSearches.InitialStatic()
#     linesearch!::L# = LineSearches.HagerZhang()
#     initial_invH::H# = initial_diaginvH
#     manifold::M
# end
# function BFGS(;
#         alphaguess::IL = InitialStatic(),
#         linesearch::L = BackTracking(),
#         initial_invH::H = initial_diaginvH,
#         manifold::M = Flat()
#     ) where {IL, L, H, M <: Manifold}
#     BFGS{M, IL, L, H}(alphaguess, linesearch, initial_invH, manifold)
# end
# function BFGS{T}(;
#         alphaguess::IL = InitialStatic{T,false}(one(T)),
#         linesearch::L = BackTracking{T}(),
#         initial_invH::H = initial_diaginvH,
#         manifold::M = Flat()
#     ) where {T, IL <: InitialStatic{T}, L <: BackTracking{T}, H, M <: Manifold}
#     BFGS{M, IL, L, H}(alphaguess, linesearch, initial_invH, manifold)
# end

# # Base.summary(::BFGS) = "BFGS"

# @generated function initial_diaginvH(v::MVector{P,T}, α::T = one(T)) where {T,P}
#     quote
#         H = Symmetric(MMatrix{$P,$P,$T}(undef))
#         ind = 0
#         @inbounds for i = 1:$P
#             for j ∈ 1:i-1
#                 H.data[ind] = zero($T)
#             end
#             H.data[ind] = α
#             for j ∈ i+1:$P
#                 H.data[ind] = zero($T)
#             end
#         end
#         H
#     end
# end
# @generated function initial_diaginvH(v::SizedSIMDVector{P,T}, α::T = one(T)) where {T,P}
#     quote
#         H = Symmetric(SizedSIMDArray(undef, Val(($P,$P)), $T))
#         ind = 0
#         @inbounds for i = 1:$P
#             for j ∈ 1:i-1
#                 H.data[j,i] = zero($T)
#             end
#             H.data[i,i] = α
#             for j ∈ i+1:$P
#                 H.data[j,i] = zero($T)
#             end
#         end
#         H
#     end
# end
# @inline function initial_diaginvH(H::SymmetricMatrix{P,T},V::SizedVector{P,T}, α = one(T)) where {T,P}
#     # println("Resseting H.")
#     # @show H
#     @inbounds for i = 1:P
#         for j ∈ 1:i-1
#             H.data[j,i] = zero(T)
#         end
#         H.data[i,i] = α
#         for j ∈ i+1:P
#             H.data[j,i] = zero(T)
#         end
#     end
#     H
# end
# @inline function initial_diaginvH(H::SymmetricMatrix{P,T}, α = one(T)) where {T,P}
#     # println("Resseting H.")
#     # @show H
#     @inbounds for i = 1:P
#         for j ∈ 1:i-1
#             H.data[j,i] = zero(T)
#         end
#         H.data[i,i] = α
#         for j ∈ i+1:P
#             H.data[j,i] = zero(T)
#         end
#     end
#     H
# end

# mutable struct BFGSState{P,T,L,SV<:SizedSIMDVector{P,T,L},SM<:SymmetricMatrix{P,T}} <: AbstractOptimizerState
#     x::SV
#     x_previous::SV
#     g_previous::SV
#     f_x_previous::T
#     dx::SV
#     dg::SV
#     u::SV
#     invH::SM
#     s::SV
#     x_ls::SV
#     alpha::T
# end

# function initial_state(method::BFGS{M}, options, d, initial_x::SizedVector{P,T}) where {M,T,P}
#     initial_x = copy(initial_x)
#     retract!(method.manifold, initial_x)

#     value_gradient!!(d, initial_x)

#     project_tangent!(method.manifold, gradient(d), initial_x)
#     # Maintain a cache for line search results
#     # Trace the history of states visited
#     BFGSState(
#         initial_x, # Maintain current state in state.x
#         similar(initial_x), # Maintain previous state in state.x_previous
#         copy(gradient(d)), # Store previous gradient in state.g_previous
#         real(T)(NaN), # Store previous f in state.f_x_previous
#         similar(initial_x), # Store changes in position in state.dx
#         similar(initial_x), # Store changes in gradient in state.dg
#         similar(initial_x), # Buffer stored in state.u
#         method.initial_invH(initial_x), # Store current invH in state.invH
#         similar(initial_x), # Store current search direction in state.s
#         similar(initial_x), # Buffer of x for line search in state.x_ls
#         real(one(T))
#     )
# end

# function reset_search_direction!(state, d, method::BFGS)
#     method.initial_invH(state.invH, state.x)
#     grad = gradient(d)
#     @inbounds for i ∈ eachindex(state.s, grad)
#         state.s[i] = - grad[i]
#     end
#     return true
# end
# function update_state!(d, state::AbstractOptimizerState, method)

#     # @show state.s, state.invH, vec(gradient(d))
#     # Set the search direction
#     # Search direction is the negative gradient divided by the approximate Hessian
#     mul!(vec(state.s), state.invH.data, vec(gradient(d)))
#     SIMDArrays.scale!(state.s, -one(eltype(state.s)))
#     # # @show state.s
#     project_tangent!(method.manifold, state.s, state.x)
#     # @show state.s, state.x

#     # Maintain a record of the previous gradient
#     copyto!(state.g_previous, gradient(d))

#     # Determine the distance of movement along the search line
#     # This call resets invH to initial_invH is the former in not positive
#     # semi-definite
#     # lssuccess = perform_linesearch!(state, method, ManifoldObjective(M, d))
#     perform_linesearch!(state, method, ManifoldObjective(method.manifold, d)) # 120 ns
#     # @show state.alpha, state.s

#     # Update current position
#     @inbounds @simd for i ∈ 1:full_length(state.dx)
#         state.dx[i] = state.alpha * state.s[i]
#         state.x[i] = state.x[i] + state.dx[i]
#     end
#     retract!(method.manifold, state.x)
#     # @show state.x

#     # lssuccess == false # break on linesearch error
#     false
# end
# function perform_linesearch!(state, method, d)
#     # Calculate search direction dphi0
#     # # @show typeof(d)
#     # # @show @which gradient(d)
#     # @show gradient(d), state.s
#     dphi_0 = real(dot(gradient(d), state.s))
#     # @show dphi_0
#     # reset the direction if it becomes corrupted
#     if dphi_0 >= zero(dphi_0) && reset_search_direction!(state, d, method)
#         dphi_0 = real(dot(gradient(d), state.s)) # update after direction reset
#     end
#     phi_0  = value(d)

#     # Guess an alpha
#     method.alphaguess!(method.linesearch!, state, phi_0, dphi_0, d)

#     # Store current x and f(x) for next iteration
#     state.f_x_previous = phi_0
#     copyto!(state.x_previous, state.x)

#     # Perform line search; catch LineSearchException to allow graceful exit
#     # try
#         # state.alpha, ϕalpha =
#     linesearch!(method.linesearch!, d, state.x, state.s, state.alpha, state.x_ls, phi_0, dphi_0)
#     return true # lssuccess = true
#     # catch ex
#     #     if isa(ex, LineSearches.LineSearchException)
#     #         state.alpha = ex.alpha
#     #         Base.warn("Linesearch failed, using alpha = $(state.alpha) and exiting optimization.\nThe linesearch exited with message:\n$(ex.message)")
#     #         return false # lssuccess = false
#     #     else
#     #         rethrow(ex)
#     #     end
#     # end
# end
# #
# # function create_quote()
# #     q = quote @fastmath @inbounds begin end end
# #     @static if VERSION > v"0.7-"
# #         qa = q.args[2].args[3].args[3].args
# #     else
# #         qa = q.args[2].args[2].args[2].args
# #     end
# #     q, qa
# # end

# function update_h!(d, state::BFGSState{P,T,L}, method) where {P,T,L}
#     # Measure the change in the gradient
#     ∇ = gradient(d)
#     @inbounds @simd for i ∈ 1:L
#         state.dg[i] =  ∇[i] - state.g_previous[i]
#     end

#     # @show gradient(d), state.g_previous, state.dg
#     # Update the inverse Hessian approximation using Sherman-Morrison
#     dx_dg = real(dot(state.dx, state.dg))
#     dx_dg == 0.0 && return true # force stop
#     # @show dx_dg, state.invH
#     # @show state.u
#     # I may not support the use of `vec` here.
#     mul!(vec(state.u), state.invH.data, vec(state.dg))

#     c2 = 1 / dx_dg
#     c1 = fma(real(dot(state.dg, state.u)), c2*c2, c2)

#     # c1v2 = ( dx_dg + real(dot(state.dg, state.u))) / (dx_dg' * dx_dg)
#     # c2v2 = 1 / dx_dg
#     # @show state.u
#     # @show state.dx
#     # @show state.dg

#     # hx = copy(state.invH)
#     # @inbounds for i ∈ 1:P
#     #     sdxi = state.dx[i]
#     #     sui = state.u[i]
#     #     @simd for j ∈ 1:L
#     #         hx[j,i] += c1 * sdxi * state.dx[j] - c2 * (sui * state.dx[j] + state.u[j] * sdxi)
#     #     end
#     # end
#     # @show hx

#     # @show state.x, state.invH
#     SIMDArrays.BFGS_update!(state.invH, state.dx, state.u, c1, c2)
#     # @show state.x, state.invH
#     false
# end

# function maxdiff(x::AbstractArray{T1}, y::AbstractArray{T2}) where {T1,T2}
#     res = real(zero(promote_type(T1,T2)))
#     @inbounds for i in 1:length(x)
#         delta = abs(x[i] - y[i])
#         if delta > res
#             res = delta
#         end
#     end
#     return res
# end
# f_abschange(d::AbstractDifferentiableObject, state) = f_abschange(value(d), state.f_x_previous)
# f_abschange(f_x::T, f_x_previous) where T = abs(f_x - f_x_previous)
# x_abschange(state) = x_abschange(state.x, state.x_previous)
# x_abschange(x, x_previous) = maxdiff(x, x_previous)
# g_residual(d::AbstractDifferentiableObject) = g_residual(gradient(d))
# # g_residual(g) = norm(g, Inf)
# @inline g_residual(g) = SIMDArrays.maximum_abs(g)
# function gradient_convergence_assessment(state::AbstractOptimizerState, d, options)
#     g_residual(gradient(d)) ≤ options.g_tol
# end
# # gradient_convergence_assessment(state::ZerothOrderState, d, options) = false

# function assess_convergence(state::BFGSState, d, options)
#     default_convergence_assessment(state, d, options)
# end
# function default_convergence_assessment(state::AbstractOptimizerState, d, options)
#     x_converged, f_converged, f_increased, g_converged = false, false, false, false

#     # TODO: Create function for x_convergence_assessment
#     if x_abschange(state.x, state.x_previous) ≤ options.x_tol
#         x_converged = true
#     end

#     # Relative Tolerance
#     # TODO: Create function for f_convergence_assessment
#     f_x = value(d)
#     if f_abschange(f_x, state.f_x_previous) ≤ options.f_tol*abs(f_x)
#         f_converged = true
#     end

#     if f_x > state.f_x_previous
#         f_increased = true
#     end
#     g_converged = gradient_convergence_assessment(state,d,options)

#     converged = x_converged || f_converged || g_converged

#     return x_converged, f_converged, g_converged, converged, f_increased
# end

function BFGS_update_quote(Mₖ,Pₖ,stride_AD,T)
    T_size = sizeof(T)
    AD_stride = stride_AD * T_size
    reg_size = VectorizationBase.REGISTER_SIZE
    W = reg_size ÷ T_size
    Q, r = divrem(stride_AD, W) #Assuming stride_AD is a multiple of W
    if Q > 0
        r == 0 || throw("Number of rows plus padding $stride_AD not a multiple of register size: $(VectorizationBase.REGISTER_SIZE).")
        L = VectorizationBase.CACHELINE_SIZE
    else
        W = r
        Q = 1
    end
    V = Vec{W,T}
    # C = VectorizationBase.CACHELINE_SIZE ÷ T_size

    common = quote
        vC1 = vbroadcast($V, c1)
        vC2 = vbroadcast($V, -c2)
        ptr_invH, ptr_S, ptr_U = pointer(invH), pointer(s), pointer(u)
    end
    if Q > 1
        return quote
            $common
            # @nexprs $Pₖ p -> @nexprs $Q q -> invH_p_q = vload($V, ptr_invH + $REGISTER_SIZE*(q-1) + $AD_stride*(p-1))

            # invH = invH + c1 * (s * s') - c2 * (u * s' + s * u')
            # @nexprs $Pₖ p -> begin
            for p ∈ 0:$T_size:$((Pₖ-1)*T_size)
                # @nexprs $Q q -> invH_q = vload($V, ptr_invH + $REGISTER_SIZE*(q-1) + $AD_stride*(p-1))
                # @nexprs $Q q -> invH_q = vload($V, ptr_invH + $REGISTER_SIZE*(q-1) + $AD_stride*p)
                # vSb = vbroadcast($V, unsafe_load(ptr_S + p-1 ))
                vSb = vbroadcast($V, VectorizationBase.load(ptr_S + p ))
                vSbc1 = vmul(vSb, vC1)
                vSbc2 = vmul(vSb, vC2)
                # vUbc2 = vbroadcast($V, unsafe_load(ptr_U + p-1 )) * vC2
                vUbc2 = vmul(vbroadcast($V, VectorizationBase.load(ptr_U + p )), vC2)
                # or
                # vSbc1 = vbroadcast($V, c1*unsafe_load(ptr_S + p-1 ))
                # vSbc2 = vbroadcast($V, c2*unsafe_load(ptr_S + p-1 ))
                # vUbc2 = vbroadcast($V, c2*unsafe_load(ptr_U + p-1 ))
                Base.Cartesian.@nexprs $Q q -> begin # I am concerned over the size of these dependency chains.
                    invH_q = vload($V, ptr_invH + $reg_size*(q-1) + $stride_AD*p)
                    vU_q = vload($V, ptr_U + $reg_size*(q-1))
                    invH_q = vmuladd(vU_q, vSbc2, invH_q)
                    vS_q = vload($V, ptr_S + $reg_size*(q-1))
                    invH_q = vmuladd(vS_q, vSbc1, vmuladd(vS_q, vUbc2, invH_q))
                    vstore!(ptr_invH + $reg_size*(q-1) + $stride_AD*p, invH_q)
                end
            end
            nothing
        end
    else
        return quote
            $common
            # invH = invH + c1 * (s * s') - c2 * (u * s' + s * u')
            Base.Cartesian.@nexprs $Pₖ p -> begin
                invH_p = vload($V, ptr_invH + $AD_stride*(p-1))
                vSb = vbroadcast($V, VectorizationBase.load(ptr_S + (p-1)*$T_size ))
                vS = vload($V, ptr_S)
                # Split up the dependency chain and reduce number of operations, when we can't save on vC1 and vC2 multiplications.
                invH_p = vmuladd(vC2, vmuladd( vload($V, ptr_U), vSb, vmul(vS, vbroadcast($V, VectorizationBase.load(ptr_U + (p-1)*$T_size )))), vmuladd( vS, vmul(vSb,vC1), invH_p ))
                vstore!(ptr_invH + $AD_stride*(p-1), invH_p)
            end
            nothing
        end
    end
end

@generated function BFGS_update!(invH::Union{Symmetric{T,<:AbstractMutableFixedSizeMatrix{P,P,T,R,L}},<:AbstractMutableFixedSizeMatrix{P,P,T,R,L}},
    s::AbstractMutableFixedSizeVector{P,T,R}, u::AbstractMutableFixedSizeVector{P,T,R}, c1::T, c2::T) where {P,T,L,R}

    BFGS_update_quote(P,P,R,T)

end


text:
export cross_validation

struct CrossValidation
    abs_errors::Matrix{Float64}
    mae::Vector{Float64}
    crps_scores::Matrix{Float64}
    mean_crps::Vector{Float64}
    function CrossValidation(n::Int, steps_ahead::Int)
        abs_errors = Matrix{Float64}(undef, steps_ahead, n)
        crps_scores = Matrix{Float64}(undef, steps_ahead, n)
        mae = Vector{Float64}(undef, steps_ahead)
        mean_crps = Vector{Float64}(undef, steps_ahead)
        return new(abs_errors, mae, crps_scores, mean_crps)
    end
end

discrete_crps_indicator_function(val::T, z::T) where {T} = val < z
function crps(val::T, scenarios::Vector{T}) where {T}
    sorted_scenarios = sort(scenarios)
    m = length(scenarios)
    crps_score = zero(T)
    for i = 1:m
        crps_score +=
            (sorted_scenarios[i] - val) *
            (m * discrete_crps_indicator_function(val, sorted_scenarios[i]) - i + 0.5)
    end
    return (2 / m^2) * crps_score
end
evaluate_abs_error(y::Vector{T}, forecast::Vector{T}) where T = abs.(y - forecast)
function evaluate_crps(y::Vector{T}, scenarios::Matrix{T}) where {T}
    crps_scores = Vector{T}(undef, length(y))
    for k = 1:length(y)
        crps_scores[k] = crps(y[k], scenarios[k, :])
    end
    return crps_scores
end

"""
TODO
"""
function cross_validation(gas::Model{<:Distribution, T}, y::Vector{T}, steps_ahead::Int, start_idx::Int;
                  S::Int = 10_000,
                  initial_params::Matrix{T} = DEFAULT_INITIAL_PARAM,
                  opt_method = NelderMead(gas, DEFAULT_NUM_SEEDS)) where T
    num_mle = length(y) - start_idx - steps_ahead
    b = CrossValidation(num_mle, steps_ahead)
    for i in 1:num_mle
        println("CrossValidation: step $i of $num_mle")
        gas_to_fit = deepcopy(gas)
        y_to_fit = y[1:start_idx - 1 + i]
        y_to_verify = y[start_idx + i:start_idx - 1 + i + steps_ahead]
        fit!(gas_to_fit, y_to_fit; initial_params=initial_params, opt_method=opt_method)
        forec = initial_params !== DEFAULT_INITIAL_PARAM ?
                    forecast(y_to_fit, gas_to_fit, steps_ahead; S=S, initial_params=initial_params) :
                    forecast(y_to_fit, gas_to_fit, steps_ahead; S=S)
        abs_errors = evaluate_abs_error(y_to_verify, forec.observation_forecast)
        crps_scores = evaluate_crps(y_to_verify, forec.observation_scenarios)
        b.abs_errors[:, i] = abs_errors
        b.crps_scores[:, i] = crps_scores
    end
    for i in 1:steps_ahead
        b.mae[i] = mean(b.abs_errors[i, :])
        b.mean_crps[i] = mean(b.crps_scores[i, :])
    end
    return b
end


text:
@testset "IdentityAxis" begin
    r = IdentityAxis(0, -5)
    @test isempty(r)
    @test length(r) == 0
    @test size(r) == (0,)
    r = IdentityAxis(0, 2)
    @test !isempty(r)
    @test length(r) == 3
    @test size(r) == (3,)
    @test axes(r) == (0:2,)
    @test step(r) == 1
    @test first(r) == 0
    @test last(r) == 2
    @test minimum(r) == 0
    @test maximum(r) == 2
    @test r[0] == 0
    @test r[1] == 1
    @test r[2] == 2
    @test_throws BoundsError r[3]
    @test_throws BoundsError r[-1]
    @test r[0:2] == IdentityAxis(0:2)
    @test r[r] == r

    @test r .+ 1 != 1:3
    # TODO @test r .+ 1 == AxisArray(1:3, axes(r))
    @test r .+ 1 === 1 .+ r
    # TODO @test r .- 1 === AxisArray(-1:1, axes(r))
    # TODO @test 1 .- r === OffsetArray(1:-1:-1, axes(r))
    # TODO @test 2 * r === 2 .* r === r * 2 === r .* 2 === OffsetArray(0:2:4, axes(r))
    k = -1
    for i in r
        @test i == (k+=1)
    end
    @test k == length(r)-1
    @test collect(r) == [0,1,2]
    # TODO @test intersect(r, IdentityAxis(-1,1)) === intersect(IdentityAxis(-1,1), r) === IdentityAxis(0,1)
    @test intersect(r, -1:5) === intersect(-1:5, r) === 0:2
    @test intersect(r, 2:5) === intersect(2:5, r) === 2:2
    # TODO @test string(r) == "IdentityAxis(0:2)"

    r = IdentityAxis(2:4)
    # TODO @test r != 2:4
    @test checkindex(Bool, r, 4)
    @test !checkindex(Bool, r, 5)
    @test checkindex(Bool, r, :)
    @test checkindex(Bool, r, 2:4)
    @test !checkindex(Bool, r, 1:5)
    @test !checkindex(Bool, r, trues(4))
    @test !checkindex(Bool, r, trues(5))
    # TODO @test convert(UnitRange, r) == 2:4
    @test convert(StepRange, r) == 2:1:4
    # TODO  @test !in(1, r)
    @test in(2, r)
    # TODO @test in(4, r)
    @test !in(5, r)
    @test issorted(r)
    @test maximum(r) == 4
    @test minimum(r) == 2
    # TODO @test sortperm(r) == r
    # TODO @test r != 2:4
    @test 1:4 == IdentityAxis(1:4) == 1:4
    @test r+r == AxisArray(4:2:8, axes(r))
    # TODO this can't be done with other AbstractUnitRange types so why here?
    # @test r-r == OffsetArray([0,0,0], axes(r))
    @test (9:2:13)-r == 7:9
    @test -r == AxisArray(-2:-1:-4, axes(r))
    @test reverse(r) == AxisArray(4:-1:2, axes(r))
    @test r / 2 == r ./ 2 == AxisArray(1:0.5:2, axes(r))
    @test 2 \ r == 2 .\ r == AxisArray(1:0.5:2, axes(r))

    r = IdentityAxis{Int16}(0, 4)
    @test length(r) === 5
    @test iterate(r) == (0,0)
    k = -1
    for i in r
        @test i == (k+=1)
    end
    @test k == length(r)-1
    #=
    x, y = promote(IdentityAxis(2,4), IdentityAxis{Int16}(3,7))
    @test x === IdentityAxis(2:4)
    @test y === IdentityAxis(3:7)
    x, y = promote(IdentityAxis(4:5), 0:7)
    @test x === 4:5
    @test y === 0:7
    @test convert(IdentityAxis{Int16}, IdentityAxis(2:5)) === IdentityAxis{Int16}(2:5)
    @test convert(IdentityAxis{Int}, IdentityAxis(2:5)) === IdentityAxis(2:5)
    @test convert(UnitRange, IdentityAxis(2:4)) === 2:4
    r = IdentityAxis(Int128(1),Int128(10))
    @test length(r) === Int128(10)
    =#
end