text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

function boykov_kolmogorov_cut(A::AbstractArray{T,N}, B::AbstractArray{T,N}, dim::Integer) where {N,T}
  @assert size(A) == size(B) "arrays must have the same size for cut"

  # size and number of voxels
  sz   = size(A)
  nvox = prod(sz)

  # source and sink terminals
  s = nvox + 1; t = nvox + 2

  # lattice graph and capacity matrix
  G = DiGraph(nvox+2)
  C = spzeros(nvox+2, nvox+2)

  # fill lattice graph with original vertices
  for d=1:N
    # loop over all indices except the borders using an increment direction
    inds = CartesianIndices(ntuple(i -> i == d ? (1:sz[d]-1) : (1:sz[i]), N))
    incr = CartesianIndex(ntuple(i -> i == d ? 1 : 0, N))

    for ind in inds
      # adjacent vertices along direction
      u = cart2lin(sz, ind)
      v = cart2lin(sz, ind + incr)

      # simple difference
      Du = abs(A[u] - B[u])
      Dv = abs(A[v] - B[v])

      # gradient along direction
      ∇Au = abs(A[v] - A[u])
      ∇Bu = abs(B[v] - B[u])

      # next vertex along direction
      w = ind + 2incr

      # repeat gradient on the border when
      # outside valid indices for arrays
      ok = all(w.I .≤ sz)
      ∇Av = ok ? abs(A[w] - A[v]) : ∇Au
      ∇Bv = ok ? abs(B[w] - B[v]) : ∇Bu

      add_edge!(G, u, v)
      add_edge!(G, v, u)

      C[u,v] = C[v,u] = (Du + Dv) / (∇Au + ∇Av + ∇Bu + ∇Bv + eps())
    end
  end

  # fill edges from source terminal to left slice
  linds = CartesianIndices(ntuple(i -> i == dim ? (1:1) : (1:sz[i]), N))
  for ind in linds
    u = cart2lin(sz, ind)
    add_edge!(G, s, u)
    C[s,u] = Inf
  end

  # fill edges from right slice to sink terminal
  rinds = CartesianIndices(ntuple(i -> i == dim ? (sz[dim]:sz[dim]) : (1:sz[i]), N))
  for ind in rinds
    v = cart2lin(sz, ind)
    add_edge!(G, v, t)
    C[v,t] = Inf
  end

  # Boykov-Kolmogorov max-flow/min-cut
  _, __, labels = maximum_flow(G, s, t, C, algorithm=BoykovKolmogorovAlgorithm())

  # remove source and sink terminals
  labels = labels[1:end-2]

  # cut mask
  M = falses(size(A))
  M[labels .== 1] .= true
  M[labels .== 0] .= true

  M
end


text:
include("solver_preamble.jl")
using COSMO
# sosdemo9: 33000 < max_iter < 34000
factory = with_optimizer(COSMO.Optimizer, verbose=false, max_iter = 34000)
config = MOI.Test.TestConfig(atol=1e-3, rtol=1e-3)
@testset "Linear" begin
    Tests.linear_test(factory, config, [
        # See https://github.com/oxfordcontrol/COSMO.jl/issues/96
        "dsos_horn"
    ])
end
@testset "SOC" begin
    Tests.soc_test(factory, config, [
        "sdsos_options_pricing",
        # See https://github.com/oxfordcontrol/COSMO.jl/issues/96
        "sdsos_horn"
    ])
end
@testset "SDP" begin
    Tests.sd_test(factory, config, [
        "sos_options_pricing",
        # Iteration limit
        "chebyshev",
        # Expression: JuMP.termination_status(model) == MOI.INFEASIBLE
        # Evaluated: MathOptInterface.OPTIMAL == MathOptInterface.INFEASIBLE
        "quartic_infeasible_lyapunov_switched_system",
        "quartic_infeasible_scaled_lyapunov_switched_system",
        "maxcut", "sos_horn", "motzkin"
    ])
end

ext: jl


text:
Base.rand(m::SequentialPointProcess, window, npoints::Integer; giveup=-1, returnanyway=false) =
  rpattern(m, window, npoints, giveup, returnanyway)
function rpattern(model, window, npoints, giveup=-1, returnanyway=false)
  xy = Vector{NTuple{2,Float64}}(undef, npoints)
  if npoints == 0
    return PointPattern(xy, window)
  end
  xy[1] = randpoint(Uniform(), xy, window, giveup)
  for i in 2:npoints
    prop = randpoint(model, view(xy, 1:i-1), window, giveup)
    if isnothing(prop) # If sampling from the sequential model fails at some point return nothing.
      # This works with ABC, since that kind of parameter configuration is too unlikely.
      # In other contexts, giveup is not used.
      if returnanyway
        return PointPattern(xy[1:i], window)
      else
        return PointPattern(empty(xy), window)
      end
    end
    xy[i] = prop
  end
  PointPattern(xy, window)
end
function acceptpoint(m::Softcore, xyprop, xyold)
  p = fk(m, xyprop, xyold)
  return rand() <= p
end
function acceptpoint(m::OverlappingDiscs, xyprop, xyold)
  k = length(xyold) + 1
  p = m.theta(countoverlaps(m, xyprop, xyold), k) / maximum(n -> m.theta(n,k), 0:length(xyold))
  return rand() <= p
end
function acceptpoint(::Uniform, xyprop, xyold)
  true
end
function randpoint(m::Mixture, pattern, window, giveup)
  if rand() < m.theta
    randpoint(m.first, pattern, window, giveup)
  else
    randpoint(m.second, pattern, window, giveup)
  end
end
function randp(window)
  randcoord(x) = window[x][1] + (window[x][2]-window[x][1])*rand()
  randx() = randcoord(:x)
  randy() = randcoord(:y)
  randx(), randy()
end
function randpoint(model, pattern, window, giveup=-1)
  tries = 0
  while true
    xyprop = randp(window)
    if acceptpoint(model, xyprop, pattern)
      return xyprop
    end
    tries += 1
    if giveup != -1 && tries > giveup
      return nothing
    end
  end
  return nothing
end

max_stars_repo_name: mikkoku/SequentialPointProcesses

max_issues_repo_name: mikkoku/SequentialPointProcesses

max_forks_repo_name: mikkoku/SequentialPointProcesses


text:
using Pkg
Pkg.activate(".")
pkg"dev .. ../CairoMakie ../GLMakie ../WGLMakie"
Pkg.instantiate()
Pkg.precompile()

using NodeJS
run(`$(npm_cmd()) install highlight.js`)
run(`$(npm_cmd()) install lunr`)
run(`$(npm_cmd()) install cheerio`)

using Franklin
using Documenter: deploydocs, deploy_folder, GitHubActions

cfg = GitHubActions() # this should pick up all details via GHA environment variables

repo = "github.com/JuliaPlots/Makie.jl.git"
push_preview = true

deploydecision = deploy_folder(cfg; repo, push_preview,
    devbranch = "master",
    devurl = "dev",
)

@info "Setting PREVIEW_FRANKLIN_WEBSITE_URL to $repo"
ENV["PREVIEW_FRANKLIN_WEBSITE_URL"] = repo
@info "Setting PREVIEW_FRANKLIN_PREPATH to $(deploydecision.subfolder)"
ENV["PREVIEW_FRANKLIN_PREPATH"] = deploydecision.subfolder

serve(single=true, cleanup=false)
lunr()
optimize(minify=false, prerender=false)

deploydocs(;
    repo,
    push_preview,
    target = "__site",
)


text:
function my_retry(f::Function, num_retries::Integer = 1)
    result = retry(f, delays=ExponentialBackOff(n=num_retries))()
    return result
end

function my_retry_suppress_exceptions(f::Function, num_retries::Integer = 1)
    result = try
        my_retry(f, num_retries)
    catch ex
        showerror(stderr, ex)
        Base.show_backtrace(stderr, catch_backtrace())
        println(stderr)
        nothing
    end
    return result
end


text:
module CommunalHelperDreamFallingBlock
using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/DreamFallingBlock" DreamFallingBlock(x::Integer, y::Integer,
	width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight,
	featherMode::Bool = false, oneUse::Bool = false, doubleRefill::Bool=false)

const placements = Ahorn.PlacementDict(
    "Dream Falling Block (Communal Helper)" => Ahorn.EntityPlacement(
        DreamFallingBlock,
		"rectangle"
    )
)

Ahorn.minimumSize(entity::DreamFallingBlock) = 8, 8
Ahorn.resizable(entity::DreamFallingBlock) = true, true

function Ahorn.selection(entity::DreamFallingBlock)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, width, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DreamFallingBlock, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

	 Ahorn.Cairo.save(ctx)

    Ahorn.set_antialias(ctx, 1)
    Ahorn.set_line_width(ctx, 1)

    fillColor = get(entity.data, "featherMode", false) ? (0.31, 0.69, 1.0, 0.4) : (0.0, 0.0, 0.0, 0.4)
	 lineColor = get(entity.data, "oneUse", false) ? (1.0, 0.0, 0.0, 1.0) : (1.0, 1.0, 1.0, 1.0)
    Ahorn.drawRectangle(ctx, 0, 0, width, height, fillColor, lineColor)

    Ahorn.restore(ctx)
end

end

max_stars_repo_name: FlynnSC/CommunalHelper

max_issues_repo_name: FlynnSC/CommunalHelper

max_forks_repo_name: FlynnSC/CommunalHelper


text:
using LinearAlgebra

"""
Result of Newton's iterations.
Contain `converged::Bool` flag, the minimum `argument::Vector`,
number of `iterations::Integer` and `calls::Integer`.
"""
struct NewtonResult{T}
    converged::Bool
    argument::Vector{T}
    iterations::Int
    calls::Int

    NewtonResult(conv, x, iters=-1, calls=-1) = new{Float64}(conv, copy(x), iters, calls)
end

"""
    backtracking_line_search(f, x₀, d, y₀[; α₀, p, buf])

Find `x = x₀ + α*d` (`x₀::AbstractVector`, `d::AbstractVector`) that `f(x)::Real`
is < `y₀ == f(x₀)` by exponentially decreasing `::Real` `α = α₀ * pⁱ` (`p::Real=0.5` must be < 1).
Optional `buf` is similar to `x₀` for storing `x₀ + α*d` trial state.

Return tuple of `α`, `f(x₀ + α*d)` and number of `f`'s calls.
"""
function backtracking_line_search(
    f::Function,
    x₀::AbstractVector{T},
    d::AbstractVector{T},
    y₀::T;
    α::T=one(T),
    p::Real=one(T)/2,
    buf::AbstractVector{T}=similar(x₀),
) where {T<:Real}
    xtry = buf
    calls = 0
    ftry = T(NaN)

    # (?) better `maxiter = 1 + ceil(Int, log(eps(α))/log(p))`
    # If use above, then `p^maxiter ≈ eps(α₀)`
    maxiter = 200
    for i in 1:maxiter
        @. xtry = x₀ + α*d
        calls += 1
        ftry = f(xtry)

        @debug "backtracking_line_search" i α p ftry y₀ ftry-y₀

        ftry < y₀ && break
        α *= p

        i == maxiter && error("Line search number of iterations ($(maxiter)) exceeded.")
    end
    return α, ftry, calls
end

"""
    newton(f, grad!, hess!, x₀[; gtol, maxiter, constrain_step])

Find constrained minimum of `f(x)::Real` using Newton's solver with line search
and hessian modified by Cholesky factorization.
`grad!(g, x)` must update and return `g`radient of `f` at `x`.
`hess!(H, x)` must update and return `H`essian matrix of `f` at `x`.

The optimization starts at `x₀` performing no more than `maxiter` Newton's steps.
The steps finish when norm of gradient is ≤ `gtol`.
New point `xnew <- x + αmax * d` is constrained in Newton's `d`irection by
`constrain_step(x, d) -> α::Real`, `αmax = min(1.0, α)`.

Return [`NewtonResult`](@ref) object.
"""
function newton(
    f::Function,
    ∇f!::Function,
    H!::Function,
    x::AbstractVector;
    gtol::Real=1e-6,
    maxiter::Integer=200,
    constrain_step::Function=(x, δ)->Inf,
)
    x = convert(Vector{Float64}, x)
    ∇ = similar(x)
    δx = similar(x)
    hess_full = Matrix{Float64}(undef, size(∇, 1), size(∇, 1))
    vec = similar(x)

    fval = f(x)
    totfcalls = 1

    """
    Returns closure for Downhill.strong_backtracking!
    Vector `v` for determine types.
    """
    function fdfclosure(v::AbstractVector)
        vecx = similar(v)
        vecg = similar(v)

        function fdf(x::AbstractVector, α::Real, d::AbstractVector)
            vecx .= x .+ α .* d
            return f(vecx), ∇f!(vecg, vecx)
        end

        return fdf
    end

    fdf = fdfclosure(x)

    # check gradient norm before iterations
    ∇ = ∇f!(∇, x)
    if norm(∇, 2) ≤ gtol
        return NewtonResult(true, x, 0, totfcalls)
    end

    for i in 1:maxiter
        # descent direction `δx`

        hess_full = H!(hess_full, x)
        hess = Downhill.mcholesky!(hess_full)

        vec .= ∇
        ldiv!(hess, vec)  # `hess \ ∇`, result in `vec`
        @. δx = -vec

        # choosing step in `δx`
        αmax = constrain_step(x, δx)
        α = Downhill.strong_backtracking!(fdf, x, δx;
            α=1.0,  # the function takes care of initial `α` and `αmax`
            αmax=αmax,
            σ=0.9,
        )
        fcalls = 0  # by current implementation they are unknown :c

        # previous code, which uses "simple" backtracking
        # αmax = min(1.0, constrain_step(x, δx))
        # α, fval, fcalls = backtracking_line_search(f, x, δx, fval; α=αmax, p=0.5, buf=vec)

        totfcalls += fcalls

        if α ≤ 0
            @error "newton: linesearch failed" α
            return NewtonResult(false, x, i, totfcalls)
        end

        # update argument `x`, function value and gradient in `x`
        @. x += α * δx
        fval = f(x)
        ∇ = ∇f!(∇, x)

        @debug "newton" i repr(δx) norm_step=norm(δx, 2) α αmax repr(x) fval fcalls norm_gprev=norm(∇f!(similar(x), x - α*δx)) norm_gnew=norm(∇) prod(diag(hess.U))

        # check for convergence
        if norm(∇, 2) ≤ gtol
            return NewtonResult(true, x, i, totfcalls)
        end
    end
    return NewtonResult(false, x, maxiter, totfcalls)
end

"""
Construct functions of proper signature for [`newton`](@ref) algorithm.
The function reuses `helmholtz_diff!` and `helmholtz_diff_grad!` which are
same as for BFGS-version of vtflash.

(!) BFGS' closures must be constructed from the same thermodynamical properties.
"""
function __vt_flash_newton_closures(
    helmholtz_diff!::Function,
    helmholtz_diff_grad!::Function,
    mix::BrusilovskyEoSMixture{T},
    nmol::AbstractVector{T},
    volume::T,
    RT::T,
) where {T<:Real}
    vecnc₊ = Vector{T}(undef, size(nmol, 1) + 1)
    thermo_buf = thermo_buffer(mix)
    hessian_buf = HessianBuffer(mix)

    # below `x` is mixtures' state
    # helmholtz
    function func(x::AbstractVector{T})
        ΔA, _ = helmholtz_diff!(x, vecnc₊)
        return ΔA
    end

    # helmholtz gradient
    function gradient!(gradient_::AbstractVector{T}, x::AbstractVector{T})
        helmholtz_diff_grad!(x, gradient_)  # yes, base func uses this signature
        return gradient_
    end

    # helmholtz hessian
    function hessian!(hessian::AbstractMatrix{T}, x::AbstractVector{T})
        __vt_flash_hessian!(hessian, x, mix, nmol, volume, RT; buf=hessian_buf)
        return hessian
    end

    return func, gradient!, hessian!
end

function vt_flash_newton(
    mix::BrusilovskyEoSMixture{T},
    nmol::AbstractVector,
    volume::Real,
    RT::Real,
    unstable_state::AbstractVector,
) where {T}
    state = copy(unstable_state)

    # bfgs closures, just to reuse later
    constrain_step, helmholtz_diff_grad!, helmholtz_diff! = vt_flash_closures(
        mix, nmol, volume, RT
    )
    # closures for newton algorithm
    newton_helmholtz_diff, newton_helmholtz_diff_grad!, newton_helmholtz_diff_hessian! =
        __vt_flash_newton_closures(
            helmholtz_diff!,
            helmholtz_diff_grad!,
            mix,
            nmol,
            volume,
            RT,
        )

    # run optimizer
    result = newton(
        newton_helmholtz_diff,
        newton_helmholtz_diff_grad!,
        newton_helmholtz_diff_hessian!,
        state;
        constrain_step=constrain_step,
        gtol=1e-3,
        maxiter=200,
    )

    return __vt_flash_two_phase_result(mix, nmol, volume, RT, result)
end

"""
    vt_flash_newton(mix, nmol, volume, RT)

Find VT-equilibrium of `mix`, at given `nmol`, `volume` and thermal energy `RT`
using Newton's minimization. Return [`VTFlashResult`](@ref).
"""
function vt_flash_newton(
    mix::BrusilovskyEoSMixture{T},
    nmol::AbstractVector,
    volume::Real,
    RT::Real,
) where {T}
    # run vt-stability to find out whether a state single phase or not
    singlephase, vt_stab_tries = vt_stability(mix, nmol, volume, RT)
    @debug "VTFlash: VTStability result" singlephase

    if singlephase
        return VTFlashResult{T}(
            converged=true,
            singlephase=true,
            RT=RT,
            nmol_1=nmol,
            V_1=volume,
            nmol_2=similar(nmol),
            V_2=0
        )
    end

    # two-phase state case
    init_found, state = __vt_flash_initial_state(
        mix, nmol, volume, RT, vt_stab_tries
    )

    @debug "VTFlash: initial state search result" found=init_found state=repr(state) ΔA=helmholtz_diff!(state, similar(state))

    if !init_found
        @error "VTFlash: Initial state was not found!" mixture=mix nmol=repr(nmol) volume=volume RT=RT
        error("VTFlash: Initial state was not found!")
    end

    return vt_flash_newton(mix, nmol, volume, RT, state)
end


text:
using Test, FluxReconstruction, KitBase

cd(@__DIR__)
include("test_triangle.jl")

let u = rand(3), u0 = rand(3)
    FR.L1_error(u, u0, 0.1)
    FR.L2_error(u, u0, 0.1)
    FR.L∞_error(u, u0, 0.1)
end

shock_detector(-Inf, 3)
shock_detector(log10(0.1), 3)
shock_detector(log10(1e5), 3)

deg = 5
ps2 = FRPSpace2D(0.0, 1.0, 20, 0.0, 1.0, 20, deg, 1, 1)

rs_jacobi(ps2.xpl, [0 0; √3 -1; √3+1 √3-1; 1 √3])
rs_jacobi(ps2.xpl, rand(3, 3, 4, 2))

ps1 = TriFRPSpace("../assets/linesource.msh", 2)
ps = FRPSpace1D(0.0, 1.0, 20, deg)

positive_limiter(ones(6), 1 / 6, ps.ll, ps.lr)
positive_limiter(ones(6, 3), 5 / 3, 1 / 6, ps.ll, ps.lr)

let u = rand(deg + 1)
    ℓ = FR.basis_norm(deg)

    modal_filter!(u, 1e-6; filter = :l2)
    modal_filter!(u, 1e-6; filter = :l2opt)
    modal_filter!(u, 1e-6, ℓ; filter = :l1)
    modal_filter!(u, ℓ; filter = :lasso)
    modal_filter!(u, 10; filter = :exp)
    modal_filter!(u, 10; filter = :houli)

    # 2D exponential filter
    FR.filter_exp(2, 10, Array(ps.V))
end

let u = rand(deg+1, deg+1)
    modal_filter!(u, 1e-6, 1e-6; filter = :l2)
    modal_filter!(u, 1e-6, 1e-6; filter = :l2opt)
end

f = randn(5, deg + 1)
fδ = randn(5, 2)
FR.interp_face!(fδ, f, ps.ll, ps.lr)

FR.standard_lagrange(ps.xpl)


text:
module TestMax
    using StreamStats
    using Distributions
    using Base.Test

    # Max of uniform draws
    for n in rand(1:1_000_000, 100)
        xs = rand(n)
        stat = StreamStats.Max()
        for x in xs
            update!(stat, x)
        end
        online_m = maximum(stat)
        online_ms = state(stat)
        online_n = nobs(stat)
        batch_m = maximum(xs)
        @test online_m == batch_m
        @test online_n == n
    end
end


text:
workspace = pwd()
f = openxl(string(workspace,"/tsp/roadtrip_quebec.xls"))

data_cities = readxl(f, "Donnees!A9:A17")

cities = map(strip, data_cities)

data_distances = readxl(f, "Donnees!C9:K17")

nb_row, nb_col = size(data_distances)
distances = zeros(Int, nb_row, nb_col)
for i = 1:nb_row
    for j = 1:nb_col
        distances[i,j] = Int(data_distances[i,j])
    end
end

save(string(workspace,"/tsp/distances_quebec.jld"),"distances",distances)
save(string(workspace,"/tsp/cities_quebec.jld"),"cities",cities)

f2 = openxl(string(workspace,"/tsp/roadtrip_usa.xls"))

data_cities2 = readxl(f2, "Feuille1!A4:A53")

cities2 = map(strip, data_cities2)

data_distances2 = readxl(f2, "Feuille1!B4:B53")
string_distances2 = map(strip, data_distances2)
n = length(string_distances2)
latitude = zeros(n)
longitude = zeros(n)
for i = 1:n
	string_distances2[i] = string_distances2[i][2:end-1]
	latlon = split(string_distances2[i], ",")
	latitude[i] = parse(Float64, latlon[1])
	longitude[i] = parse(Float64, latlon[2])
end

function distance_earth(lat1, lon1, lat2, lon2)
    # approximate radius of earth in km
    R = 6371.0

    rlat1 = 2 * pi * lat1 / 360
    rlon1 = 2 * pi * lon1 / 360
    rlat2 = 2 * pi * lat2 / 360
    rlon2 = 2 * pi * lon2 / 360

    dlon = (lon2 - lon1) * 2 * pi / 360
    dlat = (lat2 - lat1) * 2 * pi / 360

    a = sin(dlat/2)^2 + cos(rlat1)*cos(rlat2)*sin(dlon/2)^2
    c = 2 * atan(sqrt(a), sqrt(1 - a))

    return R * c
end

distances2 = zeros(n,n)
for i = 1:n
	for j = 1:n
		distances2[i,j] = distance_earth(latitude[i], longitude[i], latitude[j], longitude[j])
	end
end

save(string(workspace,"/tsp/cities_usa.jld"),"cities2",cities2)
save(string(workspace,"/tsp/distances_usa.jld"),"distances2",distances2)


text:
##############################################################################################
module TestL96
##############################################################################################
##############################################################################################
# imports and exports
using DataAssimilationBenchmarks.DeSolvers, DataAssimilationBenchmarks.L96

##############################################################################################
##############################################################################################
# Test the derivative function for known behavior with Euler method,
# the initial condiiton of zeros give back h * F in all components

function EMZerosStep()
    # step size
    h = 0.01

    # forcing parameter
    F = 8.0
    dx_params = Dict{String, Array{Float64}}("F" => [F])

    # initial conditions and arguments
    x = zeros(40)

    # parameters to test
    kwargs = Dict{String, Any}(
        "h" => h,
        "diffusion" => 0.0,
        "dx_params" => dx_params,
        "dx_dt" => L96.dx_dt,
        )

    # em_step! writes over x in place
    em_step!(x, 0.0, kwargs)

    # evaluate test pass/fail if the vector of x is equal to (f*h) in every instance
    if sum(x .== (F*h)) == 40
        true
    else
        false
    end

end

##############################################################################################
# Test the derivative function for known behavior with Euler method,
# the vector with all components equal to F is a fixed point for the system

function EMFStep()
    # step size
    h = 0.01

    # forcing parameter
    F = 8.0
    dx_params = Dict{String, Array{Float64}}("F" => [F])

    # initial conditions and arguments
    x = ones(40)
    x = x * F

    # parameters to test
    kwargs = Dict{String, Any}(
        "h" => h,
        "diffusion" => 0.0,
        "dx_params" => dx_params,
        "dx_dt" => L96.dx_dt,
        )

    # em_step! writes over x in place
    em_step!(x, 0.0, kwargs)

    # evaluate test pass/fail if the vector of x is equal to (f*h) in every instance
    if sum(x .== (F)) == 40
        true
    else
        false
    end

end

##############################################################################################
# end module

end

max_stars_repo_name: sukhreens/DataAssimilationBenchmarks

max_issues_repo_name: sukhreens/DataAssimilationBenchmarks

max_forks_repo_name: sukhreens/DataAssimilationBenchmarks


text:
#
# Convex Hull
#
immutable ConvexHull{T}
    extremepoints::Vector{Point{2, T}}

    # TODO: Add inner constructor
end

length(ch::ConvexHull) = length(ch.extremepoints)
isequal(ch1::ConvexHull, ch2::ConvexHull) = ch1.extremepoints == ch2.extremepoints

#
# Interval
#
immutable Interval
    lb
    ub

    function Interval(lb, ub)
        if lb<ub
            return new(lb, ub)
        else
            return new(ub, lb)
        end
    end
end

in(x, i::Interval) = (i.lb <= x <= i.ub) ? true : false

#
# Line Segment
#
immutable LineSegment{T}
    # Points
    p1::Point{2, T}
    p2::Point{2, T}

    # Slope
    _slope::Float64
    _norm::Float64
end

function LineSegment(p1, p2)
    x1 = p1[1]
    y1 = p1[2]
    x2 = p2[1]
    y2 = p2[2]

    _slope = (y2 - y1)/(x2 - x1)
    _norm = sqrt((y2-y1)*(y2-y1) + (x2-x1)*(x2-x1))

    return LineSegment(p1, p2, _slope, _norm)
end

isfinite(l::LineSegment) = all((isfinite(l.p1[1]), isfinite(l.p1[2]), isfinite(l.p2[1]), isfinite(l.p2[2])))

function evaluatey(l::LineSegment, x::Float64)
    x1, y1 = l.p1
    x2, y2 = l.p2

    if x1 <= x <= x2
        y = y1 + l._slope*(x-x1)
    else
        error("x is not on line segment")
    end

    return y
end

function evaluatex(l::LineSegment, y::Float64)
    y1 = l.p1[2]
    y2 = l.p2[2]

    if y1 < y < y2
        x = l.p1[1] + (y-y1)/l._slope
    else
        error("y is not on line segment")
    end

    return x
end

# TODO: This isn't quite right. It returns distance 0 for points on same line
# (but not in line segment)
function distance(l::LineSegment, p::Point)
    # Pull out points
    x0, y0 = p[1], p[2]
    x1, y1 = l.p1[1], l.p1[2]
    x2, y2 = l.p2[1], l.p2[2]

    # Compute numerator and denominator of Wikipedia expression
    numerator = abs((y2-y1)*x0 - (x2-x1)*y0 + x2*y1 - y2*x1)
    denominator = l._norm

    return numerator/denominator
end

#
# Quadrant
#
immutable Quadrant{T}
    # The origin of quadrant (x, y)
    origin::Point{2, T}

    # The end point in x and y directions (should be Inf or -Inf)
    lx::LineSegment
    ly::LineSegment

    # Functions which say whether points are in correct piece of quadrant
    dir_x::Function
    dir_y::Function

end

function Quadrant(origin::Point, _end_x, _end_y)
    # Get origin points
    x = origin[1]
    y = origin[2]

    # If end is positive infinity then want to check right of origin
    # TODO: Come back and allow this to be Inf -- The intersection
    # method for line segments can't handle this currently
    if _end_x == Inf
        dir_x(p) = p[1] >= x
        lx = LineSegment(origin, Point(1e12, y))
    elseif _end_x == -Inf
        dir_x(p) = p[1] <= x
        lx = LineSegment(Point(-1e12, y), origin)
    end

    # If end is positive infinity then want to check above origin
    if _end_y == Inf
        dir_y(p) = p[2] >= y
        ly = LineSegment(origin, Point(x, Inf))
    elseif _end_y == -Inf
        dir_y(p) = p[2] <= y
        ly = LineSegment(Point(x, -Inf), origin)
    end

    return Quadrant(origin, lx, ly, dir_x, dir_y)
end

in(x::Point, q::Quadrant) = (q.dir_x(x) && q.dir_y(x)) ? true : false


text:
import TaskDAGRecorders
const DAGr = TaskDAGRecorders

function dac(depth = 2)
    function f(n)
        if n <= 1
            sleep(0.01)
        else
            DAGr.@sync begin
                DAGr.@spawn f(n ÷ 2)
                f(n ÷ 2)
            end
        end
    end
    f(2 ^ depth)
end

function dac4(depth = 2)
    function f(n)
        if n <= 1
            sleep(0.01)
        else
            DAGr.@sync begin
                DAGr.@spawn f(n ÷ 4)
                DAGr.@spawn f(n ÷ 4)
                DAGr.@spawn f(n ÷ 4)
                f(n ÷ 4)
            end
        end
    end
    f(4 ^ depth)
end


text:
using Immerse, Gtk

x = range(0,stop=4pi,length=101)
y = sin.(x)

figure()
p = plot(x=x,y=y,Geom.point(tag=:dots),Geom.line(tag=:line))
hfig = display(p)
hit((hfig,:dots), @guarded (figtag, i, event, dist) -> if dist < 1 Immerse.circle_center(figtag, i) end)
hit((hfig,:line), @guarded (figtag, index, event, dist) -> begin
    if dist < 2
        println("You clicked on line ", index[1], " in segment ", index[2], " at $(round(Int,100*index[3]))% along the segment")
    end
end)


text:
struct Euc2DFunction
    components  ::Dict{Int, Point2D}
    rounded     ::Bool
end

struct Euc3DFunction
    components  ::Dict{Int, Point3D}
    rounded     ::Bool
end

struct Euc2DObjective<:ObjectiveFunction
    f           ::Euc2DFunction
    type        ::ObjectiveType
end

struct Euc3DObjective<:ObjectiveFunction
    f           ::Euc3DFunction
    type        ::ObjectiveType
end

Euc2DObjective(map::Dict{Int, Point2D}, rounded::Bool,type::ObjectiveType) = Euc2DObjective(Euc2DFunction(map,rounded),type)
Euc3DObjective(map::Dict{Int, Point3D}, rounded::Bool,type::ObjectiveType) = Euc3DObjective(Euc3DFunction(map,rounded),type)

"""
    evaluateDecision(f::Euc2DFunction, i1::Int, i2::Int)

Takes in a decision and returns the impact on the function

# Arguments
- `f::Euc2DFunction`: the function to be evaluated
- `i1::Int`: the index of the first point
- `i2::Int`: the index of the second point
"""
function evaluateDecision(f::Euc2DFunction, i1::Int, i2::Int)
    if f.rounded
        return euclidean2DRounded(i1::Int,i2::Int, f.components)
    else
        return euclidean2D(i1::Int,i2::Int, f.components)
    end
end

"""
    evaluateDecision(f::Euc3DFunction, i1::Int, i2::Int)

Takes in a decision and returns the impact on the function

# Arguments
- `f::Euc3DFunction`: the function to be evaluated
- `i1::Int`: the index of the first point
- `i2::Int`: the index of the second point
"""
function evaluateDecision(f::Euc3DFunction, i1::Int, i2::Int)
    if f.rounded
        return euclidean3DRounded(i1::Int,i2::Int, f.components)
    else
        return euclidean3D(i1::Int,i2::Int, f.components)
    end
end

"""
    evaluateDecision(o::Union{Euc2DObjective,Euc3DObjective}, i1::Int, i2::Int)

Takes in a decision and returns the impact on the objective

# Arguments
- `f::Union{Euc2DObjective,Euc3DObjective}`: the objective to be evaluated
- `i1::Int`: the index of the first point
- `i2::Int`: the index of the second point
"""
function evaluateDecision(o::Union{Euc2DObjective,Euc3DObjective}, i1::Int, i2::Int)
    evaluateDecision(o.f, i1, i2)
end

"""
    evaluateDecisions(o::Union{Euc2DObjective,Euc3DObjective}, decisions::Vector{T})where{T<:Int}

Takes in a sequence of decsions and returns the value of that sequence
Returns: <:Real

# Arguments
- `f::Union{Euc2DObjective,Euc3DObjective}`: the objective to be evaluated
- `decisions::Vector{T}`: the index of the first point
"""
function evaluateDecisions(o::Union{Euc2DObjective,Euc3DObjective}, decisions::Vector{T})where{T<:Int}
    value = 0
    for i in 1:length(decisions)-1
        value += evaluateDecision(o.f, decisions[i], decisions[i+1])
    end
    return value
end


text:
struct Output
    exitcode::Int
    stdout::String
    stderr::String
end

"""
    apply(fns, obj)
    apply(fns) -> Function
    apply(fn::Function, nt::NamedTuple) -> NamedTuple

Apply function `fn` or functions `fns` to object.
The functions are applied in order, unlike the behaviour of function composition.
Also defines partial function.
(For partial declarations in `Base`, see issue #35052 or `endswith(suffix)`.)
"""
apply(fns, obj) = ∘(reverse(fns)...)(obj)
apply(fns) = obj -> apply(fns, obj)

apply(fn::Function, nt::NamedTuple)::NamedTuple = (; Dict([(t[1], string(t[2])) for t in zip(keys(nt), nt)])...)

"""
    dirparent(path)::String
    dirparent(path, n)::String

Returns the parent or `n`-th parent directory for `path`, where `path` can be a file or directory.

```@example
dirparent("/a/b/c")
```
"""
dirparent(path)::String = splitdir(endswith(path, '/') ? path[1:end-1] : path)[1]
dirparent(path, n)::String = ∘(repeat([dirparent], n)...)(path)

"""
    has_duplicates(A::AbstractArray)::Bool

Returns whether `A` contains duplicates.
"""
has_duplicates(A::AbstractArray) = length(A) != length(unique(A))

"""
    map_by_df(a::Array, df::DataFrame, from::Symbol, to::Symbol; missing=nothing)::Array

Return array `A` where all elements are mapped `from` `U` `to` `V`.
Leaving all elements of `A` for which no match is found unchanged.

!! Note-to-self. This is dumb function. Use joins instead.
"""
function map_by_df(A::AbstractArray, df::DataFrame, from::Symbol, to::Symbol)::Array
    U = df[:, from]

    if eltype(A) != eltype(U)
        @warn TypeError(:map_by_df, "trying to map `A` with `from`", typeof(A), typeof(U))
    end
    function map_element(e)
        if ismissing(e)
            return missing
        else
            filtered = filter(row -> row[from] == e, df)
            nrow(filtered) == 1 ? first(filtered)[to] : e
        end
    end
    map(map_element, A)
end

rmextension(path::AbstractString) = string(splitext(path)[1])::String

"""
    rescale(a, a_l, a_u, b_l, b_u)::Number

Apply feature scaling to `a` from the range `[a_l, a_u]` to the range `[b_l, b_u]`.
"""
rescale(a, a_l, a_u, b_l, b_u) = b_l + ( ((a - a_l)*(b_u - b_l)) / (a_u - a_l) )

"""
    nrow_per_group(df::DataFrame, group::Symbol; col1="group", col2="nrow")::DataFrame

Return the group name and the number of rows per group in `df`.
"""
function nrow_per_group(df, group::Symbol; col1="group", col2="nrow")::DataFrame
	grouped = groupby(df, group)
	out = DataFrame([(col1 = first(k), col2 = nrow(g)) for (k, g) in zip(keys(grouped), grouped)])
    select(out, :col1 => col1, :col2 => col2)
end

"""
    output(f::Function) -> Output
    output(cmd::Cmd) -> Output

Evaluates `f` of type `f(out::String, err::String)::CmdRedirect` or `cmd::Cmd`.
"""
function output(f::Function)::Output
    # Don't need fancy live scrolling log because it runs in CI anyway.
    out = IOBuffer()
    err = IOBuffer()
    take_str!(io) = String(take!(io))
    exitcode = 99
    try
        process = run(f(out, err))
        exitcode = process.exitcode
    catch
        @error "Error while running function. Printing stdout and stderr:"
        out_s = take_str!(out)
        println("stdout = ")
        println.(split(out_s, '\n'))
        err_s = take_str!(err)
        println("stderr = ")
        println.(split(err_s, '\n'))
        return Output(-1, out_s, err_s)
    end

    out_s = take_str!(out)
    err_s = take_str!(err)
    Output(exitcode, out_s, err_s)
end

output(cmd::Cmd)::Output = output((out, err) -> pipeline(cmd, stdout=out, stderr=err))

nt2dict(nt::NamedTuple)::Dict = Dict(zip(keys(nt), nt))

function graded()
    q = [:P1, :P2, :P3, :P4, :P5, :M1, :M2, :M3, :M4, :M5, :C1, :Penalty]
    s = [
        0.5, # P1 [0.5]
        0.75, # P2 [0.75]
        0.5, # P3 [0.5]
        1, # P4 [1.5]
        1, # P5 [1.0]
        0.5, # M1 [0.5]
        1, # M2 [1.5]
        1, # M3 [1.0]
        0.25, # M4 [0.25]
        1.5, # M5 [1.5]
        1, # C1 [1.0]
        0, # Penalty [-0.5]
    ]
    max = [0.5, 0.75, 0.5, 1.5, 1.0, 0.5, 1.5, 1.0, 0.25, 1.5, 1.0, -0.5]

    df = DataFrame(; q, s, max)
    @show df
    println()
    @show sum(df.s)
    nothing
end

"""
    split_data(df::DataFrame, group_col, data_col, groups)::Tuple

Return a NamedTuple containing data from `data_col` for each group in `groups`.
The ordering of the data is the same as the ordering of `groups`.
"""
function split_data(df::DataFrame, group_col, data_col, groups)::Tuple
    typeassert.(groups, eltype(df[:, group_col]))
    return Tuple(filter(group_col => ==(g), df)[:, data_col] for g in groups)
end


text:
function init_constraint_struct(com, set::ComplementSet{F,S}, internals) where {F,S}
    inner_set = set.set
    fct = internals.fct
    if !(internals.fct isa SAF) && F <: SAF
        fct = get_saf(internals.fct)
    end

    constraint = get_constraint(com, fct, inner_set)
    complement_constraint = get_complement_constraint(com, constraint)
    return complement_constraint
end


text:
# Stiff 1
using ParameterizedFunctions, OrdinaryDiffEq, Plots
plotly()
rober = @ode_def Rober begin
  dy₁ = -k₁*y₁+k₃*y₂*y₃
  dy₂ =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
  dy₃ =  k₂*y₂^2
end k₁ k₂ k₃
prob = ODEProblem(rober,[1.0;0.0;0.0],(0.0,1e11),(0.04,3e7,1e4))
sol = solve(prob, Rosenbrock23())
plot(sol, xscale=:log10, tspan=(1e-6, 1e5), layout=(3,1))

#######################################
#######################################
#######################################
# Stiff 2

using OrdinaryDiffEq, Plots
function orego(du,u,p,t)
  s,q,w = p
  y1,y2,y3 = u
  du[1] = s*(y2+y1*(1-q*y1-y2))
  du[2] = (y3-(1+y1)*y2)/s
  du[3] = w*(y1-y3)
end
p = [77.27,8.375e-6,0.161]
prob = ODEProblem(orego,[1.0,2.0,3.0],(0.0,360.0),p)
sol = solve(prob, Rodas5())
plot(sol)
# phase space

#######################################
#######################################
#######################################
# Stiff SDE

using StochasticDiffEq
function orego(du,u,p,t)
  s,q,w = p
  y1,y2,y3 = u
  du[1] = s*(y2+y1*(1-q*y1-y2))
  du[2] = (y3-(1+y1)*y2)/s
  du[3] = w*(y1-y3)
end
function g(du,u,p,t)
  du[1] = 0.1u[1]
  du[2] = 0.1u[2]
  du[3] = 0.1u[3]
end
p = [77.27,8.375e-6,0.161]
prob = SDEProblem(orego,g,[1.0,2.0,3.0],(0.0,30.0),p)
sol = solve(prob,SOSRI())
plot(sol)

#######################################
#######################################
#######################################
# DAE

function f(out, da, a, p, t)
   L1, m1, L2, m2, g = p

   u1, v1, x1, y1, T1,
   u2, v2, x2, y2, T2 = a

   du1, dv1, dx1, dy1, dT1,
   du2, dv2, dx2, dy2, dT2 = da

   out[1]  = x2*T2/(m2*L2) - du2
   out[2]  = y2*T2/(m2*L2) - g - dv2
   out[3]  = u2 - dx2
   out[4]  = v2 - dy2
   out[5]  = u2^2 + v2^2 -y2*g + T2/m2

   out[6]  = x1*T1/(m1*L1) - x2*T2/(m2*L2) - du1
   out[7]  = y1*T1/(m1*L1) - g - y2*T2/(m2*L2) - dv1
   out[8]  = u1 - dx1
   out[9]  = v1 - dy1
   out[10] = u1^2 + v1^2 + T1/m1 +
                (-x1*x2 - y1*y2)/(m1*L2)*T2 - y1*g
   nothing
end

# Release pendulum from top right
u0 = zeros(10)
u0[3] = 1.0
u0[8] = 1.0
du0 = zeros(10)
du0[2] = 9.8
du0[7] = 9.8

p = [1,1,1,1,9.8]
tspan = (0.,100.)

differential_vars = [true, true, true, true, false,
                     true, true, true, true, false]
prob = DAEProblem(f, du0, u0, tspan, p, differential_vars = differential_vars)
using Sundials
sol = solve(prob, IDA())
plot(sol, vars=((x1,y1,x2,y2)->(x1, y1),3,4,8,9), lab="")
plot!(sol, vars=((x1,y1,x2,y2)->(x1+x2, y1+y2),3,4,8,9), xlims=(-2.2, 2.2), aspect_ratio=1, lab="")

#######################################
#######################################
#######################################
# DDE
using DiffEqProblemLibrary.DDEProblemLibrary, LinearAlgebra
using Plots
using DelayDiffEq
DDEProblemLibrary.importddeproblems()
prob = DDEProblemLibrary.prob_dde_RADAR5_waltman_5
sol = solve(prob, MethodOfSteps(Rosenbrock23(linsolve=LinSolveFactorize(lu))); reltol = 1e-6, abstol = [1e-18, 1e-18, 1e-18, 1e-18, 1e-6, 1e-6])
plot(sol)

#######################################
#######################################
#######################################
# Stiff PDE

using OrdinaryDiffEq, Sundials, Plots

# initial condition
function init_brusselator_2d(xyd)
    N = length(xyd)
    u = zeros(N, N, 2)
    for I in CartesianIndices((N, N))
        x = xyd[I[1]]
        y = xyd[I[2]]
        u[I,1] = 22*(y*(1-y))^(3/2)
        u[I,2] = 27*(x*(1-x))^(3/2)
    end
    u
end

N = 32

xyd_brusselator = range(0,stop=1,length=N)

u0 = vec(init_brusselator_2d(xyd_brusselator))

tspan = (0, 22.)

p = (3.4, 1., 10., xyd_brusselator)

brusselator_f(x, y, t) = ifelse((((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) &&
                                (t >= 1.1), 5., 0.)
function brusselator_2d_loop(du, u, p, t)
    A, B, α, xyd = p
    dx = step(xyd)
    N = length(xyd)
    α = α/dx^2
    limit = a -> let N=N
        a == N+1 ? 1 :
        a == 0 ? N :
        a
    end
    II = LinearIndices((N, N, 2))

    @inbounds begin
        for I in CartesianIndices((N, N))
            x = xyd[I[1]]
            y = xyd[I[2]]
            i = I[1]
            j = I[2]
            ip1 = limit(i+1)
            im1 = limit(i-1)
            jp1 = limit(j+1)
            jm1 = limit(j-1)

            ii1 = II[i,j,1]
            ii2 = II[i,j,2]

            du[II[i,j,1]] = α*(u[II[im1,j,1]] + u[II[ip1,j,1]] + u[II[i,jp1,1]] + u[II[i,jm1,1]] - 4u[ii1]) +
            B + u[ii1]^2*u[ii2] - (A + 1)*u[ii1] + brusselator_f(x, y, t)

            du[II[i,j,2]] = α*(u[II[im1,j,2]] + u[II[ip1,j,2]] + u[II[i,jp1,2]] + u[II[i,jm1,2]] - 4u[II[i,j,2]]) +
            A*u[ii1] - u[ii1]^2*u[ii2]
        end
    end
    nothing
end
gr()
function plot_sol(sol)
    off = N^2
    for t in sol.t[1]:0.1:sol.t[end]
        solt = sol(t)
        plt1 = surface(reshape(solt[1:off], N, N), zlims=(0, 5), leg=false)
        surface!(plt1, reshape(solt[off+1:end], N, N), zlims=(0, 5), leg=false)
        display(plt1)
        sleep(0.05)
    end
    nothing
end

prob2 = ODEProblem(brusselator_2d_loop, u0, tspan, p)

using Sundials
sol2 = @time solve(prob2, CVODE_BDF(linear_solver=:GMRES))

plot_sol(sol2)

max_stars_repo_name: Runfa-Zhang/MeetUpMaterials

max_issues_repo_name: Runfa-Zhang/MeetUpMaterials

max_forks_repo_name: Runfa-Zhang/MeetUpMaterials


text:
export extended_source_modeling

# Setup time-domain linear or nonlinear foward and adjoint modeling and interface to OPESCI/devito
function extended_source_modeling(model_full::Model, srcData, recGeometry, recData, weights, dm, op::Char, mode::Int64, options)

    # Load full geometry for out-of-core geometry containers
    recGeometry = Geometry(recGeometry)

    model = model_full

    # Set up Python model structure
    modelPy = devito_model(model, options; dm=dm)
    if op=='J' && mode == 1
        if modelPy.dm == 0 && options.return_array == false
            return judiVector(recGeometry, zeros(Float32, recGeometry.nt[1], length(recGeometry.xloc[1])))
        elseif modelPy.dm == 0 && options.return_array == true
            return vec(zeros(Float32, recGeometry.nt[1], length(recGeometry.xloc[1])))
        end
    end

    # Load shot record if stored on disk
    typeof(recData) == SegyIO.SeisCon && (recData = convert(Array{Float32,2}, recData[1].data))

    # Remove receivers outside the modeling domain (otherwise leads to segmentation faults)
    recGeometry, recData = remove_out_of_bounds_receivers(recGeometry, recData, model)

    isnothing(weights) ? nothing : weights = pad_array(weights, pad_sizes(model, options; so=0); mode=:zeros)
    # Devito interface
    argout = devito_interface(modelPy, model, srcData, recGeometry, recData, weights, dm, options)

    # Extend gradient back to original model size
    if op=='J' && mode==-1 && options.limit_m==true
        argout = extend_gradient(model_full, model, argout)
    end

    return argout
end

# Function instance without options
extended_source_modeling(model::Model, srcData, recGeometry, recData,  weights, perturbation, srcnum::Int64, op::Char, mode::Int64) =
    extended_source_modeling(model, srcData, recGeometry, recData, weights, perturbation, srcnum, op, mode, Options())

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
export hours

function hours{T<:DateTime}(t::T, n::Int; direction::Symbol=:forward, shift::Int=0)
 forward = isforward(direction)
 nextfun = forward ? t->t+Dates.Hour(1) : t->t-Dates.Hour(1)
 jumpstart = n>0 && forward ? t-Dates.Hour(1-shift) :
             n>0 && !forward ? t+Dates.Hour(n+shift) :
             n<=0 && forward ? t+Dates.Hour(n-1+shift) :
             t+Dates.Hour(shift) # <-- n<=0 && !forward
 ManyTimesIter(jumpstart, abs(n), nextfun, forward, Dates.Hour(1))
end

hours{T<:DateTime}(t::T, nd::Base.Dates.Day; direction::Symbol=:forward, shift::Int=0) =
    hours(t, nd.value*24, direction=direction)

hours{T<:DateTime}(t::T, nd::Base.Dates.Hour; direction::Symbol=:forward, shift::Int=0) =
    hours(t, nd.value, direction=direction)

hours{T<:Date}(d::T, n::Int; direction::Symbol=:forward, shift::Int=0) =
    hours(DateTime(d), n, direction=direction)

hours{T<:Date}(d::T, nd::Base.Dates.Day; direction::Symbol=:forward, shift::Int=0) =
    hours(DateTime(d), nd.value*24, direction=direction)

hours{T<:Date}(d::T, nd::Base.Dates.Hour; direction::Symbol=:forward, shift::Int=0) =
    hours(DateTime(d), nd.value, direction=direction)

function hours{T<:DateTime}(t1::T, t2::T; shift::Int=0)
 if t2 > t1
    tstart = t1
    direction = :forward
 else
    tstart = t2
    direction = :backward
 end
 n = abs( Int( Dates.Hour( t2 - t1 )))
 hours(tstart, n, direction=direction, shift=shift)
end

hours{ T1<:DateTime, T2<:Date }(t1::T1, d2::T2; shift::Int=0) =
    hours(t1, DateTime(d2), shift=shift)

hours{ T1<:Date, T2<:DateTime }(d1::T1, t2::T2; shift::Int=0) =
    hours(DateTime(d1), t2, shift=shift)

hours{ T1<:Date, T2<:Date }(d1::T1, d2::T2; shift::Int=0) =
    hours(DateTime(d1), DateTime(d2), shift=shift)

"""
Creating hours like days(). Except using DateTime(), not Date()
for representing a each hour in parameters.

Many.hours(DateTime(2016,08,15,06),3)

From tuples:

hours((2016,08,15,06),3)

hours((2016,08,15,06),(2016,08,15,10))

"""
hours(h1::Tuple, args...) =
 hours(DateTime(h1...), args...)

hours(h1::Tuple, h2::Tuple, args...) =
 hours(DateTime(h1...), DateTime(h2...), args...)

max_stars_repo_name: closescreen/Many

max_issues_repo_name: closescreen/Many

max_forks_repo_name: closescreen/Many


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayHttpListenerPropertiesFormat <: SwaggerModel
    frontendIPConfiguration::Any # spec type: Union{ Nothing, SubResource } # spec name: frontendIPConfiguration
    frontendPort::Any # spec type: Union{ Nothing, SubResource } # spec name: frontendPort
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol
    hostName::Any # spec type: Union{ Nothing, String } # spec name: hostName
    sslCertificate::Any # spec type: Union{ Nothing, SubResource } # spec name: sslCertificate
    requireServerNameIndication::Any # spec type: Union{ Nothing, Bool } # spec name: requireServerNameIndication
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function ApplicationGatewayHttpListenerPropertiesFormat(;frontendIPConfiguration=nothing, frontendPort=nothing, protocol=nothing, hostName=nothing, sslCertificate=nothing, requireServerNameIndication=nothing, provisioningState=nothing)
        o = new()
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("frontendIPConfiguration"), frontendIPConfiguration)
        setfield!(o, Symbol("frontendIPConfiguration"), frontendIPConfiguration)
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("frontendPort"), frontendPort)
        setfield!(o, Symbol("frontendPort"), frontendPort)
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("hostName"), hostName)
        setfield!(o, Symbol("hostName"), hostName)
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("sslCertificate"), sslCertificate)
        setfield!(o, Symbol("sslCertificate"), sslCertificate)
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("requireServerNameIndication"), requireServerNameIndication)
        setfield!(o, Symbol("requireServerNameIndication"), requireServerNameIndication)
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ApplicationGatewayHttpListenerPropertiesFormat

const _property_map_ApplicationGatewayHttpListenerPropertiesFormat = Dict{Symbol,Symbol}(Symbol("frontendIPConfiguration")=>Symbol("frontendIPConfiguration"), Symbol("frontendPort")=>Symbol("frontendPort"), Symbol("protocol")=>Symbol("protocol"), Symbol("hostName")=>Symbol("hostName"), Symbol("sslCertificate")=>Symbol("sslCertificate"), Symbol("requireServerNameIndication")=>Symbol("requireServerNameIndication"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ApplicationGatewayHttpListenerPropertiesFormat = Dict{Symbol,String}(Symbol("frontendIPConfiguration")=>"SubResource", Symbol("frontendPort")=>"SubResource", Symbol("protocol")=>"String", Symbol("hostName")=>"String", Symbol("sslCertificate")=>"SubResource", Symbol("requireServerNameIndication")=>"Bool", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ ApplicationGatewayHttpListenerPropertiesFormat }) = collect(keys(_property_map_ApplicationGatewayHttpListenerPropertiesFormat))
Swagger.property_type(::Type{ ApplicationGatewayHttpListenerPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ApplicationGatewayHttpListenerPropertiesFormat[name]))}
Swagger.field_name(::Type{ ApplicationGatewayHttpListenerPropertiesFormat }, property_name::Symbol) =  _property_map_ApplicationGatewayHttpListenerPropertiesFormat[property_name]

const _allowed_ApplicationGatewayHttpListenerPropertiesFormat_protocol = ["Http", "Https"]

function check_required(o::ApplicationGatewayHttpListenerPropertiesFormat)
    true
end

function validate_property(::Type{ ApplicationGatewayHttpListenerPropertiesFormat }, name::Symbol, val)
    if name === Symbol("protocol")
        Swagger.validate_param(name, "ApplicationGatewayHttpListenerPropertiesFormat", :enum, val, _allowed_ApplicationGatewayHttpListenerPropertiesFormat_protocol)
    end
end


text:
module WindMixing

export data, read_les_output,
       animate_prediction,
       mse,
       predict_flux, predict_NDE,
       write_metadata_NN_training, write_data_NN_training, write_data_NN,
       write_metadata_NDE_training, write_data_NDE_training,
       train_NDE, train_NN,
       oceananigans_modified_pacanowski_philander_nn,
       solve_NDE_mutating, solve_oceananigans_modified_pacanowski_philander_nn,
       NDE_profile, NDE_profile_oceananigans,
       column_model_1D_kpp,
       animate_NN, animate_profile, animate_flux, animate_profile_flux, animate_profiles, animate_local_richardson_profile,
       animate_profiles_fluxes, animate_training_data_profiles_fluxes, animate_profiles_fluxes_comparison,
       animate_training_results, animate_training_results_oceananigans,
       local_richardson,
       smoothing_filter,
       loss, loss_gradient,
       optimise_modified_pacanowski_philander,
       extract_NN, extract_parameters_modified_pacanowski_philander_optimisation,
       diurnal_fluxes,
       directories

using Flux, Plots
using Oceananigans.Grids: Center, Face
using Oceananigans: OceananigansLogger
using Oceananigans
using OceanParameterizations
import OceanTurb
using JLD2
using FileIO
using BSON
using OrdinaryDiffEq, DiffEqSensitivity
using GalacticOptim
using Statistics
using Random
using Logging
using CairoMakie
using Printf
using CUDA
using LinearAlgebra

mse(x::Tuple{Array{Float64,2}, Array{Float64,2}}) = Flux.mse(x[1], x[2])
mse(x::Tuple{Array{Float32,2}, Array{Float64,2}}) = Flux.mse(Float64.(x[1]), x[2])
mse(x::Tuple{Array{Float32,2}, Array{Float32,2}}) = Flux.mse(x[1], x[2])

include("lesbrary_data.jl")
include("data_containers.jl")
include("NDE_training.jl")
include("NN_training.jl")
include("animation.jl")
include("data_writing.jl")
include("filtering_operators.jl")
include("diffusivity_parameter_optimisation.jl")
include("data_extraction.jl")
include("NDE_oceananigans.jl")
include("training_postprocessing.jl")
include("k_profile_parameterization.jl")
include("loss.jl")
include("plotting.jl")

function __init__()
    Logging.global_logger(OceananigansLogger())
end

end


text:
abstract type Trajectory end

function Trajectory( filename :: String,
                     solute :: Selection, solvent :: Selection;
                     format :: String = "")
  if format == "dcd" || FileOperations.file_extension(filename) == "dcd"
    NamdDCD(filename,solute,solvent)
  elseif format == "PDBTraj"
    PDBTraj(filename,solute,solvent)
  else
    ChemFile(filename,solute,solvent,format=format)
  end
end

# If only one selection is provided, assume that the solute and the
# solvent are the same

function Trajectory( filename :: String,
                     solvent :: Selection;
                     format :: String = "")
  Trajectory(filename,solvent,solvent,format=format)
end


text:
abstract type BibliographyBlock <: Expanders.ExpanderPipeline end

Selectors.order(::Type{BibliographyBlock}) = 12.0  # Expand bibliography last
Selectors.matcher(::Type{BibliographyBlock}, node, page, doc) = Expanders.iscode(node, r"^@bibliography")

const tex2unicode_replacements = (
    "---" => "—", # em dash needs to go first
    "--"  => "–",
    r"\\`\{(\S)\}" => s"\1\u300", # \`{o} 	ò 	grave accent
    r"\\'\{(\S)\}" => s"\1\u301", # \'{o} 	ó 	acute accent
    r"\\\^\{(\S)\}" => s"\1\u302", # \^{o} 	ô 	circumflex
    r"\\~\{(\S)\}" => s"\1\u303", # \~{o} 	õ 	tilde
    r"\\=\{(\S)\}" => s"\1\u304", # \={o} 	ō 	macron accent (a bar over the letter)
    r"\\u\{(\S)\}" => s"\1\u306",  # \u{o} 	ŏ 	breve over the letter
    r"\\\.\{(\S)\}" => s"\1\u307", # \.{o} 	ȯ 	dot over the letter
    r"\\\\\"\{(\S)\}" => s"\1\u308", # \"{o} 	ö 	umlaut, trema or dieresis
    r"\\r\{(\S)\}" => s"\1\u30A",  # \r{a} 	å 	ring over the letter (for å there is also the special command \aa)
    r"\\H\{(\S)\}" => s"\1\u30B",  # \H{o} 	ő 	long Hungarian umlaut (double acute)
    r"\\v\{(\S)\}" => s"\1\u30C",  # \v{s} 	š 	caron/háček ("v") over the letter
    r"\\d\{(\S)\}" => s"\1\u323",  # \d{u} 	ụ 	dot under the letter
    r"\\c\{(\S)\}" => s"\1\u327",  # \c{c} 	ç 	cedilla
    r"\\k\{(\S)\}" => s"\1\u328",  # \k{a} 	ą 	ogonek
    r"\\b\{(\S)\}" => s"\1\u331",  # \b{b} 	ḇ 	bar under the letter

    # TODO:
    # \t{oo} 	o͡o 	"tie" (inverted u) over the two letters
    # \o{} 	ø 	slashed o (o with stroke)
    # {\i} 	ı 	dotless i (i without tittle)
    # \l{} 	ł 	barred l (l with stroke)
)

function tex2unicode(s)
    for replacement in tex2unicode_replacements
        s = replace(s, replacement)
    end
    return Unicode.normalize(s)
end

function Selectors.runner(::Type{BibliographyBlock}, x, page, doc)
    @info "Expanding bibliography."
    raw_bib = "
"
    for (id, entry) in doc.plugins[CitationBibliography].bib
        @info "Expanding bibliography entry: $id."

        # Add anchor that citations can link to from anywhere in the docs.
        Anchors.add!(doc.internal.headers, entry, entry.id, page.build)

        authors = xnames(entry)
        year = xyear(entry)
        link = xlink(entry)
        title = xtitle(entry) |> tex2unicode
        published_in = xin(entry)

        entry_text = """
$id

$authors ($year), $title, $published_in

"""
        raw_bib *= entry_text
    end
    raw_bib *= "\n
"

    page.mapping[x] = Documents.RawNode(:html, raw_bib)
end


text:
mutable struct LDA
    K::Int          # topics
    D::Int          # documents
    N::Vector{Int}  # observations per document modality
    V::Int          # vocab items per modality

    η::Float64              # topic hyperparameter
    λ::Matrix{Float64}      # topic variational parameters
    β::Matrix{Float64}      # topic-word distributions
    Elnβ::Matrix{Float64}   # expectation of ln(β)

    α::Float64                  # doc-topic hyperparameter
    γ::Matrix{Float64}          # doc-topic variational parameters
    θ::Matrix{Float64}          # doc-topic distributions
    Elnθ::Matrix{Float64}       # expectation of ln(θ)
    ϕ::Vector{Matrix{Float64}}  # Z variational parameters

    X::Vector{Matrix{Int}}

    converged::Bool
    elbo::Float64
    ll::Float64

    function LDA(k::Int, α::Float64, η::Float64, V::Int,
            X::Vector{Matrix{Int}})
        model = new()

        model.K = k
        model.α = α
        model.η = η
        model.X = X
        model.D = length(X)
        model.N = [sum(X[d][:, 2]) for d in 1:model.D]
        model.V = V

        model.λ = rand(1:100, model.V, model.K)
        model.β = Array{Float64}(undef, model.V, model.K)
        model.Elnβ = Array{Float64}(undef, model.V, model.K)
        update_Elnβ!(model)

        model.γ = fill(1.0, model.K, model.D)
        model.θ = Array{Float64}(undef, model.K, model.D)
        model.Elnθ = Array{Float64}(undef, model.K, model.D)
        update_Elnθ!(model)

        model.ϕ = [
            fill(1.0 / model.K, model.K, size(model.X[d], 1))
            for d in 1:model.D
        ]

        model.converged = false

        return model
    end
end

function LDA(k::Int, α::Float64, η::Float64, X::Vector{Matrix{Int}})
    D = length(X)
    V = 0
    for d in 1:D
        if size(X[d], 1) > 0
            V = max(V, maximum(X[d][:, 1]))
        end
    end

    return LDA(k, α, η, V, X)
end

function update_ϕ!(model::LDA)
    for d in 1:model.D
        model.ϕ[d] .= exp.(
            model.Elnθ[:, d] .+ model.Elnβ[model.X[d][:, 1], :]'
        )
        model.ϕ[d] ./= sum(model.ϕ[d], dims=1)
    end
end

function update_Elnθ!(model::LDA)
    model.Elnθ .= digamma.(model.γ) .- digamma.(sum(model.γ, dims=1))
end

function update_γ!(model::LDA)
    model.γ .= model.α

    for d in 1:model.D
        model.γ[:, d] .+= model.ϕ[d] * model.X[d][:, 2]
    end

    update_Elnθ!(model)
end

function update_θ!(model::LDA)
    model.θ .= model.γ ./ sum(model.γ, dims=1)
end

function update_Elnβ!(model::LDA)
    model.Elnβ .= digamma.(model.λ) .- digamma.(sum(model.λ, dims=1))
end

function update_λ!(model::LDA)
    model.λ .= model.η

    for d in 1:model.D
        model.λ[model.X[d][:, 1], :] .+= model.ϕ[d]' .* model.X[d][:, 2]
    end

    update_Elnβ!(model)
end

function update_β!(model::LDA)
    model.β .= model.λ ./ sum(model.λ, dims=1)
end

function calculate_ElnPβ(model::LDA)
    lnp = model.K * (lgamma(model.V * model.η) - model.V * lgamma(model.η))
    lnp += (model.η - 1) * sum(model.Elnβ)
    return lnp
end

function calculate_ElnPθ(model::LDA)
    lnp = model.D * (lgamma(model.K * model.α) - model.K * lgamma(model.α))
    lnp += (model.α - 1) * sum(model.Elnθ)
    return lnp
end

function calculate_ElnPZ(model::LDA)
    lnp = 0.0
    for d in 1:model.D
        lnp += sum(model.ϕ[d] .* model.Elnθ[:, d] .* model.X[d][:, 2]')
    end
    return lnp
end

function calculate_ElnPX(model::LDA)
    lnp = 0.0
    for d in 1:model.D
        lnp += sum(model.ϕ[d]' .* model.Elnβ[model.X[d][:, 1], :] .* model.X[d][:, 2])
    end
    return lnp
end

function calculate_ElnQβ(model::LDA)
    lnq = sum(lgamma.(model.λ)) - sum(lgamma.(sum(model.λ, dims=1)))
    lnq -= sum((model.λ .- 1) .* model.Elnβ)
    return lnq
end

function calculate_ElnQθ(model::LDA)
    lnq = sum(lgamma.(model.γ)) - sum(lgamma.(sum(model.γ, dims=1)))
    lnq -= sum((model.γ .- 1) .* model.Elnθ)
    return lnq
end

function calculate_ElnQZ(model::LDA)
    lnq = 0.0
    for d in 1:model.D
        lnq += sum(log.(model.ϕ[d] .^ model.ϕ[d]))
    end
    return lnq
end

function calculate_elbo(model::LDA)
    elbo = 0.0
    elbo += calculate_ElnPβ(model)
    elbo += calculate_ElnPθ(model)
    elbo += calculate_ElnPZ(model)
    elbo += calculate_ElnPX(model)
    elbo -= calculate_ElnQβ(model)
    elbo -= calculate_ElnQθ(model)
    elbo -= calculate_ElnQZ(model)
    return elbo
end

function calculate_loglikelihood(X::Vector{Matrix{Int}}, θ::Matrix{Float64},
                                 β::Matrix{Float64})
    ll = 0.0

    N = 0
    for d in 1:length(X)
        N += sum(X[d][:, 2])
        for w in 1:size(X[d], 1)
            v = X[d][w, 1]
            ll += X[d][w, 2] * log(dot(θ[:, d], β[v, :]))
        end
    end

    return ll / N
end

function calculate_loglikelihood(X::Vector{Matrix{Int}}, model::LDA)
    return calculate_loglikelihood(X, model.θ, model.β)
end

function calculate_loglikelihood(model::LDA)
    return calculate_loglikelihood(model.X, model.θ, model.β)
end

function fit!(model::LDA; maxiter=1000, tol=1e-4, verbose=true)
    ll = Float64[]

    for iter in 1:maxiter
        update_γ!(model)
        update_ϕ!(model)
        update_λ!(model)

        update_β!(model)
        update_θ!(model)

        push!(ll, calculate_loglikelihood(model))

        if verbose
            println("$iter\tLog-likelihood: ", ll[end])
        end

        if length(ll) > 10 && check_convergence(ll, tol=tol)
            model.converged = true
            break
        end
    end
    model.elbo = calculate_elbo(model)
    model.ll = ll[end]

    return ll
end

function unsmoothed_update_ϕ!(model::LDA)
    for d in 1:model.D
        model.ϕ[d] .= exp.(model.Elnθ[:, d]) .* model.β[model.X[d][:, 1], :]'
        model.ϕ[d] ./= sum(model.ϕ[d], dims=1)
    end
end

function transform(model::LDA, X::Vector{Matrix{Int}};
                   maxiter=1000, tol=1e-4, verbose=false)

    newmodel = LDA(model.K, model.α, model.η, X)
    newmodel.β .= model.β

    ll = Float64[]

    for iter in 1:maxiter
        update_γ!(newmodel)
        unsmoothed_update_ϕ!(newmodel)
        update_θ!(newmodel)

        push!(ll, calculate_loglikelihood(newmodel))

        if verbose
            println("$iter\tLog-likelihood: ", ll[end])
        end

        if length(ll) > 10 && check_convergence(ll, tol=tol)
            newmodel.converged = true
            break
        end
    end

    if !newmodel.converged
        warn("transform did not converge")
    end

    return newmodel.θ
end

function fit_heldout(Xheldout::Vector{Matrix{Int}}, model::LDA;
        maxiter=100, verbose=false)

    heldout_model = LDA(model.K, model.α, model.η, Xheldout)
    heldout_model.λ = deepcopy(model.λ)
    heldout_model.β = deepcopy(model.β)
    heldout_model.Elnβ = deepcopy(model.Elnβ)

    ll = Float64[]
    for iter in 1:maxiter
        update_γ!(heldout_model)
        update_ϕ!(heldout_model)

        update_θ!(heldout_model)

        push!(ll, calculate_loglikelihood(Xheldout, heldout_model))

        if verbose
            println("$iter\tLog-likelihood: ", ll[end])
        end

        if length(ll) > 10 && check_convergence(ll)
            heldout_model.converged = true
            break
        end
    end
    heldout_model.elbo = calculate_elbo(heldout_model)
    heldout_model.ll = ll[end]

    return heldout_model
end


text:
# interfacing with NNlib.jl

import NNlib: stride, padding, dilation, flipkernel, spatial_dims, kernel_size,
              conv!, ∇conv_filter!, ∇conv_data!,
              maxpool!, meanpool!, ∇maxpool!, ∇meanpool!,
              softmax, softmax!, ∇softmax, ∇softmax!,
              logsoftmax, logsoftmax!, ∇logsoftmax, ∇logsoftmax!

import DataStructures: DefaultDict

const CUDNNFloat = Union{Float16,Float32,Float64}

# Since CUDNN does not support 1D convolution, Conv in Flux will give a CUDNNError if the size is 1-dimensional.
fix1d(x) = x
fix1d(x::DenseCuArray{T, 3}) where T = reshape(x, size(x, 1), 1, size(x, 2), size(x, 3))
fix1d(cdims::DenseConvDims{1,K,C_in,C_out,S,P,D,F}) where {K,C_in,C_out,S,P,D,F} =
  DenseConvDims{2,(K...,1),C_in,C_out,(S...,1),(P...,0,0),(D...,1),F}((cdims.I...,1))
fix1d(pdims::PoolDims{1,K,S,P,D}) where {K,S,P,D,F} =
  PoolDims{2,(K...,1),(S...,1),(P...,0,0),(D...,1)}((pdims.I..., 1), pdims.C_in)

# We have to reshape the CuArray/PoolDims/DenseConvDims to 4D before feeding to CUDNN.
reshape4D(x::AbstractVector) = reshape(x, 1, 1, length(x), 1)
reshape4D(x::AbstractMatrix) = reshape(x, 1, 1, size(x)...)
reshape4D(x::AbstractArray{T,3}) where T = reshape(x, size(x, 1), 1, size(x, 2), size(x, 3))
reshape4D(x::AbstractArray{T}) where T = x

workspacesize(x) = min(Mem.info()[1] ÷ 16, sizeof(x) * 2)

function perfChoose(perfResults, returnedAlgoCount)::UInt32
  if perfResults[1].status != 0
    return 0
  else
    (best_algo,best_time,best_memory) = (perfResults[1].algo,perfResults[1].time,perfResults[1].memory)
    for i = 2:returnedAlgoCount
      if perfResults[i].status == 0 && perfResults[i].memory < best_memory && perfResults[i].time < best_time * 1.1
        (best_algo,best_memory) = (perfResults[i].algo,perfResults[i].memory)
      end
    end
    return best_algo
  end
end

# Softmax

# in-place for x or dy
softmax(x::DenseCuArray{T}; dims=1) where T<:CUDNNFloat =
  softmax!(x, x, dims=dims)

∇softmax(dy::DenseCuArray{T}, x::DenseCuArray{T}; dims=1) where T<:CUDNNFloat =
  ∇softmax!(dy, dy, x, dims=dims)

logsoftmax(x::DenseCuArray{T}; dims=1) where T<:CUDNNFloat =
  logsoftmax!(x, x, dims=dims)

∇logsoftmax(dy::DenseCuArray{T}, x::DenseCuArray{T}; dims=1) where T<:CUDNNFloat =
  ∇logsoftmax!(dy, dy, x, dims=dims)

function softmax!(y::DenseCuArray{T}, x::DenseCuArray{T}; dims=1) where T<:CUDNNFloat
  cudnnSoftmaxForward(reshape4D(x), reshape4D(y),
                      algo=CUDNN_SOFTMAX_FAST, mode=cudnnSoftmaxMode_t(dims-1))
  return y
end

function ∇softmax!(dx::DenseCuArray{T}, dy::DenseCuArray{T}, x::DenseCuArray{T}; dims=1) where T<:CUDNNFloat
  y = softmax(x, dims=dims)
  cudnnSoftmaxBackward(reshape4D(y), reshape4D(dy), reshape4D(dx),
                       algo=CUDNN_SOFTMAX_FAST, mode=cudnnSoftmaxMode_t(dims-1))
  return dx
end

function logsoftmax!(y::DenseCuArray{T}, x::DenseCuArray{T}; dims=1) where T<:CUDNNFloat
  cudnnSoftmaxForward(reshape4D(x), reshape4D(y),
                      algo=CUDNN_SOFTMAX_LOG, mode=cudnnSoftmaxMode_t(dims-1))
  return y
end

function ∇logsoftmax!(dx::DenseCuArray{T}, dy::DenseCuArray{T}, x::DenseCuArray{T}; dims=1) where T<:CUDNNFloat
  y = logsoftmax(x, dims=dims)
  cudnnSoftmaxBackward(reshape4D(y), reshape4D(dy), reshape4D(dx),
                       algo=CUDNN_SOFTMAX_LOG, mode=cudnnSoftmaxMode_t(dims-1))
  return dx
end

# Convolution

const conv_forward_algos = DefaultDict{Tuple, Int32}(Int32(-1))
function conv!(y::DenseCuArray{T}, x::DenseCuArray{T}, w::DenseCuArray{T}, cdims::DenseConvDims;
               alpha=1, algo=-1) where T<:CUDNNFloat
  if version() < v"6"
    all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
  end

  if algo < 0
    global conv_forward_algos
    key = (T, strides(x), strides(w), strides(y), cdims, size(x)[end])
    algo = conv_forward_algos[key]
    if algo < 0 # not in conv_forward_algos
      # algo = UInt32(cudnnGetConvolutionForwardAlgorithm(fix1d(y), fix1d(x), fix1d(w), fix1d(cdims), preference=2, workspacesize=workspacesize(x)) # will be removed in cuDNN 8
      # returnedAlgoCount, perfResults = cudnnGetConvolutionForwardAlgorithm_v7(fix1d(y), fix1d(x), fix1d(w), fix1d(cdims))
      # returnedAlgoCount, perfResults = cudnnFindConvolutionForwardAlgorithm(fix1d(y), fix1d(x), fix1d(w), fix1d(cdims))
      returnedAlgoCount, perfResults = cudnnFindConvolutionForwardAlgorithmEx(fix1d(y), fix1d(x), fix1d(w), fix1d(cdims), workspacesize=workspacesize(x))
      algo = perfChoose(perfResults, returnedAlgoCount)
      conv_forward_algos[key] = algo
    end
  end

  cudnnConvolutionForward(fix1d(y), fix1d(x), fix1d(w), fix1d(cdims), alpha=alpha, algo=algo)
  return y
end

if isdefined(NNlib, :conv_bias_act!)
function NNlib.conv_bias_act!(y::DenseCuArray{T}, x::DenseCuArray{T}, w::DenseCuArray{T}, cdims::DenseConvDims, b::DenseCuArray{T}, σ=identity;
                              z::DenseCuArray{T}=y, alpha1=1, alpha2=0, algo=-1) where T<:CUDNNFloat
  if version() < v"6"
    all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
  end

  if algo < 0
    global conv_forward_algos
    key = (T, strides(x), strides(w), strides(y), cdims, size(x)[end])
    algo = conv_forward_algos[key]
    if algo < 0 # not in conv_forward_algos
      # algo = UInt32(cudnnGetConvolutionForwardAlgorithm(fix1d(y), fix1d(x), fix1d(w), fix1d(cdims), preference=2, workspacesize=workspacesize(x))) # will be removed in cuDNN 8
      # returnedAlgoCount, perfResults = cudnnGetConvolutionForwardAlgorithm_v7(fix1d(y), fix1d(x), fix1d(w), fix1d(cdims))
      # returnedAlgoCount, perfResults = cudnnFindConvolutionForwardAlgorithm(fix1d(y), fix1d(x), fix1d(w), fix1d(cdims))
      returnedAlgoCount, perfResults = cudnnFindConvolutionForwardAlgorithmEx(fix1d(y), fix1d(x), fix1d(w), fix1d(cdims), workspacesize=workspacesize(x))
      algo = perfChoose(perfResults, returnedAlgoCount)
      conv_forward_algos[key] = algo
    end
  end

  # only relu and identity are supported
  if σ == NNlib.relu # always merge convolutions, bias, and relu, even when bias is turned off
    cudnnConvolutionBiasActivationForward(fix1d(y), fix1d(x), fix1d(w), fix1d(z), fix1d(b),
                                          fix1d(cdims), algo=algo, alpha1=alpha1, alpha2=alpha2,
                                          activationMode=CUDNN_ACTIVATION_RELU, activationCoeff=0.0)
  elseif algo == 1 && b != nothing # only merge convolution and bias if the fastest algorithm is also the only supported algorithm and the bias is not turned off
    # algo must be CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_PRECOMP_GEMM (1) when activationMode equals CUDNN_ACTIVATION_IDENTITY
    cudnnConvolutionBiasActivationForward(fix1d(y), fix1d(x), fix1d(w), fix1d(z), fix1d(b),
                                          fix1d(cdims), algo=algo, alpha1=alpha1, alpha2=alpha2,
                                          activationMode=CUDNN_ACTIVATION_IDENTITY, activationCoeff=0.0)
    σ.(y)
  else # fallback
    if b == nothing # bias is turned off
      σ.(conv!(y, x, w, cdims, alpha=alpha1, algo=algo))
    else # bias is turned on
      σ.(add_bias(conv!(y, x, w, cdims, alpha=alpha1, algo=algo), b))
    end
  end

  return y
end
end

const conv_data_algos = DefaultDict{Tuple, Int32}(Int32(-1))
function ∇conv_data!(dx::DenseCuArray{T}, dy::DenseCuArray{T}, w::DenseCuArray{T},
                     cdims::DenseConvDims; alpha=1, algo=-1) where T<:CUDNNFloat
  if version() < v"6"
    all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
  end

  if algo < 0
    global conv_data_algos
    key = (T, strides(dx), strides(w), strides(dy), cdims, size(dx)[end])
    algo = conv_data_algos[key]
    if algo < 0 # not in conv_data_algos
      # algo = UInt32(cudnnGetConvolutionBackwardDataAlgorithm(fix1d(dx), fix1d(w), fix1d(dy), fix1d(cdims), preference=2, workspacesize=workspacesize(dx))) # will be removed in cuDNN 8
      # returnedAlgoCount, perfResults = cudnnGetConvolutionBackwardDataAlgorithm_v7(fix1d(dx), fix1d(w), fix1d(dy), fix1d(cdims))
      # returnedAlgoCount, perfResults = cudnnFindConvolutionBackwardDataAlgorithm(fix1d(dx), fix1d(w), fix1d(dy), fix1d(cdims))
      returnedAlgoCount, perfResults = cudnnFindConvolutionBackwardDataAlgorithmEx(fix1d(dx), fix1d(w), fix1d(dy), fix1d(cdims), workspacesize=workspacesize(dx))
      algo = perfChoose(perfResults, returnedAlgoCount)
      conv_data_algos[key] = algo
    end
  end

  cudnnConvolutionBackwardData(fix1d(dx), fix1d(w), fix1d(dy), fix1d(cdims), alpha=alpha, algo=algo)
  return dx
end

const conv_filter_algos = DefaultDict{Tuple, Int32}(Int32(-1))
function ∇conv_filter!(dw::DenseCuArray{T}, x::DenseCuArray{T}, dy::DenseCuArray{T},
                       cdims::DenseConvDims; alpha=1, algo=-1) where T<:CUDNNFloat
  if version() < v"6"
    all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
  end

  if algo < 0
    global conv_filter_algos
    # (type, batchsize, conv descriptor)
    key = (T, strides(x), strides(dw), strides(dy), cdims, size(x)[end])
    algo = conv_filter_algos[key]
    if algo < 0 # not in conv_filter_algos
      # algo = UInt32(cudnnGetConvolutionBackwardFilterAlgorithm(fix1d(dw), fix1d(x), fix1d(dy), fix1d(cdims), preference=2, workspacesize=workspacesize(x))) # will be removed in cuDNN 8
      # returnedAlgoCount, perfResults = cudnnGetConvolutionBackwardFilterAlgorithm_v7(fix1d(dw), fix1d(x), fix1d(dy), fix1d(cdims))
      # returnedAlgoCount, perfResults = cudnnFindConvolutionBackwardFilterAlgorithm(fix1d(dw), fix1d(x), fix1d(dy), fix1d(cdims))
      returnedAlgoCount, perfResults = cudnnFindConvolutionBackwardFilterAlgorithmEx(fix1d(dw), fix1d(x), fix1d(dy), fix1d(cdims), workspacesize=workspacesize(x))
      algo = perfChoose(perfResults, returnedAlgoCount)
      conv_filter_algos[key] = algo
    end
  end

  cudnnConvolutionBackwardFilter(fix1d(dw), fix1d(x), fix1d(dy), fix1d(cdims), alpha=alpha, algo=algo)
  return dw
end

# Bias

# in-place for x (add b to x)
add_bias(x::DenseCuArray{T}, b::DenseCuArray{T})  where {T<:CUDNNFloat} =
  (cudnnAddTensor(reshape4D(x), reshape4D(b)); return x)

∇conv_bias!(db::DenseCuArray{T}, dy::DenseCuArray{T}; alpha=1, beta=0) where T<:CUDNNFloat =
  (cudnnConvolutionBackwardBias(fix1d(db), fix1d(dy), alpha=alpha, beta=beta); return db)

# Pooling

maxpool!(y::DenseCuArray{T}, x::DenseCuArray{T}, pdims::PoolDims) where T<:CUDNNFloat =
  (cudnnPoolingForward(fix1d(y), fix1d(x), fix1d(pdims); mode=0); return y)

∇maxpool!(dx::DenseCuArray{T}, dy::DenseCuArray{T}, y::DenseCuArray{T}, x::DenseCuArray{T}, pdims::PoolDims) where T<:CUDNNFloat =
  (cudnnPoolingBackward(fix1d(dx), fix1d(dy), fix1d(x), fix1d(y), fix1d(pdims), mode=0); return dx)

meanpool!(y::DenseCuArray{T}, x::DenseCuArray{T}, pdims::PoolDims) where T<:CUDNNFloat =
  (cudnnPoolingForward(fix1d(y), fix1d(x), fix1d(pdims), mode=1); return y)

∇meanpool!(dx::DenseCuArray{T}, dy::DenseCuArray{T}, y::DenseCuArray{T}, x::DenseCuArray{T}, pdims::PoolDims) where T<:CUDNNFloat =
  (cudnnPoolingBackward(fix1d(dx), fix1d(dy), fix1d(x), fix1d(y), fix1d(pdims), mode=1); return dx)

# Activation

# in-place for x
Base.broadcasted(::typeof(NNlib.σ), x::DenseCuArray{T}) where {T<:CUDNNFloat} =
  (cudnnActivationForward(reshape4D(x), mode=CUDNN_ACTIVATION_SIGMOID, coeff=0.0); return x)

Base.broadcasted(::typeof(NNlib.relu), x::DenseCuArray{T}) where {T<:CUDNNFloat} =
  (cudnnActivationForward(reshape4D(x), mode=CUDNN_ACTIVATION_RELU, coeff=0.0); return x)

Base.broadcasted(::typeof(NNlib.tanh), x::DenseCuArray{T}) where {T<:CUDNNFloat} =
  (cudnnActivationForward(reshape4D(x), mode=CUDNN_ACTIVATION_TANH, coeff=0.0); return x)

Base.broadcasted(::typeof(NNlib.relu6), x::DenseCuArray{T}) where {T<:CUDNNFloat} =
  (cudnnActivationForward(reshape4D(x), mode=CUDNN_ACTIVATION_CLIPPED_RELU, coeff=6.0); return x)

Base.broadcasted(::typeof(NNlib.elu), x::DenseCuArray{T}) where {T<:CUDNNFloat} =
  (cudnnActivationForward(reshape4D(x), mode=CUDNN_ACTIVATION_ELU, coeff=1.0); return x)

# CUDNN_ACTIVATION_IDENTITY does not work with cudnnActivationForward
Base.broadcasted(::typeof(NNlib.identity), x::DenseCuArray{T}) where {T<:CUDNNFloat} = x

Base.broadcasted(::typeof(NNlib.leakyrelu), x::DenseCuArray{T}, a=T(0.01)) where {T<:CUDNNFloat} =
  (cudnnOpTensor(CUDNN_OP_TENSOR_MAX, reshape4D(x), reshape4D(x), reshape4D(x), alpha1=a); return x)


text:
# by JMW

function unquote(e::Expr)
    @assert e.head == :quote
    return e.args[1]
end

function unquote(e::QuoteNode)
    return e.value
end

macro delegate(source::Expr, targets::Expr)
    typename = esc(source.args[1])
    fieldname = unquote(source.args[2])
    funcnames = targets.args
    n = length(funcnames)
    fdefs = Vector{Any}(undef, n)
    for i in 1:n
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typename), args...) =
                       ($funcname)(a.$fieldname, args...)
                   end
    end
    return Expr(:block, fdefs...)
end

macro delegate_return_parent(source::Expr, targets::Expr)
    typename = esc(source.args[1])
    fieldname = unquote(source.args[2])
    funcnames = targets.args
    n = length(funcnames)
    fdefs = Vector{Any}(undef, n)
    for i in 1:n
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typename), args...) =
                       (($funcname)(a.$fieldname, args...); a)
                   end
    end
    return Expr(:block, fdefs...)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
###
### Mer Iteration
###
###
### Iterator over all k-mers in a biological sequence.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# Note about the variable names:

# it.step is the distance from start of one kmer to start of next
#
# filled is the number of nucleotides in a kmer that has the correct value set.
# e.g. when moving from xxxxxxx to yyyyyyy, filled goes from K to 3.
# when moving from xxxxxxx to zzzzzzz, filled goes from K to 0.
#
# increment is how far the index jumps ahead when going to the next kmer.
# for close kmers where no jump is possible, it's 1, else it can be arbitrary high.
#
# For a kmeriterator that first emits xxxxxxx, then zzzzzzz:
#
#       |--------- step ---------|
#       xxxxxxx                  zzzzzzz
# -------------------------------------------------------------
#           yyyyyyy
#             |---- increment ---|
#
# The state returned at each iteration is the state upon return, not the state
# needed for the following iteration.

struct MerIterResult{T<:AbstractMer}
    position::Int
    fw::T
    bw::T
end

@inline function Base.iterate(x::MerIterResult, state = 1)
    if 1 ≤ state ≤ 3
        return getfield(x, state), state + 1
    else
        return nothing
    end
end

@inline position(x::MerIterResult) = x.position
@inline fwmer(x::MerIterResult) = x.fw
@inline bwmer(x::MerIterResult) = x.bw
@inline canonical(x::MerIterResult) = min(fwmer(x), bwmer(x))

function Base.show(io::IO, ::MIME"text/plain", x::MerIterResult)
    println(io, "Mer iteration result:")
    println(io, "Position: ", x.position)
    print(io, "Forward: ")
    showcompact(io, fwmer(x))
    print(io, "\nBackward: ")
    showcompact(io, bwmer(x))
    print(io, '\n')
end

abstract type AbstractMerIterator{T,S} end

struct EveryMerIterator{T<:AbstractMer,S<:BioSequence} <: AbstractMerIterator{T,S}
    seq::S
    start::Int
    stop::Int
end

struct SpacedMerIterator{T<:AbstractMer,S<:BioSequence} <: AbstractMerIterator{T,S}
    seq::S
    start::Int
    step::Int
    stop::Int
    filled::Int # This is cached for speed
    increment::Int # This is cached for speed
end

"""
    each(::Type{T}, seq::BioSequence) where {T<:AbstractMer}

Initialize an iterator over all overlapping k-mers in a sequence `seq` skipping
ambiguous nucleotides without changing the reading frame.
"""
function each(::Type{T}, seq::BioSequence) where {T<:AbstractMer}
    if eltype(seq) ∉ (DNA, RNA)
        throw(ArgumentError("element type must be either DNA or RNA nucleotide"))
    elseif !(1 ≤ ksize(T) ≤ capacity(T))
        throw(ArgumentError("k-mer length must be between 0 and $(capacity(T))"))
    end
    return EveryMerIterator{T,typeof(seq)}(seq, 1, lastindex(seq))
end

"""
    each(::Type{T}, seq::BioSequence, step::Integer) where {T<:AbstractMer}

Initialize an iterator over k-mers separated by a `step` parameter, in a
sequence `seq` skipping ambiguous nucleotides without changing the reading frame.
"""
function each(::Type{T}, seq::BioSequence, step::Integer) where {T<:AbstractMer}
    if eltype(seq) ∉ (DNA, RNA)
        throw(ArgumentError("element type must be either DNA or RNA nucleotide"))
    elseif !(1 ≤ ksize(T) ≤ capacity(T))
        throw(ArgumentError("k-mer length must be between 0 and $(capacity(T))"))
    elseif step < 1
        throw(ArgumentError("step size must be positive"))
    end
    if step == 1
        return EveryMerIterator{T,typeof(seq)}(seq, 1, lastindex(seq))
    else
        filled = max(0, ksize(T) - step)
        increment = max(1, step - ksize(T) + 1)
        return SpacedMerIterator{T,typeof(seq)}(seq, 1, step, lastindex(seq), filled, increment)
    end
end

@inline Base.eltype(::Type{<:AbstractMerIterator{T,S}}) where {T,S} = MerIterResult{T}
@inline Base.IteratorSize(::Type{<:AbstractMerIterator{T,S}}) where {T,S<:Union{ReferenceSequence,LongSequence{<:NucleicAcidAlphabet{4}}}} = Base.SizeUnknown()
@inline Base.IteratorSize(::Type{<:AbstractMerIterator{T,S}}) where {T,S<:LongSequence{<:NucleicAcidAlphabet{2}}} = Base.HasLength()

@inline function Base.length(it::AbstractMerIterator{T,S}) where {T,S<:LongSequence{<:NucleicAcidAlphabet{2}}}
    return max(0, fld(it.stop - it.start + 1 - ksize(T), step(it)) + 1)
end

Base.step(x::EveryMerIterator) = 1
Base.step(x::SpacedMerIterator) = x.step

# A nucleotide with bitvalue B has kmer-bitvalue kmerbits[B+1].
# ambiguous nucleotides have no kmervalue, here set to 0xff
const kmerbits = (0xff, 0x00, 0x01, 0xff,
                  0x02, 0xff, 0xff, 0xff,
                  0x03, 0xff, 0xff, 0xff,
                  0xff, 0xff, 0xff, 0xff)

# Initializers for two-bit nucleic acid alphabets...
@inline function Base.iterate(it::AbstractMerIterator{T,S}) where {T,S<:LongSequence{<:NucleicAcidAlphabet{2}}}
    UT = encoded_data_type(T)
    filled, i = 0, it.start
    fwkmer = rvkmer = zero(UT)

    while i ≤ it.stop
        nt = reinterpret(UInt8, inbounds_getindex(it.seq, i))
        @inbounds fbits = UT(kmerbits[nt + 1])
        rbits = ~fbits & typeof(fbits)(0x03)
        fwkmer = (fwkmer << 0x02 | fbits)
        rvkmer = (rvkmer >> 0x02) | (UT(rbits) << unsigned(offset(T, 1)))
        filled += 1
        if filled == ksize(T)
            return MerIterResult(1, T(fwkmer), T(rvkmer)), (i, fwkmer, rvkmer)
        end
        i += 1
    end
    return nothing
end

@inline function Base.iterate(it::EveryMerIterator{T,S}, state) where {T,S<:LongSequence{<:NucleicAcidAlphabet{2}}}
    UT = encoded_data_type(T)
    i, fwkmer, rvkmer = state
    i += 1

    if i > it.stop
        return nothing
    else
        nt = reinterpret(UInt8, inbounds_getindex(it.seq, i))
        @inbounds fbits = UT(kmerbits[nt + 1])
        rbits = ~fbits & typeof(fbits)(0x03)
        fwkmer = (fwkmer << 0x02 | fbits)
        rvkmer = (rvkmer >> 0x02) | (UT(rbits) << unsigned(offset(T, 1)))
        pos = i - ksize(T) + 1

        return MerIterResult(pos, T(fwkmer), T(rvkmer)), (i, fwkmer, rvkmer)
    end
end

@inline function Base.iterate(it::SpacedMerIterator{T,S}, state
    ) where {T,S<:LongSequence{<:NucleicAcidAlphabet{2}}}
    UT = encoded_data_type(T)
    i, fkmer, rkmer = state
    filled = it.filled
    i += it.increment

    while i ≤ it.stop
        nt = reinterpret(Int8, inbounds_getindex(it.seq, i))
        @inbounds val = kmerbits[nt + 1]
        rbits = ~val & typeof(val)(0x03)
        fkmer = fkmer << 2 | val
        rkmer = rkmer >> 2 | (UT(rbits) << unsigned(offset(T, 1)))
        filled += 1
        if filled == ksize(T)
            pos = i - ksize(T) + 1
            return MerIterResult(pos, T(fkmer), T(rkmer)), (i, fkmer, rkmer)
        end
        i += 1
    end
    return nothing
end

@inline function Base.iterate(it::EveryMerIterator{T,S}, state=(it.start-1,1,encoded_data_type(T)(0),encoded_data_type(T)(0))
    ) where {T,S<:Union{ReferenceSequence, LongSequence{<:NucleicAcidAlphabet{4}}}}

    UT = encoded_data_type(T)
    i, filled, fkmer, rkmer = state
    i += 1
    filled -= 1

    while i ≤ it.stop
        nt = reinterpret(Int8, inbounds_getindex(it.seq, i))
        @inbounds fbits = UT(kmerbits[nt + 1])
        rbits = ~fbits & typeof(fbits)(0x03)
        fkmer = fkmer << 0x02 | fbits
        rkmer = (rkmer >> 0x02) | (UT(rbits) << unsigned(offset(T, 1)))
        filled = ifelse(fbits == 0xff, 0, filled + 1)
        if filled == ksize(T)
            pos = i - ksize(T) + 1
            return MerIterResult(pos, T(fkmer), T(rkmer)), (i, ksize(T), fkmer, rkmer)
        end
        i += 1
    end
    return nothing
end

@inline function Base.iterate(it::SpacedMerIterator{T,S}, state=(it.start-it.increment, 1, 0, zero(encoded_data_type(T)), zero(encoded_data_type(T)))
    ) where {T,S<:Union{ReferenceSequence, LongSequence{<:NucleicAcidAlphabet{4}}}}
    UT = encoded_data_type(T)
    i, pos, filled, fwkmer, rvkmer = state
    i += it.increment

    while i ≤ it.stop
        nt = reinterpret(Int8, inbounds_getindex(it.seq, i))
        @inbounds val = kmerbits[nt + 1]
        rbits = ~val & typeof(val)(0x03)
        if val == 0xff # ambiguous
            filled = 0
            # Find the beginning of next possible kmer after i
            pos = i + it.step - Core.Intrinsics.urem_int(i-pos, it.step)
            i = pos - 1
        else
            filled += 1
            fwkmer = fwkmer << 2 | val
            rvkmer = (rvkmer >> 0x02) | UT(rbits) << unsigned(offset(T, 1))
        end
        if filled == ksize(T)
            state = (i, i - ksize(T) + 1 + it.step, it.filled, fwkmer, rvkmer)
            return MerIterResult(pos, T(fwkmer), T(rvkmer)), state
        end
        i += 1
    end
    return nothing
end


text:
function descent(p::Param, η::Real)
  function ()
    @. p.x -= η * p.Δ
    @. p.Δ = 0
  end
end

function momentum(p::Param, ρ, η)
  v = zeros(p.x)
  function ()
    @. v = ρ * v - η * p.Δ
    @. p.Δ = -v
  end
end

# Ref. https://arxiv.org/pdf/1212.0901.pdf
function nesterov(p::Param, ρ, η)
  v = zeros(p.x)
  function ()
    d = @. ρ^2 * v - (1+ρ) * η * p.Δ
    @. v = ρ*v - η*p.Δ
    @. p.Δ = -d
  end
end

function rmsprop(p::Param; η::Real = 0.001, ρ::Real = 0.9, ϵ::Real = 1e-8)
  acc  = zeros(p.x)
  function ()
    @. acc = ρ * acc + (1 - ρ) * p.Δ^2
    @. p.Δ *= η / √(acc + ϵ)
  end
end

function adagrad(p::Param; η::Real = 0.01, ϵ::Real = 1e-8)
  acc = zeros(p.x) .+ ϵ
  function ()
    @. acc += p.Δ^2
    @. p.Δ *= η / √(acc + ϵ)
  end
end

function adadelta(p::Param; ρ::Real = 0.9, ϵ::Real = 1e-8)
  acc = zeros(p.x)
  Δacc = zeros(p.x)
  function ()
    @. acc = ρ * acc + (1 - ρ) * p.Δ^2
    @. p.Δ *= √(Δacc + ϵ) / √(acc + ϵ)
    @. Δacc = ρ * Δacc + (1 - ρ) * p.Δ^2
   end
end

function adam(p::Param; η::Real = 0.001, β1::Real = 0.9, β2::Real = 0.999, ϵ::Real = 1e-8)
  mt = zeros(p.x)
  vt = zeros(p.x)
  β1p, β2p = β1, β2
  function ()
    @. mt = β1 * mt + (1 - β1) * p.Δ
    @. vt = β2 * vt + (1 - β2) * p.Δ^2
    @. p.Δ =  mt / (1 - β1p) / √(vt / (1 - β2p) + ϵ) * η
    β1p *= β1
    β2p *= β2
  end
end

function adamax(p::Param; η::Real = 0.002, β1::Real = 0.9, β2::Real = 0.999, ϵ::Real = 1e-8)
  mt = zeros(p.x)
  ut = zeros(p.x)
  β1p = β1
  function ()
    @. mt = β1 * mt + (1 - β1) * p.Δ
    @. ut = max(β2 * ut, abs(p.Δ))
    @. p.Δ = (η/(1 - β1p)) * mt/(ut + ϵ)
    β1p *= β1
  end
end

function amsgrad(p::Param; η::Real = 0.001, β1::Real = 0.9, β2::Real = 0.999, ϵ::Real = 1e-8)
  mt = zeros(p.x)
  vt = zeros(p.x) .+ ϵ
  v̂t = zeros(p.x) .+ ϵ
  function ()
    @. mt = β1 * mt + (1 - β1) * p.Δ
    @. vt = β2 * vt + (1 - β2) * p.Δ ^ 2
    @. v̂t = max.(v̂t, vt)
    @. p.Δ = η * mt / √v̂t
  end
end

function nadam(p::Param; η::Real = 0.001, β1::Real = 0.9, β2::Real = 0.999, ϵ::Real = 1e-8)
  mt = zeros(p.x)
  vt = zeros(p.x)
  β1p, β2p = β1, β2
  function ()
    @. mt = β1 * mt + (1 - β1) * p.Δ
    @. vt = β2 * vt + (1 - β2) * p.Δ^2
    @. p.Δ = (β1 * mt / (1 - β1 * β1p) + (1 - β1) * p.Δ / (1 - β1p)) / √(vt * β2 / (1 - β2p) + ϵ) * η
    β1p *= β1
    β2p *= β2
  end
end

clip(p::Param, thresh::Real) = () -> clamp!(p.Δ, -thresh, thresh)

function expdecay(p::Param, γ::Real)
  if γ != 0
    return () -> p.Δ .+= γ .* p.x
  else
    return () -> nothing
  end
end

function invdecay(p::Param, γ::Real)
  if γ != 0
    n = 0
    return () -> begin
      p.Δ .*= 1 / (1 + γ * n)
      n += 1
    end
  else
    return () -> nothing
  end
end


text:
using Test

@testset "promotion" begin
    @testset "number_promotions" begin
    @test promote_type(Int, BD{Int}) == BD{Int64}
    @test promote_type(BD{Int}, Int) == BD{Int64}

    @test promote_type(Float64, BD{Float64}) == BD{Float64}
    @test promote_type(BD{Float64}, Float64) == BD{Float64}
    end

    @testset "wrapped_number_promotions" begin
    @test promote_type(BD{Int}, BD{Float64}) == BD{Float64}
    @test promote_type(BD{Float64}, BD{Int}) == BD{Float64}
    end

    @testset "wrappedarray_number_promotions" begin
    @test promote_type(BD{Array{Int}}, Int) == Any
    @test promote_type(BD{Int}, Array{Int}) == Any

    @test promote_type(BD{Array{Float64}}, Float64) == Any
    @test promote_type(BD{Float64}, Array{Float64}) == Any

    @test promote_type(BD{Array{Float64}}, Int) == Any
    @test promote_type(Array{Float64}, BD{Int}) == Any

    @test promote_type(BD{Array{Int}}, BD{Float64}) == Any
    @test promote_type(BD{Array{Float64}}, BD{Array{Int}}) == BD{Array{Float64}}
    end

end


text:
#
# Example script to wrap PETSc
#

using Clang.cindex
using Clang.wrap_c
using Compat

import Clang.wrap_c.repr_jl

include("rewriter.jl")
#=
function repr_jl(ptr::cindex.Pointer)
  ptee = pointee_type(ptr)
  ex1 = Expr(:curly, :Ptr, repr_jl(ptee))
  ex2 = Expr(:curly, :AbstractArray, repr_jl(ptee))
  Expr(:curly, Union, ex1, ex2)
end
=#
PETSC_INCLUDE = "/home/jared/build/petsc-3.6.0/include"
#PETSC_INCLUDE = "/home/jared/.julia/v0.4/PETSc/deps/petsc-3.6.0/arch-linux2-c-opt/include"
MPI_INCLUDE = "/usr/include/mpich"

petsc_header = [joinpath(PETSC_INCLUDE, "petscksp.h")]
#petsc_header = [joinpath(PETSC_INCLUDE, "petsc.h")]
h1 = "/usr/include"
h2 = joinpath(PETSC_INCLUDE, "petsc/mpiuni")
h3 = joinpath(PETSC_INCLUDE, "petsc/private")
h4 = joinpath(PETSC_INCLUDE, "petsc/private/kernels")

# Set up include paths
clang_includes = ASCIIString[]
push!(clang_includes, PETSC_INCLUDE)
#push!(clang_includes, MPI_INCLUDE)
#push!(clang_includes, h1)
#push!(clang_includes, h2)
#push!(clang_includes, h3)
#push!(clang_includes, h4)

println("clang_includes = ", clang_includes)
# Clang arguments
clang_extraargs = ["-v"]
# clang_extraargs = ["-D", "__STDC_LIMIT_MACROS", "-D", "__STDC_CONSTANT_MACROS"]

# Callback to test if a header should actually be wrapped (for exclusion)
function wrap_header(top_hdr::ASCIIString, cursor_header::ASCIIString)
    return startswith(dirname(cursor_header), PETSC_INCLUDE)
end

lib_file(hdr::ASCIIString) = "petsc"
output_file(hdr::ASCIIString) = "PETSc.jl"

function wrap_cursor(name::ASCIIString, cursor)
    println("name = ", name)
    exc = true
    if length(name) > 0
      if name[1] == '_'
        exc = false
      end
#    exc = false
#    exc |= contains(name, "MPI")
    end
    return exc
end

const wc = wrap_c.init(;
                        headers = petsc_header,
                        output_file = "libPETSc_h.jl",
                        common_file = "libPETSc_common.jl",
                        clang_includes      = clang_includes,
                        clang_args          = clang_extraargs,
                        header_wrapped      = wrap_header,
                        header_library      = lib_file,
                        header_outputfile   = output_file,
                        cursor_wrapped      = wrap_cursor,
                        clang_diagnostics = true,
                        rewriter = petsc_rewriter)

run(wc)


text:
function extrap_prep{N,d}(::Type{Linear}, ::Val{N}, ::Val{d}, ::Val{:lo})
    coords = [Symbol("xs_", k) for k in 1:N]
    xs_d = coords[d]
    quote
        if $xs_d < lbound(etp.itp, $d)
            $xs_d = lbound(etp.itp, $d)
            return etp[$(coords...)] + gradient(etp, $(coords...))[$d] * (xs[$d] - $xs_d)
        end
    end
end
function extrap_prep{N,d}(::Type{Linear}, ::Val{N}, ::Val{d}, ::Val{:hi})
    coords = [Symbol("xs_", k) for k in 1:N]
    xs_d = coords[d]
    quote
        if $xs_d > ubound(etp, $d)
            $xs_d = ubound(etp, $d)
            return etp[$(coords...)] + gradient(etp, $(coords...))[$d] * (xs[$d] - $xs_d)
        end
    end
end

extrap_prep{N,d,l}(::Val{:gradient}, ::Type{Linear}, n::Val{N}, dim::Val{d}, lohi::Val{l}) =
    extrap_prep(Flat, n, dim, lohi)

extrap_prep{N,d}(::Val{:gradient}, ::Type{Linear}, n::Val{N}, dim::Val{d}) =
    extrap_prep(Flat, n, dim)


text:
using ForwardDiff
using LaTeXStrings
using StatsPlots
using Plots
pyplot()

let

"""Function"""

# Objective
f(x::Vector) = -(x[1]*x[2] + x[2]*x[3] + x[1]*x[3])
Ω(x::Vector) = sum((abs.(x) .- 1).^2)
h(x::Vector) = f(x) + Ω(x)

# Minimization step
function g(x::Vector, i::Integer)

	if i == 1
		# Closed-form solution
		g_i = (1 + abs(x[2] + x[3])/2)*sign(x[2] + x[3])

	elseif i == 2
		# Closed form solution
		g_i = (1 + abs(x[1] + x[3])/2)*sign(x[1] + x[3])

	elseif i == 3
		# Closed form solution
		g_i = (1 + abs(x[1] + x[2])/2)*sign(x[1] + x[2])
	end
	return g_i
end

"""Coordinate Descent"""

# Number of iterations
K = 2

# Fix step size
α = 1.0

# Initialize coordinates
ϵ = 1e-3
x = [-1. - ϵ, 1. + ϵ/2, -1. - ϵ/4]

# Preallocate result array
global xx = zeros(K*3,3)

# Initalize visualization
scatter([x[1]], [x[2]], [x[3]],
		color="red",
		markerstrokecolor=nothing,
		label="x^0",
		xlims=(-1.5, 1.5),
		ylims=(-1.5, 1.5),
		zlims=(-1.5, 1.5),
		xlabel=L"x_1",
		ylabel=L"x_2",
		zlabel=L"x_3")

# Visualize stationary points
scatter!([+1], [+1], [+1], color="green", label="min", markerstrokecolor=nothing)
scatter!([-1], [-1], [-1], color="green", label="", markerstrokecolor=nothing)

# Loop over iterations
anim = @animate for k = 1:K*3

	# Keep track of x
	println(x)
	xx[k, :] = x

	# Update index
	i = (k-1) % 3 + 1
	println(i)

	# Update coordinate
	x[i] = g(x, i)

	# Visualization
	plot!([xx[k,1], x[1]],
	 	  [xx[k,2], x[2]],
		  [xx[k,3], x[3]],
	      color="blue",
	 	  label="",
	 	  markerstrokecolor=nothing)
end

gif(anim, "powells.gif", fps=2)

end

max_stars_repo_name: wmkouw/CD-powells

max_issues_repo_name: wmkouw/CD-powells

max_forks_repo_name: wmkouw/CD-powells


text:
#!/usr/bin/env julia

using Pkg
Pkg.activate(".")
using PhysicsEngine

PhysicsEngine.simulate()

max_stars_repo_name: jonathanbuchanan/PhysicsEngine

max_issues_repo_name: jonathanbuchanan/PhysicsEngine

max_forks_repo_name: jonathanbuchanan/PhysicsEngine


text:
module SimpleTasksTests
include("utils/include.jl")
include("tasks/include.jl")
include("services/include.jl")
include("examples/include.jl")
end # module CloudTest

max_forks_repo_name: seung-lab/SimpleTasks